From a6c311f41069365cc40f252fdbdb17b37a2e2c5f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Nov 2021 12:16:32 +0100 Subject: [PATCH 1/5] commented --- src/lattice.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 4b87fccbf..637d585e5 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -548,6 +548,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & !-------------------------------------------------------------------------------------------------- ! elasticity matrix of the target phase in cube orientation if (lattice_target == 'hP') then + ! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19) + ! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48) if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) & call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target)) C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal @@ -2091,7 +2093,7 @@ function lattice_equivalent_nu(C,assumption) result(nu) K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & / 9.0_pReal elseif(IO_lc(assumption) == 'reuss') then - call math_invert(S,error,C) + call math_invert(S,error,C) ! ToDo: correct for Voigt? if(error) error stop 'matrix inversion failed' K = 1.0_pReal & / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3))) @@ -2123,7 +2125,7 @@ function lattice_equivalent_mu(C,assumption) result(mu) mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) & / 15.0_pReal elseif(IO_lc(assumption) == 'reuss') then - call math_invert(S,error,C) + call math_invert(S,error,C) ! ToDo: correct for Voigt? if(error) error stop 'matrix inversion failed' mu = 15.0_pReal & / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6))) From 87100d1dce0bc2264b5e403b82ad141bd8521a27 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Nov 2021 15:15:59 +0100 Subject: [PATCH 2/5] avoid conversion Voigt-6x6 to 3x3x3x3 --- src/lattice.f90 | 4 +-- src/rotations.f90 | 69 ++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 61 insertions(+), 12 deletions(-) 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 From f6440f7f17003b8d3adac4e8e5176d62dd3702a5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Nov 2021 15:18:48 +0100 Subject: [PATCH 3/5] need to check 'any', adjusted tolerance --- src/rotations.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index d4f226065..aeb0cb48c 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1491,20 +1491,20 @@ 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 (any(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 (any(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 (any(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))) & + if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-12_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 From a857285e345ae4ef456277214cd35dedae191505 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Nov 2021 22:09:01 +0100 Subject: [PATCH 4/5] backward transformation for Voigt --- src/math.f90 | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 2515c64d7..48c736f76 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -853,6 +853,40 @@ pure function math_66toSym3333(m66,weighted) end function math_66toSym3333 +!-------------------------------------------------------------------------------------------------- +!> @brief Convert 6 Voigt stress vector into symmetric 3x3 tensor. +!-------------------------------------------------------------------------------------------------- +pure function math_Voigt6to33_stress(sigma_tilde) result(sigma) + + real(pReal), dimension(3,3) :: sigma + real(pReal), dimension(6), intent(in) :: sigma_tilde + + + sigma = reshape([sigma_tilde(1), sigma_tilde(6), sigma_tilde(5), & + sigma_tilde(6), sigma_tilde(2), sigma_tilde(4), & + sigma_tilde(5), sigma_tilde(4), sigma_tilde(3)],[3,3]) + + +end function math_Voigt6to33_stress + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert 6 Voigt strain vector into symmetric 3x3 tensor. +!-------------------------------------------------------------------------------------------------- +pure function math_Voigt6to33_strain(epsilon_tilde) result(epsilon) + + real(pReal), dimension(3,3) :: epsilon + real(pReal), dimension(6), intent(in) :: epsilon_tilde + + + epsilon = reshape([ epsilon_tilde(1), 0.5_pReal*epsilon_tilde(6), 0.5_pReal*epsilon_tilde(5), & + 0.5_pReal*epsilon_tilde(6), epsilon_tilde(2), 0.5_pReal*epsilon_tilde(4), & + 0.5_pReal*epsilon_tilde(5), 0.5_pReal*epsilon_tilde(4), epsilon_tilde(3)],[3,3]) + + +end function math_Voigt6to33_strain + + !-------------------------------------------------------------------------------------------------- !> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. !-------------------------------------------------------------------------------------------------- @@ -864,12 +898,12 @@ pure function math_Voigt66to3333(m66) integer :: i,j - do i=1,6; do j=1, 6 + do i=1,6; do j=1,6 math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) - enddo; enddo + end do; end do end function math_Voigt66to3333 From 28ff027b388917a80e0305ab66017bfd443ec59a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 20 Nov 2021 22:36:01 +0100 Subject: [PATCH 5/5] exploit symmetry for stress calculation --- src/math.f90 | 37 ++++++++++++++++++++++++++++---- src/phase_mechanical_elastic.f90 | 6 ++++-- 2 files changed, 37 insertions(+), 6 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 48c736f76..aa81a0ed3 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -841,14 +841,14 @@ pure function math_66toSym3333(m66,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL - endif + end if do i=1,6; do j=1,6 math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(1,j),MAPNYE(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(MAPNYE(2,i),MAPNYE(1,i),MAPNYE(2,j),MAPNYE(1,j)) = w(i)*w(j)*m66(i,j) - enddo; enddo + end do; end do end function math_66toSym3333 @@ -866,7 +866,6 @@ pure function math_Voigt6to33_stress(sigma_tilde) result(sigma) sigma_tilde(6), sigma_tilde(2), sigma_tilde(4), & sigma_tilde(5), sigma_tilde(4), sigma_tilde(3)],[3,3]) - end function math_Voigt6to33_stress @@ -883,10 +882,40 @@ pure function math_Voigt6to33_strain(epsilon_tilde) result(epsilon) 0.5_pReal*epsilon_tilde(6), epsilon_tilde(2), 0.5_pReal*epsilon_tilde(4), & 0.5_pReal*epsilon_tilde(5), 0.5_pReal*epsilon_tilde(4), epsilon_tilde(3)],[3,3]) - end function math_Voigt6to33_strain +!-------------------------------------------------------------------------------------------------- +!> @brief Convert 3x3 tensor into 6 Voigt stress vector. +!-------------------------------------------------------------------------------------------------- +pure function math_33toVoigt6_stress(sigma) result(sigma_tilde) + + real(pReal), dimension(6) :: sigma_tilde + real(pReal), dimension(3,3), intent(in) :: sigma + + + sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), & + sigma(3,2), sigma(3,1), sigma(1,2)] + +end function math_33toVoigt6_stress + + +!-------------------------------------------------------------------------------------------------- +!> @brief Convert 3x3 tensor into 6 Voigt strain vector. +!-------------------------------------------------------------------------------------------------- +pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde) + + real(pReal), dimension(6) :: epsilon_tilde + real(pReal), dimension(3,3), intent(in) :: epsilon + + + epsilon_tilde = [ epsilon(1,1), epsilon(2,2), epsilon(3,3), & + 2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] + +end function math_33toVoigt6_strain + + + !-------------------------------------------------------------------------------------------------- !> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 5792ea3cd..bb0422bd1 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -148,15 +148,17 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient real(pReal), dimension(3,3) :: E + real(pReal), dimension(6,6) :: C66 real(pReal), dimension(3,3,3,3) :: C integer :: & i, j - C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)) + C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en) + C = math_Voigt66to3333(C66) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration - S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration + S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration do i =1,3; do j=1,3 dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko