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)
{
  _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)
{
  _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)
{
  _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;
  _mlawCoales->setCfTOffsetMethod(MethodNumber);
  if (_CfTOffsetMethod == 0 or _CfTOffsetMethod == 1 or _CfTOffsetMethod == 2 or _CfTOffsetMethod == 3){
    Msg::Info("mlawNonLocalPorousCoupledLaw::useCfTOffset: Offset method %d is used",_CfTOffsetMethod);
  }
  else{
    Msg::Fatal("mlawNonLocalPorousCoupledLaw::useCfTOffset: Offset method %d is not existing",_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::setFailureTolerance(const double NewFailureTol){
  if ((NewFailureTol > 1.) or (NewFailureTol < 0.)){
    Msg::Fatal("mlawNonLocalPorousCoupledLaw::setFailureTolerance: Uncorrect value of tolerance failure");
  }
  _failedTol = NewFailureTol;
  _mlawGrowth->setFailureTolerance(NewFailureTol);
  _mlawCoales->setFailureTolerance(NewFailureTol);
};

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::setStressTriaxialityFunction_kw(const scalarFunction& fct){
  mlawNonLocalPorosity::setStressTriaxialityFunction_kw(fct);
  _mlawGrowth->setStressTriaxialityFunction_kw(fct);
  _mlawCoales->setStressTriaxialityFunction_kw(fct);
};

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 = q1->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 = q1->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 = q1->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 = q1->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);
  }
}

double mlawNonLocalPorousCoupledLaw::getOnsetCriterion(const IPNonLocalPorosity* q1) const{
  return q1->getConstRefToIPCoalescence().getCrackOffsetOnCft();
}

void mlawNonLocalPorousCoupledLaw::nucleation(const double p1, const double p0, const double f0, IPNonLocalPorosity* q1, const IPNonLocalPorosity *q0) const{
  bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
  if (coalesFlag){
    _mlawCoales->nucleation(p1,p0,f0,q1,q0);
  }
  else{
    _mlawGrowth->nucleation(p1,p0,f0,q1,q0);
  }
};

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!");
  }

  // Get stress state
  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());

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

  // Coupling management
  if (_withCrackTransition){
    // If crack transition is used,
      // At the interface: coalescence is used, namely as crack insertion criterion
      // In the bulk, nothing is done (== always Gurson)
    q1Thom->getRefToYieldOffset() = 0.;
    
    
    if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS){
      // At the interface:
      // Check if transition has already occurs or not
      if (!q1Thom->getCoalescenceOnsetFlag()){
        // If Coalescence has not yet occured at least once
        // Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface

        // Get normal traction on the interface from tau stress
        double sigmax = 0.;
        static SVector3 Norm;
        Norm = q1->getConstRefToReferenceOutwardNormal();
        if (Norm.norm() > 0){ // is it usefulll ???
          Norm.normalize();
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              sigmax += Norm(i)*kCor(i,j)*Norm(j);
            }
          }
        }
        else{
          Msg::Fatal("mlawNonLocalPorousCoupledLaw::updatePlasticState: zero normal");
        }

        // Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
        // and manage offset
          // Why is the offset needed?
          // We want to ensure a continuity with the Thomason yield surface at coalescence onset
          // The offset correspond to a factor to obtain fT = 0 with the actual stress state
            // -actual stress state:  yieldThomason = sthingOfStress - Cft*R/R0
            // -correction to satidfy fT=0 with actual stress:  yieldThomason + Cft*R/R0 - CftOffset*Cft*R/R0  = 0
            // -solving for offset gives: CftOffset = (yieldThomason + Cft*R/R0)/(Cft*R/R0) = 1+yieldThomason*R0/(Cft*R)
            // NB: if _CfTOffsetMethod > 1: fT>0,
        // Compute Thomason yield surface
        double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
        
        // Offset management
        if(_CfTOffsetMethod == 2 or _CfTOffsetMethod == 3){
          // Method 2: Offset is always computed to ensure contuinity if coalescence is forced (offset can be <1)
          q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
        }

        if ((sigmax - Cft*R > _tol) and (yieldThomason > _tol)){
          // If coalescence occurs, update Cft offset value and others variables

          // Update onset coalescence variables
            // Activate coalescence onset
          q1Thom->getRefToCoalescenceOnsetFlag() = true;
            // Set onset variables to the actual values (yield porosity and geometrical parameters)
          q1Thom->getRefToPorosityAtCoalescenceOnset() = q1->getYieldPorosity();
          q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
          q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
          q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q1Thom->getShapeFactor();

          q1Thom->getRefToCoalescenceActiveFlag() = true;
          q1Thom->getRefToAccelerateRate() = 1.;

          // Offset management
          if (_CfTOffsetMethod == 1){
            // Method 1: offset computed at coalescence onset only
            q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
            Msg::Info("Coalescence occurs for the first time: Offset =%e",q1Thom->getRefToCrackOffsetOnCft());
          }
          else if (_CfTOffsetMethod == 0){
            // Method 0: Offset is not used: always equal to 1.
            q1Thom->getRefToCrackOffsetOnCft() = 1.;
          }


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

          q1Thom->getRefToCoalescenceActiveFlag() = false;
          q1Thom->getRefToAccelerateRate() = 1.;

          // Offset management
          if (_CfTOffsetMethod == 1 or _CfTOffsetMethod == 0){
            q1Thom->getRefToCrackOffsetOnCft() = 1.;
          }
        }
      }
      else{
        // Coalescence has occured at least once
        q1Thom->getRefToCoalescenceOnsetFlag() = true; // already done....
        
        // Update onset value (yield porosity and geometrical parameters)
        q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();
        q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
        q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
        q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
        
        // Cft onset - Method 0-1-2: keep it always constant once the coalescence occurs
        q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();


        // Determine which mode should be used for the next time step.
        if(!q0Thom->getCoalescenceOnsetFlag()){
          // First time that coalescence is detected: Thomason == ON
          // This case is only appearing with with selective update
          q1Thom->getRefToCoalescenceActiveFlag() = true;
          q1Thom->getRefToAccelerateRate() = 1.0;
          q1Thom->getRefToYieldOffset() = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
        }
        else{
          // Not the first time that coalescence appears
          if (q0Thom->getCoalescenceActiveFlag()){
            // If Thomason is/was used for this time step : check if Gurson is more restrictive
            double yieldGurson = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
            q1Thom->getRefToYieldOffset() = yieldGurson;
            if(yieldGurson > _tol){
              Msg::Info("Change yield surface to Gurson fG =%e",yieldGurson);
              q1Thom->getRefToCoalescenceActiveFlag() = true; // 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 more restrictive
            double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
            q1Thom->getRefToYieldOffset() = yieldThomason;
            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;
              
              if (_CfTOffsetMethod == 3){
                q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft() + q1Thom->getYieldOffset()*R0/(Cft*R);
                if (q1Thom->getRefToCrackOffsetOnCft() > 2.){Msg::Error("High predicted value of Cft Offset = %e",q1Thom->getRefToCrackOffsetOnCft());};
              }
              
            }
            else{
              q1Thom->getRefToCoalescenceActiveFlag() = false;  // Gurson will be used at the next time step
              q1Thom->getRefToAccelerateRate() = 1.0;
            }
          }
        }
      }
    }
    else
    {
      // In the bulk:
        // Coalescence never occurs
      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.;
      q1Thom->getRefToYieldOffset() = 0.;
    }
  }
  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);

      if (yieldThomason > _tol){
				q1Thom->getRefToCoalescenceCriterion() = 0.;
        // 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->getRefToCoalescenceCriterion() = yieldThomason; // save data
        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->getRefToCoalescenceCriterion() = 0.;
      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 yieldfV = q1->getYieldPorosity();
    double Chi = q1Thom->getLigamentRatio();
    double Cft = q1Thom->getConcentrationFactor();
    double CftOffest = q1Thom->getCrackOffsetOnCft();
    
    if (_CfTOffsetMethod==2 or _CfTOffsetMethod==3){
      if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffest = 1.;};
    }
    
    if (Cft*CftOffest < (1.-_failedTol)*5.){
      q1->setFailed(true);
      Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with Cft at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
    }
    else if (Chi > _failedTol){
      q1->setFailed(true);
      Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with Chi at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
    }
    else if (yieldfV > _failedTol*_mlawGrowth->getFailureYieldPorosityValue() ){
      q1->setFailed(true);
      Msg::Info("mlawNonLocalPorousCoupledLaw::checkFailed: an ip failed with fVy at Cft= %e, chi= %e, yieldfV= %e", Cft, Chi, yieldfV);
    }
  }
};

void mlawNonLocalPorousCoupledLaw::localPorosityGrowth(double& DeltafV,
                      const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
                      const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr,
                      const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1) const{
  bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
  if (coalesFlag){
    _mlawCoales->localPorosityGrowth(DeltafV,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,q0,q1);
  }
  else{
    _mlawGrowth->localPorosityGrowth(DeltafV,DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,q0,q1);
  }
};

// Local porosity growth (with der.)
void mlawNonLocalPorousCoupledLaw::localPorosityGrowth(double& DeltafV, double& DDeltafVDDeltaHatD,double& DDeltafVDDeltaHatQ, double& DDeltafVDDeltaHatP, STensor3& DDeltafVDEpr,
                  const double DeltaHatD, const double DeltaHatQ, const double DeltaHatP,
                  const STensor3& devKcorpr, const double kcorEqpr, const double pcorpr, const STensor43& DdevKcorprDEpr,
                  const IPNonLocalPorosity *q0, const IPNonLocalPorosity *q1,
                  const double* T, const STensor3* DdevKcorprDT, double* DDeltafVDT) const{
  bool coalesFlag = q1->getConstRefToIPCoalescence().getCoalescenceActiveFlag();
  if (coalesFlag){
    _mlawCoales->localPorosityGrowth(DeltafV,DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,DDeltafVDEpr,
                                    DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,DdevKcorprDEpr,q0,q1,
                                    T,DdevKcorprDT,DDeltafVDT);
  }
  else{
    _mlawGrowth->localPorosityGrowth(DeltafV,DDeltafVDDeltaHatD,DDeltafVDDeltaHatQ,DDeltafVDDeltaHatP,DDeltafVDEpr,
                                    DeltaHatD,DeltaHatQ,DeltaHatP,devKcorpr,kcorEqpr,pcorpr,DdevKcorprDEpr,q0,q1,
                                    T,DdevKcorprDT,DDeltafVDT);
  }
};


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 = q1->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 = q1->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);
  }

};