From 1df409376c884e1138dc3b2eeb0b7f9eee76481b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Jan 2021 23:32:54 +0100 Subject: [PATCH] sourceState is now damage state --- src/constitutive.f90 | 283 ++---------------------------------- src/constitutive_damage.f90 | 193 ++++++++++++++++++++++++ src/constitutive_mech.f90 | 82 ++++++++++- 3 files changed, 286 insertions(+), 272 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 6a023022f..43d9b6b3f 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -235,14 +235,22 @@ module constitutive ! == cleaned:end =================================================================================== module function integrateThermalState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt integer, intent(in) :: & el, & !< element index in element loop ip, & !< integration point index in ip loop co !< grain index in grain loop logical :: broken - end function + end function integrateThermalState + + module function integrateDamageState(dt,co,ip,el) result(broken) + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function integrateDamageState module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt @@ -395,7 +403,6 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LiAndItsTangents, & constitutive_damage_getRateAndItsTangents, & constitutive_thermal_getRateAndItsTangents, & constitutive_results, & @@ -413,7 +420,8 @@ module constitutive crystallite_push33ToRef, & constitutive_restartWrite, & constitutive_restartRead, & - integrateSourceState, & + integrateThermalState, & + integrateDamageState, & constitutive_thermal_setT, & constitutive_mech_getP, & constitutive_mech_setF, & @@ -555,173 +563,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki end function kinematics_active -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: MD: S is Mi? -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Li !< intermediate velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dS, & !< derivative of Li with respect to S - dLi_dFi - - real(pReal), dimension(3,3) :: & - my_Li, & !< intermediate velocity gradient - FiInv, & - temp_33 - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal) :: & - detFi - integer :: & - k, i, j, & - instance, of - - Li = 0.0_pReal - dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticityType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select plasticityType - - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - enddo KinematicsLoop - - FiInv = math_inv33(Fi) - detFi = math_det33(Fi) - Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration - temp_33 = matmul(FiInv,Li) - - do i = 1,3; do j = 1,3 - dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - -end subroutine constitutive_LiAndItsTangents - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - integer :: & - so !< counter in source loop - logical :: broken - - - broken = .false. - - SourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& - co, ip, el) ! correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState(co, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState(co, ip, el) - - end select sourceType - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) - - enddo SourceLoop - -end function constitutive_damage_collectDotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief for constitutive models having an instantaneous change of state -!> will return false if delta state is not needed/supported by the constitutive model -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - real(pReal), intent(in), dimension(3,3) :: & - Fe !< elastic deformation gradient - integer :: & - so, & - myOffset, & - mySize - logical :: & - broken - - - broken = .false. - - sourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & - co, ip, el) - broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) - if(.not. broken) then - myOffset = sourceState(ph)%p(so)%offsetDeltaState - mySize = sourceState(ph)%p(so)%sizeDeltaState - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) - endif - - end select sourceType - - enddo SourceLoop - -end function constitutive_damage_deltaState - - !-------------------------------------------------------------------------------------------------- !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- @@ -1030,107 +871,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33) end function crystallite_push33ToRef -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -function integrateSourceState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - - integer :: & - NiterationState, & !< number of iterations in state loop - ph, & - me, & - so - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken, converged_ - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - converged_ = .true. - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) return - - do so = 1, phase_Nsources(ph) - size_so(so) = sourceState(ph)%p(so)%sizeDotState - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) exit iteration - - do so = 1, phase_Nsources(ph) - zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & - - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) then - broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) - exit iteration - endif - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end function integrateSourceState !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 3ce614666..be47d92b6 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -214,6 +214,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module function integrateDamageState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) return + + do so = 1, phase_Nsources(ph) + size_so(so) = sourceState(ph)%p(so)%sizeDotState + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, phase_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) exit iteration + + do so = 1, phase_Nsources(ph) + zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & + - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & + - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt + sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + sourceState(ph)%p(so)%state(1:size_so(so),me), & + sourceState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) then + broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) + exit iteration + endif + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateDamageState + + !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- @@ -250,4 +355,92 @@ module subroutine damage_results(group,ph) end subroutine damage_results +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + integer :: & + so !< counter in source loop + logical :: broken + + + broken = .false. + + SourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& + co, ip, el) ! correct stress? + + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_dotState(co, ip, el) + + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_dotState(co, ip, el) + + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) + + enddo SourceLoop + +end function constitutive_damage_collectDotState + + + +!-------------------------------------------------------------------------------------------------- +!> @brief for constitutive models having an instantaneous change of state +!> will return false if delta state is not needed/supported by the constitutive model +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + real(pReal), intent(in), dimension(3,3) :: & + Fe !< elastic deformation gradient + integer :: & + so, & + myOffset, & + mySize + logical :: & + broken + + + broken = .false. + + sourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & + co, ip, el) + broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) + if(.not. broken) then + myOffset = sourceState(ph)%p(so)%offsetDeltaState + mySize = sourceState(ph)%p(so)%sizeDeltaState + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & + sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) + endif + + end select sourceType + + enddo SourceLoop + +end function constitutive_damage_deltaState + + end submodule constitutive_damage diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 9a065a829..7c403eeea 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -1653,7 +1653,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif @@ -1938,5 +1938,85 @@ module subroutine constitutive_mech_setF(F,co,ip,el) end subroutine constitutive_mech_setF + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient + FiInv, & + temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS + real(pReal) :: & + detFi + integer :: & + k, i, j, & + instance, of + + Li = 0.0_pReal + dLi_dS = 0.0_pReal + dLi_dFi = 0.0_pReal + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticityType + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + case default plasticityType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select plasticityType + + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + case (KINEMATICS_cleavage_opening_ID) kinematicsType + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_slipplane_opening_ID) kinematicsType + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_thermal_expansion_ID) kinematicsType + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) + case default kinematicsType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select kinematicsType + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + enddo KinematicsLoop + + FiInv = math_inv33(Fi) + detFi = math_det33(Fi) + Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration + temp_33 = matmul(FiInv,Li) + + do i = 1,3; do j = 1,3 + dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + enddo; enddo + +end subroutine constitutive_LiAndItsTangents + end submodule constitutive_mech