### 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

where is elastic stiffness matrix of the undamaged material and is the damage parameter. Initially, 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 . For , the stiffness completely disappears.

In the present context, the matrix represents the secant stiffness that relates the total strain to the total stress

Similarly to the theory of plasticity, a loading function 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 , describing the evolution of the damage. Physically, is a scalar measure of the largest strain level ever reached. The loading function usually has the form

where 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 (). This is expressed by the following loading/unloading conditions

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

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 and the damage parameter and reducing the effective stress according to .

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 (). In addition, the corresponding damage level will be stored. This is not necessary because damage can be always computed from corresponding . 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 and 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;

/**
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 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)
*/
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 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
*/
GaussPoint* aGaussPoint,
InternalStateType type,
TimeStep* atTime) ;

/**
Returns the mask of reduced indexes of
Internal Variable component .
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.
*/
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)
*/
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 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)
*/
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 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)
*/
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 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)
*/
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
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);

return ;

}

void
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) {
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()) {
omega = status->giveDamage();
}

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

// update gp
status-> letTempStrainVectorBe (totalStrain);
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
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);

return ;

}

void
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);

return ;

}
void
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);

return ;

}


int
GaussPoint* aGaussPoint,
InternalStateType type,
TimeStep* atTime)
{
IsotropicDamageMaterialStatus* status =
(IsotropicDamageMaterialStatus*) this -> giveStatus (aGaussPoint);
if (type == IST_DamageTensor) {
return 1;
} else if (type == IST_DamageTensorTemp) {
return 1;
}  else if (type == IST_MaxEquivalentStrainLevel) {
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
InternalStateType type,
MaterialMode mmode)
{
if ((type == IST_DamageTensor)) {
return 1;
} else
}

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

void
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
//
{

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;

if ((iores = StructuralMaterialStatus :: restoreContext(stream,obj))
!= CIO_OK) THROW_CIOERR(iores);