ZebFront example

Description

This is an example of a ZebFront model using modular components and Runge-kutta integration for the simulation mode of operation.

#include <Elasticity.h>
#include <External_parameter.h>
#include <Basic_nl_behavior.h>
#include <Flow.h>
#include <Criterion.h>
#include <Isotropic.h>
#include <Print.h>

@Class MODULAR_PLASTIC_BEHAVIOR : BASIC_NL_BEHAVIOR, BASIC_SIMULATOR {
   @Name     modular_plastic;
   @SubClass ELASTICITY elasticity;
   @SubClass CRITERION  criterion;
   @SubClass FLOW       flow;
   @SubClass ISOTROPIC_HARDENING  isotropic;

   @Coefs    C [0-N] @Factor 1./1.5;
   @Coefs    D [!C];
   @tVarInt  eel, alpha [!C];
   @sVarInt  evcum;
   @tVarAux  X [!C], evi;
   @tVarAux  Xtot;
   @tVarUtil m [!C];
   @tVarUtil Xdot;
};

@PostStep {
  evi  = eto - eel;
  sig  = *elasticity*eel;
  if (integration&LOCAL_INTEGRATION::THETA_ID) {
     SMATRIX tmp(psz,f_grad,0,0);
     if (Dtime>0.0) m_tg_matrix=*elasticity*tmp;
     else           m_tg_matrix=*elasticity;
  } else m_tg_matrix=*elasticity;
}

@Derivative {
  int c;
  ELASTICITY& E=*elasticity;
  sig = E*eel;

  Xtot = 0.0;
  for (c=0;c<!C;c++) Xtot += X[c] = C[c]*alpha[c];
  double    radius = isotropic->radius(evcum);
  TENSOR2   sigeff = deviator(sig) - Xtot;
  double    f      = criterion->yield(sigeff, radius);

  int yld = (flow->plasticity()) ? (f>=0.0) : (f>0.0);

  if (yld) {
     double miso=0.0;
     m.resize(!C);
     TENSOR2          norm     = criterion->normal();

     Xdot = 0.0;
     for (c=0;c<!C;c++) {
         m[c]  = norm - X[c]*(D[c]/C[c]);
         Xdot += C[c]*m[c];
     }

     if (flow->plasticity()) {
        devcum    = norm|(E*deto);
        @Simulation { ERROR("can't resolve multiplier in the simulation"); }
        if (!if_constant_coefs()) {
             LIST<EXTERNAL_PARAM*>& ep=EXTERNAL_PARAM::Get_EP_list();
             for (int i=0;i<!ep;i++) {
                 double dparam = ((*curr_mat_data->param_set())[i] -
                                 (*curr_mat_data->param_set_ini())[i])/Dtime;
                 devcum -= isotropic->dradius_dparam(ep[i]->name())*dparam;;
                 for (c=0;c<!C;c++) devcum -= (norm|alpha[c])*C[c].d_param(ep[i]->name()())*dparam;
             }
        }
        double          H1        = norm|(E*norm);
        double          H2        = norm|Xdot;
        double          H3        = isotropic->dradius_dflow();
                        devcum   /= H1 + H2 + H3;
     } else {
        devcum    = flow->flow_rate(evcum,f);
     }
     for (c=0;c<!C;c++) dalpha[c] = devcum*m[c];
     resolve_flux_grad(E, deel, deto, devcum*norm);

  } else resolve_flux_grad(E, deel, deto);

  @Simulation { evi = eto - eel; }
}