diff --git a/src/phase.f90 b/src/phase.f90 index 628bfc443..82b82d71d 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -181,6 +181,11 @@ module phase real(pReal) :: dot_T end function thermal_dot_T + module function damage_phi(ph,me) result(phi) + integer, intent(in) :: ph,me + real(pReal) :: phi + end function damage_phi + module subroutine phase_mechanical_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 6063f3583..8defcf3a1 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -528,4 +528,15 @@ module function phase_damage_get_phi(co,ip,el) result(phi) end function phase_damage_get_phi +module function damage_phi(ph,me) result(phi) + + integer, intent(in) :: ph, me + real(pReal) :: phi + + + phi = current(ph)%phi(me) + +end function damage_phi + + end submodule damagee diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 61151d33b..f7631a38d 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -98,12 +98,9 @@ submodule(phase) mechanics end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -583,7 +580,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) enddo LpLoop call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, co, ip, el) + S, Fi_new, ph,me) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -1303,7 +1300,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & phase_mechanical_S(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), & - co,ip,el) + ph,me) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me)) diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 3fb2d1bf2..9eb3c5767 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -16,11 +16,8 @@ submodule(phase:mechanics) eigendeformation logical, dimension(:,:), allocatable :: myKinematics end function kinematics_thermal_expansion_init - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -29,11 +26,8 @@ submodule(phase:mechanics) eigendeformation dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) end subroutine kinematics_cleavage_opening_LiAndItsTangent - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) - integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -44,7 +38,6 @@ submodule(phase:mechanics) eigendeformation module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me) integer, intent(in) :: ph, me - !< element number real(pReal), intent(out), dimension(3,3) :: & Li !< thermal velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & @@ -165,12 +158,10 @@ end function kinematics_active2 ! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) + S, Fi, ph,me) integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element + ph,me real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -190,18 +181,16 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & real(pReal) :: & detFi integer :: & - k, i, j, & - instance, of, me, ph + k, i, j Li = 0.0_pReal dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + + plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me) case default plasticType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal @@ -210,17 +199,13 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & 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))) + KinematicsLoop: do k = 1, phase_Nkinematics(ph) + kinematicsType: select case (phase_kinematics(k,ph)) case (KINEMATICS_cleavage_opening_ID) kinematicsType - print*, 'clea' - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) case (KINEMATICS_slipplane_opening_ID) kinematicsType - print*, 'slipp' - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) case (KINEMATICS_thermal_expansion_ID) kinematicsType - me = material_phaseMemberAt(co,ip,el) - ph = material_phaseAt(co,el) call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) case default kinematicsType my_Li = 0.0_pReal diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index 8878c10bd..5bb9847f4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -95,12 +95,10 @@ end function kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + ph,me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -117,9 +115,9 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, Ld = 0.0_pReal dLd_dTstar = 0.0_pReal - associate(prm => param(material_phaseAt(co,el))) + associate(prm => param(ph)) do i = 1,prm%sum_N_cl - traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal + traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) if (abs(traction_d) > traction_crit + tol_math_check) then diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index e8410a935..64988d87f 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -115,12 +115,10 @@ end function kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & - co, & !< grain number - ip, & !< integration point number - el !< element number + ph, me real(pReal), intent(in), dimension(3,3) :: & S real(pReal), intent(out), dimension(3,3) :: & @@ -135,7 +133,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt - associate(prm => param(material_phaseAt(co,el))) + associate(prm => param(ph)) Ld = 0.0_pReal dLd_dTstar = 0.0_pReal do i = 1, prm%sum_N_sl @@ -144,7 +142,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i)) traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i)) - traction_crit = prm%g_crit(i)* phase_damage_get_phi(co,ip,el) ! degrading critical load carrying capacity by damage + traction_crit = prm%g_crit(i)* damage_phi(ph,me) udotd = sign(1.0_pReal,traction_d)* prm%dot_o* ( abs(traction_d)/traction_crit & - abs(traction_d)/prm%g_crit(i))**prm%q