From 67c3eb3a2ef2e7a7b92b3744dba626fb9a9c53f2 Mon Sep 17 00:00:00 2001
From: Pratheek	Shanthraj <p.shanthraj@mpie.de>
Date: Mon, 8 Dec 2014 14:05:38 +0000
Subject: [PATCH] bounding the analytically integrated damage variable

---
 code/damage_isoBrittle.f90 |  9 ++++++---
 code/damage_phaseField.f90 | 16 +++++++++-------
 2 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/code/damage_isoBrittle.f90 b/code/damage_isoBrittle.f90
index 52687a12e..9c13ab338 100644
--- a/code/damage_isoBrittle.f90
+++ b/code/damage_isoBrittle.f90
@@ -266,6 +266,8 @@ end subroutine damage_isoBrittle_aTolState
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
 subroutine damage_isoBrittle_microstructure(C, Fe, subdt, ipc, ip, el)
+ use numerics, only: &
+   residualStiffness
  use material, only: &
    mappingConstitutive, &
    phase_damageInstance, &
@@ -303,13 +305,14 @@ subroutine damage_isoBrittle_microstructure(C, Fe, subdt, ipc, ip, el)
  stress = math_mul66x6(C,strain) 
  
  damageState(phase)%state(2,constituent) = &
-   min(damageState(phase)%state0(2,constituent), &
-       damage_isoBrittle_critStrainEnergy(instance)/sum(abs(stress*strain)))
+   max(residualStiffness, &
+       min(damageState(phase)%state0(2,constituent), &
+           damage_isoBrittle_critStrainEnergy(instance)/(2.0_pReal*sum(abs(stress*strain)))))                   !< residualStiffness < damage < damage0
        
  damageState(phase)%state(1,constituent) = &
    damageState(phase)%state(2,constituent) + &
    (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* &
-   exp(-subdt/lattice_DamageMobility(phase))
+   exp(-subdt/(damageState(phase)%state(2,constituent)*lattice_DamageMobility(phase)))
               
 end subroutine damage_isoBrittle_microstructure
  
diff --git a/code/damage_phaseField.f90 b/code/damage_phaseField.f90
index 60bd39d56..d7a565a56 100644
--- a/code/damage_phaseField.f90
+++ b/code/damage_phaseField.f90
@@ -271,6 +271,8 @@ end subroutine damage_phaseField_aTolState
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
 subroutine damage_phaseField_microstructure(C, Fe, Cv, subdt, ipc, ip, el)
+ use numerics, only: &
+   residualStiffness
  use material, only: &
    mappingConstitutive, &
    phase_damageInstance, &
@@ -311,14 +313,15 @@ subroutine damage_phaseField_microstructure(C, Fe, Cv, subdt, ipc, ip, el)
  stress = math_mul66x6(C,strain) 
  
  damageState(phase)%state(2,constituent) = &
-   min(damageState(phase)%state0(2,constituent), &
-       damage_phaseField_surfaceEnergy(instance)/ &
-       (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance))))
+   max(residualStiffness, &
+       min(damageState(phase)%state0(2,constituent), &
+           (1.0_pReal - Cv)*damage_phaseField_surfaceEnergy(instance)/ &
+           (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance)))))
 
  damageState(phase)%state(1,constituent) = &
-   damageState(phase)%state(2,constituent)*(1.0_pReal - Cv) + &
-   (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent)*(1.0_pReal - Cv))* &
-   exp(-subdt/lattice_DamageMobility(phase))
+   damageState(phase)%state(2,constituent) + &
+   (damageState(phase)%subState0(1,constituent) - damageState(phase)%state(2,constituent))* &
+   exp(-subdt/(damageState(phase)%state(2,constituent)*lattice_DamageMobility(phase)))
  
 end subroutine damage_phaseField_microstructure
  
@@ -419,7 +422,6 @@ function damage_phaseField_getDamageDiffusion33(ipc, ip, el)
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  damage_phaseField_getDamageDiffusion33 = &
-   damageState(phase)%state(2,constituent)* &
    lattice_DamageDiffusion33(1:3,1:3,phase)
     
 end function damage_phaseField_getDamageDiffusion33