From 3f7a389ff73fb4703aedd12d19741f95f0afc91d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 13 Mar 2014 06:43:49 +0000 Subject: [PATCH] changed state parsing for local models (and for delta_state) such that only the needed part of the state array (for the given material point) is used --- code/constitutive.f90 | 97 +++++++-------- code/constitutive_dislotwin.f90 | 170 +++++++++++++------------- code/constitutive_j2.f90 | 34 +++--- code/constitutive_nonlocal.f90 | 72 +++++------ code/constitutive_phenopowerlaw.f90 | 78 ++++++------ code/constitutive_titanmod.f90 | 181 ++++++++++++++-------------- 6 files changed, 316 insertions(+), 316 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 14ae66128..6d25802dc 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -192,8 +192,8 @@ subroutine constitutive_init select case(phase_plasticity(phase)) ! split per constititution case (PLASTICITY_NONE_ID) outputName = PLASTICITY_NONE_label - thisOutput => NULL() ! constitutive_none_output - thisSize => NULL() ! constitutive_none_sizePostResult + thisOutput => null() ! constitutive_none_output + thisSize => null() ! constitutive_none_sizePostResult case (PLASTICITY_J2_ID) outputName = PLASTICITY_J2_label thisOutput => constitutive_j2_output @@ -506,11 +506,11 @@ pure function constitutive_homogenizedC(ipc,ip,el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_DISLOTWIN_ID) - constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el) - + constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state(ipc,ip,el), & + ipc,ip,el) case (PLASTICITY_TITANMOD_ID) - constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el) - + constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state(ipc,ip,el), & + ipc,ip,el) case default constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase(ipc,ip,el)) @@ -550,11 +550,11 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_DISLOTWIN_ID) - call constitutive_dislotwin_microstructure(temperature,constitutive_state,ipc,ip,el) - + call constitutive_dislotwin_microstructure(temperature,constitutive_state(ipc,ip,el), & + ipc,ip,el) case (PLASTICITY_TITANMOD_ID) - call constitutive_titanmod_microstructure(temperature,constitutive_state,ipc,ip,el) - + call constitutive_titanmod_microstructure(temperature,constitutive_state(ipc,ip,el), & + ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el) @@ -610,20 +610,20 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip dLp_dTstar = math_identity2nd(9) case (PLASTICITY_J2_ID) - call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el) - + call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el) - + call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) - call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el) - + call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & + temperature,constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_TITANMOD_ID) - call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el) - + call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & + temperature,constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) - call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, constitutive_state(ipc,ip,el), ipc,ip,el) - + call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & + temperature, constitutive_state(ipc,ip,el), ipc,ip,el) end select end subroutine constitutive_LpAndItsTangent @@ -700,7 +700,8 @@ end subroutine constitutive_hooke_TandItsTangent !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, subfrac, ipc, ip, el) +subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, subdt, subfracArray,& + ipc, ip, el) use prec, only: & pLongInt use debug, only: & @@ -742,10 +743,10 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub Temperature, & subdt !< timestep real(pReal), intent(in), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - subfrac !< subfraction of timestep + subfracArray !< subfraction of timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Fe, & !< elastic deformation gradient - Fp !< plastic deformation gradient + FeArray, & !< elastic deformation gradient + FpArray !< plastic deformation gradient real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pLongInt) :: & @@ -762,20 +763,21 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub constitutive_dotState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? case (PLASTICITY_J2_ID) - constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,constitutive_state,ipc,ip,el) - + constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,& + constitutive_state(ipc,ip,el), ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_state,ipc,ip,el) - + constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,& + constitutive_state(ipc,ip,el), ipc,ip,el) case (PLASTICITY_TITANMOD_ID) - constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - + constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,& + constitutive_state(ipc,ip,el), ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) - constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - + constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,& + constitutive_state(ipc,ip,el), ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) - constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, constitutive_state, & - constitutive_state0, subdt, subfrac, ipc, ip, el) + constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, & + Temperature, constitutive_state, constitutive_state0, subdt, & + subfracArray, ipc, ip, el) end select @@ -830,7 +832,8 @@ subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el) select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_NONLOCAL_ID) - call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,ipc,ip,el) + call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),& + constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el) case default constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? @@ -853,7 +856,7 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el) +function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -888,7 +891,7 @@ function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el) real(pReal), intent(in) :: & temperature real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - Fe !< elastic deformation gradient + FeArray !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) @@ -899,20 +902,20 @@ function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el) case (PLASTICITY_NONE_ID) case (PLASTICITY_TITANMOD_ID) - constitutive_postResults = constitutive_titanmod_postResults(constitutive_state,ipc,ip,el) - + constitutive_postResults = constitutive_titanmod_postResults(& + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_J2_ID) - constitutive_postResults = constitutive_j2_postResults(Tstar_v,constitutive_state,ipc,ip,el) - + constitutive_postResults = constitutive_j2_postResults(Tstar_v,& + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) - constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_state,ipc,ip,el) - + constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,& + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) - constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - + constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,& + constitutive_state(ipc,ip,el),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) - constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, constitutive_state, & - constitutive_dotstate(ipc,ip,el), ipc, ip, el) + constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, & + constitutive_state, constitutive_dotstate(ipc,ip,el), ipc, ip, el) end select end function constitutive_postResults diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index d82dde906..b5a0b1530 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -939,7 +939,7 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state integer(pInt) :: instance,ns,nt,i,phase @@ -952,13 +952,13 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el) nt = constitutive_dislotwin_totalNtwin(instance) !* Total twin volume fraction - sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 + sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 !* Homogenized elasticity matrix constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*lattice_C66(1:6,1:6,phase) do i=1_pInt,nt constitutive_dislotwin_homogenizedC = & - constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*lattice_C66(1:6,1:6,phase) + constitutive_dislotwin_homogenizedC + state%p(3_pInt*ns+i)*lattice_C66(1:6,1:6,phase) enddo end function constitutive_dislotwin_homogenizedC @@ -992,7 +992,7 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el) el !< element real(pReal), intent(in) :: & temperature !< temperature at IP - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + type(p_vec), intent(inout) :: & state !< microstructure state integer(pInt) :: & @@ -1022,7 +1022,7 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el) !* State: 7*ns+5*nt+1 : 7*ns+6*nt twin volume !* Total twin volume fraction - sumf = sum(state(ipc,ip,el)%p((3*ns+1):(3*ns+nt))) ! safe for nt == 0 + sumf = sum(state%p((3*ns+1):(3*ns+nt))) ! safe for nt == 0 !* Stacking fault energy sfe = constitutive_dislotwin_SFE_0K(instance) + & @@ -1031,61 +1031,61 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el) !* rescaled twin volume fraction for topology forall (t = 1_pInt:nt) & fOverStacksize(t) = & - state(ipc,ip,el)%p(3_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,instance) + state%p(3_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,instance) !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(3_pInt*ns+2_pInt*nt+s) = & - sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),& + state%p(3_pInt*ns+2_pInt*nt+s) = & + sqrt(dot_product((state%p(1:ns)+state%p(ns+1_pInt:2_pInt*ns)),& constitutive_dislotwin_forestProjectionEdge(1:ns,s,instance)))/ & constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,instance) !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation !$OMP CRITICAL (evilmatmul) - state(ipc,ip,el)%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt)) = 0.0_pReal + state%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt)) = 0.0_pReal if (nt > 0_pInt .and. ns > 0_pInt) & - state(ipc,ip,el)%p((4_pInt*ns+2_pInt*nt+1):(5_pInt*ns+2_pInt*nt)) = & + state%p((4_pInt*ns+2_pInt*nt+1):(5_pInt*ns+2_pInt*nt)) = & matmul(constitutive_dislotwin_interactionMatrix_SlipTwin(1:ns,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin !$OMP CRITICAL (evilmatmul) if (nt > 0_pInt) & - state(ipc,ip,el)%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt)) = & + state%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt)) = & matmul(constitutive_dislotwin_interactionMatrix_TwinTwin(1:nt,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* mean free path between 2 obstacles seen by a moving dislocation do s = 1_pInt,ns if (nt > 0_pInt) then - state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+s) = & + state%p(5_pInt*ns+3_pInt*nt+s) = & constitutive_dislotwin_GrainSize(instance)/(1.0_pReal+constitutive_dislotwin_GrainSize(instance)*& - (state(ipc,ip,el)%p(3_pInt*ns+2_pInt*nt+s)+state(ipc,ip,el)%p(4_pInt*ns+2_pInt*nt+s))) + (state%p(3_pInt*ns+2_pInt*nt+s)+state%p(4_pInt*ns+2_pInt*nt+s))) else - state(ipc,ip,el)%p(5_pInt*ns+s) = & + state%p(5_pInt*ns+s) = & constitutive_dislotwin_GrainSize(instance)/& - (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state(ipc,ip,el)%p(3_pInt*ns+s))) + (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*(state%p(3_pInt*ns+s))) endif enddo !* mean free path between 2 obstacles seen by a growing twin forall (t = 1_pInt:nt) & - state(ipc,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = & + state%p(6_pInt*ns+3_pInt*nt+t) = & (constitutive_dislotwin_Cmfptwin(instance)*constitutive_dislotwin_GrainSize(instance))/& - (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state(ipc,ip,el)%p(5_pInt*ns+2_pInt*nt+t)) + (1.0_pReal+constitutive_dislotwin_GrainSize(instance)*state%p(5_pInt*ns+2_pInt*nt+t)) !* threshold stress for dislocation motion if(lattice_structure(phase) /= LATTICE_BCC_ID) then ! bcc value remains constant forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ & + state%p(6_pInt*ns+4_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(instance)+ & lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(s,instance)*& - sqrt(dot_product((state(ipc,ip,el)%p(1:ns)+state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),& + sqrt(dot_product((state%p(1:ns)+state%p(ns+1_pInt:2_pInt*ns)),& constitutive_dislotwin_interactionMatrix_SlipSlip(s,1:ns,instance))) endif !* threshold stress for growing twin forall (t = 1_pInt:nt) & - state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+t) = & + state%p(7_pInt*ns+4_pInt*nt+t) = & constitutive_dislotwin_Cthresholdtwin(instance)*& (sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,instance))+& 3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,instance)*lattice_mu(phase)/& @@ -1093,8 +1093,8 @@ subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el) !* final twin volume after growth forall (t = 1_pInt:nt) & - state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+t) = & - (pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,instance)*state(ipc,ip,el)%p(6*ns+3*nt+t)**(2.0_pReal) + state%p(7_pInt*ns+5_pInt*nt+t) = & + (pi/4.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,instance)*state%p(6*ns+3*nt+t)**(2.0_pReal) !* equilibrium seperation of partial dislocations do t = 1_pInt,nt @@ -1145,12 +1145,12 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat LATTICE_fcc_ID implicit none - integer(pInt), intent(in) :: ipc,ip,el - real(pReal), intent(in) :: Temperature - real(pReal), dimension(6), intent(in) :: Tstar_v - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: state - real(pReal), dimension(3,3), intent(out) :: Lp - real(pReal), dimension(9,9), intent(out) :: dLp_dTstar + integer(pInt), intent(in) :: ipc,ip,el + real(pReal), intent(in) :: Temperature + real(pReal), dimension(6), intent(in) :: Tstar_v + type(p_vec), intent(inout) :: state + real(pReal), dimension(3,3), intent(out) :: Lp + real(pReal), dimension(9,9), intent(out) :: dLp_dTstar integer(pInt) :: instance,phase,ns,nt,f,i,j,k,l,m,n,index_myFamily,s1,s2 real(pReal) :: sumf,StressRatio_p,StressRatio_pminus1,StressRatio_r,BoltzmannRatio,DotGamma0,Ndot0 @@ -1190,7 +1190,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat nt = constitutive_dislotwin_totalNtwin(instance) !* Total twin volume fraction - sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 + sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal @@ -1214,16 +1214,16 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat StressRatio_p = 0.0_pReal StressRatio_pminus1 = 0.0_pReal else - StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))& + StressRatio_p = (abs(tau_slip(j))/state%p(6*ns+4*nt+j))& **constitutive_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6*ns+4*nt+j))& + StressRatio_pminus1 = (abs(tau_slip(j))/state%p(6*ns+4*nt+j))& **(constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) endif !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) !* Initial shear rates DotGamma0 = & - state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*& + state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*& constitutive_dislotwin_v0PerSlipSystem(j,instance) !* Shear rates due to slip @@ -1234,7 +1234,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat !* Derivatives of shear rates dgdot_dtauslip(j) = & ((abs(gdot_slip(j))*BoltzmannRatio*constitutive_dislotwin_pPerSlipFamily(f,instance)& - *constitutive_dislotwin_qPerSlipFamily(f,instance))/state(ipc,ip,el)%p(6*ns+4*nt+j))*& + *constitutive_dislotwin_qPerSlipFamily(f,instance))/state%p(6*ns+4*nt+j))*& StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal) !* Plastic velocity gradient for dislocation glide @@ -1320,15 +1320,15 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat !* Stress ratios if (tau_twin(j) > tol_math_check) then - StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance) + StressRatio_r = (state%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance) !* Shear rates and their derivatives due to twin select case(lattice_structure(phase)) case (LATTICE_fcc_ID) s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) if (tau_twin(j) < constitutive_dislotwin_tau_r(j,instance)) then - Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+& - abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/& + Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+& + abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/& (constitutive_dislotwin_L0(instance)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))*& (1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& (constitutive_dislotwin_tau_r(j,instance)-tau_twin(j)))) @@ -1340,7 +1340,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat end select gdot_twin(j) = & (constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*& - state(ipc,ip,el)%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r) + state%p(7*ns+5*nt+j)*Ndot0*exp(-StressRatio_r) dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_rPerTwinFamily(f,instance))/tau_twin(j))*StressRatio_r endif @@ -1392,15 +1392,15 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e LATTICE_bcc_ID implicit none - real(pReal), dimension(6), intent(in):: & + real(pReal), dimension(6), intent(in):: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & temperature !< temperature at integration point - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_dislotwin_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & constitutive_dislotwin_dotState @@ -1421,7 +1421,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e nt = constitutive_dislotwin_totalNtwin(instance) !* Total twin volume fraction - sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 + sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 constitutive_dislotwin_dotState = 0.0_pReal @@ -1441,16 +1441,16 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e StressRatio_p = 0.0_pReal StressRatio_pminus1 = 0.0_pReal else - StressRatio_p = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_p = (abs(tau_slip(j))/state%p(6_pInt*ns+4_pInt*nt+j))**& constitutive_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = (abs(tau_slip(j))/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_pminus1 = (abs(tau_slip(j))/state%p(6_pInt*ns+4_pInt*nt+j))**& (constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) endif !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) !* Initial shear rates DotGamma0 = & - state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*& + state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)*& constitutive_dislotwin_v0PerSlipSystem(j,instance) !* Shear rates due to slip @@ -1459,7 +1459,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e !* Multiplication DotRhoMultiplication(j) = abs(gdot_slip(j))/& - (constitutive_dislotwin_burgersPerSlipSystem(j,instance)*state(ipc,ip,el)%p(5*ns+3*nt+j)) + (constitutive_dislotwin_burgersPerSlipSystem(j,instance)*state%p(5*ns+3*nt+j)) !* Dipole formation EdgeDipMinDistance = & @@ -1470,22 +1470,22 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e EdgeDipDistance(j) = & (3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/& (16.0_pReal*pi*abs(tau_slip(j))) - if (EdgeDipDistance(j)>state(ipc,ip,el)%p(5*ns+3*nt+j)) EdgeDipDistance(j)=state(ipc,ip,el)%p(5*ns+3*nt+j) + if (EdgeDipDistance(j)>state%p(5*ns+3*nt+j)) EdgeDipDistance(j)=state%p(5*ns+3*nt+j) if (EdgeDipDistance(j) tol_math_check) then - StressRatio_r = (state(ipc,ip,el)%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance) + StressRatio_r = (state%p(7*ns+4*nt+j)/tau_twin(j))**constitutive_dislotwin_rPerTwinFamily(f,instance) !* Shear rates and their derivatives due to twin select case(lattice_structure(phase)) @@ -1535,8 +1535,8 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) if (tau_twin(j) < constitutive_dislotwin_tau_r(j,instance)) then - Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+& - abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/& + Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+& + abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/& (constitutive_dislotwin_L0(instance)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))*& (1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& (constitutive_dislotwin_tau_r(j,instance)-tau_twin(j)))) @@ -1548,7 +1548,7 @@ pure function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,ipc,ip,e end select constitutive_dislotwin_dotState(3_pInt*ns+j) = & (constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*& - state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r) + state%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r) !* Dotstate for accumulated shear due to twin constitutive_dislotwin_dotState(3_pInt*ns+nt+j) = constitutive_dislotwin_dotState(3_pInt*ns+j) * & lattice_sheartwin(index_myfamily+i,phase) @@ -1592,15 +1592,15 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) LATTICE_fcc_ID implicit none - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & temperature !< temperature at integration point - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_dislotwin_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & @@ -1625,7 +1625,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) nt = constitutive_dislotwin_totalNtwin(instance) !* Total twin volume fraction - sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 + sumf = sum(state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 !* Required output c = 0_pInt @@ -1638,10 +1638,10 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) select case(constitutive_dislotwin_outputID(o,instance)) case (edge_density_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(1_pInt:ns) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state%p(1_pInt:ns) c = c + ns case (dipole_density_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state%p(ns+1_pInt:2_pInt*ns) c = c + ns case (shear_rate_slip_ID) j = 0_pInt @@ -1657,16 +1657,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) StressRatio_p = 0.0_pReal StressRatio_pminus1 = 0.0_pReal else - StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_p = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**& constitutive_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_pminus1 = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**& (constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) endif !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) !* Initial shear rates DotGamma0 = & - state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & + state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & constitutive_dislotwin_v0PerSlipSystem(j,instance) !* Shear rates due to slip @@ -1677,11 +1677,11 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) c = c + ns case (accumulated_shear_slip_ID) constitutive_dislotwin_postResults(c+1_pInt:c+ns) = & - state(ipc,ip,el)%p((2_pInt*ns+1_pInt):(3_pInt*ns)) + state%p((2_pInt*ns+1_pInt):(3_pInt*ns)) c = c + ns case (mfp_slip_ID) constitutive_dislotwin_postResults(c+1_pInt:c+ns) =& - state(ipc,ip,el)%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt)) + state%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt)) c = c + ns case (resolved_stress_slip_ID) j = 0_pInt @@ -1695,7 +1695,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) c = c + ns case (threshold_stress_slip_ID) constitutive_dislotwin_postResults(c+1_pInt:c+ns) = & - state(ipc,ip,el)%p((6_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+4_pInt*nt)) + state%p((6_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+4_pInt*nt)) c = c + ns case (edge_dipole_distance_ID) j = 0_pInt @@ -1706,8 +1706,8 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) constitutive_dislotwin_postResults(c+j) = & (3.0_pReal*lattice_mu(phase)*constitutive_dislotwin_burgersPerSlipSystem(j,instance))/& (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)))) - constitutive_dislotwin_postResults(c+j)=min(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(5*ns+3*nt+j)) - ! constitutive_dislotwin_postResults(c+j)=max(constitutive_dislotwin_postResults(c+j),state(ipc,ip,el)%p(4*ns+2*nt+j)) + constitutive_dislotwin_postResults(c+j)=min(constitutive_dislotwin_postResults(c+j),state%p(5*ns+3*nt+j)) + ! constitutive_dislotwin_postResults(c+j)=max(constitutive_dislotwin_postResults(c+j),state%p(4*ns+2*nt+j)) enddo; enddo c = c + ns case (resolved_stress_shearband_ID) @@ -1740,7 +1740,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) enddo c = c + 6_pInt case (twin_fraction_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt)) c = c + nt case (shear_rate_twin_ID) if (nt > 0_pInt) then @@ -1758,16 +1758,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) StressRatio_p = 0.0_pReal StressRatio_pminus1 = 0.0_pReal else - StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + StressRatio_p = (abs(tau)/state%p(5_pInt*ns+3_pInt*nt+j))**& constitutive_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + StressRatio_pminus1 = (abs(tau)/state%p(5_pInt*ns+3_pInt*nt+j))**& (constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) endif !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) !* Initial shear rates DotGamma0 = & - state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & + state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & constitutive_dislotwin_v0PerSlipSystem(j,instance) !* Shear rates due to slip @@ -1784,7 +1784,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) !* Resolved shear stress on twin system tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,phase)) !* Stress ratios - StressRatio_r = (state(ipc,ip,el)%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_rPerTwinFamily(f,instance) + StressRatio_r = (state%p(7_pInt*ns+4_pInt*nt+j)/tau)**constitutive_dislotwin_rPerTwinFamily(f,instance) !* Shear rates due to twin if ( tau > 0.0_pReal ) then @@ -1793,8 +1793,8 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) if (tau < constitutive_dislotwin_tau_r(j,instance)) then - Ndot0=(abs(gdot_slip(s1))*(state(ipc,ip,el)%p(s2)+state(ipc,ip,el)%p(ns+s2))+& - abs(gdot_slip(s2))*(state(ipc,ip,el)%p(s1)+state(ipc,ip,el)%p(ns+s1)))/& + Ndot0=(abs(gdot_slip(s1))*(state%p(s2)+state%p(ns+s2))+& + abs(gdot_slip(s2))*(state%p(s1)+state%p(ns+s1)))/& (constitutive_dislotwin_L0(instance)*& constitutive_dislotwin_burgersPerSlipSystem(j,instance))*& (1.0_pReal-exp(-constitutive_dislotwin_VcrossSlip(instance)/(kB*Temperature)*& @@ -1807,17 +1807,17 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) end select constitutive_dislotwin_postResults(c+j) = & (constitutive_dislotwin_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*& - state(ipc,ip,el)%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r) + state%p(7_pInt*ns+5_pInt*nt+j)*Ndot0*exp(-StressRatio_r) endif enddo ; enddo endif c = c + nt case (accumulated_shear_twin_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((3_pInt*ns+nt+1_pInt):(3_pInt*ns+2_pInt*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((3_pInt*ns+nt+1_pInt):(3_pInt*ns+2_pInt*nt)) c = c + nt case (mfp_twin_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt)) c = c + nt case (resolved_stress_twin_ID) if (nt > 0_pInt) then @@ -1831,7 +1831,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) endif c = c + nt case (threshold_stress_twin_ID) - constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(ipc,ip,el)%p((7_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+5_pInt*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state%p((7_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+5_pInt*nt)) c = c + nt case (stress_exponent_ID) j = 0_pInt @@ -1846,16 +1846,16 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) StressRatio_p = 0.0_pReal StressRatio_pminus1 = 0.0_pReal else - StressRatio_p = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_p = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**& constitutive_dislotwin_pPerSlipFamily(f,instance) - StressRatio_pminus1 = (abs(tau)/state(ipc,ip,el)%p(6_pInt*ns+4_pInt*nt+j))**& + StressRatio_pminus1 = (abs(tau)/state%p(6_pInt*ns+4_pInt*nt+j))**& (constitutive_dislotwin_pPerSlipFamily(f,instance)-1.0_pReal) endif !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(j,instance)/(kB*Temperature) !* Initial shear rates DotGamma0 = & - state(ipc,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & + state%p(j)*constitutive_dislotwin_burgersPerSlipSystem(j,instance)* & constitutive_dislotwin_v0PerSlipSystem(j,instance) !* Shear rates due to slip @@ -1866,7 +1866,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el) dgdot_dtauslip = & ((abs(gdot_slip(j))*BoltzmannRatio*& constitutive_dislotwin_pPerSlipFamily(f,instance)*constitutive_dislotwin_qPerSlipFamily(f,instance))/& - state(ipc,ip,el)%p(6*ns+4*nt+j))*StressRatio_pminus1*(1_pInt-StressRatio_p)**& + state%p(6*ns+4*nt+j))*StressRatio_pminus1*(1_pInt-StressRatio_p)**& (constitutive_dislotwin_qPerSlipFamily(f,instance)-1.0_pReal) !* Stress exponent diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 255e5da51..13e91b968 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -319,18 +319,18 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip phase_plasticityInstance implicit none - real(pReal), dimension(3,3), intent(out) :: & + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & + real(pReal), dimension(9,9), intent(out) :: & dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(3,3) :: & @@ -355,7 +355,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip dLp_dTstar99 = 0.0_pReal else gamma_dot = constitutive_j2_gdot0(instance) & - * (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state(ipc,ip,el)%p(1)) & + * (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state%p(1)) & **constitutive_j2_n(instance) Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(instance) @@ -395,13 +395,13 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el) implicit none real(pReal), dimension(1) :: & constitutive_j2_dotState - real(pReal), dimension(6), intent(in):: & + real(pReal), dimension(6), intent(in):: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(6) :: & @@ -425,7 +425,7 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el) ! strain rate gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!----------------------------------------------------------------------------------- - (constitutive_j2_fTaylor(instance) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(instance) + (constitutive_j2_fTaylor(instance) * state%p(1)) ) ** constitutive_j2_n(instance) !-------------------------------------------------------------------------------------------------- ! hardening coefficient @@ -447,8 +447,8 @@ pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el) ) endif hardening = ( constitutive_j2_h0(instance) + constitutive_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) & - * abs( 1.0_pReal - state(ipc,ip,el)%p(1)/saturation )**constitutive_j2_a(instance) & - * sign(1.0_pReal, 1.0_pReal - state(ipc,ip,el)%p(1)/saturation) + * abs( 1.0_pReal - state%p(1)/saturation )**constitutive_j2_a(instance) & + * sign(1.0_pReal, 1.0_pReal - state%p(1)/saturation) else hardening = 0.0_pReal endif @@ -476,13 +476,13 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el) phase_Noutput implicit none - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & @@ -511,13 +511,13 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el) outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_j2_outputID(o,instance)) case (flowstress_ID) - constitutive_j2_postResults(c+1_pInt) = state(ipc,ip,el)%p(1) + constitutive_j2_postResults(c+1_pInt) = state%p(1) c = c + 1_pInt case (strainrate_ID) constitutive_j2_postResults(c+1_pInt) = & constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!---------------------------------------------------------------------------------- - (constitutive_j2_fTaylor(instance) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(instance) + (constitutive_j2_fTaylor(instance) * state%p(1)) ) ** constitutive_j2_n(instance) c = c + 1_pInt end select enddo outputsLoop diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 8048abe2b..89218ab4f 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -1297,7 +1297,7 @@ end function constitutive_nonlocal_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, gr, ip, el) +subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, ipc, ip, el) use IO, only: & IO_error @@ -1347,7 +1347,7 @@ use lattice, only: & implicit none !*** input variables -integer(pInt), intent(in) :: gr, & ! current grain ID +integer(pInt), intent(in) :: ipc, & ! current grain ID ip, & ! current integration point el ! current element real(pReal), dimension(3,3), intent(in) :: & @@ -1387,7 +1387,7 @@ real(pReal), dimension(2) :: rhoExcessGradient, & rhoTotal real(pReal), dimension(3) :: rhoExcessDifferences, & normal_latticeConf -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & rhoForest, & ! forest dislocation density tauBack, & ! back stress from pileup on same slip system tauThreshold ! threshold shear stress @@ -1397,24 +1397,24 @@ real(pReal), dimension(3,3) :: invFe, & ! inverse of elast invConnections real(pReal), dimension(3,mesh_maxNipNeighbors) :: & connection_latticeConf -real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: & +real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & rhoExcess -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & rhoDip ! dipole dislocation density (edge, screw) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))), & - totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))), & + totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & myInteractionMatrix ! corrected slip interaction matrix real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: & neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: & +real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & m ! direction of dislocation motion logical inversionError -phase = material_phase(gr,ip,el) +phase = material_phase(ipc,ip,el) instance = phase_plasticityInstance(phase) ns = totalNslip(instance) @@ -1422,11 +1422,11 @@ ns = totalNslip(instance) !*** get basic states forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state(gr,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,instance)) + rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities + rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl) < significantRho(instance)) & rhoSgl = 0.0_pReal @@ -1488,7 +1488,7 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) if (neighbor_el > 0 .and. neighbor_ip > 0) then - neighbor_phase = material_phase(gr,neighbor_ip,neighbor_el) + neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) neighbor_instance = phase_plasticityInstance(neighbor_phase) neighbor_ns = totalNslip(neighbor_instance) if (.not. phase_localPlasticity(neighbor_phase) & @@ -1497,14 +1497,14 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance nRealNeighbors = nRealNeighbors + 1_pInt forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) neighbor_rhoExcess(c,s,n) = & - max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles - - max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) ! negative mobiles + max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles + - max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) ! negative mobiles neighbor_rhoTotal(c,s,n) = & - max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles - + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) & ! negative mobiles - + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads - + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) & ! negative deads - + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_instance)), 0.0_pReal) ! dipoles + max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles + + max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) & ! negative mobiles + + abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads + + abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) & ! negative deads + + max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_instance)), 0.0_pReal) ! dipoles endforall connection_latticeConf(1:3,n) = & math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & @@ -1586,17 +1586,17 @@ endif !*** set dependent states -state(gr,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest -state(gr,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold -state(gr,ip,el)%p(iTauB(1:ns,instance)) = tauBack +state(ipc,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest +state(ipc,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold +state(ipc,ip,el)%p(iTauB(1:ns,instance)) = tauBack #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == gr)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,gr + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,ipc write(6,*) write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 @@ -2004,7 +2004,7 @@ integer(pInt), intent(in) :: ipc, & ! current real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation !*** input/output variables -type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & +type(p_vec), intent(inout) :: & state ! current microstructural state !*** output variables @@ -2056,15 +2056,15 @@ ns = totalNslip(instance) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) - v(s,t) = state(ipc,ip,el)%p(iV(s,t,instance)) + rhoSgl(s,t) = max(state%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = state%p(iRhoB(s,t,instance)) + v(s,t) = state%p(iV(s,t,instance)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities - dUpperOld(s,c) = state(ipc,ip,el)%p(iD(s,c,instance)) + rhoDip(s,c) = max(state%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities + dUpperOld(s,c) = state%p(iD(s,c,instance)) endforall -tauBack = state(ipc,ip,el)%p(iTauB(1:ns,instance)) +tauBack = state%p(iTauB(1:ns,instance)) where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl) < significantRho(instance)) & @@ -2168,7 +2168,7 @@ forall (t=1_pInt:4_pInt) & !*** store new maximum dipole height in state forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - state(ipc,ip,el)%p(iD(s,c,instance)) = dUpper(s,c) + state%p(iD(s,c,instance)) = dUpper(s,c) @@ -3362,14 +3362,14 @@ pure function constitutive_nonlocal_postResults(Tstar_v,Fe,state,dotState,ipc,ip lattice_nu implicit none - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe !< elastic deformation gradient type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state !< microstructure state type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index a7c2d7adb..9c2511c22 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -630,18 +630,18 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar phase_plasticityInstance implicit none - real(pReal), dimension(3,3), intent(out) :: & + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & + real(pReal), dimension(9,9), intent(out) :: & dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state integer(pInt) :: & @@ -693,12 +693,12 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,phase) enddo gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*& + ((abs(tau_slip_pos(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))*& sign(1.0_pReal,tau_slip_pos(j)) gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* & - ((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*& + ((abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))*& sign(1.0_pReal,tau_slip_neg(j)) - Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + Lp = Lp + (1.0_pReal-state%p(index_F))*& ! 1-F (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase) !-------------------------------------------------------------------------------------------------- @@ -723,16 +723,16 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar j = 0_pInt twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family - do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family j = j+1_pInt !-------------------------------------------------------------------------------------------------- ! Calculation of Lp tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) - gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F constitutive_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& + (abs(tau_twin(j))/state%p(nSlip+j))**& constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,phase) @@ -777,13 +777,13 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) phase_plasticityInstance implicit none - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & @@ -820,17 +820,17 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(instance)*& - (1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state(ipc,ip,el)%p(index_F)**& + (1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state%p(index_F)**& constitutive_phenopowerlaw_twinB(instance)) c_SlipTwin = 0.0_pReal c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(instance)*& - state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(instance) + state%p(index_Gamma)**constitutive_phenopowerlaw_twinE(instance) c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(instance)*& - state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(instance) + state%p(index_F)**constitutive_phenopowerlaw_twinD(instance) !-------------------------------------------------------------------------------------------------- ! calculate left and right vectors and calculate dot gammas - ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state(ipc,ip,el)%p(index_F)) + ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state%p(index_F)) j = 0_pInt slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family @@ -838,10 +838,10 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) j = j+1_pInt left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part - right_SlipSlip(j) = abs(1.0_pReal-state(ipc,ip,el)%p(j) / & + right_SlipSlip(j) = abs(1.0_pReal-state%p(j) / & (constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) & **constitutive_phenopowerlaw_a_slip(instance)& - *sign(1.0_pReal,1.0_pReal-state(ipc,ip,el)%p(j) / & + *sign(1.0_pReal,1.0_pReal-state%p(j) / & (constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) right_TwinSlip(j) = 1.0_pReal ! no system-dependent part @@ -856,8 +856,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase)) enddo gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & - ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) & - +(abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))& + ((abs(tau_slip_pos(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) & + +(abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))& *sign(1.0_pReal,tau_slip_pos(j)) enddo enddo slipFamiliesLoop1 @@ -875,9 +875,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! Calculation of dot vol frac tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) - gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F constitutive_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& + (abs(tau_twin(j))/state%p(nSlip+j))**& constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) enddo enddo twinFamiliesLoop1 @@ -913,7 +913,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) c_TwinTwin * left_TwinTwin(j) * & dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,instance), & right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor - if (state(ipc,ip,el)%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 + if (state%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + & gdot_twin(j)/lattice_shearTwin(index_myFamily+i,phase) constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j)) @@ -950,13 +950,13 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) mesh_maxNips implicit none - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & @@ -987,11 +987,11 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_phenopowerlaw_outputID(o,instance)) case (resistance_slip_ID) - constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip) + constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(1:nSlip) c = c + nSlip case (accumulatedshear_slip_ID) - constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(index_accshear_slip:& + constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(index_accshear_slip:& index_accshear_slip+nSlip) c = c + nSlip @@ -1010,8 +1010,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase)) enddo constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & - ((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) & - +(abs(tau_slip_neg)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))& + ((abs(tau_slip_pos)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) & + +(abs(tau_slip_neg)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))& *sign(1.0_pReal,tau_slip_pos) enddo enddo slipFamiliesLoop1 @@ -1031,17 +1031,17 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) case (totalshear_ID) constitutive_phenopowerlaw_postResults(c+1_pInt) = & - state(ipc,ip,el)%p(index_Gamma) + state%p(index_Gamma) c = c + 1_pInt case (resistance_twin_ID) constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & - state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip) + state%p(1_pInt+nSlip:nTwin+nSlip) c = c + nTwin case (accumulatedshear_twin_ID) constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & - state(ipc,ip,el)%p(index_accshear_twin:index_accshear_twin+nTwin) + state%p(index_accshear_twin:index_accshear_twin+nTwin) c = c + nTwin case (shearrate_twin_ID) @@ -1051,9 +1051,9 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family j = j + 1_pInt tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) - constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state%p(index_F))*& ! 1-F constitutive_phenopowerlaw_gdot0_twin(instance)*& - (abs(tau)/state(ipc,ip,el)%p(j+nSlip))**& + (abs(tau)/state%p(j+nSlip))**& constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau)) enddo enddo twinFamiliesLoop1 @@ -1072,7 +1072,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) c = c + nTwin case (totalvolfrac_ID) - constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_F) + constitutive_phenopowerlaw_postResults(c+1_pInt) = state%p(index_F) c = c + 1_pInt end select diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 9cd22d707..8a9a8a008 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -1108,7 +1108,7 @@ implicit none ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & volumefraction_PerTwinSys @@ -1130,7 +1130,7 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance !-------------------------------------------------------------------------------------------------- ! total twin volume fraction do i=1_pInt,nt - volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ & + volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance) enddo sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0 @@ -1170,7 +1170,7 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el) el !< element real(pReal), intent(in) :: & temperature !< temperature at IP - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + type(p_vec), intent(inout) :: & state !< microstructure state integer(pInt) :: & @@ -1192,13 +1192,10 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! total twin volume fraction - do i=1_pInt,nt - volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ & + forall (i = 1_pInt:nt) & + volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance) - enddo - - sumf = sum(abs(volumefraction_PerTwinSys(1:nt))) ! safe for nt == 0 sfe = 0.0002_pReal*Temperature-0.0396_pReal @@ -1206,43 +1203,43 @@ subroutine constitutive_titanmod_microstructure(temperature,state,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! average segment length for edge dislocations in matrix forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ & - sqrt(dot_product(state(ipc,ip,el)%p(1:ns), & + state%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSys(s,instance)/ & + sqrt(dot_product(state%p(1:ns), & constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ & - dot_product(state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns), & + dot_product(state%p(ns+1_pInt:2_pInt*ns), & constitutive_titanmod_forestProjectionScrew(1:ns,s,instance))) !-------------------------------------------------------------------------------------------------- ! average segment length for screw dislocations in matrix forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ & - sqrt(dot_product(state(ipc,ip,el)%p(1:ns), & + state%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSys(s,instance)/ & + sqrt(dot_product(state%p(1:ns), & constitutive_titanmod_forestProjectionEdge(1:ns,s,instance))+ & - dot_product(state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns), & + dot_product(state%p(ns+1_pInt:2_pInt*ns), & constitutive_titanmod_forestProjectionScrew(1:ns,s,instance))) !-------------------------------------------------------------------------------------------------- ! threshold stress or slip resistance for edge dislocation motion forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(5_pInt*ns+nt+s) = & + state%p(5_pInt*ns+nt+s) = & lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*& - sqrt(dot_product((state(ipc,ip,el)%p(1:ns)),& + sqrt(dot_product((state%p(1:ns)),& constitutive_titanmod_interactionMatrix_ee(1:ns,s,instance))+ & - dot_product((state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),& + dot_product((state%p(ns+1_pInt:2_pInt*ns)),& constitutive_titanmod_interactionMatrix_es(1:ns,s,instance))) !-------------------------------------------------------------------------------------------------- ! threshold stress or slip resistance for screw dislocation motion forall (s = 1_pInt:ns) & - state(ipc,ip,el)%p(6_pInt*ns+nt+s) = & + state%p(6_pInt*ns+nt+s) = & lattice_mu(phase)*constitutive_titanmod_burgersPerSlipSys(s,instance)*& - sqrt(dot_product((state(ipc,ip,el)%p(1:ns)),& + sqrt(dot_product((state%p(1:ns)),& constitutive_titanmod_interactionMatrix_es(1:ns,s,instance))+ & - dot_product((state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns)),& + dot_product((state%p(ns+1_pInt:2_pInt*ns)),& constitutive_titanmod_interactionMatrix_ss(1:ns,s,instance))) !-------------------------------------------------------------------------------------------------- ! threshold stress or slip resistance for dislocation motion in twin forall (t = 1_pInt:nt) & - state(ipc,ip,el)%p(7_pInt*ns+nt+t) = & + state%p(7_pInt*ns+nt+t) = & lattice_mu(phase)*constitutive_titanmod_burgersPerTwinSys(t,instance)*& - (dot_product((abs(state(ipc,ip,el)%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),& + (dot_product((abs(state%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),& constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,instance))) end subroutine constitutive_titanmod_microstructure @@ -1278,20 +1275,20 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& phase_plasticityInstance implicit none - real(pReal), dimension(3,3), intent(out) :: & + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & + real(pReal), dimension(9,9), intent(out) :: & dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress - real(pReal), dimension(6), intent(in) :: & + real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & temperature !< temperature at IP - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + type(p_vec), intent(inout) :: & state !< microstructure state integer(pInt) :: & index_myFamily, instance,phase, & @@ -1318,7 +1315,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& nt = constitutive_titanmod_totalNtwin(instance) do i=1_pInt,nt - volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ & + volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance) enddo @@ -1346,7 +1343,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,phase)) if(lattice_structure(phase)==LATTICE_hex_ID) then ! only for prismatic and pyr systems in hex screwvelocity_prefactor=constitutive_titanmod_debyefrequency(instance)* & - state(ipc,ip,el)%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ & + state%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSys(j,instance)/ & constitutive_titanmod_kinkcriticallength_PerSlipSys(j,instance))**2 !* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress @@ -1371,7 +1368,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& else ! if the structure is not hex or the slip family is basal screwvelocity_prefactor=constitutive_titanmod_v0s_PerSlipSys(j,instance) - bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+state(ipc,ip,el)%p(6*ns+nt+j) + bottomstress_screw=constitutive_titanmod_tau0s_PerSlipSys(j,instance)+state%p(6*ns+nt+j) StressRatio_screw_p = ((abs(tau_slip(j)))/( bottomstress_screw ))**constitutive_titanmod_ps_PerSlipSys(j,instance) if((1.0_pReal-StressRatio_screw_p)>0.001_pReal) then @@ -1389,7 +1386,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& endif !* Stress ratio for edge - bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+state(ipc,ip,el)%p(5*ns+nt+j) + bottomstress_edge=constitutive_titanmod_tau0e_PerSlipSys(j,instance)+state%p(5*ns+nt+j) StressRatio_edge_p = ((abs(tau_slip(j)))/ & ( bottomstress_edge) & )**constitutive_titanmod_pe_PerSlipSys(j,instance) @@ -1415,29 +1412,29 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& constitutive_titanmod_qe_PerSlipSys(j,instance)) !* Shear rates due to edge slip - gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state(ipc,ip,el)%p(j)* & + gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(j)* & edge_velocity(j))* sign(1.0_pReal,tau_slip(j)) !* Shear rates due to screw slip - gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state(ipc,ip,el)%p(ns+j) * & + gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(state%p(ns+j) * & screw_velocity(j))* sign(1.0_pReal,tau_slip(j)) !Total shear rate gdot_slip(j) = gdot_slip_edge(j) + gdot_slip_screw(j) - state(ipc,ip,el)%p(7*ns+2*nt+j)=edge_velocity(j) - state(ipc,ip,el)%p(8*ns+2*nt+j)=screw_velocity(j) - state(ipc,ip,el)%p(9*ns+2*nt+j)=tau_slip(j) - state(ipc,ip,el)%p(10*ns+2*nt+j)=gdot_slip_edge(j) - state(ipc,ip,el)%p(11*ns+2*nt+j)=gdot_slip_screw(j) - state(ipc,ip,el)%p(12*ns+2*nt+j)=StressRatio_edge_p - state(ipc,ip,el)%p(13*ns+2*nt+j)=StressRatio_screw_p + state%p(7*ns+2*nt+j)=edge_velocity(j) + state%p(8*ns+2*nt+j)=screw_velocity(j) + state%p(9*ns+2*nt+j)=tau_slip(j) + state%p(10*ns+2*nt+j)=gdot_slip_edge(j) + state%p(11*ns+2*nt+j)=gdot_slip_screw(j) + state%p(12*ns+2*nt+j)=StressRatio_edge_p + state%p(13*ns+2*nt+j)=StressRatio_screw_p !* Derivatives of shear rates dgdot_dtauslip(j) = constitutive_titanmod_burgersPerSlipSys(j,instance)*(( & ( & ( & ( & - (edge_velocity(j)*state(ipc,ip,el)%p(j))) * & + (edge_velocity(j)*state%p(j))) * & BoltzmannRatioedge*& constitutive_titanmod_pe_PerSlipSys(j,instance)* & constitutive_titanmod_qe_PerSlipSys(j,instance) & @@ -1450,7 +1447,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& ( & ( & ( & - (state(ipc,ip,el)%p(ns+j) * screw_velocity(j)) * & + (state%p(ns+j) * screw_velocity(j)) * & BoltzmannRatioscrew* & constitutive_titanmod_ps_PerSlipSys(j,instance)* & constitutive_titanmod_qs_PerSlipSys(j,instance) & @@ -1492,20 +1489,20 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& !************************************************************************************** !* Stress ratios -! StressRatio_r = (state(ipc,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_titanmod_r(instance) +! StressRatio_r = (state%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_titanmod_r(instance) !* Shear rates and their derivatives due to twin ! if ( tau_twin(j) > 0.0_pReal ) !then ! gdot_twin(j) = 0.0_pReal!& ! (constitutive_titanmod_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,phase)*& -! state(ipc,ip,el)%p(6*ns+4*nt+j)*constitutive_titanmod_Ndot0PerTwinSys(f,instance)*exp(-StressRatio_r) +! state%p(6*ns+4*nt+j)*constitutive_titanmod_Ndot0PerTwinSys(f,instance)*exp(-StressRatio_r) ! dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_titanmod_r(instance))/tau_twin(j))*StressRatio_r ! endif !************************************************************************************** !* Stress ratio for edge twinStressRatio_p = ((abs(tau_twin(j)))/ & - ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) & + ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) & )**constitutive_titanmod_twinp_PerTwinSys(j,instance) if((1.0_pReal-twinStressRatio_p)>0.001_pReal) then @@ -1515,7 +1512,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& endif twinStressRatio_pminus1 = ((abs(tau_twin(j)))/ & - ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) & + ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) & )**(constitutive_titanmod_twinp_PerTwinSys(j,instance)-1.0_pReal) !* Boltzmann ratio @@ -1584,15 +1581,15 @@ function constitutive_titanmod_dotState(Tstar_v,temperature,state,ipc,ip,el) phase_plasticityInstance implicit none - real(pReal), dimension(6), intent(in):: & + real(pReal), dimension(6), intent(in):: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & + real(pReal), intent(in) :: & temperature !< temperature at integration point - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_titanmod_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & constitutive_titanmod_dotState @@ -1622,7 +1619,7 @@ implicit none nt = constitutive_titanmod_totalNtwin(instance) do i=1_pInt,nt - volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ & + volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance) enddo @@ -1638,13 +1635,13 @@ implicit none j = j+1_pInt DotRhoEdgeGeneration(j) = & ! multiplication of edge dislocations - state(ipc,ip,el)%p(ns+j)*state(ipc,ip,el)%p(8*ns+2*nt+j)/state(ipc,ip,el)%p(4*ns+nt+j) + state%p(ns+j)*state%p(8*ns+2*nt+j)/state%p(4*ns+nt+j) DotRhoScrewGeneration(j) = & ! multiplication of screw dislocations - state(ipc,ip,el)%p(j)*state(ipc,ip,el)%p(7*ns+2*nt+j)/state(ipc,ip,el)%p(3*ns+nt+j) - DotRhoEdgeAnnihilation(j) = -((state(ipc,ip,el)%p(j))**2)* & ! annihilation of edge dislocations - constitutive_titanmod_capre_PerSlipSys(j,instance)*state(ipc,ip,el)%p(7*ns+2*nt+j)*0.5_pReal - DotRhoScrewAnnihilation(j) = -((state(ipc,ip,el)%p(ns+j))**2)* & ! annihilation of screw dislocations - constitutive_titanmod_caprs_PerSlipSys(j,instance)*state(ipc,ip,el)%p(8*ns+2*nt+j)*0.5_pReal + state%p(j)*state%p(7*ns+2*nt+j)/state%p(3*ns+nt+j) + DotRhoEdgeAnnihilation(j) = -((state%p(j))**2)* & ! annihilation of edge dislocations + constitutive_titanmod_capre_PerSlipSys(j,instance)*state%p(7*ns+2*nt+j)*0.5_pReal + DotRhoScrewAnnihilation(j) = -((state%p(ns+j))**2)* & ! annihilation of screw dislocations + constitutive_titanmod_caprs_PerSlipSys(j,instance)*state%p(8*ns+2*nt+j)*0.5_pReal constitutive_titanmod_dotState(j) = & ! edge dislocation density rate of change DotRhoEdgeGeneration(j)+DotRhoEdgeAnnihilation(j) @@ -1652,7 +1649,7 @@ implicit none DotRhoScrewGeneration(j)+DotRhoScrewAnnihilation(j) constitutive_titanmod_dotState(2*ns+j) = & ! sum of shear due to edge and screw - state(ipc,ip,el)%p(10*ns+2*nt+j)+state(ipc,ip,el)%p(11*ns+2*nt+j) + state%p(10*ns+2*nt+j)+state%p(11*ns+2*nt+j) enddo enddo slipFamiliesLoop @@ -1668,7 +1665,7 @@ implicit none !* Stress ratio for edge twinStressRatio_p = ((abs(tau_twin(j)))/ & - ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state(ipc,ip,el)%p(7*ns+nt+j)) & + ( constitutive_titanmod_twintau0_PerTwinSys(j,instance)+state%p(7*ns+nt+j)) & )**(constitutive_titanmod_twinp_PerTwinSys(j,instance)) @@ -1708,11 +1705,11 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el) phase_Noutput implicit none - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + type(p_vec), intent(in) :: & state !< microstructure state real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & constitutive_titanmod_postResults @@ -1734,7 +1731,7 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el) nt = constitutive_titanmod_totalNtwin(instance) do i=1_pInt,nt - volumefraction_PerTwinSys(i)=state(ipc,ip,el)%p(3_pInt*ns+i)/ & + volumefraction_PerTwinSys(i)=state%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSys(i,instance) enddo @@ -1749,91 +1746,91 @@ pure function constitutive_titanmod_postResults(state,ipc,ip,el) do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_titanmod_outputID(o,instance)) case (rhoedge_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(1_pInt:ns) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(1_pInt:ns) c = c + ns case (rhoscrew_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(ns+1_pInt:2_pInt*ns) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(ns+1_pInt:2_pInt*ns) c = c + ns case (segment_edge_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt) c = c + ns case (segment_screw_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt) c = c + ns case (resistance_edge_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt) c = c + ns case (resistance_screw_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt) c = c + ns case (velocity_edge_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(7*ns+2*nt+1:8*ns+2*nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(7*ns+2*nt+1:8*ns+2*nt) c = c + ns case (velocity_screw_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(ipc,ip,el)%p(8*ns+2*nt+1:9*ns+2*nt) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state%p(8*ns+2*nt+1:9*ns+2*nt) c = c + ns case (tau_slip_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(9*ns+2*nt+1:10*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(9*ns+2*nt+1:10*ns+2*nt)) c = c + ns case (gdot_slip_edge_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(10*ns+2*nt+1:11*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt)) c = c + ns case (gdot_slip_screw_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(11*ns+2*nt+1:12*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(11*ns+2*nt+1:12*ns+2*nt)) c = c + ns case (gdot_slip_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(10*ns+2*nt+1:11*ns+2*nt)) + & - abs(state(ipc,ip,el)%p(11*ns+2*nt+1:12*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(10*ns+2*nt+1:11*ns+2*nt)) + & + abs(state%p(11*ns+2*nt+1:12*ns+2*nt)) c = c + ns case (stressratio_edge_p_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(12*ns+2*nt+1:13*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(12*ns+2*nt+1:13*ns+2*nt)) c = c + ns case (stressratio_screw_p_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(13*ns+2*nt+1:14*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(13*ns+2*nt+1:14*ns+2*nt)) c = c + ns case (shear_system_ID) - constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(ipc,ip,el)%p(2*ns+1:3*ns)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state%p(2*ns+1:3*ns)) c = c + ns case (shear_basal_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+1:2*ns+3))) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:2*ns+3))) c = c + 1_pInt case (shear_prism_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+4:2*ns+6))) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+4:2*ns+6))) c = c + 1_pInt case (shear_pyra_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+7:2*ns+12))) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+7:2*ns+12))) c = c + 1_pInt case (shear_pyrca_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+13:2*ns+24))) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+13:2*ns+24))) c = c + 1_pInt case (rhoedge_basal_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(1:3)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(1:3)) c = c + 1_pInt case (rhoedge_prism_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(4:6)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(4:6)) c = c + 1_pInt case (rhoedge_pyra_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(7:12)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(7:12)) c = c + 1_pInt case (rhoedge_pyrca_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(13:24)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(13:24)) c = c + 1_pInt case (rhoscrew_basal_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+1:ns+3)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+1:ns+3)) c = c + 1_pInt case (rhoscrew_prism_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+4:ns+6)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+4:ns+6)) c = c + 1_pInt case (rhoscrew_pyra_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+7:ns+12)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+7:ns+12)) c = c + 1_pInt case (rhoscrew_pyrca_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(ipc,ip,el)%p(ns+13:ns+24)) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state%p(ns+13:ns+24)) c = c + 1_pInt case (shear_total_ID) - constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(ipc,ip,el)%p(2*ns+1:3*ns))) + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state%p(2*ns+1:3*ns))) c = c + 1_pInt case (twin_fraction_ID) constitutive_titanmod_postResults(c+1_pInt:c+nt) = abs(volumefraction_PerTwinSys(1:nt))