HYPERELASTIC_LAW example

Description

This model is a hyper-elastic component for either the Hyperelastic model or the Hyper-viscoelastic model.

\[\ten \sigma = \tenf D_{el}:\ln(\ten U) = \tenf D_{el}:\ln(\ten B^{1/2})\]

Code Listing

#include <Behavior.h>
#include <Integration_result.h>
#include <Mechanical_behavior.h>
#include <Hyperelastic_behavior.h>
#include <Swap.h>
#include <ZMath.h>
#include <Print.h>
#include <Object_factory.h>
#include <Hyperelastic_behavior.h>
#include <Arruda_boyce_hyper.h>
#include <Tensor_operator.h>
#include <Print.h>

class LOG_HYPER_LAW : public HYPERELASTIC_LAW {
  protected :
      PTR<S_COEFFICIENT_TENSOR4> elasticity;

      TENSOR2       one;
      TENSOR2_VAUX  pk2_0;
      TENSOR4       dlnU_dB;
      TENSOR4       dB_dD;
      TENSOR4       C_star;
      bool          pull_back;
      bool          fd_log;
      bool          modify_for_jaumann;
      bool          modify_for_truesdale;
      bool          symmetrize;

  public :
      LOG_HYPER_LAW();
      virtual ~LOG_HYPER_LAW();

      virtual void initialize(ASCII_FILE&,MATERIAL_PIECE*);
      virtual int base_read(const STRING&,ASCII_FILE&);
      virtual void setup(int& flux_pos, int& grad_pos, int& vi_pos, int& va_pos);
      virtual bool calc_coef();

      virtual INTEGRATION_RESULT* integrate(const  VECTOR& delta_grad, int flags);

      TENSOR2& give_pk2_0();
      virtual BEHAVIOR::STRESS_MEASURE get_stress_measure() const;
};

DECLARE_OBJECT(HYPERELASTIC_LAW, LOG_HYPER_LAW, logarithmic);
BEHAVIOR_READER(HYPERELASTIC_BEHAVIOR,hyper_elastic_logarithmic);

LOG_HYPER_LAW::LOG_HYPER_LAW()
{ pull_back = FALSE;
  fd_log    = FALSE;
  modify_for_jaumann = FALSE;
  modify_for_truesdale = FALSE;
  symmetrize = FALSE;
}

LOG_HYPER_LAW::~LOG_HYPER_LAW()
{
}

void LOG_HYPER_LAW::initialize(ASCII_FILE& file,MATERIAL_PIECE* mp)
{ HYPERELASTIC_LAW::initialize(file,mp);
  sig.initialize(this,"sig",tsz(),1);
  tg.resize(tsz());
  one.resize(tsz());
  dlnU_dB.resize(tsz());
  dB_dD.resize(utsz());
  C_star.resize(tsz());
  one = TENSOR2::unity(tsz());
}

int LOG_HYPER_LAW::base_read(const STRING& strin,ASCII_FILE& file)
{ STRING str;
  prn("strin:"+strin);
  if (strin=="**pull_back") pull_back = TRUE;
  else if (strin=="**fd_log") fd_log = TRUE;
  else if (strin=="**modify_for_jaumann") modify_for_jaumann = TRUE;
  else if (strin=="**modify_for_truesdale") modify_for_truesdale = TRUE;
  else if (strin=="**dont_symmetrize") symmetrize = FALSE;
  else if (strin=="**symmetrize") symmetrize = TRUE;
  else if (strin=="**elasticity") {
       elasticity = S_COEFFICIENT_TENSOR4::read(file,this);
  }
  else return 0;
  return 1;
}

TENSOR2& LOG_HYPER_LAW::give_pk2_0()
{ return pk2_0;
}

BEHAVIOR::STRESS_MEASURE LOG_HYPER_LAW::get_stress_measure() const
{ if (pull_back) return BEHAVIOR::PK2_0;
  return BEHAVIOR::CAUCHY;
}


bool LOG_HYPER_LAW::calc_coef()
{ HYPERELASTIC_LAW::calc_coef();
  return TRUE;
}

void LOG_HYPER_LAW::setup(int& flux_pos, int& grad_pos, int& vi_pos, int& va_pos)
{
  if (pull_back) pk2_0.initialize(this,"pk2",tsz(),0);
  else           pk2_0.resize(tsz());

  HYPERELASTIC_LAW::setup(flux_pos,grad_pos,vi_pos,va_pos);
}

INTEGRATION_RESULT* LOG_HYPER_LAW::integrate(const VECTOR& df, int)
{ int i,c;
  set_var_aux_to_var_aux_ini();

  TENSOR2 Ft      = transpose(F);
  TENSOR2 B       = syme(F*Ft);

  //
  // All this is the calculation of ln(U) and its tangent
  //
  SMATRIX dlnU_dB(tsz());
  TENSOR2 lnU   = 0.5*log_tensor(B, dlnU_dB, TRUE);
  dlnU_dB      *= 0.5;
  dB_dD         = 2.0*Mr(to_5_9(B));

  //
  // Ok, now the model is simple
  //
  sig     = (*elasticity)*lnU;
  tg      = (*elasticity)*expand_in(expand_out(dlnU_dB)*dB_dD);

  //
  // Convert a trusdale rate to jaumann rate
  //
  if (modify_for_jaumann) {
      TENSOR2 ss      = to_5_9(sig);
      TENSOR4 sig_one = sig^one;
      C_star = expand_in(Mr(ss) + Ml(ss)) - sig_one;
      tg     = tg + C_star;
  }
  //
  // Convert a trusdale rate to jaumann rate
  //
  else if (modify_for_truesdale) {
      TENSOR2 ss      = to_5_9(sig);
      TENSOR4 sig_one = sig^one;
      C_star = expand_in(Mr(ss) + Ml(ss)) - sig_one;
      tg     = tg - C_star;
  }

  //
  // And we give the option if a total lagrange model would
  // make the user happy. this is a nice check on the meaning
  // of the tangents.
  //
  if (pull_back) {
      prn("pull back");
      double  J  = F.determin();
      TENSOR2 Fi = inverse(F);
      pk2_0 = J*rotate_tensor(sig,Fi);
      tg    = J*rotate_matrix(tg,Fi);
      sig   = pk2_0;
  }
  if (symmetrize) {
      tg = 0.5*(tg + transpose(tg));
  }

  return NULL;
}