From 29771feaaed13bd99c26f77d9b46709e92549c7c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 14:43:05 +0100 Subject: [PATCH 01/13] cleaning --- src/phase_mechanical_elastic.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 9c9a159a0..587614865 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -106,7 +106,7 @@ module function elastic_C66(ph,en) result(C66) associate(prm => param(ph)) - C66 = get_C66(ph,en) + C66 = get_C66(ph,en) C66 = math_sym3333to66(math_Voigt66to3333(C66)) ! Literature data is in Voigt notation end associate @@ -126,7 +126,7 @@ module function elastic_mu(ph,en) result(mu) mu - mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') + mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') end function elastic_mu @@ -142,7 +142,7 @@ module function elastic_nu(ph,en) result(nu) nu - nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') + nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') end function elastic_nu From 72c07cfc17b3c239f19759df80cb71ba82cebf6f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 16:25:51 +0100 Subject: [PATCH 02/13] 'present' propagates --- src/phase_mechanical_elastic.f90 | 9 +-------- src/rotations.f90 | 8 +++----- 2 files changed, 4 insertions(+), 13 deletions(-) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 587614865..4d7037830 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -100,16 +100,9 @@ module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & en - real(pReal), dimension(6,6) :: & - C66 - associate(prm => param(ph)) - - C66 = get_C66(ph,en) - C66 = math_sym3333to66(math_Voigt66to3333(C66)) ! Literature data is in Voigt notation - - end associate + C66 = math_sym3333to66(math_Voigt66to3333(get_C66(ph,en))) ! Literature data is in Voigt notation end function elastic_C66 diff --git a/src/rotations.f90 b/src/rotations.f90 index 68f9273a3..1ca114acb 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -383,11 +383,8 @@ pure function rotTensor4sym(self,T,active) result(tRot) real(pReal), intent(in), dimension(6,6) :: T logical, intent(in), optional :: active - if (present(active)) then - tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) - else - tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T))) - endif + + tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) end function rotTensor4sym @@ -400,6 +397,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other + misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) end function misorientation From 8d64a1c2f24a975dd1fc84fcca8685cd4321110d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 16:37:34 +0100 Subject: [PATCH 03/13] mark compressed notation --- src/homogenization_mechanical_RGC.f90 | 2 +- src/phase.f90 | 6 +++--- src/phase_damage.f90 | 2 +- src/phase_mechanical_elastic.f90 | 6 +++--- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 71d1542c8..f0c2900e2 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -652,7 +652,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) real(pReal), dimension(6,6) :: C - C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) + C = phase_homogenizedC66(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu diff --git a/src/phase.f90 b/src/phase.f90 index a36fc26b9..9c4212c71 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -230,10 +230,10 @@ module phase end function phase_mechanical_constitutive !ToDo: Merge all the stiffness functions - module function phase_homogenizedC(ph,en) result(C) + module function phase_homogenizedC66(ph,en) result(C) integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C - end function phase_homogenizedC + end function phase_homogenizedC66 module function phase_damage_C(C_homogenized,ph,en) result(C) real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized integer, intent(in) :: ph,en @@ -299,7 +299,7 @@ module phase public :: & phase_init, & - phase_homogenizedC, & + phase_homogenizedC66, & phase_f_phi, & phase_f_T, & phase_K_phi, & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 49b9f9528..8e4b66ae1 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -417,7 +417,7 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) + call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if (.not. broken) then myOffset = damageState(ph)%offsetDeltaState diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 4d7037830..c9e5b6490 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -166,7 +166,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j - C = math_66toSym3333(phase_homogenizedC(ph,en)) + C = math_66toSym3333(phase_homogenizedC66(ph,en)) C = phase_damage_C(C,ph,en) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration @@ -184,7 +184,7 @@ end subroutine phase_hooke_SandItsTangents !> @brief returns the homogenized elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module function phase_homogenizedC(ph,en) result(C) +module function phase_homogenizedC66(ph,en) result(C) real(pReal), dimension(6,6) :: C integer, intent(in) :: ph, en @@ -197,7 +197,7 @@ module function phase_homogenizedC(ph,en) result(C) C = elastic_C66(ph,en) end select plasticType -end function phase_homogenizedC +end function phase_homogenizedC66 end submodule elastic From f7a42bdc1ac176339714c035aa9e8401c461d489 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 16:56:36 +0100 Subject: [PATCH 04/13] avoid conversion to 3333 --- src/phase.f90 | 10 +++++----- src/phase_damage.f90 | 16 +++++++++------- src/phase_mechanical_elastic.f90 | 5 ++--- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 9c4212c71..2159dba00 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -234,11 +234,11 @@ module phase integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C end function phase_homogenizedC66 - module function phase_damage_C(C_homogenized,ph,en) result(C) - real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized - integer, intent(in) :: ph,en - real(pReal), dimension(3,3,3,3) :: C - end function phase_damage_C + module function phase_damage_C66(C66_homogenized,ph,en) result(C66) + real(pReal), dimension(6,6), intent(in) :: C66_homogenized + integer, intent(in) :: ph,en + real(pReal), dimension(6,6) :: C66 + end function phase_damage_C66 module function phase_f_phi(phi,co,ce) result(f) integer, intent(in) :: ce,co diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 8e4b66ae1..78322d7f4 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -139,6 +139,7 @@ module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) integer :: & ph, en + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) @@ -150,20 +151,21 @@ end function phase_damage_constitutive !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- -module function phase_damage_C(C_homogenized,ph,en) result(C) +module function phase_damage_C66(C66_homogenized,ph,en) result(C66) + + real(pReal), dimension(6,6), intent(in) :: C66_homogenized + integer, intent(in) :: ph,en + real(pReal), dimension(6,6) :: C66 - real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized - integer, intent(in) :: ph,en - real(pReal), dimension(3,3,3,3) :: C damageType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) damageType - C = C_homogenized * damage_phi(ph,en)**2 + C66 = C66_homogenized * damage_phi(ph,en)**2 case default damageType - C = C_homogenized + C66 = C66_homogenized end select damageType -end function phase_damage_C +end function phase_damage_C66 !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index c9e5b6490..d7b8d5d1c 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -100,7 +100,7 @@ module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & en - + real(pReal), dimension(6,6) :: C66 C66 = math_sym3333to66(math_Voigt66to3333(get_C66(ph,en))) ! Literature data is in Voigt notation @@ -166,8 +166,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j - C = math_66toSym3333(phase_homogenizedC66(ph,en)) - C = phase_damage_C(C,ph,en) + C = math_66toSym3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)) 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 From 4d29393cedf21e47f51febcc2120c356034d1243 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 16:59:23 +0100 Subject: [PATCH 05/13] simplified --- src/phase_mechanical_elastic.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index d7b8d5d1c..afde0b087 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -102,6 +102,7 @@ module function elastic_C66(ph,en) result(C66) en real(pReal), dimension(6,6) :: C66 + C66 = math_sym3333to66(math_Voigt66to3333(get_C66(ph,en))) ! Literature data is in Voigt notation end function elastic_C66 @@ -119,7 +120,7 @@ module function elastic_mu(ph,en) result(mu) mu - mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') + mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') end function elastic_mu @@ -135,7 +136,7 @@ module function elastic_nu(ph,en) result(nu) nu - nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') + nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') end function elastic_nu From 8a8fdfc93cc8b7c52e686f75240db25aafac45e3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 17:33:08 +0100 Subject: [PATCH 06/13] use Voigt notation there is no advantage of using the symmetrized conversions --- src/phase_mechanical_elastic.f90 | 57 ++++++++-------------- src/phase_mechanical_plastic_dislotwin.f90 | 2 +- 2 files changed, 22 insertions(+), 37 deletions(-) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index afde0b087..55aa4cefb 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -64,38 +64,7 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- !> @brief Return 6x6 elasticity tensor. !-------------------------------------------------------------------------------------------------- -function get_C66(ph,en) - - integer, intent(in) :: & - ph, & - en - real(pReal), dimension(6,6) :: get_C66 - - - associate(prm => param(ph)) - get_C66 = 0.0_pReal - get_C66(1,1) = prm%C_11 - get_C66(1,2) = prm%C_12 - get_C66(4,4) = prm%C_44 - - if (any(phase_lattice(ph) == ['hP','tI'])) then - get_C66(1,3) = prm%C_13 - get_C66(3,3) = prm%C_33 - end if - - if (phase_lattice(ph) == 'tI') get_C66(6,6) = prm%C_66 - - get_C66 = lattice_symmetrize_C66(get_C66,phase_lattice(ph)) - - end associate - -end function get_C66 - - -!-------------------------------------------------------------------------------------------------- -!> @brief Return 6x6 elasticity tensor. -!-------------------------------------------------------------------------------------------------- -module function elastic_C66(ph,en) result(C66) +function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & @@ -103,7 +72,22 @@ module function elastic_C66(ph,en) result(C66) real(pReal), dimension(6,6) :: C66 - C66 = math_sym3333to66(math_Voigt66to3333(get_C66(ph,en))) ! Literature data is in Voigt notation + associate(prm => param(ph)) + C66 = 0.0_pReal + C66(1,1) = prm%C_11 + C66(1,2) = prm%C_12 + C66(4,4) = prm%C_44 + + if (any(phase_lattice(ph) == ['hP','tI'])) then + C66(1,3) = prm%C_13 + C66(3,3) = prm%C_33 + end if + + if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66 + + C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) + + end associate end function elastic_C66 @@ -120,10 +104,11 @@ module function elastic_mu(ph,en) result(mu) mu - mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') + mu = lattice_equivalent_mu(elastic_C66(ph,en),'voigt') end function elastic_mu + !-------------------------------------------------------------------------------------------------- !> @brief Return Poisson ratio. !-------------------------------------------------------------------------------------------------- @@ -136,7 +121,7 @@ module function elastic_nu(ph,en) result(nu) nu - nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') + nu = lattice_equivalent_nu(elastic_C66(ph,en),'voigt') end function elastic_nu @@ -194,7 +179,7 @@ module function phase_homogenizedC66(ph,en) result(C) case (PLASTICITY_DISLOTWIN_ID) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType - C = elastic_C66(ph,en) + C = math_sym3333to66(math_Voigt66to3333(elastic_C66(ph,en))) end select plasticType end function phase_homogenizedC66 diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 62119c819..934773ad7 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -484,7 +484,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) associate(prm => param(ph), stt => state(ph)) - C66 = elastic_C66(ph,en) + C66 = math_sym3333to66(math_Voigt66to3333(elastic_C66(ph,en))) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & From 020ef64d7d2a732ac4fe891bcad33bda69d8e08a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 17:44:06 +0100 Subject: [PATCH 07/13] explicit conversions --- src/lattice.f90 | 5 +++-- src/phase_mechanical_elastic.f90 | 2 +- src/rotations.f90 | 19 ------------------- 3 files changed, 4 insertions(+), 22 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 11d1dca0f..7fa6b8ab5 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -505,6 +505,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA) type(rotation) :: R integer :: i + select case(lattice) case('cF') coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& @@ -521,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) = R%rotTensor4sym(C66) + lattice_C66_twin(1:6,1:6,i) = math_sym3333to66(R%rotTensor4(math_66toSym3333(C66))) enddo end function lattice_C66_twin @@ -580,7 +581,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) = R%rotTensor4sym(C_target_unrotated66) + lattice_C66_trans(1:6,1:6,i) = math_sym3333to66(R%rotTensor4(math_66toSym3333(C_target_unrotated66))) enddo end function lattice_C66_trans diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 55aa4cefb..8ab6884f9 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -64,7 +64,7 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- !> @brief Return 6x6 elasticity tensor. !-------------------------------------------------------------------------------------------------- -function elastic_C66(ph,en) result(C66) +module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & diff --git a/src/rotations.f90 b/src/rotations.f90 index 1ca114acb..e73fb2782 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -75,7 +75,6 @@ module rotations procedure, public :: rotVector procedure, public :: rotTensor2 procedure, public :: rotTensor4 - procedure, public :: rotTensor4sym procedure, public :: misorientation procedure, public :: standardize end type rotation @@ -371,24 +370,6 @@ pure function rotTensor4(self,T,active) result(tRot) end function rotTensor4 -!--------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief rotate a symmetric rank-4 tensor stored as (6,6) passively (default) or actively -!! ToDo: Need to check active/passive !!! -!--------------------------------------------------------------------------------------------------- -pure function rotTensor4sym(self,T,active) result(tRot) - - real(pReal), dimension(6,6) :: tRot - class(rotation), intent(in) :: self - real(pReal), intent(in), dimension(6,6) :: T - logical, intent(in), optional :: active - - - tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) - -end function rotTensor4sym - - !--------------------------------------------------------------------------------------------------- !> @brief misorientation !--------------------------------------------------------------------------------------------------- From dfe6d0a1957f516fe2c1d8d38329224ede53dcd9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 21:06:38 +0100 Subject: [PATCH 08/13] more support for Voigt notation --- src/math.f90 | 31 ++++++++++++++++++++++++++++++- src/phase_mechanical_elastic.f90 | 1 - 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 31bd3b952..2515c64d7 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -854,12 +854,13 @@ end function math_66toSym3333 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 Voigt matrix into symmetric 3x3x3x3 matrix +!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. !-------------------------------------------------------------------------------------------------- pure function math_Voigt66to3333(m66) real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix + integer :: i,j @@ -873,6 +874,31 @@ pure function math_Voigt66to3333(m66) end function math_Voigt66to3333 +!-------------------------------------------------------------------------------------------------- +!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix. +!-------------------------------------------------------------------------------------------------- +pure function math_3333toVoigt66(m3333) + + real(pReal), dimension(6,6) :: math_3333toVoigt66 + real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check) + + integer :: i,j + + +#ifndef __INTEL_COMPILER + do concurrent(i=1:6, j=1:6) + math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + end do +#else + do i=1,6; do j=1,6 + math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + end do; end do +#endif + +end function math_3333toVoigt66 + + + !-------------------------------------------------------------------------------------------------- !> @brief draw a random sample from Gauss variable !-------------------------------------------------------------------------------------------------- @@ -1261,6 +1287,9 @@ subroutine selfTest if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' + if(any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + error stop 'math_3333toVoigt66/math_Voigt66to3333' + call random_number(v6) if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & error stop 'math_symmetric33' diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 8ab6884f9..0f7394e86 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -167,7 +167,6 @@ end subroutine phase_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenized elasticity matrix -!> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- module function phase_homogenizedC66(ph,en) result(C) From 038dd1fc40dadd20ffb52f5289d89f69167faf18 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 21:31:08 +0100 Subject: [PATCH 09/13] correct names whether C66 is homogenized or not is a decision of the caller --- src/phase.f90 | 6 +++--- src/phase_damage.f90 | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index 2159dba00..dc8525939 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -234,10 +234,10 @@ module phase integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C end function phase_homogenizedC66 - module function phase_damage_C66(C66_homogenized,ph,en) result(C66) - real(pReal), dimension(6,6), intent(in) :: C66_homogenized + module function phase_damage_C66(C66,ph,en) + real(pReal), dimension(6,6), intent(in) :: C66 integer, intent(in) :: ph,en - real(pReal), dimension(6,6) :: C66 + real(pReal), dimension(6,6) :: phase_damage_C66 end function phase_damage_C66 module function phase_f_phi(phi,co,ce) result(f) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 78322d7f4..5382d65e3 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -151,18 +151,18 @@ end function phase_damage_constitutive !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- -module function phase_damage_C66(C66_homogenized,ph,en) result(C66) +module function phase_damage_C66(C66,ph,en) - real(pReal), dimension(6,6), intent(in) :: C66_homogenized + real(pReal), dimension(6,6), intent(in) :: C66 integer, intent(in) :: ph,en - real(pReal), dimension(6,6) :: C66 + real(pReal), dimension(6,6) :: phase_damage_C66 damageType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) damageType - C66 = C66_homogenized * damage_phi(ph,en)**2 + phase_damage_C66 = C66 * damage_phi(ph,en)**2 case default damageType - C66 = C66_homogenized + phase_damage_C66 = C66 end select damageType end function phase_damage_C66 From fa8218124ace56b392ec1b639379106e7738fce7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Nov 2021 21:59:09 +0100 Subject: [PATCH 10/13] avoid conversions --- src/homogenization_mechanical_RGC.f90 | 8 ++++---- src/phase.f90 | 4 ++-- src/phase_damage.f90 | 8 ++++---- src/phase_damage_isobrittle.f90 | 6 +++++- src/phase_mechanical_elastic.f90 | 8 ++++---- 5 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index f0c2900e2..36a335370 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -643,17 +643,17 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) !------------------------------------------------------------------------------------------------- !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !------------------------------------------------------------------------------------------------- - real(pReal) function equivalentMu(grainID,ce) + real(pReal) function equivalentMu(co,ce) integer, intent(in) :: & - grainID,& + co,& ce real(pReal), dimension(6,6) :: C - C = phase_homogenizedC66(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) - equivalentMu = lattice_equivalent_mu(C,'voigt') + C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) + equivalentMu = lattice_equivalent_mu(math_sym3333to66(math_Voigt66to3333(C)),'voigt') !ToDo: Bug, should be Voigt not sym end function equivalentMu diff --git a/src/phase.f90 b/src/phase.f90 index dc8525939..22c35416b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -234,10 +234,10 @@ module phase integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C end function phase_homogenizedC66 - module function phase_damage_C66(C66,ph,en) + module function phase_damage_C66(C66,ph,en) result(C66_degraded) real(pReal), dimension(6,6), intent(in) :: C66 integer, intent(in) :: ph,en - real(pReal), dimension(6,6) :: phase_damage_C66 + real(pReal), dimension(6,6) :: C66_degraded end function phase_damage_C66 module function phase_f_phi(phi,co,ce) result(f) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 5382d65e3..3b3ace8ab 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -151,18 +151,18 @@ end function phase_damage_constitutive !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- -module function phase_damage_C66(C66,ph,en) +module function phase_damage_C66(C66,ph,en) result(C66_degraded) real(pReal), dimension(6,6), intent(in) :: C66 integer, intent(in) :: ph,en - real(pReal), dimension(6,6) :: phase_damage_C66 + real(pReal), dimension(6,6) :: C66_degraded damageType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) damageType - phase_damage_C66 = C66 * damage_phi(ph,en)**2 + C66_degraded = C66 * damage_phi(ph,en)**2 case default damageType - phase_damage_C66 = C66 + C66_degraded = C66 end select damageType end function phase_damage_C66 diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 25ce1d3e9..9d5b4f508 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -96,6 +96,7 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state +! ToDo: Use Voigt directly !-------------------------------------------------------------------------------------------------- module subroutine isobrittle_deltaState(C, Fe, ph,en) @@ -109,13 +110,16 @@ module subroutine isobrittle_deltaState(C, Fe, ph,en) epsilon real(pReal) :: & r_W + real(pReal), dimension(6,6) :: & + C_sym + C_sym = math_sym3333to66(math_Voigt66to3333(C)) epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) - r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit + r_W = 2.0_pReal*dot_product(epsilon,matmul(C_sym,epsilon))/prm%W_crit dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en)) end associate diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 0f7394e86..38bca44b4 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -152,12 +152,12 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j - C = math_66toSym3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)) + C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)) 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 - do i =1, 3;do j=1,3 + 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 dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn end do; end do @@ -176,9 +176,9 @@ module function phase_homogenizedC66(ph,en) result(C) plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_DISLOTWIN_ID) plasticType - C = plastic_dislotwin_homogenizedC(ph,en) + C = math_3333toVoigt66(math_66toSym3333(plastic_dislotwin_homogenizedC(ph,en))) case default plasticType - C = math_sym3333to66(math_Voigt66to3333(elastic_C66(ph,en))) + C = elastic_C66(ph,en) end select plasticType end function phase_homogenizedC66 From ff9fa1d4f7eb9b318b5b35cc804ba2d0add5f4ff Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Nov 2021 07:29:40 +0100 Subject: [PATCH 11/13] using Voigt notation instead of proprietary scaled 6x6 notation Note: This results in a change of behavior for the transformation systems of dislotwin. I assume that this fixes a bug, but still need to confirm where the equations in lattice_C66_trans come from --- src/lattice.f90 | 8 ++++---- src/phase_mechanical_elastic.f90 | 5 +++-- src/phase_mechanical_plastic_dislotwin.f90 | 18 ++++++++++-------- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 7fa6b8ab5..4b87fccbf 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_sym3333to66(R%rotTensor4(math_66toSym3333(C66))) + lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66))) enddo end function lattice_C66_twin @@ -572,16 +572,16 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) endif - do i = 1, 6 + do i = 1,6 if (abs(C_target_unrotated66(i,i)) @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law +! ToDo: Use Voigt matrix directly !-------------------------------------------------------------------------------------------------- module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi, ph, en) @@ -166,7 +167,7 @@ end subroutine phase_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenized elasticity matrix +!> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- module function phase_homogenizedC66(ph,en) result(C) @@ -176,7 +177,7 @@ module function phase_homogenizedC66(ph,en) result(C) plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_DISLOTWIN_ID) plasticType - C = math_3333toVoigt66(math_66toSym3333(plastic_dislotwin_homogenizedC(ph,en))) + C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType C = elastic_C66(ph,en) end select plasticType diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 934773ad7..de027ae44 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -473,27 +473,29 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) integer, intent(in) :: & ph, en real(pReal), dimension(6,6) :: & - homogenizedC, & - C66 + homogenizedC + + real(pReal), dimension(6,6) :: & + C real(pReal), dimension(:,:,:), allocatable :: & C66_tw, & C66_tr - integer :: i real(pReal) :: f_unrotated - associate(prm => param(ph), stt => state(ph)) - C66 = math_sym3333to66(math_Voigt66to3333(elastic_C66(ph,en))) + C = elastic_C66(ph,en) + + associate(prm => param(ph), stt => state(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * C66 + homogenizedC = f_unrotated * C twinActive: if (prm%sum_N_tw > 0) then - C66_tw = lattice_C66_twin(prm%N_tw,C66,phase_lattice(ph),phase_cOverA(ph)) + C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) do i=1,prm%sum_N_tw homogenizedC = homogenizedC & + stt%f_tw(i,en)*C66_tw(1:6,1:6,i) @@ -501,7 +503,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) end if twinActive transActive: if (prm%sum_N_tr > 0) then - C66_tr = lattice_C66_trans(prm%N_tr,C66,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) + C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) do i=1,prm%sum_N_tr homogenizedC = homogenizedC & + stt%f_tr(i,en)*C66_tr(1:6,1:6,i) From 021d614dafd57a25b0635cd2f6802a678c32ca42 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Nov 2021 09:57:40 +0100 Subject: [PATCH 12/13] using Voigt notation This is a bugfix with a change of behavior --- PRIVATE | 2 +- src/homogenization_mechanical_RGC.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/PRIVATE b/PRIVATE index 277b4a693..76bb51348 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 277b4a693a57c1e4583bcf8ea2205f9b5d147198 +Subproject commit 76bb51348de75207d483d369628670e5ae51dca9 diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 36a335370..76676726a 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -652,8 +652,8 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) real(pReal), dimension(6,6) :: C - C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) - equivalentMu = lattice_equivalent_mu(math_sym3333to66(math_Voigt66to3333(C)),'voigt') !ToDo: Bug, should be Voigt not sym + C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) ! damage not included! + equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu From da23c916cadd4e0bdb18d2d086635c28464c546d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 21 Nov 2021 15:49:04 -0500 Subject: [PATCH 13/13] polish --- src/lattice.f90 | 80 +++++++++--------- src/math.f90 | 94 +++++++++++----------- src/phase_mechanical_elastic.f90 | 10 +-- src/phase_mechanical_plastic_dislotwin.f90 | 28 +++---- 4 files changed, 105 insertions(+), 107 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 4b87fccbf..0e24d1b18 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -405,7 +405,7 @@ module lattice contains !-------------------------------------------------------------------------------------------------- -!> @brief Module initialization +!> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init @@ -417,7 +417,7 @@ end subroutine lattice_init !-------------------------------------------------------------------------------------------------- -!> @brief Characteristic shear for twinning +!> @brief characteristic shear for twinning !-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) @@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for twinning in 66-vector notation +!> @brief rotated elasticity matrices for twinning in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,lattice,CoverA) @@ -529,7 +529,7 @@ end function lattice_C66_twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for transformation in 66-vector notation +!> @brief rotated elasticity matrices for transformation in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & cOverA_trans,a_bcc,a_fcc) @@ -588,7 +588,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & !-------------------------------------------------------------------------------------------------- -!> @brief Non-schmid projections for bcc with up to 6 coefficients +!> @brief non-Schmid projections for bcc with up to 6 coefficients ! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1 !-------------------------------------------------------------------------------------------------- @@ -635,7 +635,7 @@ end function lattice_nonSchmidMatrix !-------------------------------------------------------------------------------------------------- -!> @brief Slip-slip interaction matrix +!> @brief slip-slip interaction matrix !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) @@ -883,7 +883,7 @@ end function lattice_interaction_SlipBySlip !-------------------------------------------------------------------------------------------------- -!> @brief Twin-twin interaction matrix +!> @brief twin-twin interaction matrix !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) @@ -981,7 +981,7 @@ end function lattice_interaction_TwinByTwin !-------------------------------------------------------------------------------------------------- -!> @brief Trans-trans interaction matrix +!> @brief trans-trans interaction matrix !> details only active trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) @@ -1010,7 +1010,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc - if(lattice == 'cF') then + if (lattice == 'cF') then interactionTypes = FCC_INTERACTIONTRANSTRANS NtransMax = FCC_NTRANSSYSTEM else @@ -1023,7 +1023,7 @@ end function lattice_interaction_TransByTrans !-------------------------------------------------------------------------------------------------- -!> @brief Slip-twin interaction matrix +!> @brief slip-twin interaction matrix !> details only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) @@ -1186,7 +1186,7 @@ end function lattice_interaction_SlipByTwin !-------------------------------------------------------------------------------------------------- -!> @brief Slip-trans interaction matrix +!> @brief slip-trans interaction matrix !> details only active slip and trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) @@ -1239,7 +1239,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) !-------------------------------------------------------------------------------------------------- -!> @brief Twin-slip interaction matrix +!> @brief twin-slip interaction matrix !> details only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) @@ -1411,7 +1411,7 @@ end function lattice_SchmidMatrix_twin !-------------------------------------------------------------------------------------------------- -!> @brief Schmid matrix for twinning +!> @brief Schmid matrix for transformation !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) @@ -1483,7 +1483,7 @@ end function lattice_SchmidMatrix_cleavage !-------------------------------------------------------------------------------------------------- -!> @brief Slip direction of slip systems (|| b) +!> @brief slip direction of slip systems (|| b) !-------------------------------------------------------------------------------------------------- function lattice_slip_direction(Nslip,lattice,cOverA) result(d) @@ -1501,7 +1501,7 @@ end function lattice_slip_direction !-------------------------------------------------------------------------------------------------- -!> @brief Normal direction of slip systems (|| n) +!> @brief normal direction of slip systems (|| n) !-------------------------------------------------------------------------------------------------- function lattice_slip_normal(Nslip,lattice,cOverA) result(n) @@ -1519,7 +1519,7 @@ end function lattice_slip_normal !-------------------------------------------------------------------------------------------------- -!> @brief Transverse direction of slip systems (|| t = b x n) +!> @brief transverse direction of slip systems (|| t = b x n) !-------------------------------------------------------------------------------------------------- function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) @@ -1537,7 +1537,7 @@ end function lattice_slip_transverse !-------------------------------------------------------------------------------------------------- -!> @brief Labels for slip systems +!> @brief labels of slip systems !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_labels_slip(Nslip,lattice) result(labels) @@ -1578,7 +1578,7 @@ end function lattice_labels_slip !-------------------------------------------------------------------------------------------------- -!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice +!> @brief return 3x3 tensor with symmetry according to given Bravais lattice !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrize_33(T,lattice) result(T_sym) @@ -1605,7 +1605,7 @@ end function lattice_symmetrize_33 !-------------------------------------------------------------------------------------------------- -!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice +!> @brief return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) @@ -1651,7 +1651,7 @@ end function lattice_symmetrize_C66 !-------------------------------------------------------------------------------------------------- -!> @brief Labels for twin systems +!> @brief labels of twin systems !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_labels_twin(Ntwin,lattice) result(labels) @@ -1689,7 +1689,7 @@ end function lattice_labels_twin !-------------------------------------------------------------------------------------------------- -!> @brief Projection of the transverse direction onto the slip plane +!> @brief projection of the transverse direction onto the slip plane !> @details: This projection is used to calculate forest hardening for edge dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_transverse(Nslip,lattice,cOverA) result(projection) @@ -1713,7 +1713,7 @@ end function slipProjection_transverse !-------------------------------------------------------------------------------------------------- -!> @brief Projection of the slip direction onto the slip plane +!> @brief projection of the slip direction onto the slip plane !> @details: This projection is used to calculate forest hardening for screw dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_direction(Nslip,lattice,cOverA) result(projection) @@ -1779,7 +1779,7 @@ end function coordinateSystem_slip !-------------------------------------------------------------------------------------------------- -!> @brief Populate reduced interaction matrix +!> @brief populate reduced interaction matrix !-------------------------------------------------------------------------------------------------- function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix) @@ -1822,7 +1822,7 @@ end function buildInteraction !-------------------------------------------------------------------------------------------------- -!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems +!> @brief build a local coordinate system on slip, twin, trans, cleavage systems !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- function buildCoordinateSystem(active,potential,system,lattice,cOverA) @@ -1889,7 +1889,7 @@ end function buildCoordinateSystem !-------------------------------------------------------------------------------------------------- -!> @brief Helper function to define transformation systems +!> @brief helper function to define transformation systems ! Needed to calculate Schmid matrix and rotated stiffness matrices. ! @details: set c/a = 0.0 for fcc -> bcc transformation ! set a_Xcc = 0.0 for fcc -> hex transformation @@ -2073,7 +2073,7 @@ end function getlabels !-------------------------------------------------------------------------------------------------- -!> @brief Equivalent Poisson's ratio (ν) +!> @brief equivalent Poisson's ratio (ν) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- function lattice_equivalent_nu(C,assumption) result(nu) @@ -2087,12 +2087,12 @@ function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then 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 + elseif (IO_lc(assumption) == 'reuss') then call math_invert(S,error,C) - if(error) error stop 'matrix inversion failed' + 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))) else @@ -2100,13 +2100,13 @@ function lattice_equivalent_nu(C,assumption) result(nu) endif mu = lattice_equivalent_mu(C,assumption) - nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu) + nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu) end function lattice_equivalent_nu !-------------------------------------------------------------------------------------------------- -!> @brief Equivalent shear modulus (μ) +!> @brief equivalent shear modulus (μ) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- function lattice_equivalent_mu(C,assumption) result(mu) @@ -2119,12 +2119,12 @@ function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then 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 + elseif (IO_lc(assumption) == 'reuss') then call math_invert(S,error,C) - if(error) error stop 'matrix inversion failed' + 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))) else @@ -2135,7 +2135,7 @@ end function lattice_equivalent_mu !-------------------------------------------------------------------------------------------------- -!> @brief Check correctness of some lattice functions. +!> @brief check correctness of some lattice functions !-------------------------------------------------------------------------------------------------- subroutine selfTest @@ -2153,7 +2153,7 @@ subroutine selfTest system = reshape([1.0_pReal+r(1),0.0_pReal,0.0_pReal, 0.0_pReal,1.0_pReal+r(2),0.0_pReal],[6,1]) CoSy = buildCoordinateSystem([1],[1],system,'cF',0.0_pReal) - if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' + if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' do i = 1, 10 call random_number(C) @@ -2199,13 +2199,13 @@ subroutine selfTest C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) C = lattice_symmetrize_C66(C,'cI') - if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' - if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' lambda = C(1,2) - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt' - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss' end subroutine selfTest diff --git a/src/math.f90 b/src/math.f90 index 2515c64d7..5d2609237 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -15,15 +15,15 @@ module math implicit none public #if __INTEL_COMPILER >= 1900 - ! do not make use associated entities available to other modules + ! do not make use of associated entities available to other modules private :: & IO, & config #endif real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian + real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree + real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) real(pReal), dimension(3,3), parameter :: & @@ -132,19 +132,19 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, intent(in),optional :: istart,iend, sortDim integer :: ipivot,s,e,d - if(present(istart)) then + if (present(istart)) then s = istart else s = lbound(a,2) endif - if(present(iend)) then + if (present(iend)) then e = iend else e = ubound(a,2) endif - if(present(sortDim)) then + if (present(sortDim)) then d = sortDim else d = 1 @@ -467,7 +467,7 @@ pure function math_inv33(A) logical :: error call math_invert33(math_inv33,DetA,error,A) - if(error) math_inv33 = 0.0_pReal + if (error) math_inv33 = 0.0_pReal end function math_inv33 @@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -802,7 +802,7 @@ pure function math_sym3333to66(m3333,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -822,7 +822,7 @@ end function math_sym3333to66 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix +!> @brief convert 6x6 matrix into symmetric 3x3x3x3 matrix !> @details Weighted conversion (default) rearranges according to Nye and weights shear ! components according to Mandel. Advisable for matrix operations. ! Unweighted conversion only rearranges order according to Nye @@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -854,7 +854,7 @@ end function math_66toSym3333 !-------------------------------------------------------------------------------------------------- -!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. +!> @brief convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix !-------------------------------------------------------------------------------------------------- pure function math_Voigt66to3333(m66) @@ -875,7 +875,7 @@ end function math_Voigt66to3333 !-------------------------------------------------------------------------------------------------- -!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix. +!> @brief convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix !-------------------------------------------------------------------------------------------------- pure function math_3333toVoigt66(m3333) @@ -984,7 +984,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(1), & (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)] norm = norm2(v(1:3, 1)) - fallback1: if(norm < threshold) then + fallback1: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback1 v(1:3,1) = v(1:3, 1) / norm @@ -992,7 +992,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(2), & (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)] norm = norm2(v(1:3, 2)) - fallback2: if(norm < threshold) then + fallback2: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback2 v(1:3,2) = v(1:3, 2) / norm @@ -1029,7 +1029,7 @@ pure function math_rotationalPart(F) result(R) I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))] x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal) - if(dNeq0(x)) then + if (dNeq0(x)) then Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal)) lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) & *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) @@ -1091,7 +1091,7 @@ function math_eigvalsh33(m) - 2.0_pReal/27.0_pReal*I(1)**3.0_pReal & - I(3) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) - if(all(abs([P,Q]) < TOL)) then + if (all(abs([P,Q]) < TOL)) then math_eigvalsh33 = math_eigvalsh(m) else rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal @@ -1214,7 +1214,7 @@ real(pReal) pure elemental function math_clip(a, left, right) if (present(left)) math_clip = max(left,math_clip) if (present(right)) math_clip = min(right,math_clip) if (present(left) .and. present(right)) then - if(left>right) error stop 'left > right' + if (left>right) error stop 'left > right' endif end function math_clip @@ -1260,38 +1260,38 @@ subroutine selfTest error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' call math_sort(sort_in_,1,3,2) - if(any(sort_in_ /= sort_out_)) & + if (any(sort_in_ /= sort_out_)) & error stop 'math_sort' - if(any(math_range(5) /= range_out_)) & + if (any(math_range(5) /= range_out_)) & error stop 'math_range' - if(any(dNeq(math_exp33(math_I3,0),math_I3))) & + if (any(dNeq(math_exp33(math_I3,0),math_I3))) & error stop 'math_exp33(math_I3,1)' - if(any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & + if (any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & error stop 'math_exp33(math_I3,128)' call random_number(v9) - if(any(dNeq(math_33to9(math_9to33(v9)),v9))) & + if (any(dNeq(math_33to9(math_9to33(v9)),v9))) & error stop 'math_33to9/math_9to33' call random_number(t99) - if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & + if (any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & error stop 'math_3333to99/math_99to3333' call random_number(v6) - if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & + if (any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & error stop 'math_sym33to6/math_6toSym33' call random_number(t66) - if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' - if(any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_3333toVoigt66/math_Voigt66to3333' call random_number(v6) - if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & + if (any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & error stop 'math_symmetric33' call random_number(v3_1) @@ -1299,34 +1299,34 @@ subroutine selfTest call random_number(v3_3) call random_number(v3_4) - if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & + if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & error stop 'math_volTetrahedron' call random_number(t33) - if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & + if (dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & error stop 'math_det33/math_detSym33' - if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & + if (any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & error stop 'math_mul3333xx33/math_identity4th' - if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & + if (any(dNeq0(math_eye(3),math_inv33(math_I3)))) & error stop 'math_inv33(math_I3)' do while(abs(math_det33(t33))<1.0e-9_pReal) call random_number(t33) enddo - if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & + if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & error stop 'math_inv33' call math_invert33(t33_2,det,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert33: T:T^-1 != I' - if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & + if (dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & error stop 'math_invert33 (determinant)' call math_invert(t33_2,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert t33' do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1 @@ -1334,7 +1334,7 @@ subroutine selfTest enddo t33_2 = math_rotationalPart(transpose(t33)) t33 = math_rotationalPart(t33) - if(any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & + if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & error stop 'math_rotationalPart' call random_number(r) @@ -1342,33 +1342,33 @@ subroutine selfTest txx = math_eye(d) allocate(txx_2(d,d)) call math_invert(txx_2,e,txx) - if(any(dNeq0(txx_2,txx) .or. e)) & + if (any(dNeq0(txx_2,txx) .or. e)) & error stop 'math_invert(txx)/math_eye' call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix - if(any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert(t99)' - if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & + if (any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & error stop 'math_clip' - if(math_factorial(10) /= 3628800) & + if (math_factorial(10) /= 3628800) & error stop 'math_factorial' - if(math_binomial(49,6) /= 13983816) & + if (math_binomial(49,6) /= 13983816) & error stop 'math_binomial' - if(math_multinomial([1,2,3,4]) /= 12600) & + if (math_multinomial([1,2,3,4]) /= 12600) & error stop 'math_multinomial' ijk = cshift([1,2,3],int(r*1.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & error stop 'math_LeviCivita(even)' ijk = cshift([3,2,1],int(r*2.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & error stop 'math_LeviCivita(odd)' ijk = cshift([2,2,1],int(r*2.0e2_pReal)) - if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & + if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & error stop 'math_LeviCivita' end subroutine selfTest diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 5792ea3cd..1b39f51ee 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -15,7 +15,7 @@ submodule(phase:mechanical) elastic contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialize elasticity. +!> @brief initialize elasticity !-------------------------------------------------------------------------------------------------- module subroutine elastic_init(phases) @@ -62,7 +62,7 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- -!> @brief Return 6x6 elasticity tensor. +!> @brief return 6x6 elasticity tensor !-------------------------------------------------------------------------------------------------- module function elastic_C66(ph,en) result(C66) @@ -93,7 +93,7 @@ end function elastic_C66 !-------------------------------------------------------------------------------------------------- -!> @brief Return shear modulus. +!> @brief return shear modulus !-------------------------------------------------------------------------------------------------- module function elastic_mu(ph,en) result(mu) @@ -110,7 +110,7 @@ end function elastic_mu !-------------------------------------------------------------------------------------------------- -!> @brief Return Poisson ratio. +!> @brief return Poisson ratio !-------------------------------------------------------------------------------------------------- module function elastic_nu(ph,en) result(nu) @@ -128,7 +128,7 @@ end function elastic_nu !-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> @brief return the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law ! ToDo: Use Voigt matrix directly !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index de027ae44..10522737c 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) @@ -173,7 +173,7 @@ module function plastic_dislotwin_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) @@ -368,7 +368,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! parameters required for several mechanisms and their interactions - if(prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & + if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & prm%D = pl%get_asFloat('D') if (prm%sum_N_tw + prm%sum_N_tr > 0) & @@ -425,7 +425,7 @@ module function plastic_dislotwin_init() result(myPlasticity) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw @@ -473,9 +473,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) integer, intent(in) :: & ph, en real(pReal), dimension(6,6) :: & - homogenizedC - - real(pReal), dimension(6,6) :: & + homogenizedC, & C real(pReal), dimension(:,:,:), allocatable :: & C66_tw, & @@ -500,8 +498,8 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) homogenizedC = homogenizedC & + stt%f_tw(i,en)*C66_tw(1:6,1:6,i) end do - end if twinActive - + end if twinActive + transActive: if (prm%sum_N_tr > 0) then C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) do i=1,prm%sum_N_tr @@ -597,7 +595,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - shearBandingContribution: if(dNeq0(prm%v_sb)) then + shearBandingContribution: if (dNeq0(prm%v_sb)) then E_kB_T = prm%E_sb/(kB*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? @@ -848,7 +846,7 @@ module subroutine plastic_dislotwin_results(ph,group) 'threshold stress for twinning','Pa',prm%systems_tw) case('f_tr') - if(prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & + if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & 'martensite volume fraction','m³/m³') end select @@ -931,8 +929,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & end associate - if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau - if(present(tau_sl)) tau_sl = tau + if (present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau + if (present(tau_sl)) tau_sl = tau end subroutine kinetics_sl @@ -1002,7 +1000,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau end subroutine kinetics_tw @@ -1071,7 +1069,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau end subroutine kinetics_tr