diff --git a/code/constitutive.f90 b/code/constitutive.f90
index 09a5093bb..beb2c4605 100644
--- a/code/constitutive.f90
+++ b/code/constitutive.f90
@@ -29,6 +29,7 @@ module constitutive
    constitutive_LpAndItsTangent, &
    constitutive_LiAndItsTangent, &
    constitutive_getFi, &
+   constitutive_getUndamagedFi, &
    constitutive_putFi, &
    constitutive_getFi0, &
    constitutive_getPartionedFi0, &
@@ -519,10 +520,13 @@ function constitutive_damagedC(ipc,ip,el)
  use material, only: &
    material_phase, &
    LOCAL_DAMAGE_isoBrittle_ID, &
+   LOCAL_DAMAGE_isoDuctile_ID, &
    LOCAL_DAMAGE_phaseField_ID, &   
    phase_damage
  use damage_isoBrittle, only: &
-   damage_isoBrittle_getDamagedC66  
+   damage_isoBrittle_getDamagedC66
+ use damage_isoDuctile, only: &
+   damage_isoDuctile_getDamagedC66  
  use damage_phaseField, only: &
    damage_phaseField_getDamagedC66  
 
@@ -537,6 +541,11 @@ function constitutive_damagedC(ipc,ip,el)
    case (LOCAL_DAMAGE_isoBrittle_ID)
      constitutive_damagedC = damage_isoBrittle_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
                                                              ipc,ip,el)
+                                                             
+   case (LOCAL_DAMAGE_isoDuctile_ID)
+     constitutive_damagedC = damage_isoDuctile_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
+                                                             ipc,ip,el)
+                                                             
    case (LOCAL_DAMAGE_phaseField_ID)
      constitutive_damagedC = damage_phaseField_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
                                                              ipc,ip,el)
@@ -847,6 +856,41 @@ pure function constitutive_getFi(ipc, ip, el)
  
 end function constitutive_getFi
 
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the undamaged 
+!>  intermediate deformation gradient  
+!--------------------------------------------------------------------------------------------------
+pure function constitutive_getUndamagedFi(ipc, ip, el)
+ use prec, only: &
+   pReal 
+ use math, only: &
+   math_I3, &
+   math_mul33x33
+ use material, only: &
+   phase_thermal, &
+   material_phase, &
+   LOCAL_THERMAL_adiabatic_ID
+ use thermal_adiabatic, only: &
+   thermal_adiabatic_getFT
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   dimension(3,3) :: &
+   constitutive_getUndamagedFi                                                                              !< intermediate deformation gradient
+
+ constitutive_getUndamagedFi = math_I3
+
+ select case (phase_thermal(material_phase(ipc,ip,el)))
+   case (LOCAL_THERMAL_adiabatic_ID)
+     constitutive_getUndamagedFi = math_mul33x33(constitutive_getUndamagedFi,thermal_adiabatic_getFT (ipc, ip, el))
+
+ end select
+ 
+end function constitutive_getUndamagedFi
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the intermediate deformation gradient  
diff --git a/code/damage_isoDuctile.f90 b/code/damage_isoDuctile.f90
index 124a0447c..e483d64bd 100644
--- a/code/damage_isoDuctile.f90
+++ b/code/damage_isoDuctile.f90
@@ -47,6 +47,7 @@ module damage_isoDuctile
    damage_isoDuctile_getSlipDamage, &
    damage_isoDuctile_putLocalDamage, &
    damage_isoDuctile_getLocalDamage, &
+   damage_isoDuctile_getDamagedC66, &
    damage_isoDuctile_postResults
 
 contains
@@ -330,7 +331,7 @@ function damage_isoDuctile_getDamage(ipc, ip, el)
     
    case (FIELD_DAMAGE_NONLOCAL_ID)
     damage_isoDuctile_getDamage =    fieldDamage(material_homog(ip,el))% &
-      field(1,mappingHomogenization(1,ip,el))                                                     ! Taylor type 
+      field(1,mappingHomogenization(1,ip,el))                                                       ! Taylor type 
 
  end select
  
@@ -392,6 +393,35 @@ function damage_isoDuctile_getLocalDamage(ipc, ip, el)
    damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
  
 end function damage_isoDuctile_getLocalDamage
+!--------------------------------------------------------------------------------------------------
+!> @brief returns ductile damaged stiffness tensor 
+!--------------------------------------------------------------------------------------------------
+function damage_isoDuctile_getDamagedC66(C, ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6,6) :: &
+   C
+ real(pReal),              dimension(6,6) :: &
+   damage_isoDuctile_getDamagedC66
+ integer(pInt) :: &
+   phase, constituent
+ real(pReal) :: &
+   damage
+ 
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage = damage_isoDuctile_getDamage(ipc, ip, el)
+ damage_isoDuctile_getDamagedC66 = &
+   damage*damage*C
+    
+end function damage_isoDuctile_getDamagedC66
 
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results