avoid conversion Voigt-6x6 to 3x3x3x3

This commit is contained in:
Martin Diehl 2021-11-20 15:15:59 +01:00
parent a6c311f410
commit 87100d1dce
2 changed files with 61 additions and 12 deletions

View File

@ -522,7 +522,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66)))
lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
enddo
end function lattice_C66_twin
@ -583,7 +583,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C_target_unrotated66)))
lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
enddo
end function lattice_C66_trans

View File

@ -75,6 +75,7 @@ module rotations
procedure, public :: rotVector
procedure, public :: rotTensor2
procedure, public :: rotTensor4
procedure, public :: rotStiffness
procedure, public :: misorientation
procedure, public :: standardize
end type rotation
@ -339,8 +340,7 @@ end function rotTensor2
!---------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief rotate a rank-4 tensor passively (default) or actively
!> @brief Rotate a rank-4 tensor passively (default) or actively.
!> @details: rotation is based on rotation matrix
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
@ -354,6 +354,7 @@ pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3) :: R
integer :: i,j,k,l,m,n,o,p
if (present(active)) then
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else
@ -371,7 +372,47 @@ end function rotTensor4
!---------------------------------------------------------------------------------------------------
!> @brief misorientation
!> @brief Rotate a rank-4 tensor in Voigt 6x6 notation passively (default) or actively.
!> @details: https://scicomp.stackexchange.com/questions/35600
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
pure function rotStiffness(self,C,active) result(cRot)
real(pReal), dimension(6,6) :: cRot
class(rotation), intent(in) :: self
real(pReal), intent(in), dimension(6,6) :: C
logical, intent(in), optional :: active
real(pReal), dimension(3,3) :: R
real(pReal), dimension(6,6) :: M
if (present(active)) then
R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else
R = self%asMatrix()
endif
M = reshape([R(1,1)**2.0_pReal, R(2,1)**2.0_pReal, R(3,1)**2.0_pReal, &
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
R(1,2)**2.0_pReal, R(2,2)**2.0_pReal, R(3,2)**2.0_pReal, &
R(2,2)*R(3,2), R(1,2)*R(3,2), R(1,2)*R(2,2), &
R(1,3)**2.0_pReal, R(2,3)**2.0_pReal, R(3,3)**2.0_pReal, &
R(2,3)*R(3,3), R(1,3)*R(3,3), R(1,3)*R(2,3), &
2.0_pReal*R(1,2)*R(1,3), 2.0_pReal*R(2,2)*R(2,3), 2.0_pReal*R(3,2)*R(3,3), &
R(2,2)*R(3,3)+R(2,3)*R(3,2), R(1,2)*R(3,3)+R(1,3)*R(3,2), R(1,2)*R(2,3)+R(1,3)*R(2,2), &
2.0_pReal*R(1,3)*R(1,1), 2.0_pReal*R(2,3)*R(2,1), 2.0_pReal*R(3,3)*R(3,1), &
R(2,3)*R(3,1)+R(2,1)*R(3,3), R(1,3)*R(3,1)+R(1,1)*R(3,3), R(1,3)*R(2,1)+R(1,1)*R(2,3), &
2.0_pReal*R(1,1)*R(1,2), 2.0_pReal*R(2,1)*R(2,2), 2.0_pReal*R(3,1)*R(3,2), &
R(2,1)*R(3,2)+R(2,2)*R(3,1), R(1,1)*R(3,2)+R(1,2)*R(3,1), R(1,1)*R(2,2)+R(1,2)*R(2,1)],[6,6])
cRot = matmul(M,matmul(C,transpose(M)))
end function rotStiffness
!---------------------------------------------------------------------------------------------------
!> @brief Misorientation.
!---------------------------------------------------------------------------------------------------
pure elemental function misorientation(self,other)
@ -386,7 +427,7 @@ end function misorientation
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to rotation matrix
!> @brief Convert unit quaternion to rotation matrix.
!---------------------------------------------------------------------------------------------------
pure function qu2om(qu) result(om)
@ -395,8 +436,8 @@ pure function qu2om(qu) result(om)
real(pReal) :: qq
qq = qu(1)**2-sum(qu(2:4)**2)
qq = qu(1)**2-sum(qu(2:4)**2)
om(1,1) = qq+2.0_pReal*qu(2)**2
om(2,2) = qq+2.0_pReal*qu(3)**2
@ -416,7 +457,7 @@ end function qu2om
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to Euler angles
!> @brief Convert unit quaternion to Euler angles.
!---------------------------------------------------------------------------------------------------
pure function qu2eu(qu) result(eu)
@ -425,6 +466,7 @@ pure function qu2eu(qu) result(eu)
real(pReal) :: q12, q03, chi
q03 = qu(1)**2+qu(4)**2
q12 = qu(2)**2+qu(3)**2
chi = sqrt(q03*q12)
@ -1379,6 +1421,7 @@ subroutine selfTest()
real(pReal), dimension(3) :: x, eu, ho, v3
real(pReal), dimension(3,3) :: om, t33
real(pReal), dimension(3,3,3,3) :: t3333
real(pReal), dimension(6,6) :: C
real :: A,B
integer :: i
@ -1412,6 +1455,7 @@ subroutine selfTest()
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om'
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu'
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax'
@ -1447,20 +1491,25 @@ subroutine selfTest()
call R%fromMatrix(om)
call random_number(v3)
if(all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
if (all(dNeq(R%rotVector(R%rotVector(v3),active=.true.),v3,1.0e-12_pReal))) &
error stop 'rotVector'
call random_number(t33)
if(all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
if (all(dNeq(R%rotTensor2(R%rotTensor2(t33),active=.true.),t33,1.0e-12_pReal))) &
error stop 'rotTensor2'
call random_number(t3333)
if(all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
if (all(dNeq(R%rotTensor4(R%rotTensor4(t3333),active=.true.),t3333,1.0e-12_pReal))) &
error stop 'rotTensor4'
call random_number(C)
C = C+transpose(C)
if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-14_pReal))) &
error stop 'rotStiffness'
call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML
enddo
end do
contains