diff --git a/src/phase.f90 b/src/phase.f90 index 54116f3c9..b65bf334f 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -135,18 +135,18 @@ module phase end subroutine mechanical_restartRead - module function mechanical_S(ph,me) result(S) - integer, intent(in) :: ph,me + module function mechanical_S(ph,en) result(S) + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: S end function mechanical_S - module function mechanical_L_p(ph,me) result(L_p) - integer, intent(in) :: ph,me + module function mechanical_L_p(ph,en) result(L_p) + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: L_p end function mechanical_L_p - module function mechanical_F_e(ph,me) result(F_e) - integer, intent(in) :: ph,me + module function mechanical_F_e(ph,en) result(F_e) + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: F_e end function mechanical_F_e @@ -239,8 +239,8 @@ module phase logical :: converged_ end function crystallite_stress - module function phase_homogenizedC(ph,me) result(C) - integer, intent(in) :: ph, me + module function phase_homogenizedC(ph,en) result(C) + integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C end function phase_homogenizedC diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 9f2bc4df2..6bcfc9a75 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -60,7 +60,7 @@ submodule(phase) mechanical module subroutine plastic_init end subroutine plastic_init - module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -69,34 +69,34 @@ submodule(phase) mechanical Mi !< Mandel stress integer, intent(in) :: & ph, & - me + en end subroutine plastic_isotropic_LiAndItsTangent - module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) + module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el, & !< element ph, & - me + en real(pReal), intent(in) :: & subdt !< timestep logical :: broken end function plastic_dotState - module function plastic_deltaState(ph, me) result(broken) + module function plastic_deltaState(ph, en) result(broken) integer, intent(in) :: & ph, & - me + en logical :: & broken end function plastic_deltaState module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, ph,me) + S, Fi, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -111,9 +111,9 @@ submodule(phase) mechanical module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ph,me) + S, Fi, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -155,9 +155,9 @@ submodule(phase) mechanical character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results - module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) + module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) real(pReal), dimension(6,6) :: homogenizedC - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en end function plastic_dislotwin_homogenizedC @@ -189,7 +189,7 @@ module subroutine mechanical_init(materials,phases) co, & ce, & ph, & - me, & + en, & stiffDegradationCtr, & Nmembers class(tNode), pointer :: & @@ -286,21 +286,21 @@ module subroutine mechanical_init(materials,phases) constituent => constituents%get(co) ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) - call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) + call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) & - / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) - phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3 - phase_mechanical_F0(ph)%data(1:3,1:3,me) = math_I3 + phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & + / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) + phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3 + phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 - phase_mechanical_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,me), & - phase_mechanical_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration - phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) - phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) - phase_mechanical_F(ph)%data(1:3,1:3,me) = phase_mechanical_F0(ph)%data(1:3,1:3,me) + phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), & + phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) + phase_mechanical_F(ph)%data(1:3,1:3,en) = phase_mechanical_F0(ph)%data(1:3,1:3,en) enddo enddo; enddo @@ -353,11 +353,11 @@ end subroutine mechanical_init !> the elastic and intermediate deformation gradients using Hooke's law !-------------------------------------------------------------------------------------------------- subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, ph, me) + Fe, Fi, ph, en) integer, intent(in) :: & ph, & - me + en real(pReal), intent(in), dimension(3,3) :: & Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient @@ -373,12 +373,12 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & d, & !< counter in degradation loop i, j - C = math_66toSym3333(phase_homogenizedC(ph,me)) + C = math_66toSym3333(phase_homogenizedC(ph,en)) DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) degradationType: select case(phase_stiffnessDegradation(d,ph)) case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage_phi(ph,me)**2 + C = C * damage_phi(ph,en)**2 end select degradationType enddo DegradationLoop @@ -486,7 +486,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) o, & p, & ph, & - me, & + en, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update logical :: error,broken @@ -495,12 +495,12 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) broken = .true. ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) call plastic_dependentState(co,ip,el) - Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,me) ! take as first guess - Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,me) ! take as first guess + Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess + Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess call math_invert33(invFp_current,devNull,error,subFp0) if (error) return ! error @@ -535,10 +535,10 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) B = math_I3 - Delta_t*Lpguess Fe = matmul(matmul(A,B), invFi_new) call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, ph, me) + Fe, Fi_new, ph, en) call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, ph,me) + S, Fi_new, ph,en) !* update current residuum and check for convergence of loop atol_Lp = max(num%rtol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -579,7 +579,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, ph,me) + S, Fi_new, ph,en) !* 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 @@ -629,13 +629,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - phase_mechanical_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) - phase_mechanical_S(ph)%data(1:3,1:3,me) = S - phase_mechanical_Lp(ph)%data(1:3,1:3,me) = Lpguess - phase_mechanical_Li(ph)%data(1:3,1:3,me) = Liguess - phase_mechanical_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - phase_mechanical_Fi(ph)%data(1:3,1:3,me) = Fi_new - phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) + phase_mechanical_P(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + phase_mechanical_S(ph)%data(1:3,1:3,en) = S + phase_mechanical_Lp(ph)%data(1:3,1:3,en) = Lpguess + phase_mechanical_Li(ph)%data(1:3,1:3,en) = Liguess + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = Fi_new + phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(matmul(F,invFp_new),invFi_new) broken = .false. end function integrateStress @@ -660,7 +660,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul integer :: & NiterationState, & !< number of iterations in state loop ph, & - me, & + en, & sizeDotState real(pReal) :: & zeta @@ -671,38 +671,37 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) - broken = plastic_dotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,me) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - dotState(1:sizeDotState,2) = 0.0_pReal + plasticState(ph)%state(1:sizeDotState,en) = subState0 & + + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t iteration: do NiterationState = 1, num%nState - if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1) - dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me) + dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1) + dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration - broken = plastic_dotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) exit iteration - zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),& + zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),& dotState(1:sizeDotState,2)) - plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta & + plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta & + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) - r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) & + r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) & - subState0 & - - plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t - plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) & + - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) & - r(1:sizeDotState) - if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then - broken = plastic_deltaState(ph,me) + if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then + broken = plastic_deltaState(ph,en) exit iteration endif @@ -714,20 +713,19 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul !-------------------------------------------------------------------------------------------------- !> @brief calculate the damping for correction of state and dot state !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) + real(pReal) pure function damper(omega_0,omega_1,omega_2) - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 + real(pReal), dimension(:), intent(in) :: & + omega_0, omega_1, omega_2 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 + dot_prod12 = dot_product(omega_0 - omega_1, omega_1 - omega_2) + dot_prod22 = dot_product(omega_1 - omega_2, omega_1 - omega_2) + + damper = merge(0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22), & + 1.0_pReal, & + (min(dot_product(omega_0,omega_1), dot_prod12) < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) end function damper @@ -751,21 +749,21 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res integer :: & ph, & - me, & + en, & sizeDotState ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) - broken = plastic_dotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,me) = subState0 & - + plasticState(ph)%dotState(1:sizeDotState,me) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = subState0 & + + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t - broken = plastic_deltaState(ph,me) + broken = plastic_deltaState(ph,en) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -790,34 +788,34 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip integer :: & ph, & - me, & + en, & sizeDotState real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) - broken = plastic_dotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return sizeDotState = plasticState(ph)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * Delta_t - plasticState(ph)%state(1:sizeDotState,me) = subState0 & - + plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t + residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = subState0 & + + plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t - broken = plastic_deltaState(ph,me) + broken = plastic_deltaState(ph,en) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) return - broken = plastic_dotState(Delta_t, co,ip,el,ph,me) + broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return - broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, & - plasticState(ph)%state(1:sizeDotState,me), & + broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, & + plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState)) end function integrateStateAdaptiveEuler @@ -908,55 +906,55 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D stage, & ! stage index in integration stage loop n, & ph, & - me, & + en, & sizeDotState real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) - broken = plastic_dotState(Delta_t,co,ip,el,ph,me) + broken = plastic_dotState(Delta_t,co,ip,el,ph,en) if(broken) return sizeDotState = plasticState(ph)%sizeDotState do stage = 1, size(A,1) - plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me) - plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en) + plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do n = 2, stage - plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) & + plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) & + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) enddo - plasticState(ph)%state(1:sizeDotState,me) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = subState0 & + + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) if(broken) exit - broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,me) + broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en) if(broken) exit enddo if(broken) return - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me) - plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(ph)%state(1:sizeDotState,me) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en) + plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(ph)%state(1:sizeDotState,en) = subState0 & + + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t if(present(DB)) & broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & - plasticState(ph)%state(1:sizeDotState,me), & + plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState)) if(broken) return - broken = plastic_deltaState(ph,me) + broken = plastic_deltaState(ph,en) if(broken) return broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) @@ -1087,14 +1085,14 @@ end subroutine mechanical_forward !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent !-------------------------------------------------------------------------------------------------- -module function phase_homogenizedC(ph,me) result(C) +module function phase_homogenizedC(ph,en) result(C) real(pReal), dimension(6,6) :: C - integer, intent(in) :: ph, me + integer, intent(in) :: ph, en plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_DISLOTWIN_ID) plasticType - C = plastic_dislotwin_homogenizedC(ph,me) + C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType C = lattice_C66(1:6,1:6,ph) end select plasticType @@ -1117,7 +1115,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - ph, me, sizeDotState + ph, en, sizeDotState logical :: todo real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & @@ -1131,19 +1129,19 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) + en = material_phaseMemberAt(co,ip,el) sizeDotState = plasticState(ph)%sizeDotState - subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,me) - subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) - subState0 = plasticState(ph)%State0(:,me) + subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) + subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) + subState0 = plasticState(ph)%State0(:,en) if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) + damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en) - subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) - subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) - subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) + subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) + subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) + subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. @@ -1161,29 +1159,29 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then subF0 = subF - subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,me) - subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,me) - subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) - subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) - subState0 = plasticState(ph)%state(:,me) + subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,en) + subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,en) + subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en) + subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en) + subState0 = plasticState(ph)%state(:,en) if (damageState(ph)%sizeState > 0) & - damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) + damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en) endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) else subStep = num%subStepSizeCryst * subStep - phase_mechanical_Fp(ph)%data(1:3,1:3,me) = subFp0 - phase_mechanical_Fi(ph)%data(1:3,1:3,me) = subFi0 - phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use? + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0 + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0 + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use? if (subStep < 1.0_pReal) then ! actual (not initial) cutback - phase_mechanical_Lp(ph)%data(1:3,1:3,me) = subLp0 - phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 + phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0 + phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0 endif - plasticState(ph)%state(:,me) = subState0 + plasticState(ph)%state(:,en) = subState0 if (damageState(ph)%sizeState > 0) & - damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) + damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif @@ -1192,7 +1190,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! prepare for integration if (todo) then subF = subF0 & - + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me)) + + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en)) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) endif @@ -1212,22 +1210,22 @@ module subroutine mechanical_restore(ce,includeL) includeL !< protect agains fake cutback integer :: & - co, ph, me + co, ph, en do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) + en = material_phaseEntry(co,ce) if (includeL) then - phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) - phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) + phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) + phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en) endif ! maybe protecting everything from overwriting makes more sense - phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) - phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) - phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) + phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) - plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) + plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en) enddo end subroutine mechanical_restore @@ -1245,7 +1243,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) integer :: & o, & - p, ph, me + p, ph, en real(pReal), dimension(3,3) :: devNull, & invSubFp0,invSubFi0,invFp,invFi, & temp_33_1, temp_33_2, temp_33_3 @@ -1266,20 +1264,20 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) ph = material_phaseID(co,ce) - me = material_phaseEntry(co,ce) + en = material_phaseEntry(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - phase_mechanical_Fe(ph)%data(1:3,1:3,me), & - phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) + phase_mechanical_Fe(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) 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), & - ph,me) + phase_mechanical_S(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en), & + ph,en) - 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)) - invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me)) - invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,me)) + invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) + invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en)) + invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en)) + invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,en)) if (sum(abs(dLidS)) < tol_math_check) then dFidS = 0.0_pReal @@ -1305,15 +1303,15 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) endif call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - phase_mechanical_S(ph)%data(1:3,1:3,me), & - phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me) + phase_mechanical_S(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- ! calculate dSdF temp_33_1 = transpose(matmul(invFp,invFi)) - temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invSubFp0) - temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) + temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invSubFp0) + temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invFp), invSubFi0) do o=1,3; do p=1,3 rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) @@ -1341,9 +1339,9 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) !-------------------------------------------------------------------------------------------------- ! assemble dPdF - temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,me),transpose(invFp)) - temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp) - temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,me)) + temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,en),transpose(invFp)) + temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),invFp) + temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,en)) dPdF = 0.0_pReal do p=1,3 @@ -1351,7 +1349,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) enddo do o=1,3; do p=1,3 dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & - + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) enddo; enddo @@ -1396,13 +1394,13 @@ end subroutine mechanical_restartRead !---------------------------------------------------------------------------------------------- !< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) !---------------------------------------------------------------------------------------------- -module function mechanical_S(ph,me) result(S) +module function mechanical_S(ph,en) result(S) - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: S - S = phase_mechanical_S(ph)%data(1:3,1:3,me) + S = phase_mechanical_S(ph)%data(1:3,1:3,en) end function mechanical_S @@ -1410,13 +1408,13 @@ end function mechanical_S !---------------------------------------------------------------------------------------------- !< @brief Get plastic velocity gradient (for use by non-mech physics) !---------------------------------------------------------------------------------------------- -module function mechanical_L_p(ph,me) result(L_p) +module function mechanical_L_p(ph,en) result(L_p) - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: L_p - L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,me) + L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,en) end function mechanical_L_p @@ -1424,13 +1422,13 @@ end function mechanical_L_p !---------------------------------------------------------------------------------------------- !< @brief Get elastic deformation gradient (for use by non-mech physics) !---------------------------------------------------------------------------------------------- -module function mechanical_F_e(ph,me) result(F_e) +module function mechanical_F_e(ph,en) result(F_e) - integer, intent(in) :: ph,me + integer, intent(in) :: ph,en real(pReal), dimension(3,3) :: F_e - F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,me) + F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,en) end function mechanical_F_e diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 7c4e70c1c..8587b2bc8 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -145,10 +145,10 @@ end function kinematics_active2 ! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, ph,me) + S, Fi, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & @@ -179,7 +179,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_isotropic_ID) plasticType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,me) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. @@ -188,7 +188,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & KinematicsLoop: do k = 1, Nmodels(ph) kinematicsType: select case (model(k,ph)) case (KINEMATICS_thermal_expansion_ID) kinematicsType - call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me) + call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. @@ -197,12 +197,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & select case (model_damage(ph)) case (KINEMATICS_cleavage_opening_ID) - call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. case (KINEMATICS_slipplane_opening_ID) - call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 3e9772a87..3ec275795 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -240,9 +240,9 @@ end subroutine plastic_init ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ph,me) + S, Fi, ph,en) integer, intent(in) :: & - ph,me + ph,en real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola-Kirchhoff stress Fi !< intermediate deformation gradient @@ -250,7 +250,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & Lp !< plastic velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLp_dS, & - dLp_dFi !< derivative me Lp with respect to Fi + dLp_dFi !< derivative en Lp with respect to Fi real(pReal), dimension(3,3,3,3) :: & dLp_dMp !< derivative of Lp with respect to Mandel stress @@ -270,22 +270,22 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & dLp_dMp = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticType - call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTICITY_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTICITY_KINEHARDENING_ID) plasticType - call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) + call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) + call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) + call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) + call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,en),ph,en) end select plasticType @@ -301,14 +301,14 @@ end subroutine plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) +module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point el, & !< element ph, & - me + en real(pReal), intent(in) :: & subdt !< timestep real(pReal), dimension(3,3) :: & @@ -316,30 +316,30 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) logical :: broken - Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),& - phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me)) + Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& + phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_ISOTROPIC_ID) plasticType - call isotropic_dotState(Mp,ph,me) + call isotropic_dotState(Mp,ph,en) case (PLASTICITY_PHENOPOWERLAW_ID) plasticType - call phenopowerlaw_dotState(Mp,ph,me) + call phenopowerlaw_dotState(Mp,ph,en) case (PLASTICITY_KINEHARDENING_ID) plasticType - call plastic_kinehardening_dotState(Mp,ph,me) + call plastic_kinehardening_dotState(Mp,ph,en) case (PLASTICITY_DISLOTWIN_ID) plasticType - call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) + call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType - call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) + call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) case (PLASTICITY_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) + call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) end select plasticType - broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) + broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en))) end function plastic_dotState @@ -383,11 +383,11 @@ end subroutine plastic_dependentState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -module function plastic_deltaState(ph, me) result(broken) +module function plastic_deltaState(ph, en) result(broken) integer, intent(in) :: & ph, & - me + en logical :: & broken @@ -398,18 +398,18 @@ module function plastic_deltaState(ph, me) result(broken) mySize - Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),& - phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me)) + Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& + phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en)) plasticType: select case (phase_plasticity(ph)) case (PLASTICITY_KINEHARDENING_ID) plasticType - call plastic_kinehardening_deltaState(Mp,ph,me) - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) + call plastic_kinehardening_deltaState(Mp,ph,en) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) case (PLASTICITY_NONLOCAL_ID) plasticType - call plastic_nonlocal_deltaState(Mp,ph,me) - broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) + call plastic_nonlocal_deltaState(Mp,ph,en) + broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en))) case default broken = .false. @@ -422,8 +422,8 @@ module function plastic_deltaState(ph, me) result(broken) myOffset = plasticState(ph)%offsetDeltaState mySize = plasticState(ph)%sizeDeltaState - plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) = & - plasticState(ph)%state(myOffset + 1:myOffset + mySize,me) + plasticState(ph)%deltaState(1:mySize,me) + plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = & + plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en) end select endif diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 43be614ac..b04f28f0e 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -483,10 +483,10 @@ end function plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- !> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- -module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) +module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) integer, intent(in) :: & - ph, me + ph, en real(pReal), dimension(6,6) :: & homogenizedC @@ -498,17 +498,17 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) stt => state(ph)) f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,me)) & - - sum(stt%f_tr(1:prm%sum_N_tr,me)) + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) homogenizedC = f_unrotated * prm%C66 do i=1,prm%sum_N_tw homogenizedC = homogenizedC & - + stt%f_tw(i,me)*prm%C66_tw(1:6,1:6,i) + + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) enddo do i=1,prm%sum_N_tr homogenizedC = homogenizedC & - + stt%f_tr(i,me)*prm%C66_tr(1:6,1:6,i) + + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) enddo end associate diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 548e2fd75..d1f191f28 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -208,7 +208,7 @@ end subroutine isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate inelastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- -module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) +module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -219,7 +219,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) Mi !< Mandel stress integer, intent(in) :: & ph, & - me + en real(pReal) :: & tr !< trace of spherical part of Mandel stress (= 3 x pressure) @@ -232,7 +232,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero Li = math_I3 & - * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) & + * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) & * tr * abs(tr)**(prm%n-1.0_pReal) forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n) else