BEHAVIOR example

Description

#include <Error_messager.h>
#include <File.h>

#include <Clock.h>
#include <Coefficient.h>
#include <Coefficient_T4.h>
#include <Integration_method.h>
#include <Integration_result.h>
#include <Mat_data.h>
#include <Mechanical_behavior.h>
#include <Verbose.h>

//----------------------------------------------------------------------------
class DAMAGE_ELASTIC : public M_TLE_B_SD {
     SCALAR_VAUX   damage;
     SCALAR_VINT   y_max;
     SMATRIX       tgmat;
     COEFF         Y0, alpha;

  public:
     DAMAGE_ELASTIC();
     virtual void initialize(ASCII_FILE& file,int dim, LOCAL_INTEGRATION*);

     INTEGRATION_RESULT* integrate( MAT_DATA&     mdat,
                                    const VECTOR& delta_grad,
                                    MATRIX*&      tg_matrix,
                                    int           flags);
};
//----------------------------------------------------------------------------

BEHAVIOR_READER(DAMAGE_ELASTIC,damage_elasticity)

DAMAGE_ELASTIC::DAMAGE_ELASTIC() {}

void DAMAGE_ELASTIC::initialize(ASCII_FILE& file,int dim, LOCAL_INTEGRATION*)
{ M_TLE_B_SD::initialize(file, dim,NULL);
  damage.initialize(this,"damage",1);
  y_max.initialize(this,"y_max",1);

   for (;;) {
       STRING str=file.getSTRING();
       if (str.start_with("***") ) break;
       else if (str=="**Y0")     Y0.read(str,file,this);
       else if (str=="**alpha")  alpha.read(str,file,this);
       else if (!base_read(str,file)) INPUT_ERROR("Unknown coefficient: "+str);
   } file.back();

   tgmat.resize(tsz());
}

INTEGRATION_RESULT* DAMAGE_ELASTIC::integrate(MAT_DATA& mdat,
         const VECTOR&, MATRIX*& tg_matrix, int flags_in)
{ int step=0;
  tg_matrix=&tgmat;
  if (!curr_ext_param)
     curr_ext_param = *mdat.param_set();
  attach_all(mdat);
  calc_local_coefs();

  TENSOR2 eel = eto;
  if (thermal_strain.if_not_null() && mdat.param_set())
     eel -= thermal_strain->compute_strain();

  double  y_bar   = 0.5*(eel|(*elasticity*eel));
  if (y_bar>y_max) { y_max=y_bar; step=1; }
  double  ym_root = sqrt(y_max);
  double  yz_root = sqrt(Y0);

  if (ym_root <= yz_root) damage = 0.0;
  else                    damage = alpha()*(ym_root - yz_root);
  if (damage >= 0.99)     damage = 0.99;


  if (flags_in&CALC_TG_MATRIX) {
     tgmat=(1.0-damage)*(*elasticity);
     sig = tgmat*eel;
     if (step) tgmat-=(0.5*alpha/ym_root/(1.0-damage)/(1.0-damage))*(sig^sig);
  } else sig=*elasticity * ((1.0-damage)*eel);

  return NULL;
}