223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

getting rid of conversions with unclear behavior

This commit is contained in:
Martin Diehl 2019-09-20 07:08:21 -07:00
parent 76eaa9855f
commit 42fba28fa1
2 changed files with 36 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
src/lattice.f90
View File

@ -11,6 +11,7 @@ module lattice
use IO use IO
use config use config
use math use math
use rotations
implicit none implicit none
private private
@ -916,8 +917,8 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
real(pReal),dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem
real(pReal), dimension(3,3) :: R type(rotation) :: R
integer :: i integer :: i
if (len_trim(structure) /= 3) & if (len_trim(structure) /= 3) &
@ -938,9 +939,10 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
end select end select
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
R = math_axisAngleToR(coordinateSystem(1:3,2,i), PI) ! ToDo: Why always 180 deg? call R%fromAxisAnglePair([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = math_sym3333to66(math_rotate_forward3333(math_66toSym3333(C66),R)) lattice_C66_twin(1:6,1:6,i) = math_sym3333to66(math_rotate_forward3333(math_66toSym3333(C66),R%asRotationMatrix()))
enddo enddo
end function lattice_C66_twin end function lattice_C66_twin
@ -1000,6 +1002,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
do i = 1, sum(Ntrans) do i = 1, sum(Ntrans)
lattice_C66_trans(1:6,1:6,i) = math_sym3333to66(math_rotate_forward3333(C_target_unrotated,Q(1:3,1:3,i))) lattice_C66_trans(1:6,1:6,i) = math_sym3333to66(math_rotate_forward3333(C_target_unrotated,Q(1:3,1:3,i)))
enddo enddo
end function lattice_C66_trans end function lattice_C66_trans
@ -1016,7 +1019,8 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system
real(pReal), dimension(:), allocatable :: direction, normal, np real(pReal), dimension(3) :: direction, normal, np
type(rotation) :: R
integer :: i integer :: i
if (abs(sense) /= 1) call IO_error(0,ext_msg='lattice_nonSchmidMatrix') if (abs(sense) /= 1) call IO_error(0,ext_msg='lattice_nonSchmidMatrix')
@ -1029,7 +1033,9 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
do i = 1,sum(Nslip) do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i) direction = coordinateSystem(1:3,1,i)
normal = coordinateSystem(1:3,2,i) normal = coordinateSystem(1:3,2,i)
np = matmul(math_axisAngleToR(direction,60.0_pReal*INRAD), normal) call R%fromAxisAnglePair([direction,60.0_pReal],degrees=.true.,P=1)
np = R%rotVector(normal)
if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
+ nonSchmidCoefficients(1) * math_outer(direction, np) + nonSchmidCoefficients(1) * math_outer(direction, np)
if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
@ -1044,6 +1050,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
+ nonSchmidCoefficients(6) * math_outer(direction, direction) + nonSchmidCoefficients(6) * math_outer(direction, direction)
enddo enddo
end function lattice_nonSchmidMatrix end function lattice_nonSchmidMatrix
@ -2143,10 +2150,11 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
a_bcc, & !< lattice parameter a for target bcc structure a_bcc, & !< lattice parameter a for target bcc structure
a_fcc !< lattice parameter a for parent fcc structure a_fcc !< lattice parameter a for parent fcc structure
real(pReal), dimension(3,3) :: & type(rotation) :: &
R, & !< Pitsch rotation R, & !< Pitsch rotation
B !< Rotation of fcc to Bain coordinate system
real(pReal), dimension(3,3) :: &
U, & !< Bain deformation U, & !< Bain deformation
B, & !< Rotation of fcc to Bain coordinate system
ss, sd ss, sd
real(pReal), dimension(3) :: & real(pReal), dimension(3) :: &
x, y, z x, y, z
@ -2170,17 +2178,17 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter :: & real(pReal), dimension(4,LATTICE_fcc_Ntrans), parameter :: &
LATTICE_FCCTOBCC_SYSTEMTRANS = reshape([& LATTICE_FCCTOBCC_SYSTEMTRANS = reshape([&
0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3) 0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3)
0.0, 1.0, 0.0, -10.26, & 0.0,-1.0, 0.0, 10.26, &
0.0, 0.0, 1.0, 10.26, & 0.0, 0.0, 1.0, 10.26, &
0.0, 0.0, 1.0, -10.26, & 0.0, 0.0,-1.0, 10.26, &
1.0, 0.0, 0.0, 10.26, & 1.0, 0.0, 0.0, 10.26, &
1.0, 0.0, 0.0, -10.26, & -1.0, 0.0, 0.0, 10.26, &
0.0, 0.0, 1.0, 10.26, & 0.0, 0.0, 1.0, 10.26, &
0.0, 0.0, 1.0, -10.26, & 0.0, 0.0,-1.0, 10.26, &
1.0, 0.0, 0.0, 10.26, & 1.0, 0.0, 0.0, 10.26, &
1.0, 0.0, 0.0, -10.26, & -1.0, 0.0, 0.0, 10.26, &
0.0, 1.0, 0.0, 10.26, & 0.0, 1.0, 0.0, 10.26, &
0.0, 1.0, 0.0, -10.26 & 0.0,-1.0, 0.0, 10.26 &
],shape(LATTICE_FCCTOBCC_SYSTEMTRANS)) ],shape(LATTICE_FCCTOBCC_SYSTEMTRANS))
integer, dimension(9,LATTICE_fcc_Ntrans), parameter :: & integer, dimension(9,LATTICE_fcc_Ntrans), parameter :: &
@ -2220,19 +2228,17 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
if (a_bcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation if (a_bcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
R = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & call R%fromAxisAnglePair(LATTICE_FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1)
lattice_fccTobcc_systemTrans(4,i)*INRAD) call B%fromAxisAnglePair(LATTICE_FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1)
B = math_axisAngleToR(lattice_fccTobcc_bainRot(1:3,i), & x = real(LATTICE_FCCTOBCC_BAINVARIANT(1:3,i),pReal)
lattice_fccTobcc_bainRot(4,i)*INRAD) y = real(LATTICE_FCCTOBCC_BAINVARIANT(4:6,i),pReal)
x = real(LATTICE_fccTobcc_bainVariant(1:3,i),pReal) z = real(LATTICE_FCCTOBCC_BAINVARIANT(7:9,i),pReal)
y = real(LATTICE_fccTobcc_bainVariant(4:6,i),pReal)
z = real(LATTICE_fccTobcc_bainVariant(7:9,i),pReal)
U = (a_bcc/a_fcc)*math_outer(x,x) & U = (a_bcc/a_fcc)*math_outer(x,x) &
+ (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) & + (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) &
+ (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal) + (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal)
Q(1:3,1:3,i) = matmul(R,B) Q(1:3,1:3,i) = matmul(R%asRotationMatrix(),B%asRotationMatrix())
S(1:3,1:3,i) = matmul(R,U) - MATH_I3 S(1:3,1:3,i) = matmul(R%asRotationMatrix(),U) - MATH_I3
enddo enddo
elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex transformation elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex transformation
ss = MATH_I3 ss = MATH_I3
@ -2242,8 +2248,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
sd(3,3) = cOverA/sqrt(8.0_pReal/3.0_pReal) sd(3,3) = cOverA/sqrt(8.0_pReal/3.0_pReal)
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
x = lattice_fccTohex_systemTrans(1:3,i)/norm2(lattice_fccTohex_systemTrans(1:3,i)) x = LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(1:3,i))
z = lattice_fccTohex_systemTrans(4:6,i)/norm2(lattice_fccTohex_systemTrans(4:6,i)) z = LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i)/norm2(LATTICE_FCCTOHEX_SYSTEMTRANS(4:6,i))
y = -math_cross(x,z) y = -math_cross(x,z)
Q(1:3,1,i) = x Q(1:3,1,i) = x
Q(1:3,2,i) = y Q(1:3,2,i) = y

47
src/math.f90
View File

@ -968,10 +968,7 @@ end function math_qRot
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians) !> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians)
!> @details rotation matrix is meant to represent a PASSIVE rotation, composed of INTRINSIC !> @details deprecated
!> @details rotations around the axes of the details rotating reference frame.
!> @details similar to eu2om from "D Rowenhorst et al. Consistent representations of and conversions
!> @details between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)", but R is transposed
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_EulerToR(Euler) pure function math_EulerToR(Euler)
@ -1003,48 +1000,6 @@ pure function math_EulerToR(Euler)
end function math_EulerToR end function math_EulerToR
!--------------------------------------------------------------------------------------------------
!> @brief rotation matrix from axis and angle (in radians)
!> @details rotation matrix is meant to represent a ACTIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!> @details formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html
!> @details equivalent to eu2om (P=-1) from "D Rowenhorst et al. Consistent representations of and
!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)"
!--------------------------------------------------------------------------------------------------
pure function math_axisAngleToR(axis,omega)
real(pReal), dimension(3,3) :: math_axisAngleToR
real(pReal), dimension(3), intent(in) :: axis
real(pReal), intent(in) :: omega
real(pReal), dimension(3) :: n
real(pReal) :: norm,s,c,c1
norm = norm2(axis)
wellDefined: if (norm > 1.0e-8_pReal) then
n = axis/norm ! normalize axis to be sure
s = sin(omega)
c = cos(omega)
c1 = 1.0_pReal - c
math_axisAngleToR(1,1) = c + c1*n(1)**2.0_pReal
math_axisAngleToR(1,2) = c1*n(1)*n(2) - s*n(3)
math_axisAngleToR(1,3) = c1*n(1)*n(3) + s*n(2)
math_axisAngleToR(2,1) = c1*n(1)*n(2) + s*n(3)
math_axisAngleToR(2,2) = c + c1*n(2)**2.0_pReal
math_axisAngleToR(2,3) = c1*n(2)*n(3) - s*n(1)
math_axisAngleToR(3,1) = c1*n(1)*n(3) - s*n(2)
math_axisAngleToR(3,2) = c1*n(2)*n(3) + s*n(1)
math_axisAngleToR(3,3) = c + c1*n(3)**2.0_pReal
else wellDefined
math_axisAngleToR = math_I3
endif wellDefined
end function math_axisAngleToR
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Gauss variable !> @brief draw a random sample from Gauss variable
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------