From ba9bd9120ea92175cf39f4009d5557affae4ecc8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Feb 2020 07:25:36 +0100 Subject: [PATCH] unifying code style --- src/source_damage_anisoBrittle.f90 | 78 +++++++++++++++--------------- src/source_damage_anisoDuctile.f90 | 49 ++++++++++--------- src/source_damage_isoBrittle.f90 | 43 +++++++++------- src/source_damage_isoDuctile.f90 | 34 +++++++------ 4 files changed, 110 insertions(+), 94 deletions(-) diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index c8224176a..a532e3e4e 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -9,8 +9,8 @@ module source_damage_anisoBrittle use debug use IO use math - use material use discretization + use material use config use lattice use results @@ -58,7 +58,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoBrittle_init - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + integer :: Ninstance,source,sourceOffset,NofMyPhase,p character(len=pStringLen) :: & extmsg = '' @@ -72,7 +72,6 @@ subroutine source_damage_anisoBrittle_init allocate(source_damage_anisoBrittle_instance(size(config_phase)), source=0) allocate(param(Ninstance)) - do p = 1, size(config_phase) source_damage_anisoBrittle_instance(p) = count(phase_source(:,1:p) == SOURCE_DAMAGE_ANISOBRITTLE_ID) do source = 1, phase_Nsources(p) @@ -84,42 +83,43 @@ subroutine source_damage_anisoBrittle_init associate(prm => param(source_damage_anisoBrittle_instance(p)), & config => config_phase(p)) + prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) + prm%totalNcleavage = sum(prm%Ncleavage) + prm%aTol = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('anisobrittle_ratesensitivity') prm%sdot_0 = config%getFloat('anisobrittle_sdot0') - prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray) - prm%totalNcleavage = sum(prm%Ncleavage) + prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage)) prm%critLoad = config%getFloats('anisobrittle_criticalload', requiredSize=size(prm%Ncleavage)) - prm%cleavage_systems = lattice_SchmidMatrix_cleavage (prm%Ncleavage,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%cleavage_systems = lattice_SchmidMatrix_cleavage(prm%Ncleavage,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) ! expand: family => system - prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) - prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) + prm%critDisp = math_expand(prm%critDisp, prm%Ncleavage) + prm%critLoad = math_expand(prm%critLoad, prm%Ncleavage) - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' - if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' - if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' - if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' + ! sanity checks + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity' + if (prm%sdot_0 <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0' + if (any(prm%critLoad < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticalload' + if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') + if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')') !-------------------------------------------------------------------------------------------------- ! output pararameters prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) - NofMyPhase=count(material_phaseAt==p) * discretization_nIP - instance = source_damage_anisoBrittle_instance(p) + NofMyPhase = count(material_phaseAt==p) * discretization_nIP sourceOffset = source_damage_anisoBrittle_offset(p) call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) - sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate enddo @@ -138,44 +138,43 @@ subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) el !< element real(pReal), intent(in), dimension(3,3) :: & S + integer :: & phase, & constituent, & - instance, & sourceOffset, & damageOffset, & homog, & - index + i real(pReal) :: & traction_d, traction_t, traction_n, traction_crit phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_anisoBrittle_instance(phase) sourceOffset = source_damage_anisoBrittle_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) + associate(prm => param(source_damage_anisoBrittle_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal + do i = 1, prm%totalNcleavage - do index = 1, param(instance)%totalNcleavage + traction_d = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,1,i)) + traction_t = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,2,i)) + traction_n = math_mul33xx33(S,prm%cleavage_systems(1:3,1:3,3,i)) - traction_d = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,1,index)) - traction_t = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,2,index)) - traction_n = math_mul33xx33(S,param(instance)%cleavage_systems(1:3,1:3,3,index)) + traction_crit = prm%critLoad(i)*damage(homog)%p(damageOffset)**2.0_pReal - traction_crit = param(instance)%critLoad(index)* & - damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset) - - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - param(instance)%sdot_0* & - ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + & - (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ & - param(instance)%critDisp(index) + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + + prm%sdot_0* & + ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%N + & + (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%N)/ & + prm%critDisp(i) enddo + end associate end subroutine source_damage_anisoBrittle_dotState @@ -193,16 +192,17 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi + integer :: & sourceOffset sourceOffset = source_damage_anisoBrittle_offset(phase) - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi + end subroutine source_damage_anisoBrittle_getRateAndItsTangent diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index a460272ea..169543dbd 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -21,7 +21,7 @@ module source_damage_anisoDuctile source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type, private :: tParameters !< container type for internal constitutive parameters real(pReal) :: & aTol, & N @@ -35,7 +35,7 @@ module source_damage_anisoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) public :: & @@ -53,7 +53,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_anisoDuctile_init - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + integer :: Ninstance,source,sourceOffset,NofMyPhase,p character(len=pStringLen) :: & extmsg = '' @@ -78,18 +78,21 @@ subroutine source_damage_anisoDuctile_init associate(prm => param(source_damage_anisoDuctile_instance(p)), & config => config_phase(p)) + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + prm%aTol = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal) prm%N = config%getFloat('anisoductile_ratesensitivity') - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip)) - prm%totalNslip = sum(prm%Nslip) ! expand: family => system - prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) + ! expand: family => system + prm%critPlasticStrain = math_expand(prm%critPlasticStrain, prm%Nslip) ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity' - if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity' + if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -100,11 +103,10 @@ subroutine source_damage_anisoDuctile_init prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyPhase=count(material_phaseAt==p) * discretization_nIP - instance = source_damage_anisoDuctile_instance(p) sourceOffset = source_damage_anisoDuctile_offset(p) call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) - sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate enddo @@ -121,28 +123,28 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element + integer :: & phase, & constituent, & sourceOffset, & - homog, damageOffset, & - instance, & + damageOffset, & + homog, & i phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_anisoDuctile_instance(phase) sourceOffset = source_damage_anisoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) - - do i = 1, param(instance)%totalNslip - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + & - plasticState(phase)%slipRate(i,constituent)/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(i) + associate(prm => param(source_damage_anisoDuctile_instance(phase))) + do i = 1, prm%totalNslip + sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + + plasticState(phase)%slipRate(i,constituent)/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain(i) enddo + end associate end subroutine source_damage_anisoDuctile_dotState @@ -160,16 +162,17 @@ subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalph real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi + integer :: & sourceOffset sourceOffset = source_damage_anisoDuctile_offset(phase) - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi + end subroutine source_damage_anisoDuctile_getRateAndItsTangent diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 30f76d38a..d8ae7768c 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -30,7 +30,7 @@ module source_damage_isoBrittle output end type tParameters - type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) public :: & @@ -48,7 +48,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoBrittle_init - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + integer :: Ninstance,source,sourceOffset,NofMyPhase,p character(len=pStringLen) :: & extmsg = '' @@ -78,9 +78,9 @@ subroutine source_damage_isoBrittle_init prm%critStrainEnergy = config%getFloat('isobrittle_criticalstrainenergy') ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol' - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n' - if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n' + if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -92,17 +92,17 @@ subroutine source_damage_isoBrittle_init prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyPhase = count(material_phaseAt==p) * discretization_nIP - instance = source_damage_isoBrittle_instance(p) sourceOffset = source_damage_isoBrittle_offset(p) call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,1) - sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate enddo end subroutine source_damage_isoBrittle_init + !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- @@ -116,22 +116,25 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) Fe real(pReal), intent(in), dimension(6,6) :: & C + integer :: & - phase, constituent, instance, sourceOffset + phase, & + constituent, & + sourceOffset + real(pReal), dimension(6) :: & + strain real(pReal) :: & - strain(6), & strainenergy phase = material_phaseAt(ipc,el) !< phase ID at ipc,ip,el constituent = material_phasememberAt(ipc,ip,el) !< state array offset for phase ID at ipc,ip,el ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on! - instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source sourceOffset = source_damage_isoBrittle_offset(phase) - strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) - strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy + associate(prm => param(source_damage_isoBrittle_instance(phase))) + strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%critStrainEnergy ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then @@ -142,9 +145,11 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & sourceState(phase)%p(sourceOffset)%state(1,constituent) endif + end associate end subroutine source_damage_isoBrittle_deltaState + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- @@ -158,17 +163,19 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi - integer :: & - instance, sourceOffset - instance = source_damage_isoBrittle_instance(phase) + integer :: & + sourceOffset + sourceOffset = source_damage_isoBrittle_offset(phase) - localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - & + associate(prm => param(source_damage_isoBrittle_instance(phase))) + localphiDot = (1.0_pReal - phi)**(prm%n - 1.0_pReal) - & phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* & - (1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) & + dLocalphiDot_dPhi = - (prm%n - 1.0_pReal)* & + (1.0_pReal - phi)**max(0.0_pReal,prm%n - 2.0_pReal) & - sourceState(phase)%p(sourceOffset)%state(1,constituent) + end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index de7efc90a..736ba1b8d 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -15,6 +15,7 @@ module source_damage_isoDuctile implicit none private + integer, dimension(:), allocatable :: & source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism @@ -46,7 +47,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine source_damage_isoDuctile_init - integer :: Ninstance,instance,source,sourceOffset,NofMyPhase,p + integer :: Ninstance,source,sourceOffset,NofMyPhase,p character(len=pStringLen) :: & extmsg = '' @@ -76,9 +77,9 @@ subroutine source_damage_isoDuctile_init prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain') ! sanity checks - if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol' - if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity' - if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' + if (prm%aTol < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol' + if (prm%N <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity' + if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' !-------------------------------------------------------------------------------------------------- ! exit if any parameter is out of range @@ -90,11 +91,10 @@ subroutine source_damage_isoDuctile_init prm%output = config%getStrings('(output)',defaultVal=emptyStringArray) NofMyPhase=count(material_phaseAt==p) * discretization_nIP - instance = source_damage_isoDuctile_instance(p) sourceOffset = source_damage_isoDuctile_offset(p) call material_allocateSourceState(p,sourceOffset,NofMyPhase,1,1,0) - sourceState(p)%p(sourceOffset)%aTolState=param(instance)%aTol + sourceState(p)%p(sourceOffset)%aTolState=prm%aTol end associate enddo @@ -111,23 +111,28 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element + integer :: & - phase, constituent, instance, homog, sourceOffset, damageOffset + phase, & + constituent, & + sourceOffset, & + damageOffset, & + homog phase = material_phaseAt(ipc,el) constituent = material_phasememberAt(ipc,ip,el) - instance = source_damage_isoDuctile_instance(phase) sourceOffset = source_damage_isoDuctile_offset(phase) homog = material_homogenizationAt(el) damageOffset = damageMapping(homog)%p(ip,el) + associate(prm => param(source_damage_isoDuctile_instance(phase))) sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & - sum(plasticState(phase)%slipRate(:,constituent))/ & - ((damage(homog)%p(damageOffset))**param(instance)%N)/ & - param(instance)%critPlasticStrain + sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%N)/prm%critPlasticStrain + end associate end subroutine source_damage_isoDuctile_dotState + !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- @@ -141,16 +146,17 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD real(pReal), intent(out) :: & localphiDot, & dLocalphiDot_dPhi + integer :: & sourceOffset sourceOffset = source_damage_isoDuctile_offset(phase) - localphiDot = 1.0_pReal & - - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + localphiDot = 1.0_pReal & + + dLocalphiDot_dPhi*phi + end subroutine source_damage_isoDuctile_getRateAndItsTangent