@DSL Default;
@Author Tran Thang DANG;
@Date 05/02/2016;
@Behaviour PhaseFieldDisplacement;
// @ModellingHypothesis {PlaneStrain,Tridimensional};
@Description{
  "A damage model coupled with a phase "
  "field approach"
}

@MaterialProperty stress yg;
yg.setGlossaryName("YoungModulus") ;
@MaterialProperty real nu ;
nu.setGlossaryName("PoissonRatio") ;
@Parameter real kk = 1e-6 ;

@StateVariable real H;
H.setEntryName("DamageDrivingForce");

@StateVariable real de;
de.setGlossaryName("Damage");

@ProvidesSymmetricTangentOperator ;
@Integrator{
  constexpr const strain emin = 1.e-12;
  // positive part
  const auto f  = [](const real x){return x>0 ? x : 0;};
  // derivative of the positive part
  const auto df = [&emin](const real x)
    {return std::abs(x)<emin ? 0.5 : ((x<0) ? 0 : 1);};
  // update the damage
  const auto d = de+dde;
  // lame coefficients
  const auto lambda = computeLambda(yg,nu);
  const auto mu     = computeMu(yg,nu);
  // computation of the stress, positive energy density and consistent
  // tangent operator
  const StrainStensor e_ = eto + deto;
  const auto fdf     = e_.computeIsotropicFunctionAndDerivative(f,df,emin*0.1);
  const auto& ep     = fdf.first;  // positive part of e_
  const auto& dep_de = fdf.second; // derivative of the positive part of e_
  const StrainStensor en = e_-ep;  // negative part of e_
  // energy density
  const strain tr  = trace(e_);
  const strain tr_p = max(tr,real(0));
  const strain tr_n = min(tr,real(0));
  H = max(H,(((lambda/2)*(tr_p)*(tr_p))+(mu*(ep|ep))));
  // stress
  const auto deff = ((1-d)*(1-d))+kk;
  sig = 2*mu*(deff*ep+en)+lambda*(deff*tr_p+tr_n)*StrainStensor::Id();
  // consistent tangent operator (secant one here)
  if(computeTangentOperator_){
    if(tr>=0){
      Dt=deff*(lambda*Stensor4::IxI()+2*mu*dep_de)+(2*mu*(Stensor4::Id()-dep_de));
    } else {
      Dt=deff*2*mu*dep_de+(lambda*Stensor4::IxI()+2*mu*(Stensor4::Id()-dep_de));
    }
  }
  static_cast<void>(smt);
} // end of @Integrator

@APosterioriTimeStepScalingFactor {
  // pas d'acroissement de l'endommagement sur le pas de temps
  if (dde<1.e-4){
    return {true,1.e-2/(max(dde,1e-4))};
  }
  return {true,1};
}