mlawNonLocalPorousCoupled.cpp

//
// C++ Interface: material law
//
// Description: Porous material law with Gurson void growth and coalescence
//
// Author:  <J. Leclerc>, (C) 2017
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include "mlawNonLocalPorousCoupled.h"
#include "nonLinearMechSolver.h"


mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const int num,const double E,const double nu, const double rho,
                const double q1, const double q2, const double q3,
                const double fVinitial, const double lambda0, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
                const double tol, const bool matrixbyPerturbation, const double pert) :
                    mlawNonLocalPorosity(num, E, nu, rho, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert),
                    _withCrackTransition(false),_CfTOffsetMethod(1)
{
  _mlawGrowth = new mlawNonLocalDamageGurson(num, E, nu, rho, q1, q2, q3, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert);

  _mlawCoales = new mlawNonLocalPorousThomasonLaw(num, E, nu, rho, fVinitial, lambda0, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
  // modify current coalescence law
  mlawNonLocalPorosity::setCoalescenceLaw(*_mlawCoales->getCoalescenceLaw());

  _failedTol = 0.99;
};

mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const int num,const double E,const double nu, const double rho,
                const double q1, const double q2, const double q3,
                const double fVinitial, const double lambda0, const double kappa, const J2IsotropicHardening &j2IH, const CLengthLaw &cLLaw,
                const double tol, const bool matrixbyPerturbation, const double pert) :
                    mlawNonLocalPorosity(num, E, nu, rho, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert),
                    _withCrackTransition(false),_CfTOffsetMethod(1)
{
  _mlawGrowth = new mlawNonLocalDamageGurson(num, E, nu, rho, q1, q2, q3, fVinitial, j2IH, cLLaw, tol, matrixbyPerturbation, pert);

  _mlawCoales = new mlawNonLocalPorousThomasonLaw(num, E, nu, rho, fVinitial, lambda0, kappa, j2IH, cLLaw, tol, matrixbyPerturbation, pert);
  // modify current coalescence law
  mlawNonLocalPorosity::setCoalescenceLaw(*_mlawCoales->getCoalescenceLaw());

  _failedTol = 0.99;
};



mlawNonLocalPorousCoupledLaw::mlawNonLocalPorousCoupledLaw(const mlawNonLocalPorousCoupledLaw &source) :
                    mlawNonLocalPorosity(source),_withCrackTransition(source._withCrackTransition),
                    _CfTOffsetMethod(source._CfTOffsetMethod)
{
  _mlawGrowth = NULL;
  if(source._mlawGrowth != NULL) _mlawGrowth = dynamic_cast<mlawNonLocalDamageGurson*>(source._mlawGrowth->clone());

  _mlawCoales = NULL;
  if(source._mlawCoales != NULL) _mlawCoales = dynamic_cast<mlawNonLocalPorousThomasonLaw*>(source._mlawCoales->clone());
};


void mlawNonLocalPorousCoupledLaw::setNonLocalMethod(const int i){
  mlawNonLocalPorosity::setNonLocalMethod(i);
  _mlawGrowth->setNonLocalMethod(i);
  _mlawCoales->setNonLocalMethod(i);
};

// Options settings
void mlawNonLocalPorousCoupledLaw::setOrderForLogExp(const int newOrder)
{
  mlawNonLocalPorosity::setOrderForLogExp(newOrder);
  _mlawGrowth->setOrderForLogExp(newOrder);
  _mlawCoales->setOrderForLogExp(newOrder);
  Msg::Info("mlawNonLocalPorousCoupledLaw :: order = %d is used to approximate log and exp",newOrder);
};


void mlawNonLocalPorousCoupledLaw::setViscosityLaw(const viscosityLaw& vico)
{
  mlawNonLocalPorosity::setViscosityLaw(vico);
  _mlawGrowth->setViscosityLaw(vico);
  _mlawCoales->setViscosityLaw(vico);
};

void mlawNonLocalPorousCoupledLaw::setNucleationLaw(const NucleationLaw& added_gdnLaw)
{
  mlawNonLocalPorosity::setNucleationLaw(added_gdnLaw);
  _mlawGrowth->setNucleationLaw(added_gdnLaw);
  _mlawCoales->setNucleationLaw(added_gdnLaw);
};

void mlawNonLocalPorousCoupledLaw::setScatterredInitialPorosity(double f0min, double f0max)
{
  mlawNonLocalPorosity::setScatterredInitialPorosity(f0min,f0max);
  _mlawGrowth->setScatterredInitialPorosity(f0min, f0max);
  _mlawCoales->setScatterredInitialPorosity(f0min, f0max);

  Msg::Info("mlawNonLocalPorousCoupledLaw :: Add scatter in initial porosity");
  _randMinbound = f0min;
  if(f0min < 0.){Msg::Fatal("mlawNonLocalPorousCoupledLaw :: Wrong scatter parameter fV0_min");}
  _randMaxbound = f0max;
  if(f0max < f0min){Msg::Fatal("mlawNonLocalPorousCoupledLaw :: Wrong scatter parameter fV0_max");}
};

void mlawNonLocalPorousCoupledLaw::setCoalescenceLaw(const CoalescenceLaw& added_coalsLaw)
{
  mlawNonLocalPorosity::setCoalescenceLaw(added_coalsLaw);
  _mlawCoales->setCoalescenceLaw(added_coalsLaw);
};


void  mlawNonLocalPorousCoupledLaw::setYieldSurfaceExponent(const double newN)
{
  _mlawCoales->setYieldSurfaceExponent(newN);
};

void mlawNonLocalPorousCoupledLaw::setRoundedYieldSurfaceMethod(const int method){
  _mlawCoales->setRoundedYieldSurfaceMethod(method);
};
void mlawNonLocalPorousCoupledLaw::setOnsetTriaxialityForRoundedYieldSurface(const double newTc){
  _mlawCoales->setOnsetTriaxialityForRoundedYieldSurface(newTc);
};

void mlawNonLocalPorousCoupledLaw::setCfTOffsetMethod(const int MethodNumber){
  _CfTOffsetMethod = MethodNumber;
  Msg::Info("mlawNonLocalPorousCoupledLaw::useCfTOffset: Offset method %d is used",_CfTOffsetMethod);
};

void mlawNonLocalPorousCoupledLaw::setCrackTransition(const bool fl){
  _withCrackTransition = fl;
};

void mlawNonLocalPorousCoupledLaw::setBalanceEquationType(const int type){
  _mlawCoales->setBalanceEquationType(type);
};

void mlawNonLocalPorousCoupledLaw::setShearPorosityGrowthFactor(const double k) {
  mlawNonLocalPorosity::setShearPorosityGrowthFactor(k);
  _mlawGrowth->setShearPorosityGrowthFactor(k);
  _mlawCoales->setShearPorosityGrowthFactor(k);
};
void mlawNonLocalPorousCoupledLaw::setLocalRegularizedFunction(const scalarFunction& fct){
  mlawNonLocalPorosity::setLocalRegularizedFunction(fct);
  _mlawGrowth->setLocalRegularizedFunction(fct);
  _mlawCoales->setLocalRegularizedFunction(fct);
};
void mlawNonLocalPorousCoupledLaw::setCorrectedRegularizedFunction(const scalarFunction& fct){
  mlawNonLocalPorosity::setCorrectedRegularizedFunction(fct);
  _mlawGrowth->setCorrectedRegularizedFunction(fct);
  _mlawCoales->setCorrectedRegularizedFunction(fct);
};
void mlawNonLocalPorousCoupledLaw::setPostBlockedDissipationBehavior(const int method){
  mlawNonLocalPorosity::setPostBlockedDissipationBehavior(method);
  _mlawGrowth->setPostBlockedDissipationBehavior(method);
  _mlawCoales->setPostBlockedDissipationBehavior(method);
};

void mlawNonLocalPorousCoupledLaw::clearCLengthLaw(){
  mlawNonLocalPorosity::clearCLengthLaw();
  _mlawGrowth->clearCLengthLaw();
  _mlawCoales->clearCLengthLaw();
};

void mlawNonLocalPorousCoupledLaw::setCLengthLaw(const CLengthLaw& clength){
  mlawNonLocalPorosity::setCLengthLaw(clength);
  _mlawGrowth->setCLengthLaw(clength);
  _mlawCoales->setCLengthLaw(clength);
};
//-------------------added
void mlawNonLocalPorousCoupledLaw::setReferenceT(const double referenceT){
  mlawNonLocalPorosity::setReferenceT(referenceT);
  _mlawGrowth->setReferenceT(referenceT);
  _mlawCoales->setReferenceT(referenceT);
};
void mlawNonLocalPorousCoupledLaw::setThermalExpansionCoefficient(const double alp){
  mlawNonLocalPorosity::setThermalExpansionCoefficient(alp);
  _mlawGrowth->setThermalExpansionCoefficient(alp);
  _mlawCoales->setThermalExpansionCoefficient(alp);
};
void mlawNonLocalPorousCoupledLaw::setSubStepping(const bool fl, const int maxNumStep){
  mlawNonLocalPorosity::setSubStepping(fl,maxNumStep);
  _mlawGrowth->setSubStepping(fl,maxNumStep);
  _mlawCoales->setSubStepping(fl,maxNumStep);
};

void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_rho(const scalarFunction& Tfunc){
  mlawNonLocalPorosity::setTemperatureFunction_rho(Tfunc);
  _mlawGrowth->setTemperatureFunction_rho(Tfunc);
  _mlawCoales->setTemperatureFunction_rho(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_E(const scalarFunction& Tfunc){
  mlawNonLocalPorosity::setTemperatureFunction_E(Tfunc);
  _mlawGrowth->setTemperatureFunction_E(Tfunc);
  _mlawCoales->setTemperatureFunction_E(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_nu(const scalarFunction& Tfunc){
  mlawNonLocalPorosity::setTemperatureFunction_nu(Tfunc);
  _mlawGrowth->setTemperatureFunction_nu(Tfunc);
  _mlawCoales->setTemperatureFunction_nu(Tfunc);
};
void mlawNonLocalPorousCoupledLaw::setTemperatureFunction_alp(const scalarFunction& Tfunc){
  mlawNonLocalPorosity::setTemperatureFunction_alp(Tfunc);
  _mlawGrowth->setTemperatureFunction_alp(Tfunc);
  _mlawCoales->setTemperatureFunction_alp(Tfunc);
};

void mlawNonLocalPorousCoupledLaw::createIPState(IPStateBase* &ips,const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
{
  //bool inter=true;
  //const MInterfaceElement *iele = static_cast<const MInterfaceElement*>(ele);
  //if(iele==NULL) inter=false;
  double fvInit = getInitialPorosity();
  IPVariable* ipvi = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
  IPVariable* ipv1 = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
  IPVariable* ipv2 = new IPNonLocalPorosity(fvInit,_j2IH,&_cLLaw,&_gdnLawContainer,_coalescenceLaw);
  if(ips != NULL) delete ips;
  ips = new IP3State(state_,ipvi,ipv1,ipv2);

};


double mlawNonLocalPorousCoupledLaw::yieldFunction(const double kcorEq, const double pcor, const double R, const double yieldfV,
                                const IPNonLocalPorosity* q0, const IPNonLocalPorosity* q1, const double* T) const{
  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();

  if (coalesFlag){
    return _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
  }
  else return _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
};


void mlawNonLocalPorousCoupledLaw::computeResidual(SVector3& res, STensor3& J,
                const double kcorEqpr, const double pcorpr, const double R, const double H,
                const double yieldfV,  const double DyieldfVDfV,
                const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
                const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
                const double DeltafV, const double DDeltafVDDeltaHatD, const double DDeltafVDDeltaHatQ, const double DDeltafVDDeltaHatP,
                const double* T
              ) const{

  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();

  if (coalesFlag){
    _mlawCoales->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
                                DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
                                DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
  }
  else {
    _mlawGrowth->computeResidual(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
                                DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,
                                DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,T);
  }
}

void mlawNonLocalPorousCoupledLaw::computeDResidual(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
                                  double &dres0dtildefV, double &dres1dtildefV, double &dres2dtildefV,
                                  const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
                                  const double yieldfV, const double DyieldfVDfV, const double DyieldfVDtildefV,
                                  const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
                                  const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
                                  const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
                                  const double* T, const double* DRDT,
                                  const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
                                  double *dres0dT, double *dres1dT, double *dres2dT
                                  ) const{

  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();

  if (coalesFlag){
    _mlawCoales->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
                                  dres0dtildefV,dres1dtildefV,dres2dtildefV,
                                  devKcorpr,kcorEqpr,pcorpr,R,
                                  yieldfV,DyieldfVDfV,DyieldfVDtildefV,
                                  DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD,DeltaHatQ,DeltaHatP,
                                  T,DRDT,
                                  DdevKcorprDT ,DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }
  else {
    _mlawGrowth->computeDResidual(dres0dEpr,dres1dEpr,dres2dEpr,
                                  dres0dtildefV,dres1dtildefV,dres2dtildefV,
                                  devKcorpr,kcorEqpr,pcorpr,R,
                                  yieldfV,DyieldfVDfV,DyieldfVDtildefV,
                                  DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD,DeltaHatQ,DeltaHatP,
                                  T,DRDT,
                                  DdevKcorprDT ,DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }

};


void mlawNonLocalPorousCoupledLaw::computeDResidualLocal(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
                                  const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
                                  const double yieldfV, const double DyieldfVDfV,
                                  const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
                                  const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
                                  const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
                                  const double* T, const double* DRDT,
                                  const STensor3* DdevKcorprDT , const double* DpcorprDT, const double* DDeltafVDT,
                                  double *dres0dT, double *dres1dT, double *dres2dT
                                  ) const{
  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();

  if (coalesFlag){
    _mlawCoales->computeDResidualLocal(dres0dEpr, dres1dEpr, dres2dEpr,
                                  devKcorpr, kcorEqpr, pcorpr,R,
                                  yieldfV, DyieldfVDfV,
                                  DDeltafVDEpr, DdevKcorprDEpr, DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD, DeltaHatQ,DeltaHatP,
                                  T, DRDT,
                                  DdevKcorprDT , DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }
  else{
    _mlawGrowth->computeDResidualLocal(dres0dEpr, dres1dEpr, dres2dEpr,
                                  devKcorpr, kcorEqpr, pcorpr,R,
                                  yieldfV, DyieldfVDfV,
                                  DDeltafVDEpr, DdevKcorprDEpr, DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD, DeltaHatQ,DeltaHatP,
                                  T, DRDT,
                                  DdevKcorprDT , DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }
}

void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0, const double* T) const{
   // Get ipvs
  const IPThomasonCoalescence* q0Thom = dynamic_cast<const IPThomasonCoalescence*>(&q0->getConstRefToIPCoalescence());
  IPThomasonCoalescence* q1Thom = dynamic_cast<IPThomasonCoalescence*>(&q1->getRefToIPCoalescence());
  if (q1Thom == NULL){
    Msg::Fatal("mlawNonLocalPorousThomasonLaw::computeResidual: IPThomasonCoalescence must be used!");
  }

  const STensor3& kCor = q1->getConstRefToCorotationalKirchhoffStress();
  static STensor3 devKcor;
  static double kcorEq, pcor;
  STensorOperation::decomposeDevTr(kCor,devKcor,pcor);
  pcor /= 3.;
  kcorEq = sqrt(1.5*devKcor.dotprod());

  const double R = q1->getCurrentViscoplasticYieldStress();

  // get needed data for testing criterion
  double yieldfV = q1->getYieldPorosity();
  double R0 = _j2IH->getYield0();
  double Cft = q1Thom->getConcentrationFactor();


  if (_withCrackTransition){
    // If crack transition is used, coalescence is used only at interface as crack insertion criterion.
    // In the bulk, nothing is done (== always Gurson)
    if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS)
    {
      // Compute normal components on the interface from kCor
      double sigmax = 0.;
      static SVector3 Norm;
      //Norm = q1->getConstRefToCurrentOutwardNormal();
      Norm = q1->getConstRefToReferenceOutwardNormal();
      if (Norm.norm() > 0){
        Norm.normalize();
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            sigmax += Norm(i)*kCor(i,j)*Norm(j);
          }
        }
        // update maximal stress
        q1Thom->getRefToMaximalTractionStress()= q0Thom->getMaximalTractionStress();
        if (sigmax > q1Thom->getMaximalTractionStress()){
          q1Thom->getRefToMaximalTractionStress() = sigmax;
        }
      }
      else{
        Msg::Fatal("mlawNonLocalPorousCoupledLaw::updatePlasticState: zero normal");
      }

      // - Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface

      if (!q0Thom->getCoalescenceOnsetFlag())
      {
        // Check Thomason criterion:
        double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
        /* check this value */
        if ((sigmax - Cft*R > _tol) and (yieldThomason > _tol)){
          // If coalescence occurs, update Cft offset value and others variables
          //
          // why offset? because:
          // we want to follow the continuity of the Thomason yeild surface after satisfying coalescence
          // the value of offset can be larger than 1
          // this value is found by resolving the equation
          // before coalescence yieldThomason = sthingOfStress - Cft*R/R0
          // at coalescence  yieldThomason + Cft*R/R0 - CftOffset*Cft*R/R0  = 0
          // leading to offset = (yieldThomason + Cft*R/R0)/(Cft*R/R0) = 1+yieldThomason*R0/(Cft*R)

          // Compute Cft Offset
          if (_CfTOffsetMethod == 1){
            q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
            Msg::Info("coalescence occurs for the first time: Offset =%e",q1Thom->getRefToCrackOffsetOnCft());
          }
          else if (_CfTOffsetMethod == 0){
            q1Thom->getRefToCrackOffsetOnCft() = 1.;
          }

          q1Thom->getRefToCoalescenceOnsetFlag() = true;
          // Set onset variables to their actual value
          q1Thom->getRefToPorosityAtCoalescenceOnset() = q1->getYieldPorosity();
          q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
          q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
          q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q1Thom->getShapeFactor();

          q1Thom->getRefToCoalescenceActiveFlag() = true;
          q1Thom->getRefToAccelerateRate() = 1.;
        }
        else{
          // Coalescence never occurs and is not occuring at the next step
          q1Thom->getRefToCoalescenceOnsetFlag() = false;
          // Reset onset variables to their default value
            // Porosity
          q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
            // Geometrical parameters
          q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
          q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
          q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;

          q1Thom->getRefToCoalescenceActiveFlag() = false;

          q1Thom->getRefToAccelerateRate() = 1.;
          q1Thom->getRefToCrackOffsetOnCft() = 1.; // offset equal to 1 by default
        }
      }
      else{
        // Coalescence has occured at least once
        q1Thom->getRefToCoalescenceOnsetFlag() = true;
        // Update onset value
          // Porosity
        q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();
          // Geometrical parameters
        q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
        q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
        q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
          // Cft onset
        q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();


        // Determine which mode should be used for the next time step.
        if (q0Thom->getCoalescenceActiveFlag())
        {
          // If Thomason is/was used for this time step : check if Gurson is better
          double yieldGurson = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
          if(yieldGurson > _tol){
            Msg::Info("Change yield surface to Gurson fG =%e",yieldGurson);
            q1Thom->getRefToCoalescenceActiveFlag() = false; // Gurson will be used at the next time step
            q1Thom->getRefToAccelerateRate() = 1.0;
          }
          else{
            q1Thom->getRefToCoalescenceActiveFlag() = true;  // Thomason will be used at the next time step
            q1Thom->getRefToAccelerateRate() = 1.0;
          }

        }
        else
        {
          // If Gurson was/is used for this time step : check if Thomason is better
          double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
          if(yieldThomason > _tol ){
            Msg::Info("Change yield surface to Thomason fT = %e",yieldThomason);
            q1Thom->getRefToCoalescenceActiveFlag() = true; // Thomason will be used at the next time step
            q1Thom->getRefToAccelerateRate() = 1.0;
          }
          else{
            q1Thom->getRefToCoalescenceActiveFlag() = false;  // Gurson will be used at the next time step
            q1Thom->getRefToAccelerateRate() = 1.0;
          }
        }
      }
    }
    else
    {
      // Coalescence never occurs in the bulk
      q1Thom->getRefToCoalescenceOnsetFlag() = false;
      // Reset onset variables to their default value
        // Porosity
      q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
        // Geometrical parameters
      q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
      q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
      q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;

      q1Thom->getRefToCoalescenceActiveFlag() = false;
      q1Thom->getRefToCrackOffsetOnCft() = 1.;
      q1Thom->getRefToAccelerateRate() = 1.;
    }
  }
  else{
    // If crack transition is not used: no distinction between bulk and interface ipvs
    // if Thomason yield surface is first used, it always is used
    // onset and active are the same

    if (!q0Thom->getCoalescenceOnsetFlag()){
      double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
      // offset

      if (yieldThomason > _tol){
        // update offset
        if (_CfTOffsetMethod == 1){
          q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
          Msg::Info("coalescence occurs, use Thomason yield surface in next step with offset = %e",q1Thom->getCrackOffsetOnCft());
        }
        else if(_CfTOffsetMethod == 0){
          q1Thom->getRefToCrackOffsetOnCft() = 1.;
        }

        q1Thom->getRefToCoalescenceActiveFlag() = true;
        q1Thom->getRefToCoalescenceOnsetFlag() = true;
        //Msg::Info("coalescene yield surface is used");
        // initial coalesence onset value
        q1Thom->getRefToPorosityAtCoalescenceOnset() = yieldfV;
        q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
        q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
        q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q1Thom->getShapeFactor();

        q1Thom->getRefToAccelerateRate() = 1.;
      }
      else{
        q1Thom->getRefToCoalescenceActiveFlag() = false;
        q1Thom->getRefToCoalescenceOnsetFlag() = false;

        // reset to defaut value
        q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
        q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
        q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
        q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;

        q1Thom->getRefToCrackOffsetOnCft() = 1.;
        q1Thom->getRefToAccelerateRate() = 1.;
      }
    }
    else{
      q1Thom->getRefToCoalescenceActiveFlag() = true;
      q1Thom->getRefToCoalescenceOnsetFlag() = true;

      q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
      q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
      q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
      q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();

      q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
      q1Thom->getRefToAccelerateRate() = q0Thom->getAccelerateRate();
    }
  }
};

void mlawNonLocalPorousCoupledLaw::checkFailed(IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
   // Get ipvs
  const IPThomasonCoalescence* q0Thom = dynamic_cast<const IPThomasonCoalescence*>(&q0->getConstRefToIPCoalescence());
  IPThomasonCoalescence* q1Thom = dynamic_cast<IPThomasonCoalescence*>(&q1->getRefToIPCoalescence());
  if (q1Thom == NULL){
    Msg::Fatal("mlawNonLocalPorousCoupledLaw::checkFailed: IPThomasonCoalescence must be used!");
  }

  // Manage failure flag:
  // If it's previously failed, it always be failed. Otherwise, the criterion will depends on the active plastic mode
  if (q0->isFailed()){
    q1->setFailed(true);
  }
  else
  {
    q1->setFailed(false);
    // If coalescence plastic mode is currently used: criterion based on Cft
    double Cft = q1Thom->getConcentrationFactor();
    double CftOffest = q1Thom->getCrackOffsetOnCft();
    if (Cft*CftOffest < (1.-_failedTol)*5.){
      q1->setFailed(true);
      Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip has failed at Cft = %e", Cft);
    }

    // If coalescence plastic mode is currently used: criterion based on yield fv
    double yieldfV = q1->getYieldPorosity();
    if (yieldfV > _failedTol* _mlawGrowth->getFailureYieldPorosityValue()){
      q1->setFailed(true);
      Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip has failed at yieldfV = %e", yieldfV);
    }
  }
};


void mlawNonLocalPorousCoupledLaw::computeResidual_multipleNonLocalVariables(SVector3& res, STensor3& J,
                const double kcorEqpr, const double pcorpr, const double R, const double H,
                const double yieldfV,  const double DyieldfVDfV,
                const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
                const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP, // 3 unknonwn
                const double DeltafV, const double* T
              ) const{
  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
  if (coalesFlag){
    _mlawCoales->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
                                DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
  }
  else {
    _mlawGrowth->computeResidual_multipleNonLocalVariables(res,J,kcorEqpr,pcorpr,R,H,yieldfV,DyieldfVDfV,q0,q1,
                                DeltaHatD,DeltaHatQ,DeltaHatP,DeltafV,T);
  }
};

void mlawNonLocalPorousCoupledLaw::computeDResidual_multipleNonLocalVariables(STensor3 &dres0dEpr, STensor3 &dres1dEpr, STensor3 &dres2dEpr,
                              std::vector<double> &dres0dtildefV, std::vector<double> &dres1dtildefV, std::vector<double> &dres2dtildefV,
                              const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const double R,
                              const double yieldfV, const double DyieldfVDfV, const std::vector<double>& DyieldfVDtildefV,
                              const STensor3& DDeltafVDEpr, const STensor43& DdevKcorprDEpr, const STensor3& DpcorprDEpr,
                              const IPNonLocalPorosity *q0, const IPNonLocalPorosity* q1,
                              const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
                              const double* T, const double* DRDT,
                              const STensor3* DdevKcorprDT, const double* DpcorprDT, const double* DDeltafVDT,
                              double *dres0dT, double *dres1dT, double *dres2dT
                              ) const{
  bool coalesFlag = q0->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
  if (coalesFlag){
    _mlawCoales->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
                                  dres0dtildefV,dres1dtildefV,dres2dtildefV,
                                  devKcorpr,kcorEqpr,pcorpr,R,
                                  yieldfV,DyieldfVDfV,DyieldfVDtildefV,
                                  DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD,DeltaHatQ,DeltaHatP,
                                  T,DRDT,
                                  DdevKcorprDT ,DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }
  else {
    _mlawGrowth->computeDResidual_multipleNonLocalVariables(dres0dEpr,dres1dEpr,dres2dEpr,
                                  dres0dtildefV,dres1dtildefV,dres2dtildefV,
                                  devKcorpr,kcorEqpr,pcorpr,R,
                                  yieldfV,DyieldfVDfV,DyieldfVDtildefV,
                                  DDeltafVDEpr,DdevKcorprDEpr,DpcorprDEpr,
                                  q0, q1,
                                  DeltaHatD,DeltaHatQ,DeltaHatP,
                                  T,DRDT,
                                  DdevKcorprDT ,DpcorprDT,DDeltafVDT,
                                  dres0dT,dres1dT,dres2dT);
  }

};