From 76442942b1b6ab9c37d923fdc4bc6ab900360ed7 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Mon, 3 Nov 2014 10:43:36 +0000
Subject: [PATCH] damage def grad updated in stress integration loop and not as
a dot state
---
code/constitutive.f90 | 36 +++++++-
code/crystallite.f90 | 47 +++++-----
code/damage_anisoBrittle.f90 | 162 +++++++++++++++++++++++++----------
3 files changed, 176 insertions(+), 69 deletions(-)
diff --git a/code/constitutive.f90 b/code/constitutive.f90
index 5a9ab361a..0cc0e0f13 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -29,6 +29,7 @@ module constitutive
constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, &
constitutive_getFi, &
+ constitutive_putFi, &
constitutive_getFi0, &
constitutive_getPartionedFi0, &
constitutive_TandItsTangent, &
@@ -539,7 +540,6 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
PLASTICITY_nonlocal_ID, &
LOCAL_DAMAGE_isoBrittle_ID, &
LOCAL_DAMAGE_isoDuctile_ID, &
- LOCAL_DAMAGE_anisoBrittle_ID, &
LOCAL_DAMAGE_gurson_ID
use constitutive_titanmod, only: &
@@ -563,7 +563,7 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
ipc, & !< grain number
ip, & !< integration point number
el !< element number
- real(pReal), intent(in), dimension(6) :: &
+ real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
@@ -805,6 +805,38 @@ pure function constitutive_getFi(ipc, ip, el)
end function constitutive_getFi
+!--------------------------------------------------------------------------------------------------
+!> @brief contains the constitutive equation for calculating the intermediate deformation gradient
+!--------------------------------------------------------------------------------------------------
+subroutine constitutive_putFi(Tstar_v, dt, ipc, ip, el)
+ use prec, only: &
+ pReal
+ use material, only: &
+ phase_damage, &
+ material_phase, &
+ LOCAL_DAMAGE_anisoBrittle_ID
+ use damage_anisoBrittle, only: &
+ damage_anisoBrittle_putFd
+
+ 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
+
+ select case (phase_damage(material_phase(ipc,ip,el)))
+ case (LOCAL_DAMAGE_anisoBrittle_ID)
+ call damage_anisoBrittle_putFd (Tstar_v, dt, ipc, ip, el)
+
+ end select
+
+end subroutine constitutive_putFi
+
+
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the intermediate deformation gradient
!--------------------------------------------------------------------------------------------------
diff --git a/code/crystallite.f90 b/code/crystallite.f90
index d3c8e5ddf..33e488547 100644
--- a/code/crystallite.f90
+++ b/code/crystallite.f90
@@ -585,10 +585,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
convergenceFlag_backup
! local variables used for calculating analytic Jacobian
real(pReal) :: detInvFi
- real(pReal), dimension(3,3) :: junk, &
+ real(pReal), dimension(3,3) :: temp_33, &
Fi, &
invFi, &
- Fi0, &
invFi0
real(pReal), dimension(3,3,3,3) :: dSdFe, &
dSdF, &
@@ -599,7 +598,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
rhs_3333, &
lhs_3333, &
temp_3333
- real(pReal), dimension(9,9):: temp_99
+ real(pReal), dimension(9,9):: temp_99
logical :: error
@@ -1104,39 +1103,38 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --- ANALYTIC JACOBIAN ---
- !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_99,junk,dLidS,Fi,invFi,Fi0,invFi0,detInvFi,temp_3333,myNgrains)
+ !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_99,temp_33,dLidS,Fi,invFi,invFi0,detInvFi,temp_3333,myNgrains)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1_pInt,myNgrains
- call constitutive_TandItsTangent(junk,dSdFe,crystallite_Fe(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
- call constitutive_LpAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
+ call constitutive_TandItsTangent(temp_33,dSdFe,crystallite_Fe(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
+ call constitutive_LpAndItsTangent(temp_33,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
g,i,e)
dLpdS = reshape(temp_99,shape=[3,3,3,3])
- call constitutive_LiAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
+ call constitutive_LiAndItsTangent(temp_33,temp_99,crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Lp(1:3,1:3,g,i,e), g,i,e)
dLidS = reshape(temp_99,shape=[3,3,3,3])
Fi = constitutive_getFi(g,i,e)
invFi = math_inv33(Fi)
- Fi0 = constitutive_getFi0(g,i,e)
- invFi0 = math_inv33(Fi0)
+ invFi0 = math_inv33(constitutive_getFi0(g,i,e))
detInvFi = math_det33(invFi)
- junk = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),invFi))
+ temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),invFi))
rhs_3333 = 0.0_pReal
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
- rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),junk)
+ rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33)
temp_3333 = 0.0_pReal
- junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
+ temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
- temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(junk,dLpdS(1:3,1:3,p,o)),invFi)
+ temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)),invFi)
- junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
+ temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e)), invFi0)
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
- temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(junk,dLidS(1:3,1:3,p,o))
+ temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o))
lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333)
@@ -1152,31 +1150,31 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
math_mul33x33(temp_3333(1:3,1:3,p,o),invFi))
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal
- junk = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), &
+ temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), &
math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))))/detInvFi
forall(p=1_pInt:3_pInt) &
- crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(junk)
+ crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(temp_33)
- junk = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
+ temp_33 = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), &
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))/detInvFi
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
- math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),dFpinvdF(1:3,1:3,p,o)),junk)
+ math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33)
- junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
+ temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e))/detInvFi
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
- math_mul33x33(math_mul33x33(junk,dSdF(1:3,1:3,p,o)), &
+ math_mul33x33(math_mul33x33(temp_33,dSdF(1:3,1:3,p,o)), &
math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))
- junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
+ temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), &
crystallite_invFp(1:3,1:3,g,i,e)), &
math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)))/detInvFi
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) &
crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + &
- math_mul33x33(junk,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
+ math_mul33x33(temp_33,math_transpose33(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
enddo elementLooping6
@@ -3345,6 +3343,7 @@ logical function crystallite_integrateStress(&
use constitutive, only: constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, &
constitutive_getFi0, &
+ constitutive_putFi, &
constitutive_TandItsTangent
use math, only: math_mul33x33, &
math_mul33xx33, &
@@ -3769,7 +3768,7 @@ logical function crystallite_integrateStress(&
crystallite_Fp(1:3,1:3,g,i,e) = Fp_new
crystallite_Fe(1:3,1:3,g,i,e) = Fe_new
crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
-
+ call constitutive_putFi(Tstar_v, dt, g, i, e)
!* set return flag to true
diff --git a/code/damage_anisoBrittle.f90 b/code/damage_anisoBrittle.f90
index f8494e95f..1f12318d0 100644
--- a/code/damage_anisoBrittle.f90
+++ b/code/damage_anisoBrittle.f90
@@ -57,6 +57,7 @@ module damage_anisoBrittle
damage_anisoBrittle_dotState, &
damage_anisoBrittle_LdAndItsTangent, &
damage_anisoBrittle_getFd, &
+ damage_anisoBrittle_putFd, &
damage_anisoBrittle_getFd0, &
damage_anisoBrittle_getPartionedFd0, &
damage_anisoBrittle_getDamage, &
@@ -234,9 +235,9 @@ subroutine damage_anisoBrittle_init(fileUnit)
enddo outputsLoop
! Determine size of state array
sizeDotState = 2_pInt + & ! viscous and non-viscous damage values
- 3_pInt*damage_anisoBrittle_totalNcleavage(instance) + & ! opening on each damage system
+ 3_pInt*damage_anisoBrittle_totalNcleavage(instance) ! opening on each damage system
+ sizeState = sizeDotState + &
9_pInt ! Fd
- sizeState = sizeDotState
damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState
@@ -285,7 +286,8 @@ subroutine damage_anisoBrittle_stateInit(phase,instance)
tempState(2) = 1.0_pReal
tempState(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal
tempState(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11) = reshape(math_I3, shape=[9])
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11) = &
+ 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
@@ -306,9 +308,11 @@ subroutine damage_anisoBrittle_aTolState(phase,instance)
tempTol(1) = damage_anisoBrittle_aTol_damage(instance)
tempTol(2) = damage_anisoBrittle_aTol_damage(instance)
- tempTol(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp(instance)
+ tempTol(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = &
+ damage_anisoBrittle_aTol_disp(instance)
tempTol(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11) = damage_anisoBrittle_aTol_disp(instance)
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11) = &
+ damage_anisoBrittle_aTol_disp(instance)
damageState(phase)%aTolState = tempTol
end subroutine damage_anisoBrittle_aTolState
@@ -321,7 +325,6 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
phase_damageInstance, &
damageState
use math, only: &
- math_mul33x33, &
math_norm3
use lattice, only: &
lattice_Scleavage, &
@@ -337,8 +340,6 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
- real(pReal), dimension(3,3) :: &
- Ld, Fd
integer(pInt) :: &
phase, &
constituent, &
@@ -354,24 +355,23 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
instance = phase_damageInstance(phase)
damageState(phase)%dotState(1,constituent) = &
- (1.0_pReal/lattice_DamageMobility(phase))* &
- (damageState(phase)%state(2,constituent) - &
- damageState(phase)%state(1,constituent))
+ (max(0.0_pReal,damageState(phase)%state(2,constituent)) - &
+ damageState(phase)%state(1,constituent))/lattice_DamageMobility(phase)
damageState(phase)%dotState(2,constituent) = 0.0_pReal
nonLocalFactor = 1.0_pReal + &
(damage_anisoBrittle_getDamage (ipc, ip, el) - &
damage_anisoBrittle_getLocalDamage(ipc, ip, el))
- Ld = 0.0_pReal
index = 2_pInt
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
- opening = math_norm3(damageState(phase)%state(index+1:index+3,constituent))
+ opening = max(math_norm3(damageState(phase)%state0(index+1:index+3,constituent)), &
+ math_norm3(damageState(phase)%state (index+1:index+3,constituent)))
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
- traction_crit = nonLocalFactor*damage_anisoBrittle_critLoad(f,instance)* &
- 1.0_pReal/(1.0_pReal + opening/damage_anisoBrittle_critDisp(f,instance))
+ traction_crit = nonLocalFactor*damage_anisoBrittle_critLoad(f,instance)/ &
+ (1.0_pReal + opening/damage_anisoBrittle_critDisp(f,instance))
damageState(phase)%dotState(index+1,constituent) = &
sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* &
@@ -381,34 +381,20 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el)
damage_anisoBrittle_sdot_0(instance)* &
(abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance)
damageState(phase)%dotState(index+3,constituent) = &
- sign(1.0_pReal,traction_n)* &
damage_anisoBrittle_sdot_0(instance)* &
- (abs(traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
+ (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
damageState(phase)%dotState(2,constituent) = &
damageState(phase)%dotState(2,constituent) - &
(traction_d*damageState(phase)%dotState(index+1,constituent) + &
traction_t*damageState(phase)%dotState(index+2,constituent) + &
traction_n*damageState(phase)%dotState(index+3,constituent))/ &
(0.5_pReal*damage_anisoBrittle_critDisp(f,instance)*damage_anisoBrittle_critLoad(f,instance))
- Ld = Ld + &
- damageState(phase)%dotState(index+1,constituent)* &
- lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) + &
- damageState(phase)%dotState(index+2,constituent)* &
- lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) + &
- damageState(phase)%dotState(index+3,constituent)* &
- lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
+
index = index + 3_pInt
enddo
enddo
- Fd = reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
- shape=[3,3])
- damageState(phase)%dotState(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent) = &
- reshape(math_mul33x33(Ld,Fd), shape=[9])
end subroutine damage_anisoBrittle_dotState
-
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
@@ -464,12 +450,13 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
do f = 1_pInt,lattice_maxNcleavageFamily
index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
- opening = math_norm3(damageState(phase)%state(index+1:index+3,constituent))
+ opening = max(math_norm3(damageState(phase)%state0(index+1:index+3,constituent)), &
+ math_norm3(damageState(phase)%state (index+1:index+3,constituent)))
traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
- traction_crit = nonLocalFactor*damage_anisoBrittle_critLoad(f,instance)* &
- 1.0_pReal/(1.0_pReal + opening/damage_anisoBrittle_critDisp(f,instance))
+ traction_crit = nonLocalFactor*damage_anisoBrittle_critLoad(f,instance)/ &
+ (1.0_pReal + opening/damage_anisoBrittle_critDisp(f,instance))
udotd = &
sign(1.0_pReal,traction_d)* &
damage_anisoBrittle_sdot_0(instance)* &
@@ -497,9 +484,8 @@ subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip,
endif
udotn = &
- sign(1.0_pReal,traction_n)* &
damage_anisoBrittle_sdot_0(instance)* &
- (abs(traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
+ (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
if (udotn /= 0.0_pReal) then
Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
dudotn_dt = udotn*damage_anisoBrittle_N(instance)/traction_n
@@ -542,12 +528,100 @@ pure function damage_anisoBrittle_getFd(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getFd = &
- reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
+ reshape(damageState(phase)% &
+ state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
shape=[3,3])
end function damage_anisoBrittle_getFd
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine damage_anisoBrittle_putFd(Tstar_v, dt, ipc, ip, el)
+ use material, only: &
+ mappingConstitutive, &
+ phase_damageInstance, &
+ damageState
+ use math, only: &
+ math_mul33x33, &
+ math_norm3, &
+ math_inv33, &
+ math_I3
+ use lattice, only: &
+ lattice_Scleavage, &
+ lattice_Scleavage_v, &
+ lattice_maxNcleavageFamily, &
+ lattice_NcleavageSystem
+
+ 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
+ integer(pInt) :: &
+ phase, &
+ constituent, &
+ instance, &
+ f, i, index, index_myFamily
+ real(pReal) :: &
+ traction_d, traction_t, traction_n, traction_crit, &
+ udotd, udott, udotn, &
+ opening, &
+ nonLocalFactor
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_damageInstance(phase)
+
+ Ld = 0.0_pReal
+ nonLocalFactor = 1.0_pReal + &
+ (damage_anisoBrittle_getDamage (ipc, ip, el) - &
+ damage_anisoBrittle_getLocalDamage(ipc, ip, el))
+ index = 2_pInt
+ do f = 1_pInt,lattice_maxNcleavageFamily
+ index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family
+ do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family
+ opening = max(math_norm3(damageState(phase)%state0(index+1:index+3,constituent)), &
+ math_norm3(damageState(phase)%state (index+1:index+3,constituent)))
+ traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
+ traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
+ traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
+ traction_crit = nonLocalFactor*damage_anisoBrittle_critLoad(f,instance)/ &
+ (1.0_pReal + opening/damage_anisoBrittle_critDisp(f,instance))
+ udotd = &
+ sign(1.0_pReal,traction_d)* &
+ damage_anisoBrittle_sdot_0(instance)* &
+ (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance)
+ Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
+
+ udott = &
+ sign(1.0_pReal,traction_t)* &
+ damage_anisoBrittle_sdot_0(instance)* &
+ (abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance)
+ Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
+
+ udotn = &
+ damage_anisoBrittle_sdot_0(instance)* &
+ (max(0.0_pReal,traction_n)/traction_crit)**damage_anisoBrittle_N(instance)
+ Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
+
+ index = index + 3_pInt
+ enddo
+ enddo
+ damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent) = &
+ reshape(math_mul33x33(math_inv33(math_I3 - dt*Ld), &
+ damage_anisoBrittle_getFd0(ipc, ip, el)), shape=[9])
+
+end subroutine damage_anisoBrittle_putFd
+
!--------------------------------------------------------------------------------------------------
!> @brief returns local damage deformation gradient
!--------------------------------------------------------------------------------------------------
@@ -574,8 +648,9 @@ pure function damage_anisoBrittle_getFd0(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getFd0 = &
- reshape(damageState(phase)%subState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
+ reshape(damageState(phase)% &
+ subState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
shape=[3,3])
end function damage_anisoBrittle_getFd0
@@ -606,8 +681,9 @@ pure function damage_anisoBrittle_getPartionedFd0(ipc, ip, el)
instance = phase_damageInstance(phase)
damage_anisoBrittle_getPartionedFd0 = &
- reshape(damageState(phase)%partionedState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
- 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
+ reshape(damageState(phase)% &
+ partionedState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: &
+ 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), &
shape=[3,3])
end function damage_anisoBrittle_getPartionedFd0
@@ -683,7 +759,7 @@ function damage_anisoBrittle_getLocalDamage(ipc, ip, el)
damage_anisoBrittle_getLocalDamage
damage_anisoBrittle_getLocalDamage = &
- max(0.0_pReal,damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))
+ damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
end function damage_anisoBrittle_getLocalDamage