Isotropic Damage Model

In this section, the implementation of ann isotropic damage model will be described. To cover the various models based on isotropic damage concept, a base class IsotropicDamageMaterial is defined first, declaring the necessary services and providing the implementation of them, which are general. The derived classes then only implement a particular damage-evolution law.

The isotropic damage models are based on the simplifying assumption that the stiffness degradation is isotropic, i.e., stiffness moduli corresponding to different directions decrease proportionally and independently of direction of loading. Consequently, the damaged stiffness matrix is expressed as

$$\mbox{\boldmath $D$} = (1-\omega)\mbox{\boldmath $D$}_e,$$

where $\mbox{\boldmath$D$}_e$ is elastic stiffness matrix of the undamaged material and $\omega$ is the damage parameter. Initially, $\omega$ is set to zero, representing the virgin undamaged material, and the response is linear-elastic. As the material undergoes the deformation, the initiation and propagation of microdefects decreases the stiffness, which is represented by the growth of the damage parameter $\omega$. For $\omega = 1$, the stiffness completely disappears.

In the present context, the $\mbox{\boldmath$D$}$ matrix represents the secant stiffness that relates the total strain to the total stress

$$\mbox{\boldmath $\sigma$}=\mbox{\boldmath $D$}\mbox{\boldma... ...-\omega)\mbox{\boldmath $D$}_e\mbox{\boldmath $\varepsilon$}.$$

Similarly to the theory of plasticity, a loading function $f$ is introduced. In the damage theory, it is natural to work in the strain space and therefore the loading function is depending on the strain and on an additional parameter $\kappa$, describing the evolution of the damage. Physically, $\kappa$ is a scalar measure of the largest strain level ever reached. The loading function usually has the form

$$f(\mbox{\boldmath $\varepsilon$}, \kappa) = \tilde\varepsilon(\mbox{\boldmath $\varepsilon$}) - \kappa,$$

where $\tilde\varepsilon$ is the equivalent strain, i.e., the scalar measure of the strain level. Damage can grow only if current state reaches the boundary of elastic domain ($f=0$). This is expressed by the following loading/unloading conditions

$$f \le 0,\;\;\dot\kappa \ge0,\;\;\dot\kappa f = 0.$$

It remains to link the variable $\kappa$ to the damage parameter $\omega$. As both $\kappa$ and $\omega$ grow monotonically, it is convenient to postulate an explicit evolution law

$$\omega = g(\kappa).$$

The important advantage of this explicit formulation is that the stress corresponding to the given strain can be evaluated directly, without the need to solve the nonlinear system of equations. For the given strain, the corresponding stress is computed simply by evaluating the current equivalent strain, updating the maximum previously reached equivalent strain value $\kappa$ and the damage parameter and reducing the effective stress according to $\mbox{\boldmath$\sigma$} = (1-\omega)\mbox{\boldmath$D$}_e \mbox{\boldmath$\varepsilon$}$.

This general framework for computing stresses and stiffness matrix is common for all material models of this type. Therefore, it is natural to introduce the base class for all isotropic-based damage models which provides the general implementation for the stress and stiffness matrix evaluation algorithms. The particular models then only provide their equivalent strain and damage evolution law definitions. The base class only declares the virtual services for computing equivalent strain and corresponding damage. The implementation of common services uses these virtual functions, but they are only declared at IsotropicDamageMaterial class level and have to be implemented by the derived classes.

Together with the material model, the corresponding status has to be defined, containing all necessary history variables. For the isotropic-based damage models, the only history variable is the value of the largest strain level ever reached ($\kappa$). In addition, the corresponding damage level $\omega$ will be stored. This is not necessary because damage can be always computed from corresponding $\kappa$. The IsotropicDamageMaterialStatus class is derived from StructuralMaterialStatus class. The base class represents the base material status class for all structural statuses. At StructuralMaterialStatus level, the attributes common to all ``structural analysis'' material models - the strain and stress vectors (both the temporary and non-temporary) are introduced. The corresponding services for accessing, setting, initializing, and updating these attributes are provided. Therefore, only the $\kappa$ and $\omega$ parameters are introduced (both the temporary and non-temporary). The corresponding services for manipulating these attributes are added and services for context initialization, update, and store/restore operations are overloaded, to handle the history parameters properly.

 class IsotropicDamageMaterialStatus : public StructuralMaterialStatus
 {

 protected:
   /// scalar measure of the largest strain level ever reached in material
   double kappa;
   /// non-equilibrated scalar measure of the largest strain level
   double tempKappa;
   /// damage level of material
   double damage;
   /// non-equilibrated damage level of material
   double tempDamage;

 public: 
   /// Constructor
   IsotropicDamageMaterialStatus (int n, Domain*d, GaussPoint* g) ;
   /// Destructor
   ~IsotropicDamageMaterialStatus ();

   /// Prints the receiver state to stream
   void   printOutputAt (FILE *file, TimeStep* tStep) ;

   /// Returns the last equilibrated scalar measure 
   /// of the largest strain level
   double giveKappa () {return kappa;}
   /// Returns the temp. scalar measure of the 
   /// largest strain level
   double giveTempKappa () {return tempKappa;}
   /// Sets the temp scalar measure of the largest 
   /// strain level to given value
   void   setTempKappa (double newKappa) { tempKappa = newKappa;}
   /// Returns the last equilibrated damage level
   double giveDamage () {return damage;}
   /// Returns the temp. damage level
   double giveTempDamage () {return tempDamage;}
   /// Sets the temp damage level to given value
   void   setTempDamage (double newDamage) { tempDamage = newDamage;}


   // definition
   char*                 giveClassName (char* s) const
     { return strcpy(s,"IsotropicDamageMaterialModelStatus") ;}
   classType             giveClassID () const 
     { return IsotropicDamageMaterialStatusClass; }

   /**
     Initializes the temporary internal variables, 
     describing the current state according to 
     previously reached equilibrium internal variables.
   */
   virtual void initTempStatus ();
   /**
     Update equilibrium history variables 
     according to temp-variables.
     Invoked, after new equilibrium state has been reached.
   */
   virtual void updateYourself(TimeStep*); 

   // saves current context(state) into stream
   /**
     Stores context of receiver into given stream. 
     Only non-temp internal history variables are stored.
     @param stream stream where to write data
     @param obj pointer to integration point, which invokes this method
     @return nonzero if o.k.
   */
   contextIOResultType    saveContext (FILE* stream, void *obj = NULL);
   /**
     Restores context of receiver from given stream. 
     @param stream stream where to read data
     @param obj pointer to integration point, which invokes this method
     @return nonzero if o.k.
   */
   contextIOResultType    restoreContext(FILE* stream, void *obj = NULL);


 };

The base IsotropicDamageMaterial class is derived from the StructuralMaterial class.

 /**
   Base class representing general isotropic damage model.
   It is based on isotropic damage concept, 
   assuming that damage evolution law
   is postulated in explicit form, 
   relating damage parameter (omega) to scalar measure 
   of the largest strain level ever reached in material (kappa).
 */
 class IsotropicDamageMaterial : public StructuralMaterial
 {

 protected :
   /// Coefficient of thermal dilatation
   double tempDillatCoeff;

   /// Reference to bulk (undamaged) material
   LinearElasticMaterial *linearElasticMaterial;  

   /**
     variable controling type of loading/unloading law, 
     default set to idm_strainLevel
     defines the two possibilities:
     - idm_strainLevelCR the unloaing takes place, 
       when strain level is smaller than the 
       largest level ever reached;
     - idm_damageLevelCR the unloaing takes place, 
       when damage level is smaller than the 
       largest damage ever  reached;
   */
   enum loaUnloCriterium {
     idm_strainLevelCR, 
     idm_damageLevelCR 
   } llcriteria;

 public :
   /// Constructor
   IsotropicDamageMaterial (int n,Domain* d) ;

   /// Destructor
   ~IsotropicDamageMaterial () ;

   /// Returns nonzero indicating that receiver is nonlinear
   int hasNonLinearBehaviour ()   { return 1 ;}

   /**
     Tests, if material supports material mode.
     @param mode required material mode
     @return nonzero if supported, zero otherwise
   */
   int hasMaterialModeCapability (MaterialMode mode);

   char*    giveClassName (char* s) const 
     { return strcpy(s,"IsotropicDamageMaterial") ;}

   classType giveClassID ()         const 
     { return IsotropicDamageMaterialClass;}

   /// Returns reference to undamaged (bulk) material
   LinearElasticMaterial* giveLinearElasticMaterial () 
     {return linearElasticMaterial;}

   /**
     Computes full 3d material stiffness matrix at given 
     integration point, time, respecting load history 
     in integration point.
     @param answer computed results
     @param form material response form
     @param mode material response mode
     @param gp integration point
     @param atTime time step (most models are able to 
      respond only when atTime is current time step)
   */
   virtual void  give3dMaterialStiffnessMatrix (FloatMatrix& answer,
                                                MatResponseForm form, 
                                                MatResponseMode mode, 
                                                GaussPoint* gp,
                                                TimeStep* atTime);

   /**
     Computes the real stress vector for given strain 
     increment and integration point. 
     Temp vars are updated accordingly
     @param answer contains result
     @param form material response form
     @param gp integration point 
     @param reducedStrain strain  vector in reduced form
     @param tStep current time step (most models are able to 
      respond only when tStep is current time step)
   */
   void giveRealStressVector (FloatArray& answer,  MatResponseForm form, 
                              GaussPoint* gp, const FloatArray &reducedStrain, 
                              TimeStep* tStep);


   /**
     Returns the integration point corresponding 
     value in Reduced form.
     @param answer contain corresponding ip value, 
      zero sized if not available
     @param aGaussPoint integration point
     @param type determines the type of internal variable
     @returns nonzero if ok, zero if var not supported
   */
   virtual int giveIPValue (FloatArray& answer, 
                            GaussPoint* aGaussPoint, 
                            InternalStateType type, 
                            TimeStep* atTime) ;


   /**
     Returns the mask of reduced indexes of 
     Internal Variable component .
     @param answer mask of Full VectorSize, with components 
      being the indexes to reduced form vectors.
     @param type determines the internal variable requested
     @returns nonzero if ok or error is generated for unknown mat mode.
   */
   virtual int giveIntVarCompFullIndx (IntArray& answer, 
                                       InternalStateType type, 
                                       MaterialMode mmode);

   /**
     Returns the type of internal variable (scalar, vector, ...)
     @param type determines the type of internal variable
     @returns type of internal variable
   */
   virtual 
   InternalStateValueType giveIPValueType (InternalStateType type);

   /**
     Returns the corresponding integration point  
     value size in Reduced form.
     @param type determines the type of internal variable
     @returns var size, zero if var not supported
   */
   virtual int giveIPValueSize (InternalStateType type, 
                                GaussPoint* aGaussPoint) ;

   /**
     Returns a vector of coefficients of thermal dilatation in direction
     of each material principal (local) axis.
     @param answer vector of thermal dilatation coefficients
     @param gp integration point
     @param tStep time step (most models are able to 
      respond only when atTime is current time step)
   */
   virtual void giveThermalDilatationVector (FloatArray& answer, 
                                             GaussPoint*, TimeStep*) ;



   /**
     Computes the equivalent strain measure from 
     given strain vector (full form).
     @param kappa return param, comtaining the 
      corresponding equivalent strain
     @param strain total strain vector in full form
     @param gp integration point
     @param atTime timestep
   */
   virtual void computeEquivalentStrain (double& kappa, 
                                         const FloatArray& strain, 
                                         GaussPoint* gp, 
                                         TimeStep* atTime) = 0;

   /**
     computes the value of damage parameter omega, 
     based on given value of equivalent strain
     @param omega contains result
     @param kappa equivalent strain measure
   */
   virtual void computeDamageParam (double& omega, double kappa, 
                                    const FloatArray& strain, 
                                    GaussPoint* gp) =0;

  /**
    Instanciates the receiver from input record
  */
  IRResultType initializeFrom (InputRecord* ir);

 protected: 

   /** 
     Abstract service allowing to perfom some 
     initialization, when damage first appear
     @param kappa scalar measure of strain level
     @param totalStrainVector current total strain vector
     @param gp integration point
   */
   virtual void initDamaged (double kappa, 
                             FloatArray& totalStrainVector, 
                             GaussPoint* gp) {}

   /// Creates corresponding material status 
   MaterialStatus* CreateStatus (GaussPoint *gp) 
     {return new IsotropicDamageMaterialStatus (1,domain,gp);}


   // Overloaded to use specialized versions of these 
   // services possibly implemented by linearElastic member
   /**
     Method for computing plane stress stifness matrix of receiver.
     Default implementation overloaded to use direct implementation of 
     corresponding service at bulk material model level.
     @param answer stifness matrix
     @param form material response form
     @param mode material response mode
     @param gp integration point, which load history is used
     @param atTime time step (most models are able to 
      respond only when atTime is current time step)
   */
   void givePlaneStressStiffMtrx (FloatMatrix& answer, 
                                  MatResponseForm form,
                                  MatResponseMode,
                                  GaussPoint* gp,
                                  TimeStep* atTime);


   /**
     Method for computing plane strain stifness matrix of receiver.
     Default implementation overloaded to use direct implementation of 
     corresponding service at bulk material model level.
     @param answer stifness matrix
     @param form material response form
     @param mode material response mode
     @param gp integration point, which load history is used
     @param atTime time step (most models are able to 
      respond only when atTime is current time step)
   */
   void givePlaneStrainStiffMtrx (FloatMatrix& answer, 
                                  MatResponseForm form,
                                  MatResponseMode,
                                  GaussPoint* gp, 
                                  TimeStep* atTime);


   /**
     Method for computing 1d stifness matrix of receiver.
     Default implementation overloaded to use direct implementation of 
     corresponding service at bulk material model level.
     @param answer stifness matrix
     @param form material response form
     @param mode material response mode
     @param gp integration point, which load history is used
     @param atTime time step (most models are able to 
      respond only when atTime is current time step)
    */
   void give1dStressStiffMtrx (FloatMatrix& answer, 
                               MatResponseForm form,
                               MatResponseMode,
                               GaussPoint* gp,
                               TimeStep* atTime);
 } ;

#include "isodamagemodel.h"
#include "gausspnt.h"
#include "flotmtrx.h"
#include "flotarry.h"
#include "structuralcrosssection.h"
#include "mathfem.h"

IsotropicDamageMaterial :: IsotropicDamageMaterial (int n, Domain *d)
: StructuralMaterial (n,d)
//
// constructor
//
{
  linearElasticMaterial = NULL;
  llcriteria = idm_strainLevelCR;
}


IsotropicDamageMaterial :: ~IsotropicDamageMaterial ()
//
// destructor
//
{
  delete linearElasticMaterial;
}

int
IsotropicDamageMaterial :: hasMaterialModeCapability (MaterialMode mode)
//
// returns whether receiver supports given mode
//
{
  if ((mode == _3dMat) || (mode == _PlaneStress) || 
     (mode == _PlaneStrain) || (mode == _1dMat)) return 1;
  return 0;
}

Next, the implementation for core services is presented. They compute material stiffness matrices for various material modes and compute the stress vector. The later service computes the real stress corresponding to the previous history and the given strain vector. It computes a new locally consistent (equilibrated) state, which is stored in temporary variables of the corresponding status (attribute of the integration point).

void
IsotropicDamageMaterial :: give3dMaterialStiffnessMatrix (FloatMatrix& answer, 
                                                          MatResponseForm form,
                                                          MatResponseMode mode,
                                                          GaussPoint* gp,
                                                          TimeStep* atTime)
//
// computes full constitutive matrix for case of gp stress-strain state.
//
{
  IsotropicDamageMaterialStatus *status =
    (IsotropicDamageMaterialStatus*) this -> giveStatus (gp);
  double om = status->giveTempDamage();
  om = min (om, 0.999999);

  this->giveLinearElasticMaterial()->
        give3dMaterialStiffnessMatrix (answer, form, mode, gp, atTime);
  answer.times (1.0-om);

  return ;
  
}


void
IsotropicDamageMaterial :: giveRealStressVector (FloatArray& answer, 
                                                 MatResponseForm form, 
                                                 GaussPoint* gp, 
                                                 const FloatArray& totalStrain, 
                                                 TimeStep* atTime)
//
// returns real stress vector in 3d stress space of receiver according to 
// previous load history and given strain.
// After the new local point equilibrium is reached all the
// temp history variables are updated to new local equilibrium state.
// 
//
{
  IsotropicDamageMaterialStatus *status = 
    (IsotropicDamageMaterialStatus*) this -> giveStatus (gp);

  LinearElasticMaterial* lmat = this->giveLinearElasticMaterial();
  FloatArray strainVector, reducedTotalStrainVector;
  FloatMatrix de;
  double f, equivStrain, tempKappa, omega; 

  this->initGpForNewStep(gp);

  // substract stress independent part
  // note: eigenStrains (tepmerature) is not contained in mechanical
  // strain stored in gp therefore it is necessary to substract 
  // always the total eigen strain value
  this->giveStressDependentPartOfStrainVector(reducedTotalStrainVector,
                                              gp, totalStrain,atTime, TotalMode);

  // compute equivalent strain
  this->computeEquivalentStrain (equivStrain, reducedTotalStrainVector, gp, atTime);

  if (llcriteria == idm_strainLevelCR) {
    // compute value of loading function if strainLevel crit apply
    f = equivStrain - status->giveKappa();

    if (f <= 0.0) {
      // damage do not grow
      tempKappa = status->giveKappa();
      omega     = status->giveDamage();
    } else {
      // damage grow
      tempKappa = equivStrain;
      this->initDamaged (tempKappa, reducedTotalStrainVector, gp);
      // evaluate damage parameter
      this->computeDamageParam (omega, tempKappa, reducedTotalStrainVector, gp);
    }
  } else if (llcriteria == idm_damageLevelCR) {

    // evaluate damage parameter first
    tempKappa = equivStrain;
    this->initDamaged (tempKappa, reducedTotalStrainVector, gp);
    this->computeDamageParam (omega, tempKappa, reducedTotalStrainVector, gp);
    if (omega < status->giveDamage()) {
      // unloading takes place
      omega = status->giveDamage();
    }
  } else _error ("giveRealStressVector: unsupported loading/uloading criteria");


  lmat->giveCharacteristicMatrix(de, ReducedForm, SecantStiffness, gp, atTime);
  de.times(1.0-omega);
  answer.beProductOf (de, reducedTotalStrainVector);

  // update gp
  status-> letTempStrainVectorBe (totalStrain);
  status-> letTempStressVectorBe (answer);
  status-> setTempKappa (tempKappa);
  status-> setTempDamage(omega);
  return ;
}



/*
The following routines are not necessary, since the base
StructuralMaterial class implements general algorithm how
to derive the stiffness and compliance matrices for
specific modes from full 3d stiffness/compliance.

However, these general services are not efficient, since
they require matrix inversion. The better and efficient way is
to overload these service and provide efficient implementation,
if available.
*/

void 
IsotropicDamageMaterial::givePlaneStressStiffMtrx (FloatMatrix& answer, 
                                                   MatResponseForm form, 
                                                   MatResponseMode mode,
                                                   GaussPoint* gp, 
                                                   TimeStep* atTime)
{
  IsotropicDamageMaterialStatus *status = 
    (IsotropicDamageMaterialStatus*) this -> giveStatus (gp);
  double om = status->giveTempDamage();
  om = min (om, 0.999999);

  this->giveLinearElasticMaterial()->
    giveCharacteristicMatrix (answer, form, mode, gp, atTime);
  answer.times (1.0-om);

  return ;
  
}


void 
IsotropicDamageMaterial::givePlaneStrainStiffMtrx (FloatMatrix& answer, 
                                                   MatResponseForm form,
                                                   MatResponseMode mode,
                                                   GaussPoint* gp, 
                                                   TimeStep* atTime)
{
  IsotropicDamageMaterialStatus *status = 
    (IsotropicDamageMaterialStatus*) this -> giveStatus (gp);
  double om = status->giveTempDamage();
  om = min (om, 0.999999);

  this->giveLinearElasticMaterial()->
    giveCharacteristicMatrix (answer, form, mode, gp, atTime);
  answer.times (1.0-om);

  return ;

}
void 
IsotropicDamageMaterial::give1dStressStiffMtrx (FloatMatrix& answer, 
                                                MatResponseForm form,
                                                MatResponseMode mode,
                                                GaussPoint* gp, 
                                                TimeStep* atTime)
{
  IsotropicDamageMaterialStatus *status = 
    (IsotropicDamageMaterialStatus*) this -> giveStatus (gp);
  double om = status->giveTempDamage();
  om = min (om, 0.999999);

  this->giveLinearElasticMaterial()->
    giveCharacteristicMatrix (answer, form, mode, gp, atTime);
  answer.times (1.0-om);

  return ;

}

int 
IsotropicDamageMaterial::giveIPValue (FloatArray& answer, 
                                      GaussPoint* aGaussPoint, 
                                      InternalStateType type, 
                                      TimeStep* atTime)
{
  IsotropicDamageMaterialStatus* status = 
    (IsotropicDamageMaterialStatus*) this -> giveStatus (aGaussPoint);
  if (type == IST_DamageTensor) {
    answer.resize(1);
    answer.at(1) = status->giveDamage();
    return 1;
  } else if (type == IST_DamageTensorTemp) {
    answer.resize(1);
    answer.at(1) = status->giveTempDamage();
    return 1;
  }  else if (type == IST_MaxEquivalentStrainLevel) {
    answer.resize(1);
    answer.at(1) = status->giveKappa ();
    return 1;
  }else return StructuralMaterial::giveIPValue (answer, aGaussPoint, type, atTime);
}
    



InternalStateValueType 
IsotropicDamageMaterial::giveIPValueType (InternalStateType type)
{
  if ((type == IST_DamageTensor)) return ISVT_TENSOR;
  else return StructuralMaterial::giveIPValueType (type);
}


int 
IsotropicDamageMaterial::giveIntVarCompFullIndx (IntArray& answer, 
                                                 InternalStateType type, 
                                                 MaterialMode mmode)
{
  if ((type == IST_DamageTensor)) {
    answer.resize (9);
    answer.at(1) = 1;
    return 1;
  } else 
    return StructuralMaterial::giveIntVarCompFullIndx (answer, type, mmode);
}


int
IsotropicDamageMaterial::giveIPValueSize (InternalStateType type, 
                                          GaussPoint* aGaussPoint)
{
  if ((type == IST_DamageTensor)) return 1;
  else return StructuralMaterial::giveIPValueSize (type, aGaussPoint);
}


void
IsotropicDamageMaterial :: giveThermalDilatationVector (FloatArray& answer, 
                                                        GaussPoint * gp,  
                                                        TimeStep* tStep)
//
// returns a FloatArray(6) of initial strain vector
// eps_0 = {exx_0, eyy_0, ezz_0, gyz_0, gxz_0, gxy_0}^T
// caused by unit temperature in direction of 
// gp (element) local axes
//
{
  answer.resize (6); answer.zero();
  answer.at(1) = this->tempDillatCoeff;
  answer.at(2) = this->tempDillatCoeff;
  answer.at(3) = this->tempDillatCoeff;
  
  return ;
}


IRResultType
IsotropicDamageMaterial :: initializeFrom (InputRecord* ir)
{
 const char *__proc = "initializeFrom"; // Required by IR_GIVE_FIELD macro
 IRResultType result;                   // Required by IR_GIVE_FIELD macro

 IR_GIVE_FIELD (ir, tempDillatCoeff, IFT_IsotropicDamageMaterial_talpha, "talpha"); // Macro
 return StructuralMaterial::initializeFrom (ir);
}

In the last part, the implementation of the IsotropicDamageMaterialStatus class services follows. The initTempStatus and updateYourself services are used to initialize temporary history variables according to the previous equilibrium state or to update the non-temporary variables, when a new equilibrium is reached. Also, the implementation of storeContext and restoreContext services is shown. The service printYourself is also overloaded. It allows to print history variables into output file. Note that all overloaded methods call the corresponding method of the base class (to ensure that processing corresponding to the base class is performed, to guarantee consistency) and then the required functionality is added.

IsotropicDamageMaterialStatus::IsotropicDamageMaterialStatus (int n, 
                                                              Domain* d, 
                                                              GaussPoint* g)
:StructuralMaterialStatus(n,d,g)

{
  kappa = tempKappa = 0.0;
  damage = tempDamage = 0.0;
}


IsotropicDamageMaterialStatus::~IsotropicDamageMaterialStatus () 
{}


void 
IsotropicDamageMaterialStatus :: printOutputAt  (FILE *file, TimeStep* tStep)
{
  
  StructuralMaterialStatus :: printOutputAt (file, tStep);
  fprintf (file,"status { ");
  if (this->damage > 0.0) {
    fprintf (file,"kappa %f, damage %f ",this->kappa, this->damage);
  }
  
  fprintf (file,"}\n");
}
 

void 
IsotropicDamageMaterialStatus::initTempStatus ()
{
  StructuralMaterialStatus :: initTempStatus();
  this->tempKappa = this->kappa;
  this->tempDamage= this->damage;
}

void 
IsotropicDamageMaterialStatus::updateYourself(TimeStep* atTime)
{
  StructuralMaterialStatus::updateYourself(atTime);
  this->kappa = this->tempKappa;
  this->damage= this->tempDamage;
}


contextIOResultType
IsotropicDamageMaterialStatus::saveContext (FILE* stream, void *obj)
{
  contextIOResultType iores;

  // save parent class status
  if ((iores = StructuralMaterialStatus :: saveContext (stream, obj)) 
       != CIO_OK) THROW_CIOERR(iores);

  // write a raw data
  if (fwrite(&kappa,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR);
  if (fwrite(&damage,sizeof(double),1,stream)!= 1) THROW_CIOERR(CIO_IOERR);

  return CIO_OK;
}

contextIOResultType
IsotropicDamageMaterialStatus::restoreContext(FILE* stream, void *obj)
{
  contextIOResultType iores;

  // read parent class status
  if ((iores = StructuralMaterialStatus :: restoreContext(stream,obj)) 
    != CIO_OK) THROW_CIOERR(iores);

  // read raw data 
  if (fread (&kappa,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR);
  if (fread (&damage,sizeof(double),1,stream) != 1) THROW_CIOERR(CIO_IOERR);

  return CIO_OK;
}

The above general class provides the general framework. The derived classes have to implement only services computeEquivalentStrain and computeDamageParam, which are only related to particular material model under consideration.