diff --git a/src/lattice.f90 b/src/lattice.f90 index 637d585e5..d90215c7e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -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 diff --git a/src/rotations.f90 b/src/rotations.f90 index e73fb2782..d4f226065 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -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