From 9570fb894ae8a371b8d7c55d3ad106fc918181d2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 16 Sep 2018 22:16:06 +0200 Subject: [PATCH] correct names and no superflous conversions anymore --- src/constitutive.f90 | 32 ++++++++----- src/plastic_phenopowerlaw.f90 | 88 +++++++++++++++-------------------- 2 files changed, 57 insertions(+), 63 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index dba1463a7..0c40b0dd7 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -440,7 +440,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e math_Mandel33to6, & math_Plain99to3333 use material, only: & + phasememberAt, & phase_plasticity, & + phase_plasticityInstance, & material_phase, & material_homog, & temperature, & @@ -490,7 +492,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e ho, & !< homogenization tme !< thermal member position integer(pInt) :: & - i, j + i, j, instance, of ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -509,8 +511,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) - dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) @@ -825,6 +828,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac mesh_NcpElems, & mesh_maxNips use material, only: & + phasememberAt, & + phase_plasticityInstance, & phase_plasticity, & phase_source, & phase_Nsources, & @@ -882,38 +887,41 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac real(pReal), intent(in), dimension(6) :: & S6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & - Mstar + Mp integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position - s !< counter in source loop + s, & !< counter in source loop + instance, of ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) + Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el) + call plastic_isotropic_dotState (math_Mandel33to6(Mp),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el) + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_phenopowerlaw_dotState(Mp,ipc,ip,el) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el) + call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & + call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & + call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), & subdt,subfracArray,ip,el) end select plasticityType diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 5e3c8a2ee..da3acd439 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -502,10 +502,7 @@ end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar99,Mstar_v,ipc,ip,el) - use math, only: & - math_Mandel6to33, & - math_Plain3333to99 +subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ipc,ip,el) use material, only: & phasememberAt, & material_phase, & @@ -514,23 +511,19 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar99,Mstar_v,ipc,ip, implicit none real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & - dLp_dMstar99 !< derivative of Lp with respect to the Mandel stress + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), dimension(6), intent(in) :: & - Mstar_v !< Mandel stress + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt) :: & i,k,l,m,n, & of - real(pReal), dimension(3,3) :: & - S !< Second-Piola Kirchhoff stress - real(pReal), dimension(3,3,3,3) :: & - dLp_dS !< derivative of Lp with respect to Mstar as 4th order tensor real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNslip) :: & dgdot_dtauslip_pos,dgdot_dtauslip_neg, & gdot_slip_pos,gdot_slip_neg @@ -542,37 +535,35 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar99,Mstar_v,ipc,ip, ! BEGIN DEPRECATED of = phasememberAt(ipc,ip,el) - S = math_Mandel6to33(Mstar_v) associate(prm => param(phase_plasticityInstance(material_phase(ipc,ip,el))),& stt => state(phase_plasticityInstance(material_phase(ipc,ip,el)))) ! END DEPRECATED Lp = 0.0_pReal - dLp_dS = 0.0_pReal + dLp_dMp = 0.0_pReal - call kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1_pInt, prm%totalNslip Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dS(k,l,m,n) = dLp_dS(k,l,m,n) & - + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) & - *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_pos(m,n,:,i))) + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) & + *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_pos(m,n,:,i))) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dS(k,l,m,n) = dLp_dS(k,l,m,n) & - + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) & - *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_neg(m,n,:,i))) + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) & + *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_neg(m,n,:,i))) enddo slipSystems - call kinetics_twin(prm,stt,of,S,gdot_twin,dgdot_dtautwin) + call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1_pInt, prm%totalNtwin Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dS(k,l,m,n) = dLp_dS(k,l,m,n) & - + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) enddo twinSystems end associate - dLp_dMstar99 = math_Plain3333to99(dLp_dS) ! DEPRECATED end subroutine plastic_phenopowerlaw_LpAndItsTangent @@ -580,17 +571,15 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_dotState(Mstar6,ipc,ip,el) - use math, only: & - math_Mandel6to33 +subroutine plastic_phenopowerlaw_dotState(Mp,ipc,ip,el) use material, only: & material_phase, & phasememberAt, & phase_plasticityInstance implicit none - real(pReal), dimension(6), intent(in) :: & - Mstar6 !< Mandel stress + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point @@ -604,8 +593,6 @@ subroutine plastic_phenopowerlaw_dotState(Mstar6,ipc,ip,el) c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset - real(pReal), dimension(3,3) :: & - S !< Second-Piola Kirchhoff stress real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNslip) :: & left_SlipSlip,right_SlipSlip, & gdot_slip_pos,gdot_slip_neg @@ -619,7 +606,6 @@ subroutine plastic_phenopowerlaw_dotState(Mstar6,ipc,ip,el) dot => dotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) dot%whole(:,of) = 0.0_pReal - S = math_Mandel6to33(Mstar6) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices @@ -636,10 +622,10 @@ subroutine plastic_phenopowerlaw_dotState(Mstar6,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg) + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) dot%sumGamma(of) = sum(dot%gamma_slip(:,of)) - call kinetics_twin(prm,stt,of,S,dot%gamma_twin(:,of)) + call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) if (stt%sumF(of) < 0.98_pReal) dot%sumF(of) = sum(dot%gamma_twin(:,of)/prm%gamma_twin_char) !-------------------------------------------------------------------------------------------------- @@ -671,7 +657,7 @@ end subroutine plastic_phenopowerlaw_dotState !> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the !> result (i.e. intent(out)) variables are the last to have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -subroutine kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg, & +subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) use prec, only: & dNeq0 @@ -692,7 +678,7 @@ subroutine kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg, & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg real(pReal), dimension(3,3), intent(in) :: & - S + Mp real(pReal), dimension(prm%totalNslip) :: & tau_slip_pos, & @@ -701,11 +687,11 @@ subroutine kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg, & integer(pInt) :: i, j do i = 1_pInt, prm%totalNslip - tau_slip_pos(i) = math_mul33xx33(S,prm%Schmid_slip(1:3,1:3,i)) + tau_slip_pos(i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) tau_slip_neg(i) = tau_slip_pos(i) do j = 1,size(prm%nonSchmidCoeff) - tau_slip_pos(i) = tau_slip_pos(i) + math_mul33xx33(S,prm%nonSchmid_pos(1:3,1:3,j,i)) - tau_slip_neg(i) = tau_slip_neg(i) + math_mul33xx33(S,prm%nonSchmid_neg(1:3,1:3,j,i)) + tau_slip_pos(i) = tau_slip_pos(i) + math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j,i)) + tau_slip_neg(i) = tau_slip_neg(i) + math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j,i)) enddo enddo @@ -737,7 +723,7 @@ end subroutine kinetics_slip !> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the !> result (i.e. intent(out)) variables are the last to have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -subroutine kinetics_twin(prm,stt,of,S,gdot_twin,dgdot_dtau_twin) +subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) use prec, only: & dNeq0 use math, only: & @@ -751,7 +737,7 @@ subroutine kinetics_twin(prm,stt,of,S,gdot_twin,dgdot_dtau_twin) integer(pInt), intent(in) :: & of real(pReal), dimension(3,3), intent(in) :: & - S + Mp real(pReal), dimension(prm%totalNtwin), intent(out) :: & gdot_twin real(pReal), dimension(prm%totalNtwin), optional, intent(out) :: & @@ -762,7 +748,7 @@ subroutine kinetics_twin(prm,stt,of,S,gdot_twin,dgdot_dtau_twin) integer(pInt) :: i do i = 1_pInt, prm%totalNtwin - tau_twin(i) = math_mul33xx33(S,prm%Schmid_twin(1:3,1:3,i)) + tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo gdot_twin = merge((1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin, & 0.0_pReal, tau_twin>0.0_pReal) @@ -781,7 +767,7 @@ end subroutine kinetics_twin !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el) result(postResults) +function plastic_phenopowerlaw_postResults(Mp6,ipc,ip,el) result(postResults) use material, only: & material_phase, & plasticState, & @@ -793,14 +779,14 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el) result(postResults) implicit none real(pReal), dimension(6), intent(in) :: & - Mstar6 !< Mandel stress + Mp6 !< Mandel stress integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), dimension(3,3) :: & - S !< Second-Piola Kirchhoff stress + Mp !< Second-Piola Kirchhoff stress real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: & postResults @@ -821,7 +807,7 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el) result(postResults) postResults = 0.0_pReal c = 0_pInt - S = math_Mandel6to33(Mstar6) !DEPRECATED + Mp = math_Mandel6to33(Mp6) !DEPRECATED outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) @@ -833,12 +819,12 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el) result(postResults) postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of) c = c + prm%totalNslip case (shearrate_slip_ID) - call kinetics_slip(prm,stt,of,S,gdot_slip_pos,gdot_slip_neg) + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg c = c + prm%totalNslip case (resolvedstress_slip_ID) do i = 1_pInt, prm%totalNslip - tau_slip_pos = math_mul33xx33(S,prm%Schmid_slip(1:3,1:3,i)) + tau_slip_pos = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) tau_slip_neg = tau_slip_pos !do j = 1,size(prm%nonSchmidCoeff) ! tau_slip_pos = tau_slip_pos + math_mul33xx33(S,prm%nonSchmid_pos(1:3,1:3,j,i)) @@ -855,11 +841,11 @@ function plastic_phenopowerlaw_postResults(Mstar6,ipc,ip,el) result(postResults) postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of) c = c + prm%totalNtwin case (shearrate_twin_ID) - call kinetics_twin(prm,stt,of,S,postResults(c+1_pInt:c+prm%totalNtwin)) + call kinetics_twin(prm,stt,of,Mp,postResults(c+1_pInt:c+prm%totalNtwin)) c = c + prm%totalNtwin case (resolvedstress_twin_ID) do i = 1_pInt, prm%totalNtwin - postResults(c+i) = math_mul33xx33(S,prm%Schmid_twin(1:3,1:3,i)) + postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) enddo c = c + prm%totalNtwin