diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 1747b7080..40ee68c61 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -795,9 +795,12 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el phase_thermal, & material_phase, & LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_THERMAL_adiabatic_ID use damage_anisoBrittle, only: & damage_anisoBrittle_LdAndItsTangent + use damage_anisoDuctile, only: & + damage_anisoDuctile_LdAndItsTangent use thermal_adiabatic, only: & thermal_adiabatic_LTAndItsTangent @@ -827,6 +830,11 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el call damage_anisoBrittle_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el) Li = Li + Li_temp dLi_dTstar = dLi_dTstar + dLi_dTstar_temp + + case (LOCAL_DAMAGE_anisoDuctile_ID) + call damage_anisoDuctile_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el) + Li = Li + Li_temp + dLi_dTstar = dLi_dTstar + dLi_dTstar_temp end select @@ -855,9 +863,12 @@ pure function constitutive_getFi(ipc, ip, el) phase_thermal, & material_phase, & LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_THERMAL_adiabatic_ID use damage_anisoBrittle, only: & damage_anisoBrittle_getFd + use damage_anisoDuctile, only: & + damage_anisoDuctile_getFd use thermal_adiabatic, only: & thermal_adiabatic_getFT @@ -874,6 +885,9 @@ pure function constitutive_getFi(ipc, ip, el) select case (phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_anisoBrittle_ID) constitutive_getFi = math_mul33x33(constitutive_getFi,damage_anisoBrittle_getFd (ipc, ip, el)) + + case (LOCAL_DAMAGE_anisoDuctile_ID) + constitutive_getFi = math_mul33x33(constitutive_getFi,damage_anisoDuctile_getFd (ipc, ip, el)) end select @@ -896,9 +910,12 @@ subroutine constitutive_putFi(Tstar_v, Lp, dt, ipc, ip, el) phase_thermal, & material_phase, & LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_THERMAL_adiabatic_ID use damage_anisoBrittle, only: & damage_anisoBrittle_putFd + use damage_anisoDuctile, only: & + damage_anisoDuctile_putFd use thermal_adiabatic, only: & thermal_adiabatic_putFT @@ -917,6 +934,9 @@ subroutine constitutive_putFi(Tstar_v, Lp, dt, ipc, ip, el) select case (phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_anisoBrittle_ID) call damage_anisoBrittle_putFd (Tstar_v, dt, ipc, ip, el) + + case (LOCAL_DAMAGE_anisoDuctile_ID) + call damage_anisoDuctile_putFd (Tstar_v, dt, ipc, ip, el) end select @@ -943,9 +963,12 @@ pure function constitutive_getFi0(ipc, ip, el) phase_thermal, & material_phase, & LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_THERMAL_adiabatic_ID use damage_anisoBrittle, only: & damage_anisoBrittle_getFd0 + use damage_anisoDuctile, only: & + damage_anisoDuctile_getFd0 use thermal_adiabatic, only: & thermal_adiabatic_getFT0 @@ -962,6 +985,9 @@ pure function constitutive_getFi0(ipc, ip, el) select case (phase_damage(material_phase(ipc,ip,el))) case (LOCAL_DAMAGE_anisoBrittle_ID) constitutive_getFi0 = math_mul33x33(constitutive_getFi0,damage_anisoBrittle_getFd0 (ipc, ip, el)) + + case (LOCAL_DAMAGE_anisoDuctile_ID) + constitutive_getFi0 = math_mul33x33(constitutive_getFi0,damage_anisoDuctile_getFd0 (ipc, ip, el)) end select @@ -988,9 +1014,12 @@ pure function constitutive_getPartionedFi0(ipc, ip, el) phase_thermal, & material_phase, & LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_DAMAGE_anisoDuctile_ID, & LOCAL_THERMAL_adiabatic_ID use damage_anisoBrittle, only: & damage_anisoBrittle_getPartionedFd0 + use damage_anisoDuctile, only: & + damage_anisoDuctile_getPartionedFd0 use thermal_adiabatic, only: & thermal_adiabatic_getPartionedFT0 @@ -1008,6 +1037,10 @@ pure function constitutive_getPartionedFi0(ipc, ip, el) case (LOCAL_DAMAGE_anisoBrittle_ID) constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, & damage_anisoBrittle_getPartionedFd0(ipc, ip, el)) + + case (LOCAL_DAMAGE_anisoDuctile_ID) + constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, & + damage_anisoDuctile_getPartionedFd0(ipc, ip, el)) end select diff --git a/code/damage_anisoDuctile.f90 b/code/damage_anisoDuctile.f90 index b48645450..cb4fa0d89 100644 --- a/code/damage_anisoDuctile.f90 +++ b/code/damage_anisoDuctile.f90 @@ -37,10 +37,17 @@ module damage_anisoDuctile real(pReal), dimension(:,:), allocatable, private :: & damage_anisoDuctile_critPlasticStrain + real(pReal), dimension(:), allocatable, private :: & + damage_anisoDuctile_sdot_0, & + damage_anisoDuctile_N + + real(pReal), dimension(:,:), allocatable, private :: & + damage_anisoDuctile_critLoad + enum, bind(c) enumerator :: undefined_ID, & local_damage_ID - end enum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo + end enum integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & damage_anisoDuctile_outputID !< ID of each post result output @@ -51,6 +58,11 @@ module damage_anisoDuctile damage_anisoDuctile_stateInit, & damage_anisoDuctile_aTolState, & damage_anisoDuctile_microstructure, & + damage_anisoDuctile_LdAndItsTangent, & + damage_anisoDuctile_getFd, & + damage_anisoDuctile_putFd, & + damage_anisoDuctile_getFd0, & + damage_anisoDuctile_getPartionedFd0, & damage_anisoDuctile_getDamage, & damage_anisoDuctile_putLocalDamage, & damage_anisoDuctile_getLocalDamage, & @@ -131,9 +143,12 @@ subroutine damage_anisoDuctile_init(fileUnit) damage_anisoDuctile_output = '' allocate(damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) allocate(damage_anisoDuctile_Noutput(maxNinstance), source=0_pInt) + allocate(damage_anisoDuctile_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) allocate(damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) allocate(damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) allocate(damage_anisoDuctile_totalNslip(maxNinstance), source=0_pInt) + allocate(damage_anisoDuctile_N(maxNinstance), source=0.0_pReal) + allocate(damage_anisoDuctile_sdot_0(maxNinstance), source=0.0_pReal) allocate(damage_anisoDuctile_aTol_damage(maxNinstance), source=0.0_pReal) rewind(fileUnit) @@ -176,11 +191,22 @@ subroutine damage_anisoDuctile_init(fileUnit) damage_anisoDuctile_Nslip(j,instance) = IO_intValue(line,positions,1_pInt+j) enddo + case ('sdot0') + damage_anisoDuctile_sdot_0(instance) = IO_floatValue(line,positions,2_pInt) + case ('criticalplasticstrain') do j = 1_pInt, Nchunks_SlipFamilies damage_anisoDuctile_critPlasticStrain(j,instance) = IO_floatValue(line,positions,1_pInt+j) enddo + + case ('damageratesensitivity') + damage_anisoDuctile_N(instance) = IO_floatValue(line,positions,2_pInt) + case ('criticalload') + do j = 1_pInt, Nchunks_SlipFamilies + damage_anisoDuctile_critLoad(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo + end select endif; endif enddo parsingFile @@ -196,8 +222,12 @@ subroutine damage_anisoDuctile_init(fileUnit) damage_anisoDuctile_totalNslip(instance) = sum(damage_anisoDuctile_Nslip(:,instance)) if (damage_anisoDuctile_aTol_damage(instance) < 0.0_pReal) & damage_anisoDuctile_aTol_damage(instance) = 1.0e-3_pReal ! default absolute tolerance 1e-3 + if (damage_anisoDuctile_sdot_0(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') if (any(damage_anisoDuctile_critPlasticStrain(:,instance) < 0.0_pReal)) & call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') + if (damage_anisoDuctile_N(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity_damage ('//LOCAL_DAMAGE_anisoDuctile_LABEL//')') endif myPhase enddo sanityChecks @@ -223,9 +253,9 @@ subroutine damage_anisoDuctile_init(fileUnit) ! Determine size of state array sizeDotState = 0_pInt sizeState = sizeDotState + & - 1_pInt + & ! non-local damage - damage_anisoDuctile_totalNslip(instance) - + 1_pInt + & ! time regularised damage + damage_anisoDuctile_totalNslip(instance) + & ! slip system damages + 9 ! Fd damageState(phase)%sizeState = sizeState damageState(phase)%sizeDotState = sizeDotState damageState(phase)%sizePostResults = damage_anisoDuctile_sizePostResults(instance) @@ -247,8 +277,7 @@ subroutine damage_anisoDuctile_init(fileUnit) allocate(damageState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) if (any(numerics_integrator == 5_pInt)) & allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - - call damage_anisoDuctile_stateInit(phase) + call damage_anisoDuctile_stateInit(phase,instance) call damage_anisoDuctile_aTolState(phase,instance) endif @@ -258,19 +287,27 @@ end subroutine damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief sets the relevant state values for a given instance of this damage !-------------------------------------------------------------------------------------------------- -subroutine damage_anisoDuctile_stateInit(phase) +subroutine damage_anisoDuctile_stateInit(phase, instance) use material, only: & damageState - + use math, only: & + math_I3 implicit none - integer(pInt), intent(in) :: phase !< number specifying the phase of the damage + integer(pInt), intent(in) :: phase , instance !< number specifying the phase of the damage real(pReal), dimension(damageState(phase)%sizeState) :: tempState - tempState = 1.0_pReal + tempState(1) = 1.0_pReal + tempState(2 : & + 1 + damage_anisoDuctile_totalNslip(instance)) = 1.0_pReal + tempState(damage_anisoDuctile_totalNslip(instance)+2: & + damage_anisoDuctile_totalNslip(instance)+10) = reshape(math_I3, shape=[9]) + damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) damageState(phase)%state0 = damageState(phase)%state damageState(phase)%partionedState0 = damageState(phase)%state + + end subroutine damage_anisoDuctile_stateInit !-------------------------------------------------------------------------------------------------- @@ -339,7 +376,7 @@ subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el) drivingForce = plasticState(phase)%accumulatedSlip(index,constituent)/damage_anisoDuctile_critPlasticStrain(f,instance) damageState(phase)%state(index+1,constituent) = & min(damageState(phase)%state0(index+1,constituent), & - 1.0_pReal/drivingForce) + 1.0_pReal/drivingForce) ! irreversibility endif index = index + 1_pInt enddo @@ -354,6 +391,338 @@ subroutine damage_anisoDuctile_microstructure(subdt, ipc, ip, el) end subroutine damage_anisoDuctile_microstructure +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +subroutine damage_anisoDuctile_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) + use numerics, only: & + residualStiffness + use lattice, only: & + lattice_maxNslipFamily, & + lattice_NslipSystem, & + lattice_sd, & + lattice_st, & + lattice_sn + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use math, only: & + math_Plain3333to99, & + math_I3, & + math_identity4th, & + math_symmetric33, & + math_Mandel33to6, & + math_tensorproduct + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(9,9) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (2nd-order tensor) + real(pReal), dimension(3,3,3,3) :: & + dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor) + real(pReal), dimension(3,3) :: & + projection_d, projection_t, projection_n !< projection modes 3x3 tensor + real(pReal), dimension(6) :: & + projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector + integer(pInt) :: & + phase, & + constituent, & + instance, & + f, i, index, index_myFamily, k, l, m, n + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + Ld = 0.0_pReal + dLd_dTstar3333 = 0.0_pReal + + + index = 2_pInt + do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family + do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family + + projection_d = math_tensorproduct(lattice_sd(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + projection_t = math_tensorproduct(lattice_st(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + projection_n = math_tensorproduct(lattice_sn(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + + projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3))) + projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3))) + projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3))) + + traction_d = dot_product(Tstar_v,projection_d_v(1:6)) + traction_t = dot_product(Tstar_v,projection_t_v(1:6)) + traction_n = dot_product(Tstar_v,projection_n_v(1:6)) + + traction_crit = damage_anisoDuctile_critLoad(f,instance)* & + (damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage + + udotd = & + sign(1.0_pReal,traction_d)* & + damage_anisoDuctile_sdot_0(instance)* & + (abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance) + if (udotd /= 0.0_pReal) then + Ld = Ld + udotd*projection_d + dudotd_dt = udotd*damage_anisoDuctile_N(instance)/traction_d + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudotd_dt*projection_d(k,l)* projection_d(m,n) + endif + + udott = & + sign(1.0_pReal,traction_t)* & + damage_anisoDuctile_sdot_0(instance)* & + (abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance) + if (udott /= 0.0_pReal) then + Ld = Ld + udott*projection_t + dudott_dt = udott*damage_anisoDuctile_N(instance)/traction_t + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudott_dt*projection_t(k,l)*projection_t(m,n) + endif + udotn = & + damage_anisoDuctile_sdot_0(instance)* & + (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance) + if (udotn /= 0.0_pReal) then + Ld = Ld + udotn*projection_n(1:3,1:3) + dudotn_dt = udotn*damage_anisoDuctile_N(instance)/traction_n + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudotn_dt*projection_n(k,l)* projection_n(m,n) + endif + index = index + 1_pInt + enddo + enddo + + dLd_dTstar = math_Plain3333to99(dLd_dTstar3333) + +end subroutine damage_anisoDuctile_LdAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoDuctile_getFd(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoDuctile_getFd + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoDuctile_getFd = & + reshape(damageState(phase)% & + state(damage_anisoDuctile_totalNslip(instance)+2: & + damage_anisoDuctile_totalNslip(instance)+10,constituent), & + shape=[3,3]) + +end function damage_anisoDuctile_getFd + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine damage_anisoDuctile_putFd(Tstar_v, dt, ipc, ip, el) + use numerics, only: & + residualStiffness + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use math, only: & + math_mul33x33, & + math_inv33, & + math_Plain3333to99, & + math_I3, & + math_identity4th, & + math_symmetric33, & + math_Mandel33to6, & + math_tensorproduct + use lattice, only: & + lattice_maxNslipFamily, & + lattice_NslipSystem, & + lattice_sd, & + lattice_st, & + lattice_sn + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in) :: & + dt + real(pReal), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), dimension(3,3) :: & + projection_d, projection_t, projection_n !< projection modes 3x3 tensor + real(pReal), dimension(6) :: & + projection_d_v, projection_t_v, projection_n_v !< projection modes 3x3 vector + integer(pInt) :: & + phase, & + constituent, & + instance, & + f, i, index, index_myFamily + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, udott, udotn + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + index = 2_pInt + Ld = 0.0_pReal + do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family + do i = 1_pInt,damage_anisoDuctile_Nslip(f,instance) ! process each (active) slip system in family + + projection_d = math_tensorproduct(lattice_sd(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + projection_t = math_tensorproduct(lattice_st(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + projection_n = math_tensorproduct(lattice_sn(1:3,index_myFamily+i,phase),& + lattice_sn(1:3,index_myFamily+i,phase)) + + projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3))) + projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3))) + projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3))) + + traction_d = dot_product(Tstar_v,projection_d_v(1:6)) + traction_t = dot_product(Tstar_v,projection_t_v(1:6)) + traction_n = dot_product(Tstar_v,projection_n_v(1:6)) + + traction_crit = damage_anisoDuctile_critLoad(f,instance)* & + (damageState(phase)%state0(index,constituent) + residualStiffness) ! degrading critical load carrying capacity by damage + + udotd = & + sign(1.0_pReal,traction_d)* & + damage_anisoDuctile_sdot_0(instance)* & + (abs(traction_d)/traction_crit)**damage_anisoDuctile_N(instance) + if (udotd /= 0.0_pReal) then + Ld = Ld + udotd*projection_d + endif + + udott = & + sign(1.0_pReal,traction_t)* & + damage_anisoDuctile_sdot_0(instance)* & + (abs(traction_t)/traction_crit)**damage_anisoDuctile_N(instance) + if (udott /= 0.0_pReal) then + Ld = Ld + udott*projection_t + endif + + udotn = & + damage_anisoDuctile_sdot_0(instance)* & + (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoDuctile_N(instance) + if (udotn /= 0.0_pReal) then + Ld = Ld + udotn*projection_n(1:3,1:3) + endif + + index = index + 1_pInt + enddo + enddo + + damageState(phase)%state(damage_anisoDuctile_totalNslip(instance)+2: & + damage_anisoDuctile_totalNslip(instance)+10,constituent) = & + reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), & + damage_anisoDuctile_getFd0(ipc, ip, el)), shape=[9]) + +end subroutine damage_anisoDuctile_putFd + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoDuctile_getFd0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoDuctile_getFd0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoDuctile_getFd0 = & + reshape(damageState(phase)% & + subState0(damage_anisoDuctile_totalNslip(instance)+2: & + damage_anisoDuctile_totalNslip(instance)+10,constituent), & + shape=[3,3]) + +end function damage_anisoDuctile_getFd0 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoDuctile_getPartionedFd0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoDuctile_getPartionedFd0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoDuctile_getPartionedFd0 = & + reshape(damageState(phase)% & + partionedState0(damage_anisoDuctile_totalNslip(instance)+2: & + damage_anisoDuctile_totalNslip(instance)+10,constituent), & + shape=[3,3]) + +end function damage_anisoDuctile_getPartionedFd0 + !-------------------------------------------------------------------------------------------------- !> @brief returns damage !-------------------------------------------------------------------------------------------------- @@ -427,7 +796,7 @@ function damage_anisoDuctile_getLocalDamage(ipc, ip, el) damage_anisoDuctile_getLocalDamage = & damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)) - + end function damage_anisoDuctile_getLocalDamage !--------------------------------------------------------------------------------------------------