diff --git a/PRIVATE b/PRIVATE index 2a0f30ec4..571b06ac4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2a0f30ec47473f42e0da922055b72b527730378d +Subproject commit 571b06ac4b42ed217eacfa5c71281e571bdab1df diff --git a/src/phase.f90 b/src/phase.f90 index 1d08f680d..8ea1fa5f7 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -70,7 +70,7 @@ module phase type(tPlasticState), allocatable, dimension(:), public :: & plasticState - type(tSourceState), allocatable, dimension(:), public :: & + type(tState), allocatable, dimension(:), public :: & damageState @@ -350,14 +350,12 @@ subroutine phase_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%state = plasticState(ph)%state0 - forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 - end forall - - phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, & - maxval(damageState(ph)%p%sizeDotState)) + plasticState(ph)%state = plasticState(ph)%state0 + if(phase_Nsources(ph) > 0) & + damageState(ph)%state = damageState(ph)%state0 enddo PhaseLoop2 + + phase_source_maxSizeDotState = maxval(damageState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) end subroutine phase_init @@ -404,15 +402,13 @@ subroutine phase_restore(ce,includeL) integer, intent(in) :: ce integer :: & - co, & !< constituent number - so + co do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - do so = 1, phase_Nsources(material_phaseAt2(co,ce)) - damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) - enddo + if (phase_Nsources(material_phaseAt2(co,ce)) > 0) & + damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -426,16 +422,16 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph, so + integer :: ph call mechanical_forward() call thermal_forward() do ph = 1, size(damageState) - do so = 1,phase_Nsources(ph) - damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state - enddo; enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo end subroutine phase_forward @@ -531,9 +527,8 @@ subroutine crystallite_init() phases => config_material%get('phase') do ph = 1, phases%length - do so = 1, phase_Nsources(ph) - allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack - enddo + if (phase_Nsources(ph) > 0) & + allocate(damageState(ph)%subState0,source=damageState(ph)%state0) ! ToDo: hack enddo print'(a42,1x,i10)', ' # of elements: ', eMax @@ -579,9 +574,9 @@ subroutine phase_windForward(ip,el) call mechanical_windForward(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) + enddo diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 8defcf3a1..5f323c3a2 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -176,8 +176,8 @@ module subroutine damage_init phase => phases%get(ph) sources => phase%get('damage',defaultVal=emptyList) - phase_Nsources(ph) = sources%length - allocate(damageState(ph)%p(phase_Nsources(ph))) + if (sources%length > 1) error stop + enddo allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) @@ -267,57 +267,48 @@ module function integrateDamageState(dt,co,ip,el) result(broken) NiterationState, & !< number of iterations in state loop ph, & me, & - so - integer, dimension(maxval(phase_Nsources)) :: & size_so real(pReal) :: & zeta real(pReal), dimension(phase_source_maxSizeDotState) :: & r ! state residuum - real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + real(pReal), dimension(phase_source_maxSizeDotState,2) :: source_dotState logical :: & converged_ - ph = material_phaseAt(co,el) + if (phase_Nsources(ph) == 0) return me = material_phaseMemberAt(co,ip,el) converged_ = .true. broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) return - do so = 1, phase_Nsources(ph) - size_so(so) = damageState(ph)%p(so)%sizeDotState - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & - + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo + size_so = damageState(ph)%sizeDotState + damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) & + + damageState(ph)%dotState (1:size_so,me) * dt + source_dotState(1:size_so,2) = 0.0_pReal iteration: do NiterationState = 1, num%nState - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) - enddo + if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1) + source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me) broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration - do so = 1, phase_Nsources(ph) - zeta = damper(damageState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & - - damageState(ph)%p(so)%subState0(1:size_so(so),me) & - - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt - damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - damageState(ph)%p(so)%state(1:size_so(so),me), & - damageState(ph)%p(so)%atol(1:size_so(so))) - enddo + + zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2)) + damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta & + + source_dotState(1:size_so,1)* (1.0_pReal - zeta) + r(1:size_so) = damageState(ph)%state (1:size_so,me) & + - damageState(ph)%subState0(1:size_so,me) & + - damageState(ph)%dotState (1:size_so,me) * dt + damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so) + converged_ = converged_ .and. converged(r(1:size_so), & + damageState(ph)%state(1:size_so,me), & + damageState(ph)%atol(1:size_so)) + if(converged_) then broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me) @@ -423,7 +414,7 @@ function phase_damage_collectDotState(co,ip,el,ph,me) result(broken) end select sourceType - broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,me))) + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me))) enddo SourceLoop @@ -443,7 +434,6 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer :: & - so, & myOffset, & mySize logical :: & @@ -452,23 +442,22 @@ function phase_damage_deltaState(Fe, ph, me) result(broken) broken = .false. - sourceLoop: do so = 1, phase_Nsources(ph) + if (phase_Nsources(ph) == 0) return - sourceType: select case (phase_source(so,ph)) + sourceType: select case (phase_source(1,ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me) - broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) + broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me))) if(.not. broken) then - myOffset = damageState(ph)%p(so)%offsetDeltaState - mySize = damageState(ph)%p(so)%sizeDeltaState - damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) = & - damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%p(so)%deltaState(1:mySize,me) + myOffset = damageState(ph)%offsetDeltaState + mySize = damageState(ph)%sizeDeltaState + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = & + damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me) endif end select sourceType - enddo SourceLoop end function phase_damage_deltaState diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index ad7f80604..422b2a91a 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -6,8 +6,7 @@ !-------------------------------------------------------------------------------------------------- submodule (phase:damagee) anisobrittle - integer, dimension(:), allocatable :: & - source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? + integer, dimension(:), allocatable :: & source_damage_anisoBrittle_instance !< instance of source mechanism type :: tParameters !< container type for internal constitutive parameters @@ -58,7 +57,6 @@ module function anisobrittle_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_anisoBrittle_offset (phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -68,7 +66,6 @@ module function anisobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_anisoBrittle_instance(p))) src => sources%get(sourceOffset) @@ -101,9 +98,9 @@ module function anisobrittle_init(source_length) result(mySources) if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -139,14 +136,13 @@ module subroutine anisobrittle_dotState(S, co, ip, el) real(pReal) :: & traction_d, traction_t, traction_n, traction_crit + ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoBrittle_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) + associate(prm => param(source_damage_anisoBrittle_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = 0.0_pReal + damageState(ph)%dotState(1,me) = 0.0_pReal do i = 1, prm%sum_N_cl traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) @@ -154,8 +150,8 @@ module subroutine anisobrittle_dotState(S, co, ip, el) traction_crit = prm%g_crit(i)*phase_damage_get_phi(co,ip,el)**2.0_pReal - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = damageState(ph)%p(sourceOffset)%dotState(1,me) & + damageState(ph)%dotState(1,me) & + = damageState(ph)%dotState(1,me) & + prm%dot_o / prm%s_crit(i) & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & @@ -180,12 +176,8 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoBrittle_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -204,7 +196,7 @@ module subroutine anisobrittle_results(phase,group) integer :: o associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) + stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index b22f5b3d6..bafcf80a2 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) anisoductile integer, dimension(:), allocatable :: & - source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism type :: tParameters !< container type for internal constitutive parameters @@ -53,7 +52,6 @@ module function anisoductile_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_anisoDuctile_offset (phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -65,7 +63,6 @@ module function anisoductile_init(source_length) result(mySources) sources => phase%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_anisoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_anisoDuctile_instance(p))) src => sources%get(sourceOffset) @@ -87,9 +84,9 @@ module function anisoductile_init(source_length) result(mySources) if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -116,20 +113,14 @@ module subroutine anisoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_anisoDuctile_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) + associate(prm => param(source_damage_anisoDuctile_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) & - = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit) end associate end subroutine anisoductile_dotState @@ -149,12 +140,8 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_anisoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -173,7 +160,7 @@ module subroutine anisoductile_results(phase,group) integer :: o associate(prm => param(source_damage_anisoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) + stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index a2ae7dfb7..0a87493a6 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isobrittle integer, dimension(:), allocatable :: & - source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance type :: tParameters !< container type for internal constitutive parameters @@ -48,7 +47,6 @@ module function isobrittle_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_isoBrittle_offset (phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0) do p = 1, phases%length @@ -58,7 +56,6 @@ module function isobrittle_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoBrittle_offset(p) = sourceOffset associate(prm => param(source_damage_isoBrittle_instance(p))) src => sources%get(sourceOffset) @@ -74,9 +71,9 @@ module function isobrittle_init(source_length) result(mySources) if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + call phase_allocateState(damageState(p),Nconstituents,1,1,1) + damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -102,14 +99,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) real(pReal), intent(in), dimension(6,6) :: & C - integer :: & - sourceOffset real(pReal), dimension(6) :: & strain real(pReal) :: & strainenergy - sourceOffset = source_damage_isoBrittle_offset(ph) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) @@ -117,14 +111,11 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, ph,me) strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit - if (strainenergy > damageState(ph)%p(sourceOffset)%subState0(1,me)) then - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - strainenergy - damageState(ph)%p(sourceOffset)%state(1,me) - else - damageState(ph)%p(sourceOffset)%deltaState(1,me) = & - damageState(ph)%p(sourceOffset)%subState0(1,me) - & - damageState(ph)%p(sourceOffset)%state(1,me) - endif + if (strainenergy > damageState(ph)%subState0(1,me)) then + damageState(ph)%deltaState(1,me) = strainenergy - damageState(ph)%state(1,me) + else + damageState(ph)%deltaState(1,me) = damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me) + endif end associate end subroutine source_damage_isoBrittle_deltaState @@ -144,15 +135,11 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - - sourceOffset = source_damage_isoBrittle_offset(phase) associate(prm => param(source_damage_isoBrittle_instance(phase))) - localphiDot = 1.0_pReal & - - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + - phi*damageState(phase)%state(1,constituent) + dLocalphiDot_dPhi = - damageState(phase)%state(1,constituent) end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -169,7 +156,7 @@ module subroutine isobrittle_results(phase,group) integer :: o associate(prm => param(source_damage_isoBrittle_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) + stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index f378b7ad2..657f9da39 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -7,7 +7,6 @@ submodule(phase:damagee) isoductile integer, dimension(:), allocatable :: & - source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism type:: tParameters !< container type for internal constitutive parameters @@ -50,7 +49,6 @@ module function isoductile_init(source_length) result(mySources) phases => config_material%get('phase') allocate(param(Ninstances)) - allocate(source_damage_isoDuctile_offset (phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0) do p = 1, phases%length @@ -60,7 +58,6 @@ module function isoductile_init(source_length) result(mySources) sources => phase%get('damage') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then - source_damage_isoDuctile_offset(p) = sourceOffset associate(prm => param(source_damage_isoDuctile_instance(p))) src => sources%get(sourceOffset) @@ -78,9 +75,9 @@ module function isoductile_init(source_length) result(mySources) if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) - damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + call phase_allocateState(damageState(p),Nconstituents,1,1,0) + damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -107,20 +104,14 @@ module subroutine isoductile_dotState(co, ip, el) integer :: & ph, & - me, & - sourceOffset, & - damageOffset, & - homog + me ph = material_phaseAt(co,el) me = material_phasememberAt(co,ip,el) - sourceOffset = source_damage_isoDuctile_offset(ph) - homog = material_homogenizationAt(el) - damageOffset = material_homogenizationMemberAt(ip,el) + associate(prm => param(source_damage_isoDuctile_instance(ph))) - damageState(ph)%p(sourceOffset)%dotState(1,me) = & - sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit + damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me))/(phase_damage_get_phi(co,ip,el)**prm%q)/prm%gamma_crit end associate end subroutine isoductile_dotState @@ -140,12 +131,8 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo localphiDot, & dLocalphiDot_dPhi - integer :: & - sourceOffset - sourceOffset = source_damage_isoDuctile_offset(phase) - - dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -164,7 +151,7 @@ module subroutine isoductile_results(phase,group) integer :: o associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) + stt => damageState(phase)%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index f7631a38d..bc4b2cd47 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1138,7 +1138,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal) :: & formerSubStep integer :: & - so, ph, me, sizeDotState + ph, me, sizeDotState logical :: todo real(pReal) :: subFrac,subStep real(pReal), dimension(3,3) :: & @@ -1159,10 +1159,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%State0(:,me) + if (phase_Nsources(ph) > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) - enddo 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) @@ -1188,9 +1187,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) 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) - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me) + endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) @@ -1204,9 +1203,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 endif plasticState(ph)%state(:,me) = subState0 - do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) - enddo + if (phase_Nsources(ph) > 0) & + damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me) todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) endif