Skip to content
Snippets Groups Projects
Select Git revision
  • 88d91b7992485482273861c082f31a4c29e6bc98
  • master default protected
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

nodalBasis.h

Blame
  • mlawHyperelastic.cpp 78.77 KiB
    //
    // C++ Interface: material law
    //
    // Author:  <Van Dung Nguyen>, (C) 2014
    //
    // Copyright: See COPYING file that comes with this distribution
    //
    //
    #include "mlawHyperelastic.h"
    #include "STensorOperations.h"
    #include "nonLinearMechSolver.h"
    
    void mlawHyperViscoElastic::setStrainOrder(const int order){
      _order = order;
      Msg::Info("order = %d is used to approximate log and exp");
    };
    
    void mlawHyperViscoElastic::setViscoelasticMethod(const int method){
      _viscoMethod = (viscoelasticType)method;
      if (_viscoMethod == Maxwell)
        Msg::Info("generalized maxwell model is used for viscoelasticity");
      else if (_viscoMethod == KelvinVoight)
        Msg::Info("generalized Voigt-Kelvin model is used for viscoelasticity");
      else
        Msg::Error("this method has not been implemented");
    };
    
    void mlawHyperViscoElastic::setViscoElasticNumberOfElement(const int N){
      _N = N;
      Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
      _Ki.clear(); _ki.clear();
      _Gi.clear(); _gi.clear();
      _Ki.resize(_N,0.);
      _ki.resize(_N,0.);
      _Gi.resize(_N,0.);
      _gi.resize(_N,0.);
    };
    
    void mlawHyperViscoElastic::setViscoElasticData(const int i, const double Ei, const double taui){
      if (i> _N or i<1)
        Msg::Error("This setting is invalid %d > %d",i,_N);
      else{
        double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
        double KK = Ei/(3.*(1.-2.*_nu));
        double GG = Ei/(2.*(1.+_nu));
        _Ki[i-1] = KK;
        _ki[i-1] = taui;
        _Gi[i-1] = GG;
        _gi[i-1] = taui;
    
        Msg::Info("setting: element=%d Ki= %e ki = %e, Gi=%e, gi=%e",i-1,KK,taui,GG,taui);
      }
    };
    
    void mlawHyperViscoElastic::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
      if (i> _N or i<1)
        Msg::Error("This setting is invalid %d > %d",i,_N);
      else{
        _Ki[i-1] = Ki;
        _ki[i-1] = ki;
      //  Msg::Info("setting: element=%d Ki= %e ki = %e, ",i-1,Ki,ki);
      }
    };
    void mlawHyperViscoElastic::setViscoElasticData_Shear(const int i, const double Gi, const double gi){
      if (i> _N or i<1)
        Msg::Error("This setting is invalid %d > %d",i,_N);
      else{
        _Gi[i-1] = Gi;
        _gi[i-1] = gi;
    
     //   Msg::Info("setting: element=%d: Gi=%e, gi=%e",i-1,Gi,gi);
      }
    };
    
    void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
      FILE* file = fopen(filename.c_str(),"r");
      if (file == NULL) Msg::Error("file: %s does not exist",filename.c_str());
      _Ki.clear(); _ki.clear();
      _Gi.clear(); _gi.clear();
      _N = 0;
      Msg::Info("reading viscoelastic input data");
    
      while (!feof(file)){
        int ok = fscanf(file,"%d",&_N);
        Msg::Info("Numer of Maxwell elements: %d",_N);
        double KK, k, GG, g;
    
        for (int i=0; i<_N; i++){
          ok = fscanf(file,"%lf %lf %lf %lf",&KK,&k,&GG,&g);
          _Ki.push_back(KK);
          _ki.push_back(k);
          _Gi.push_back(GG);
          _gi.push_back(g);
          Msg::Info("Maxwell element %d: K[%d] = %e, k[%d] = %e, G[%d] = %e, g[%d] = %e",
                    i,i,KK,i,k,i,GG,i,g);
        }
      }
      fclose(file);
    };
    
    void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) const
    {
      int method =1;
      if(method ==0)
      {
        A_v = getVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection());
      
        dA_vdE = getVolumeCorrection()*getZetaVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection()-1.)*exp(getXiVolumeCorrection()/3.*tr*tr) *getXiVolumeCorrection()*2./3.*tr;
      
      
        B_d = getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection());
        STensorOperation::zero(dB_vddev);
    
     
        dB_vddev=dev;
        dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection()-1.)* exp(getThetaDevCorrection()*dev.dotprod())*getThetaDevCorrection();
      }
      else
      {
        A_v = getVolumeCorrection()*(tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())+tanh(getZetaVolumeCorrection()));  
        dA_vdE = getVolumeCorrection()*getXiVolumeCorrection()*2./3.*(1.-tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())*tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection()));
      
      
        B_d = getDevCorrection()*(tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())+tanh(getPiDevCorrection()));
        STensorOperation::zero(dB_vddev);
     
        dB_vddev=dev;
        dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*(1.-tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())*tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection()) );
        
      }
    }
    
    
    double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) const
    {
      double Psy = 0.;
      if (_viscoMethod == Maxwell)
      {
        double trEe;
        static STensor3 devEe;
        STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
        
        double Av, dAv,B_d;
        static STensor3 dB_d;
        evaluatePhiPCorrection(trEe, devEe, Av, dAv, B_d, dB_d);
        Psy = getUpdatedBulkModulus(&q)*Av*(0.5*trEe*trEe)+getUpdatedShearModulus(&q)*B_d*STensorOperation::doubledot(devEe,devEe); //this is not correct we should do an integral
        
    
        for (int i=0; i<_N; i++){
          static STensor3 devEebranch;
          Psy += _Ki[i]*0.5*(q._B[i])*(q._B[i])+_Gi[i]*STensorOperation::doubledot(q._A[i],q._A[i]);
        }
      }
      else if (_viscoMethod == KelvinVoight)
      {
        // need to be complete
        Psy = 0.5*STensorOperation::doubledot(q._Ee,q._kirchhoff);
        for (int i=0; i<_N; i++){
          static STensor3 Ev;
          q.getViscoElasticStrain(i+1,Ev);
          Psy += 0.5*STensorOperation::doubledot(Ev,q._kirchhoff);
        }
      }
      else{
        Msg::Error("visco elastic method %d has not been defined");
      }
      return Psy;;
    }
    
    void mlawHyperViscoElastic::dDeformationEnergydF(const IPHyperViscoElastic& q,  const STensor3 &P, STensor3 &dElasticEnergydF) const
    {
       STensorOperation::zero(dElasticEnergydF);
       dElasticEnergydF=P;
    }
    
    double mlawHyperViscoElastic::viscousEnergy(const IPHyperViscoElastic& q0,const IPHyperViscoElastic& q) const
    {
      double DViscous = 0.;
      if (_viscoMethod == Maxwell)
      {
        double trEe, trEe0;
        static STensor3 devEe, devEe0;
        STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
        STensorOperation::decomposeDevTr(q0._Ee,devEe0,trEe0);
        for (int i=0; i<_N; i++){
          static STensor3 dqdev;
          dqdev = devEe;
          dqdev -= q._A[i];
          dqdev -= devEe0;
          dqdev += q0._A[i];
          dqdev*=2.*_Gi[i];
          double dtrq= trEe-q._B[i]-trEe0+q0._B[i];
          dtrq*=_Ki[i]*3.;
          DViscous+=  STensorOperation::doubledot(q._A[i],dqdev);
          DViscous+= q._B[i]*dtrq/3.;
        }
      }
      //else
      //{
      //  Msg::Error("visco elastic method %d has not been defined");
      //}
      return DViscous;
    
    }
    
    
    
    void mlawHyperViscoElastic::dViscousEnergydF(const IPHyperViscoElastic& q0, const IPHyperViscoElastic& q,  const STensor43 &dlnCdC, const STensor3 &Fe, STensor3 &dViscousEnergydF) const
    {
      static STensor3 dDeltaViscousdEve;
      STensorOperation::zero(dDeltaViscousdEve);
      STensorOperation::zero(dViscousEnergydF);
      if (_viscoMethod == Maxwell)
      {
        double trEe, trEe0;
        static STensor3 devEe, devEe0;
        STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
        STensorOperation::decomposeDevTr(q0._Ee,devEe0,trEe0);
        double dt=this->getTimeStep();
        for (int i=0; i<_N; i++)
        {
          double dtg = dt/(_gi[i]);
          double ratiog=exp(-dtg/2.);
          double dtk = dt/(_ki[i]);
          double ratiok=exp(-dtk/2.);
    
          static STensor3 dqdev, GA2;
          dqdev = devEe;
          dqdev -= q._A[i];
          dqdev -= devEe0;
          dqdev += q0._A[i];
          dqdev*=2.*_Gi[i];
          double dtrq= trEe-q._B[i]-trEe0+q0._B[i];
          dtrq*=_Ki[i]*3.;
          GA2=q._A[i];
          GA2*=2.*_Gi[i];
          dDeltaViscousdEve=dqdev;
          dDeltaViscousdEve-=GA2;
          dDeltaViscousdEve*=ratiog;
          dDeltaViscousdEve+=GA2;
          double ppart=(dtrq/3.-q._B[i]*_Ki[i])*ratiok+q._B[i]*_Ki[i];
          for (int k=0; k<3; k++){
             dDeltaViscousdEve(k,k)+=ppart;
          }
        }
        //HERE WE NEED TO GET WITH RESPECT TO DF -> dlog
      }
      //else
      //{
      //  Msg::Error("visco elastic method %d has not been defined");
      //}
      for(int i=0.; i<3.;i++)
      {
        for(int J=0.; J<3.;J++)
        {
          for(int K=0.; K<3.;K++)
          {
            for(int L=0.; L<3.;L++)
            {
              for(int N=0.; N<3.;N++)
              {
                dViscousEnergydF(i,J)+=0.5*dDeltaViscousdEve(K,L)*(dlnCdC(K,L,J,N)*Fe(i,N)+dlnCdC(K,L,N,J)*Fe(i,N));
              }
            }
          }
        }
      }
    }
    
    
    mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
                              const bool matrixbyPerturbation, const double pert):
        materialLaw(num,true), _rho(rho),_tangentByPerturbation(matrixbyPerturbation),_perturbationfactor(pert),
        _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.), 
        _volCorrection(0.), _xivolCorrection(1.), _zetavolCorrection(1.), _devCorrection(0.), _thetadevCorrection(1.), _pidevCorrection(1.) 
    {
       _Idev = _I4;
       STensor3 mIon3(-1./3);
       STensor43 mIIon3;
       tensprod(_I,mIon3, mIIon3);
       _Idev += mIIon3;
    
      double _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
    
      _mu = 0.5*E/(1.+nu);
      _K = E/3./(1.-2.*nu);
    
      Cel=0.;
      Cel(0,0,0,0) = _lambda + 2.*_mu;
      Cel(1,1,0,0) = _lambda;
      Cel(2,2,0,0) = _lambda;
      Cel(0,0,1,1) = _lambda;
      Cel(1,1,1,1) = _lambda + 2.*_mu;
      Cel(2,2,1,1) = _lambda;
      Cel(0,0,2,2) = _lambda;
      Cel(1,1,2,2) = _lambda;
      Cel(2,2,2,2) = _lambda + 2.*_mu;
    
      Cel(1,0,1,0) = _mu;
      Cel(2,0,2,0) = _mu;
      Cel(0,1,0,1) = _mu;
      Cel(2,1,2,1) = _mu;
      Cel(0,2,0,2) = _mu;
      Cel(1,2,1,2) = _mu;
    
      Cel(0,1,1,0) = _mu;
      Cel(0,2,2,0) = _mu;
      Cel(1,0,0,1) = _mu;
      Cel(1,2,2,1) = _mu;
      Cel(2,0,0,2) = _mu;
      Cel(2,1,1,2) = _mu;
    };
    mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): materialLaw(src),
          _rho(src._rho),
          _mu(src._mu),_K(src._K),Cel(src.Cel),
          _order(src._order),
          _tangentByPerturbation(src._tangentByPerturbation),_perturbationfactor(src._perturbationfactor),
          _N(src._N),_viscoMethod(src._viscoMethod),_Ki(src._Ki),_ki(src._ki),_Gi(src._Gi),_gi(src._gi),
          _I4(src._I4), _I(src._I), _Idev(src._Idev), _volCorrection(src._volCorrection), _xivolCorrection(src._xivolCorrection), 
          _zetavolCorrection(src._zetavolCorrection), _devCorrection(src._devCorrection), _thetadevCorrection(src._thetadevCorrection), 
          _pidevCorrection(src._pidevCorrection)
    {
    
    };
    
    mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &source)
    {
      materialLaw::operator=(source);
      const mlawHyperViscoElastic* src =dynamic_cast<const mlawHyperViscoElastic*>(&source);
      if(src != NULL)
      {
          _rho = src->_rho,
          _mu = src->_mu; _K = src->_K; Cel = src->Cel;
          _order = src->_order;
          _tangentByPerturbation = src->_tangentByPerturbation; _perturbationfactor = src->_perturbationfactor;
          _N = src->_N; _viscoMethod = src->_viscoMethod; _Ki = src->_Ki; _ki = src->_ki; _Gi = src->_Gi; _gi = src->_gi;
          _zetavolCorrection=src->_zetavolCorrection; _devCorrection=src->_devCorrection;
           _thetadevCorrection=src->_thetadevCorrection; _pidevCorrection=src->_pidevCorrection;
      }
      return *this;
    }
    
    double mlawHyperViscoElastic::soundSpeed() const
    {
      double _E=9*_K*_mu/(3*_K+_mu);
      double _nu = (3.*_K-2.*_mu)/2./(3.*_K+_mu);
      double factornu = (1.-_nu)/((1.+_nu)*(1.-2.*_nu));
      return sqrt(_E*factornu/_rho);
    }
    
    void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1, double& Ke, double & Ge) const{
      if ((_Ki.size() > 0) or (_Gi.size() > 0)){
    
        double dt = this->getTimeStep();
        if (_viscoMethod == Maxwell){
          static STensor3 DE, devDE;
          static double trDE;
    
          DE = q1->_Ee;
          DE -= q0->_Ee;
          STensorOperation::decomposeDevTr(DE,devDE,trDE);
    
          // maxwell
          Ge = getUpdatedShearModulus(q1);
          Ke = getUpdatedBulkModulus(q1);
    
          for (int i=0; i<_Gi.size(); i++){
            double dtg = dt/_gi[i];
            double expmdtg = exp(-dtg);
            double ztag = exp(-dtg/2.);
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                q1->_A[i](k,l) = expmdtg*q0->_A[i](k,l) + ztag*devDE(k,l);
              }
            }
            Ge += _Gi[i]*ztag;
          }
          for (int i=0; i<_Ki.size(); i++){
            double dtk = dt/_ki[i];
            double expmdtk = exp(-dtk);
            double ztak = exp(-dtk/2.);
            q1->_B[i] = q0->_B[i]*expmdtk +ztak*trDE;
            Ke += _Ki[i]*ztak;
          }
        }
        else if (_viscoMethod == KelvinVoight){
          static STensor3 DK, devDK;
          static double trDK;
    
          DK = q1->_kirchhoff;
          DK -= q0->_kirchhoff;
    
          STensorOperation::decomposeDevTr(DK,devDK,trDK);
    
          /*update internal variable from stress increment*/
          for (int i=0; i< _Gi.size(); i++)
          {
            STensorOperation::zero(q1->_A[i]);
            double dtg = dt/(_gi[i]);
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                q1->_A[i](k,l) += exp(-dtg)*q0->_A[i](k,l) + exp(-dtg/2.)*devDK(k,l)/(2.*_Gi[i]);
              }
            }
          }
    
          for (int i=0; i< _Ki.size(); i++){
            q1->_B[i] = 0.;
            double dtk = dt/(_ki[i]);
            q1->_B[i] += exp(-dtk)*q0->_B[i] + exp(-dtk/2.)*trDK/(3.*_Ki[i]);
          }
        }
        else{
          Msg::Error("visco elastic method %d has not been defined");
        }
      }
    };
    void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
              const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
              double& Ke, double& Ge,     double& A_v,  double& B_d, double& dApr, STensor3& dB_vddev) const{
      if ((_Ki.size() > 0) or (_Gi.size() > 0)){
        static STensor3 DE, devDE;
        static double trDE;
        DE =  Ee;
        DE -= Ee0;
        STensorOperation::decomposeDevTr(DE,devDE,trDE);
    
        double dt = this->getTimeStep();
        if (_viscoMethod == Maxwell){
          Ge = getUpdatedShearModulus(q1);
          Ke = getUpdatedBulkModulus(q1);
          for (int i=0; i<_Gi.size(); i++){
            double dtg = dt/_gi[i];
            double expmdtg = exp(-dtg);
            double ztag = exp(-dtg/2.);
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                q1->_A[i](k,l) = expmdtg*q0->_A[i](k,l) + ztag*devDE(k,l);
              }
            }
            Ge += _Gi[i]*ztag;
          }
          for (int i=0; i<_Ki.size(); i++){
            double dtk = dt/_ki[i];
            double expmdtk = exp(-dtk);
            double ztak = exp(-dtk/2.);
            q1->_B[i] = q0->_B[i]*expmdtk +ztak*trDE;
            Ke += _Ki[i]*ztak;
          }
    
          static STensor3 devK, devE;
          STensorOperation::zero(devK);
          STensorOperation::zero(devE);
          
          double p, trEe;
          STensorOperation::decomposeDevTr(Ee,devE,trEe);
    
          evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
    
          devK = devE;
          devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d));  // deviatoric part
          
          
          p  = trEe;
          p *= getUpdatedBulkModulus(q1)*(1.+A_v); // pressure
    
          for (int i=0; i<_Gi.size(); i++){
            devK.daxpy(q1->_A[i],2.*_Gi[i]);
          }
          for (int i=0; i<_Ki.size(); i++){
            p += q1->_B[i]*_Ki[i];
          }
          
          q1->_kirchhoff = devK;
          q1->_kirchhoff(0,0) += p;
          q1->_kirchhoff(1,1) += p;
          q1->_kirchhoff(2,2) += p;
          
        }
        
        else if (_viscoMethod == KelvinVoight){
          double invGe = 1./getUpdatedShearModulus(q1);
          STensor3 D(0.);
          for (int i=0; i<_Gi.size(); i++){
            double dtg = dt/(_gi[i]);
            invGe += (1.-exp(-dtg/2.))/_Gi[i];
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                D(k,l) += q0->_A[i](k,l)*(exp(-dtg)-1.);
              }
            }
          }
          Ge = 1./invGe;
    
          double invKe = 1./_K;
          double V= 0.;
          for (int i=0; i<_Ki.size(); i++){
            double dtk = dt/(_ki[i]);
            invKe += (1.-exp(-dtk/2))/_Ki[i];
            V += q0->_B[i]*(exp(-dtk)-1.);
          }
          Ke = 1./invKe;
    
          // stress increment
          static STensor3 DdevK;
          DdevK = devDE; // dev corotational kirchoff stress predictor
          DdevK += D;
          DdevK *= (2.*Ge);
          double Dp = Ke*(trDE+ V); // pressure predictor
    
          /*update internal variable from stress increment*/
          for (int i=0; i< _Gi.size(); i++){
            STensorOperation::zero(q1->_A[i]);
            double dtg = dt/(_gi[i]);
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                q1->_A[i](k,l) += exp(-dtg)*q0->_A[i](k,l) + exp(-dtg/2.)*DdevK(k,l)/(2.*_Gi[i]);
              }
            }
          }
    
          for (int i=0; i< _Ki.size(); i++){
            q1->_B[i] = 0.;
            double dtk = dt/(_ki[i]);
            q1->_B[i] += exp(-dtk)*q0->_B[i] + exp(-dtk/2.)*Dp/(_Ki[i]);
          }
    
           // corotational Kirchhoff stress tenor from previous and increment
          STensor3& corK = q1->_kirchhoff;
          corK = q0->_kirchhoff;
    
          corK += DdevK;
          corK(0,0) += Dp;
          corK(1,1) += Dp;
          corK(2,2) += Dp;
        }
        else{
          Msg::Error("visco elastic method %d has not been defined");
        }
      }
      else{
        static STensor3 devK, devE;
        double p, trEe;
        STensorOperation::zero(devK);
        STensorOperation::zero(devE);
          
        STensorOperation::decomposeDevTr(Ee,devE,trEe);
    
        STensorOperation::zero(dB_vddev);
        evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
                
        devK = devE;
        devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d));  // deviatoric part
          
          
        p  = trEe;
        p *= getUpdatedBulkModulus(q1)*(1.+A_v); // pressure
        
        Ke = getUpdatedBulkModulus(q1);
        Ge = getUpdatedShearModulus(q1);
    
    
        q1->_kirchhoff = devK;
        q1->_kirchhoff(0,0) += p;
        q1->_kirchhoff(1,1) += p;
        q1->_kirchhoff(2,2) += p;
      }
    };
    
    
    void mlawHyperViscoElastic::isotropicHookTensor(const double K, const double G, STensor43& L) const{
      double lambda = K - 2.*G/3.;
      static STensor3 I(1.);
      for (int i=0; i<3; i++){
        for (int j=0; j<3; j++){
          for (int k=0; k<3; k++){
            for (int l=0; l<3; l++){
              L(i,j,k,l) = lambda*I(i,j)*I(k,l)+ G*(I(i,k)*I(j,l)+I(i,l)*I(j,k));
            }
          }
        }
      }
    };
    
    void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
                                      const bool stiff, STensor43& Tangent) const{
      static STensor3 C;
      static STensor43 dlnCdC;
      static STensor63 ddlnCddC;
    
      STensor3& E = q1->getRefToElasticStrain();
      STensorOperation::multSTensor3FirstTranspose(F,F,C);
      if (_order == 1){
        STensorOperation::logSTensor3(C,_order,E,&dlnCdC);  // as ddlogCddC = 0
      }
      else{
        STensorOperation::logSTensor3(C,_order,E,&dlnCdC,&ddlnCddC);
      }
      E *= 0.5; // strain
      static STensor3 devE;
      double trE;
      STensorOperation::decomposeDevTr(E,devE,trE);
      
      double Ke, Ge;
      double A_v_pr=0., B_d_pr = 0., dApr=0.;
      static STensor3 dB_devpr;
      STensorOperation::zero(dB_devpr);
      viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
    
      const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
      static STensor3 secondPK;
      STensorOperation::multSTensor3STensor43(corKir,dlnCdC,secondPK);
      STensorOperation::multSTensor3(F,secondPK,P);
      // first PK
      
      
      
      q1->_elasticEnergy=deformationEnergy(*q1);
      q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
      if (stiff){
        static STensor43 DsecondPKdC;
        static STensor43 DcorKDE;
        isotropicHookTensor(Ke+getUpdatedBulkModulus(q1)*A_v_pr,Ge+getUpdatedShearModulus(q1)*B_d_pr,DcorKDE);
        
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                DsecondPKdC(i,j,k,l) = 0.;
                for (int p=0; p<3; p++){
                  for (int q=0; q<3; q++){
                    if (_order != 1){
                      DsecondPKdC(i,j,k,l) += corKir(p,q)*ddlnCddC(p,q,i,j,k,l);
                    }
                    for (int r=0; r<3; r++){
                      for (int s=0; s<3; s++){
                        DsecondPKdC(i,j,k,l) += 0.5*(DcorKDE(p,q,r,s)+2.*getUpdatedShearModulus(q1)*devE(p,q)*dB_devpr(r,s)+ getUpdatedBulkModulus(q1)*dApr*trE*_I(p,q)*_I(r,s))*dlnCdC(p,q,i,j)*dlnCdC(r,s,k,l); 
                      }
                    }
                  }
                }
              }
            }
          }
        }
        STensorOperation::zero(Tangent);
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int p=0; p<3; p++){
              Tangent(i,j,i,p) += secondPK(p,j);
              for (int q=0; q<3; q++){
                for (int s=0; s<3; s++){
                  for (int k=0; k<3; k++){
                    Tangent(i,j,p,q) += 2.*F(i,k)*DsecondPKdC(k,j,q,s)*F(p,s);
                  }
                }
              }
            }
          }
        }
        dDeformationEnergydF(*q1,  P, q1->getRefToDElasticEnergyPartdF());
        dViscousEnergydF(*q0,*q1,dlnCdC, F, q1->getRefToDViscousEnergyPartdF());
      };
    
    };
    
    void mlawHyperViscoElastic::constitutive(const STensor3& F0,
                                             const STensor3& F, STensor3&P,
                                             const IPVariable *q0i,
                                             IPVariable *q1i,
                                             STensor43& Tangent,
                                             const bool stiff,
                                             STensor43* elasticTangent,
                                             const bool dTangent,
                                             STensor63* dCalgdeps) const{
    
      const IPHyperViscoElastic *q0=dynamic_cast<const IPHyperViscoElastic *>(q0i);
      IPHyperViscoElastic *q1 = dynamic_cast<IPHyperViscoElastic *>(q1i);
      if(elasticTangent!=NULL) Msg::Error("mlawHyperViscoElastic elasticTangent not defined");
    
      if (!_tangentByPerturbation){
        this->predictorCorrector_ViscoElastic(F0,F,P,q0,q1,stiff,Tangent);
      }
      else{
        this->predictorCorrector_ViscoElastic(F0,F,P,q0,q1,false,Tangent);
        if (stiff){
          static STensor3 Fplus, Pplus;
          static STensor43 Ttemp;
          static IPHyperViscoElastic qtemp(*q1);
          for (int i= 0; i<3; i++){
            for (int j=0; j<3; j++){
              Fplus = F;
              Fplus(i,j) += _perturbationfactor;
              this->predictorCorrector_ViscoElastic(F0,Fplus,Pplus,q0,&qtemp,false,Ttemp);
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  Tangent(k,l,i,j) = (Pplus(k,l)-P(k,l))/_perturbationfactor;
                }
              }
              q1->getRefToDElasticEnergyPartdF()(i,j)=(qtemp.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
              q1->getRefToDViscousEnergyPartdF()(i,j)=(qtemp.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
    
            }
          }
        }
      }
    };
    
    
    void mlawPowerYieldHyper::setCompressionHardening(const J2IsotropicHardening& comp){
      if (_compression) delete _compression;
      _compression = comp.clone();
    };
    void mlawPowerYieldHyper::setTractionHardening(const J2IsotropicHardening& trac){
      if (_traction) delete _traction;
      _traction = trac.clone();
    };
    void mlawPowerYieldHyper::setKinematicHardening(const kinematicHardening& kin){
      if (_kinematic) delete _kinematic;
      _kinematic = kin.clone();
    };
    void mlawPowerYieldHyper::setViscosityEffect(const viscosityLaw& v, const double p){
      _p = p;
      if (_viscosity) delete _viscosity;
      _viscosity = v.clone();
    };
    
    void mlawPowerYieldHyper::hardening(const IPHyperViscoElastoPlastic* q0, IPHyperViscoElastoPlastic* q) const{
      //Msg::Error("epspCompression = %e, epspTRaction = %e, epspShear = %e",q->_epspCompression,q->_epspTraction,q->_epspShear);
      if (_compression != NULL && q->_ipCompression != NULL){
        _compression->hardening(q0->_epspCompression, *q0->_ipCompression, q->_epspCompression,*q->_ipCompression);
      }
    
      if (_traction!= NULL && q->_ipTraction != NULL){
        _traction->hardening(q0->_epspTraction,*q0->_ipTraction, q->_epspTraction,*q->_ipTraction);
      }
    
      if (_kinematic!= NULL && q->_ipKinematic != NULL)
        _kinematic->hardening(q0->_epspbarre,*q0->_ipKinematic,q->_epspbarre,*q->_ipKinematic);
    };
    
    void mlawPowerYieldHyper::tangent_full_perturbation(
                                 STensor43 &T_,
                                 STensor43& dFedF,
                                 STensor43& dFpdF,
                                 const STensor3 &P,
                                 const STensor3 &F,
                                 const IPHyperViscoElastoPlastic* q0,
                                 IPHyperViscoElastoPlastic* q1
                               ) const{
    
      static STensor43 tmpSTensor43;
      static STensor3 Fplus, Pplus;
      static IPHyperViscoElastoPlastic q11(*q0);
    
      for (int i=0; i<3; i++){
        for (int j=0; j<3; j++){
          Fplus = F;
          Fplus(i,j)+=_perturbationfactor;
          this->predictorCorrector(Fplus,q0,&q11,Pplus,false,tmpSTensor43,tmpSTensor43,tmpSTensor43,NULL);
          q1->_DgammaDF(i,j) = (q11._epspbarre - q1->_epspbarre)/_perturbationfactor;
          q1->_DirreversibleEnergyDF(i,j) = (q11._irreversibleEnergy - q1->_irreversibleEnergy)/_perturbationfactor;
          for (int k=0; k<3; k++){
            for (int l=0; l<3; l++){
              T_(k,l,i,j) = (Pplus(k,l)-P(k,l))/_perturbationfactor;
              dFpdF(k,l,i,j) = (q11._Fp(k,l)-q1->_Fp(k,l))/_perturbationfactor;
              dFedF(k,l,i,j) = (q11._Fe(k,l)-q1->_Fe(k,l))/_perturbationfactor;
            }
          }
          q1->getRefToDElasticEnergyPartdF()(i,j)=(q11.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
          q1->getRefToDViscousEnergyPartdF()(i,j)=(q11.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
          q1->getRefToDPlasticEnergyPartdF()(i,j)=(q11.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/_perturbationfactor;
    
        }
      }
    };
    
    
    void mlawPowerYieldHyper::constitutive(
      const STensor3& F0,         // initial deformation gradient (input @ time n)
      const STensor3& Fn,         // updated deformation gradient (input @ time n+1)
      STensor3 &P,                // updated 1st Piola-Kirchhoff stress tensor (output)
      const IPVariable *q0i,       // array of initial internal variable
      IPVariable *q1i,             // updated array of internal variable (in ipvcur on output),
      STensor43 &Tangent,         // constitutive tangents (output)
      const bool stiff,          // if true compute the tangents
      STensor43* elasticTangent,
      const bool dTangent,
      STensor63* dCalgdeps) const{
    
      const IPHyperViscoElastoPlastic *q0=dynamic_cast<const IPHyperViscoElastoPlastic *>(q0i);
      IPHyperViscoElastoPlastic *q1 = dynamic_cast<IPHyperViscoElastoPlastic *>(q1i);
    
      static STensor43 dFedF, dFpdF;
    
      if (_tangentByPerturbation){
        this->predictorCorrector(Fn,q0,q1,P,false,Tangent,dFedF,dFpdF,elasticTangent);
        if (stiff)
          this->tangent_full_perturbation(Tangent,dFedF,dFpdF,P,Fn,q0,q1);
      }
      else{
        this->predictorCorrector(Fn,q0,q1,P,stiff,Tangent,dFedF,dFpdF,elasticTangent);
      }
    };
    
    void mlawPowerYieldHyper::updateEqPlasticDeformation(IPHyperViscoElastoPlastic *q1, const IPHyperViscoElastoPlastic *q0,
                                                const double& nup, const double& Dgamma) const{
      q1->_epspbarre = q0->_epspbarre+ Dgamma;
      q1->_epspCompression = q0->_epspCompression+ Dgamma;
      q1->_epspTraction = q0->_epspTraction+ Dgamma;
      double k = 1./(sqrt(1.+2.*nup*nup));
      q1->_epspShear = q0->_epspShear+ Dgamma/(k*sqrt(2.));
    };
    
    
    void mlawPowerYieldHyper::getYieldCoefficients(const IPHyperViscoElastoPlastic *q, fullVector<double>& coeffs) const{
      double sigc = q->_ipCompression->getR();
      double sigt = q->_ipTraction->getR();
    
      double m = sigt/sigc;
    
      //Msg::Error("m = %e, _n = %e",m,_n);
      //Msg::Error("sigc = %e, sigt = %e",sigc,sigt);
    
      coeffs.resize(3);
      coeffs(2) = pow(sigc,-_n);
      coeffs(1) = 3.*(pow(m,_n)-1.)/(m+1.)/sigc;
      coeffs(0) = (pow(m,_n)+m)/(m+1);
    };
    
    void mlawPowerYieldHyper::getYieldCoefficientDerivatives(const IPHyperViscoElastoPlastic *q, const double& nup, fullVector<double>& Dcoeffs) const{
      double sigc(0.), Hc(0.);
      sigc = q->_ipCompression->getR();
      Hc = q->_ipCompression->getDR();
    
      double sigt(0.), Ht(0.);
      sigt = q->_ipTraction->getR();
      Ht = q->_ipTraction->getDR();
    
    
      Dcoeffs.resize(3);
      double m = sigt/sigc;
      double Dm = (Ht*sigc- sigt*Hc)/(sigc*sigc);
      double Da1Dm = 3./sigc*(_n*pow(m,_n-1.)/(m+1.) - (pow(m,_n)-1.)/(m+1.)/(m+1.));
    
      Dcoeffs(2) = -_n*pow(sigc,-_n-1.)*Hc;
      Dcoeffs(1) = Da1Dm*Dm -3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*Hc;
      Dcoeffs(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*Dm;
    
    };
    
    void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{
      /* compute elastic predictor */
      STensor3& Fp1 = q1->_Fp;
      const STensor3& Fp0 = q0->_Fp;
    
      Fp1 = Fp0; // plastic deformation tensor
      q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
      q1->_epspCompression = q0->_epspCompression;
      q1->_epspTraction = q0->_epspTraction;
      q1->_epspShear = q0->_epspShear;
      q1->_backsig = q0->_backsig; // backstress tensor
      q1->_DgammaDt = 0.;
    
    
      static STensor3 Fpinv, Ce, Fepr;
      STensorOperation::inverseSTensor3(Fp1,Fpinv);
      STensorOperation::multSTensor3(F,Fpinv,Fepr);
      STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
    
      static STensor3 invFp0; // plastic predictor
      invFp0= Fpinv;
      STensor3& Fe = q1->_Fe;
      Fe = Fepr;
    
      static STensor43 DlnDCepr, DlnDCe;
      static STensor63 DDlnDDCe;
      static STensor43 dexpAdA; // estimation of dexpA/dA
    
      STensor3& Ee = q1->_Ee;
      STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
      Ee *= 0.5;
      DlnDCe = DlnDCepr;
    
      // update A, B
      double Ge, Ke, trEepr;
      static STensor3 devEepr;
      STensorOperation::zero(devEepr);
      
      
      
      
      
      
      
      double A_v_pr, B_d_pr, dApr;
      static STensor3 dB_devpr;
      STensorOperation::zero(dB_devpr);
      viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
      
      STensorOperation::decomposeDevTr(Ee,devEepr,trEepr); 
      //std::cout << "trEepr = "<<trEepr <<std::endl;
      static STensor3 PhiPr, PhiPr_Vol, PhiPr_Dev;
      PhiPr = q1->_kirchhoff;
      PhiPr -= q1->_backsig;
      
      STensorOperation::zero(PhiPr_Dev);
      STensorOperation::zero(PhiPr_Vol);
      
      static STensor3 devPhipr,devPhi, K_dev_pr, K_dev; // effective dev stress predictor
      STensorOperation::zero(K_dev_pr);
      STensorOperation::zero(K_dev);
      
      double ptildepr,ptilde, ptildepr_Vol, ptilde_Vol;
      STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
      ptildepr/= 3.;
      
      
      
      
      
      ptildepr_Vol = getUpdatedBulkModulus(q1)*A_v_pr*trEepr;
      K_dev_pr = devEepr;
      K_dev_pr *= 2.*getUpdatedShearModulus(q1)*B_d_pr;
      
      
      
      ptilde = ptildepr; // current effective pressure
      ptilde_Vol = ptildepr_Vol; // current effective volumetric pressure  
      devPhi =devPhipr;
      K_dev = K_dev_pr;
    
      double PhiEqpr2 = 1.5*devPhipr.dotprod();
      double PhiEqpr = sqrt(PhiEqpr2);
    
      // plastic poisson ratio
      q1->_nup = (9.-2.*_b)/(18.+2.*_b);
      double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
    
      double Gamma = 0.; //  // flow rule parameter
      double PhiEq = PhiEqpr;
    
     
       // hardening
      this->hardening(q0,q1);
      static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
      this->getYieldCoefficients(q1,a);
    
      double Hb =0.;
      if (q1->_ipKinematic != NULL)
        Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
    
      double Gt= Ge + kk*Hb/2.;
      double Kt = Ke + kk*Hb/3.;
      //a.print("a init");
    
      static STensor3 devN; // dev part of yield normal
      double trN = 0.; // trace part of yield normal
      static STensor3 N; // yield normal
      STensorOperation::zero(devN);
      STensorOperation::zero(N);
    
      double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0));
    
    
      double DfDGamma = 0.;
      double dfdDgamma = 0.;
      double u = 1.;
      double v = 1.;
    
      double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
    
      double dDgammaDGamma = 0.;
      double Dgamma = 0.; // eqplastic strain
      double Dptilde_vol = 0.;
      double dptilde_VoldGamma= 0.;
      double PhiEq_cor = 0.;
      
      
      double A_v=0.;
      double B_d = 0.;
      double dAdtrEe=0.;
      static STensor3 dBddevEe;
      STensorOperation::zero(dBddevEe);
      
      double trEe = 0.;
      double J = 0.;
      trEe = trEepr-2.*_b*Gamma*((ptildepr+Dptilde_vol)/(v));  
      
      static STensor3 DK_dev, devPhipr_Dev_cor, dK_devdGamma;
      STensorOperation::zero(DK_dev);  
      static STensor3 J_dev, devEe;
      STensorOperation::zero(devEe);
      for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              devEe(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
            }
      }
      STensorOperation::zero(J_dev);
      
      
      
    
      STensorOperation::zero(devPhipr_Dev_cor);
      STensorOperation::zero(dK_devdGamma);
      
      
      for (int i=0; i<3; i++){
      	for (int j=0; j<3; j++){
      		devPhipr_Dev_cor(i,j) = devPhipr(i,j)+DK_dev(i,j);
      	}
      }
      PhiEq_cor = sqrt(1.5*devPhipr_Dev_cor.dotprod());
      
      if (q1->dissipationIsBlocked()){
        q1->getRefToDissipationActive() = false;
      }
      else{
        if (f>_tol){
          q1->getRefToDissipationActive() = true;
           // plasticity
          int ite = 0;
          int maxite = 100; // maximal number of iters
    
    
          //Msg::Error("plasticity occurs f = %e",f);
    
          //double f0 = fabs(f);
    
          while (fabs(f) >_tol or ite <1){
            double eta(0.),Deta(0.);
            if (_viscosity != NULL)
              _viscosity->get(q1->_epspbarre,eta,Deta);
            double etaOverDt = eta/this->getTimeStep();
            
            double dAdGamma = 1./(2.*A)*(-(72.*Gt*PhiEq*PhiEq)/(u)-(16.*_b*_b*_b*Kt*ptilde*ptilde/(3*v))+((8.*_b*_b*ptilde)/(3*v))*dptilde_VoldGamma);
            for (int i=0; i<3; i++){
                for (int j=0; j<3; j++){
                   dAdGamma+= 1./(2.*A)*(18./(u*u)*(devPhipr_Dev_cor(i,j)*dK_devdGamma(j,i)));
                }
            }
            dDgammaDGamma = kk*(A+Gamma*dAdGamma);
    
            this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
            dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0); //OK
            if (Gamma>0 and etaOverDt>0)
              dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p);
    
            DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/(v)-((a(1)/v)*dptilde_VoldGamma); 
            for (int i=0; i<3; i++){
                for (int j=0; j<3; j++){
                   DfDGamma+= 3.*_n*a(2)/(2.*u*u)*pow(PhiEq,_n-2.) *(devPhipr_Dev_cor(i,j)*dK_devdGamma(j,i));
                }
            }
            
            if (Gamma>0 and etaOverDt>0)
              DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
    
            double dGamma = -f/DfDGamma;
            if (Gamma + dGamma <=0.){
                Gamma /= 2.;
            }
            else
              Gamma += dGamma;
    
            //Msg::Error("Gamma = %e",Gamma);
    
            u = 1.+6.*Gt*Gamma;
            v = 1.+2.*_b*Kt*Gamma;
                    
            double dtrEedptilde = -2.*_b*Gamma/v;
            J=getUpdatedBulkModulus(q1);
            int ite_vol = 0;
            int maxite_vol = 100; // maximal number of iters
            double dJdptilde_Vol =0.;
            double Np_vol;
      
            while (1){
    
            	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
            	J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
            	dJdptilde_Vol = getUpdatedBulkModulus(q1)*( dAdtrEe*dtrEedptilde*trEe+A_v*dtrEedptilde)-1.;
           
            	ptilde_Vol = ptilde_Vol-J/dJdptilde_Vol; 
           
            	Dptilde_vol = ptilde_Vol-ptildepr_Vol;
    
            	trEe = trEepr-2.*_b*Gamma*(ptildepr+Dptilde_vol)/(v);
            	J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
           
            
            	if ((fabs(J)/getUpdatedBulkModulus(q1) <_tol and ite_vol >2) or ite_vol > maxite_vol)
            	{
                      dptilde_VoldGamma = -(2.*_b*getUpdatedBulkModulus(q1)/(v*v)*(ptildepr+Dptilde_vol)*(dAdtrEe*trEe+A_v))/(1.+2.*_b*Gamma*getUpdatedBulkModulus(q1)/(v)*(dAdtrEe*trEe+A_v));
            	  if (ite_vol > maxite_vol)
            	    Msg::Error("No convergence for phi_p_vol in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, J = %e!!",ite,J);
            	  break;
                    }
            	ite_vol++;
            }
      
     
            //################################### DEVIATORIC TERMS
      
            static STensor43 ddevEedK_dev;
            STensorOperation::zero(ddevEedK_dev);
      
      
     
            double J_dev_tol=getUpdatedShearModulus(q1);
            int ite_dev = 0;
            int maxite_dev = 100; // maximal number of iters
      
            static STensor43 dJ_dev, dJ_dev_inv;
            static STensor3 dTemp, NK_dev;
            STensorOperation::zero(dJ_dev);
            STensorOperation::zero(dJ_dev_inv);
            STensorOperation::zero(dTemp);
            STensorOperation::zero(NK_dev);
              
            while (1){
      
              	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
    
                    for (int i=0; i<3; i++){
            	        for (int j=0; j<3; j++){
            		      J_dev(i,j)=2.*getUpdatedShearModulus(q1)*B_d*devEe(i,j)-K_dev(i,j) ;
            	        }
                    }
    
            	ddevEedK_dev = _I4;
            	ddevEedK_dev *= -3.*Gamma/u;
            	STensorOperation::zero(dTemp);
            	STensorOperation::multSTensor3STensor43(dBddevEe,ddevEedK_dev,dTemp);
            	for (int i=0; i<3; i++){
            		for (int j=0; j<3; j++){
            			for (int k=0; k<3; k++){
            				for (int l=0; l<3; l++){
            					dJ_dev(i,j,k,l) = 2.*getUpdatedShearModulus(q1)*devEe(i,j)*dTemp(k,l)+ 2.*getUpdatedShearModulus(q1)*B_d*ddevEedK_dev(i,j,k,l)-_I4(i,j,k,l);
            				}
            			}
            		}
            	}
    
            	STensorOperation::inverseSTensor43(dJ_dev, dJ_dev_inv);
            	STensorOperation::zero(NK_dev);
            	STensorOperation::multSTensor43STensor3(dJ_dev_inv,J_dev,NK_dev); 
    
            
            	for (int i=0; i<3; i++){
            		for (int j=0; j<3; j++){
            			K_dev(i,j) = K_dev(i,j)-NK_dev(i,j); // tau_c
            			DK_dev(i,j)=K_dev(i,j)-K_dev_pr(i,j); //Delta Tau_c =tau_c - tau_c_pr
            			devEe(i,j) = devEepr(i,j)-3.*Gamma*(devPhipr(i,j)+DK_dev(i,j))/(u);
            			devPhipr_Dev_cor(i,j) = devPhipr(i,j)+DK_dev(i,j); //Phipr+Delta tau_c
             		        J_dev(i,j)=2.*getUpdatedShearModulus(q1)*B_d*devEe(i,j)-K_dev(i,j) ;
            	        }
                    }
            	J_dev_tol=J_dev.norm0();
    
            	if ((fabs(J_dev_tol)/getUpdatedShearModulus(q1)<_tol and ite_dev >2)  or ite_dev > maxite_dev)
            	{
            	        //#############dK_devdGamma
                           double M_d = 0.;
                           static STensor3 T_d, U_d;
     
                           for (int i=0; i<3; i++){
            	            for (int j=0; j<3; j++){
            		        T_d(i,j) = -3./(u*u)*(devPhipr_Dev_cor(i,j));
            		        M_d+=dBddevEe(i,j)*T_d(j,i);
            	            }
                           } 
                           for (int i=0; i<3; i++){
            	          for (int j=0; j<3; j++){
            		     U_d(i,j)=2.*getUpdatedShearModulus(q1)*devEe(i,j)*M_d+2.*getUpdatedShearModulus(q1)*B_d*T_d(i,j);
                              }
                           }
                           STensorOperation::multSTensor43STensor3(dJ_dev_inv,U_d,dK_devdGamma, -1.);
    
                           //#############dK_devdGamma
                           if(ite_dev > maxite_dev)
            		  Msg::Error("No convergence for J_dev in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, J_dev = %e!!",ite_dev,J_dev);
            	       break;
            	   
            	}
            	ite_dev++;
    
            }
            
      
      
            PhiEq_cor = sqrt(1.5*devPhipr_Dev_cor.dotprod());
            PhiEq = PhiEq_cor/u;
            ptilde = (ptildepr+Dptilde_vol)/v;
            A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
            Dgamma = kk*Gamma*A;
    
            //Msg::Error("it = %d, u=%e, v=%e, Dgamma=%e",ite,u,v,Dgamma);
    
            updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
            hardening(q0,q1);
            getYieldCoefficients(q1,a);
    
            f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0)); 
            
            double viscoTerm = etaOverDt*Gamma;
            if (Gamma>0 and etaOverDt>0) f-= pow(viscoTerm,_p);
    
            ite++;
            //if (ite> maxite-5)
             //Msg::Error("it = %d, DfDGamma = %e error = %e dGamma = %e, Gamma = %e",ite,DfDGamma,f,dGamma,Gamma);
             
            //std::cout << "NR F: f = "<<f <<" Gamma = "<<Gamma<<" dfdDgamma = "<<dfdDgamma<< " dDgammaDGamma = "<< dDgammaDGamma << " dptilde_VoldGamma = " << dptilde_VoldGamma<< " Dptilde_vol"<< Dptilde_vol<< " ptilde = "<<ptilde  << " dGamma = "<<dGamma<< " DfDGamma = "<< DfDGamma <<"  Iteration = "<<ite<<" ## Vc = " << getVolumeCorrection()<<" Xi = " << getXiVolumeCorrection()<<" Zeta = " << getZetaVolumeCorrection()<<" Dc = " << getDevCorrection()<<" Theta = " << getThetaDevCorrection()<<" Pi = " << getPiDevCorrection()<<std::endl;
            if (fabs(f) <_tol) break;
    
            if(ite > maxite){
              Msg::Error("No convergence for plastic correction in mlawPowerYieldHyper nonAssociatedFlow Maxwell iter = %d, f = %e!!",ite,f);
              P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
              return;
            }
          };
    
          q1->_DgammaDt = Dgamma/this->getTimeStep();
    
          // update effective stress tensor
          
          devPhi =devPhipr_Dev_cor;
          devPhi *= (1./u);
          devN = devPhi;
          devN *=  3.;
          trN =  2.*_b*ptilde;
          N = devN;
          N(0,0) += trN/3.;
          N(1,1) += trN/3.;
          N(2,2) += trN/3.;
    
          // estimate exp(GammaN)
          static STensor3 expGN;
          static STensor3 GammaN;
          GammaN = N;
          GammaN *= Gamma;
          STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
    
          // update plastic deformation tensor
          STensorOperation::multSTensor3(expGN,Fp0,Fp1);
          // update IP
          updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
    
          // update elastic deformation tensor, corotational stress
          STensorOperation::inverseSTensor3(Fp1,Fpinv);
          STensorOperation::multSTensor3(F,Fpinv,Fe);
          STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
          STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
          Ee *= 0.5;
          // update A, B
          updateViscoElasticFlow(q0,q1,Ke,Ge);
    
          // backstress
          static STensor3 DB; // increment
          DB = (N); // increment
          DB *= (kk*Hb*Gamma);
          q1->_backsig += DB; // update
    
          // corotationaal Kirchhoff stress
          q1->_kirchhoff = devPhi;
          q1->_kirchhoff += q1->_backsig;
    
          q1->_kirchhoff(0,0) += (ptilde);
          q1->_kirchhoff(1,1) += (ptilde);
          q1->_kirchhoff(2,2) += (ptilde);
        }
        else{
          q1->getRefToDissipationActive() = false;
        }
      }
    
    
      const STensor3& KS = q1->_kirchhoff;
      // second Piola Kirchhoff stress
      static STensor3 S;
      STensorOperation::multSTensor3STensor43(KS,DlnDCe,S); 
      for(int i=0; i<3; i++)
        for(int j=0; j<3; j++){
          P(i,j) = 0.;
          for(int k=0; k<3; k++)
            for(int l=0; l<3; l++)
              P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
        }
    
    
      // defo energy
      q1->getRefToElasticEnergy()=deformationEnergy(*q1);
      q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
      q1->getRefToPlasticEnergy() = q0->plasticEnergy();
      if (Gamma > 0){
        double dotKSN = dot(KS,N);
        q1->getRefToPlasticEnergy() += Gamma*dotKSN;
      }
    
      if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
        q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
      }
      else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
               (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
        q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
      }
      else{
        q1->getRefToIrreversibleEnergy() = 0.;
      }
    
    
      if (stiff){
        static STensor3 DpprDCepr;
        static STensor43 DdevKprDCepr;
        
        
        
        STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);   
        DpprDCepr*= (0.5*Ke)+ 0.5*getUpdatedBulkModulus(q1)*(dApr*trEepr+A_v_pr);
        
        
     
        static STensor43 DKpr_dev_DCepr, DKpr_dev_DCepr_Temp, Temp1, Temp2, Temp4;
        STensorOperation::zero(DKpr_dev_DCepr);
        STensorOperation::zero(Temp1);
        STensorOperation::zero(Temp2);
        STensorOperation::zero(Temp4);
    
        STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKprDCepr);
        
        DdevKprDCepr*= Ge;//+getUpdatedShearModulus(q1)*B_d_pr;
        
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                Temp1(i,j,k,l) += getUpdatedShearModulus(q1)*(devEepr(i,j)*dB_devpr(k,l)+B_d_pr*_Idev(i,j,k,l));
              }
            }
          }
        }
        
        STensorOperation::multSTensor43(Temp1,_Idev,Temp2);
        STensorOperation::multSTensor43(Temp2,DlnDCepr,Temp4);
        
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                DdevKprDCepr(i,j,k,l) += Temp4(i,j,k,l);
              }
            }
          }
        }
        
        static STensor3 DpDCepr;
        static STensor43 DdevKDCepr;
        DpDCepr = DpprDCepr;
        DdevKDCepr = DdevKprDCepr;
    
        static STensor43 dFpDCepr;
        static STensor3 DgamaDCepr;
        static STensor3 DtrNDCepr;
        static STensor43 DdevNDCepr;
        static STensor3 DGDCepr;
    
    
        if (Gamma >0){
          // plastic
          static STensor3 dAdCepr, dfDCepr, dptilde_VoldCepr, dgamadCepr, DGDCepr, dptildeOldpr_VoldCepr;
          static STensor43 dK_devdCepr, DdevKOldprDCepr;
          STensorOperation::zero(dAdCepr);
          STensorOperation::zero(dfDCepr);
          STensorOperation::zero(dptilde_VoldCepr);
          STensorOperation::zero(dgamadCepr);
          STensorOperation::zero(DGDCepr);
          STensorOperation::zero(dptildeOldpr_VoldCepr);
          STensorOperation::zero(dK_devdCepr);
          STensorOperation::zero(DdevKOldprDCepr);
                
          for (int i=0; i<3; i++){
          	for (int j=0; j<3; j++){
          		dptilde_VoldCepr(i,j) = ( getUpdatedBulkModulus(q1)*(trEe*dAdtrEe+A_v)*(_I(i,j)-2.*_b*Gamma/v*Ke*_I(i,j)) )/(1.+2.*_b*Gamma*getUpdatedBulkModulus(q1)/v*(trEe*dAdtrEe+A_v));
          	}
          }
          
          
          //########################################dK_devdEepr
          static STensor43 H_d;
          static STensor43 S_dev, S_dev_inv;
          STensorOperation::zero(H_d);
          STensorOperation::zero(S_dev);
          STensorOperation::zero(S_dev_inv);
      
      
          for (int i=0; i<3; i++){
          	for (int j=0; j<3; j++){
          		for (int k=0; k<3; k++){
          			for (int l=0; l<3; l++){
          				for (int m=0; m<3; m++){
          					for (int n=0; n<3; n++){
          						H_d(i,j,k,l) += 2.*getUpdatedShearModulus(q1)*(1.-6.*Ge*Gamma/u)*(devEe(i,j)*dBddevEe(m,n)*_Idev(n,m,k,l));
          						S_dev(i,j,k,l) += 6.*getUpdatedShearModulus(q1)*Gamma/u*(devEe(i,j)*dBddevEe(m,n)*_Idev(n,m,k,l));  				
      					}
      				}
      				H_d(i,j,k,l) += 2.*getUpdatedShearModulus(q1)*(1.-6.*Ge*Gamma/u)*B_d*_Idev(i,j,k,l);
      				S_dev(i,j,k,l) += (1.+6.*getUpdatedShearModulus(q1)*Gamma/u*B_d)*_I4(i,j,k,l);
      			}
      		}
          	}
          }
        
          STensorOperation::inverseSTensor43(S_dev, S_dev_inv);
          STensorOperation::zero(dK_devdCepr);
          STensorOperation::multSTensor43(S_dev_inv,H_d,dK_devdCepr);
      
          //########################################dK_devdEepr
    
          // d/dEpr-> d/dCpr
          STensorOperation::multSTensor43(dK_devdCepr, DlnDCepr, dK_devdCepr);
          dK_devdCepr*=0.5;
          STensorOperation::multSTensor3STensor43(dptilde_VoldCepr,DlnDCepr,dptilde_VoldCepr);
          dptilde_VoldCepr*=0.5;
          
          STensorOperation::multSTensor3STensor43(_I,DlnDCepr,dptildeOldpr_VoldCepr);
          dptildeOldpr_VoldCepr*=0.5*Ke;
          STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKOldprDCepr);
          DdevKOldprDCepr*= Ge;
          
          double fact = 1.5*a(2)*_n*pow(PhiEq,_n-2.)/(u*u); 
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              dAdCepr(i,j) = (4.*_b*_b*ptilde/(A*3.*v))*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldCepr(i,j));
              dfDCepr(i,j) =  -(a(1)/v)*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldCepr(i,j)); 
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  dAdCepr(i,j) += (9./(A*u*u))*devPhipr_Dev_cor(k,l)*(DdevKOldprDCepr(k,l,i,j)+dK_devdCepr(k,l,i,j));
                  dfDCepr(i,j) += fact*devPhipr_Dev_cor(k,l)*(DdevKOldprDCepr(k,l,i,j)+dK_devdCepr(k,l,i,j));  
                }
              }
            }
          }
    
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              DGDCepr(i,j) = (-dfDCepr(i,j)-dfdDgamma*kk*Gamma*dAdCepr(i,j))/DfDGamma;
            }
          }
    
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              DgamaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ dDgammaDGamma*DGDCepr(i,j); // remove the extra kk
            }
          }
    
    
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++)
              DtrNDCepr(i,j) = 2.*_b/v*(dptildeOldpr_VoldCepr(i,j)+dptilde_VoldGamma*DGDCepr(i,j)+dptilde_VoldCepr(i,j)) - 2.*_b*ptilde*(2.*_b*Kt)/(v)*DGDCepr(i,j);
    
    
          for (int i=0; i<3; i++)
              for (int j=0; j<3; j++)
                for (int k=0; k<3; k++)
                  for (int l=0; l<3; l++){
                  	DdevNDCepr(i,j,k,l)  = (3./u)*(DdevKOldprDCepr(i,j,k,l)+dK_devdCepr(i,j,k,l))+(3./u*dK_devdGamma(i,j)-18.*Gt/(u*u)*devPhipr_Dev_cor(i,j))*DGDCepr(k,l);
                  }
    
          static STensor43 temp1;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++)
              for (int k=0; k<3; k++)
                for (int l=0; l<3; l++)
                  temp1(i,j,k,l) = N(i,j)*DGDCepr(k,l)+ Gamma*DdevNDCepr(i,j,k,l)+ Gamma/3.*_I(i,j)*DtrNDCepr(k,l);
    
          static STensor43 EprFp0;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++)
              for (int k=0; k<3; k++)
                for (int l=0; l<3; l++){
                    EprFp0(i,j,k,l) = 0.;
                    for (int s=0; s<3; s++){
                      EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
                    }
                  }
          
          STensorOperation::multSTensor43(EprFp0,temp1,dFpDCepr);
    
          DpDCepr = dptildeOldpr_VoldCepr; 
          // update
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              DpDCepr(i,j) -= Ke*(DGDCepr(i,j)*trN+Gamma*DtrNDCepr(i,j))-dptilde_VoldGamma*DGDCepr(i,j)-dptilde_VoldCepr(i,j); 
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  DdevKDCepr(i,j,k,l) =  DdevKOldprDCepr(i,j,k,l)-2.*Ge*(DGDCepr(k,l)*devN(i,j)+Gamma*DdevNDCepr(i,j,k,l))+dK_devdCepr(i,j,k,l)+DGDCepr(k,l)*dK_devdGamma(i,j);
                }
              }
            }
          }
        }
        else{
          // elastic
          STensorOperation::zero(DgamaDCepr);
          STensorOperation::zero(dFpDCepr);
          STensorOperation::zero(DtrNDCepr);
          STensorOperation::zero(DdevNDCepr);
          STensorOperation::zero(DGDCepr);
        }
    
        static STensor43 dKcorDcepr;
        dKcorDcepr = DdevKDCepr;
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                dKcorDcepr(i,j,k,l) += _I(i,j)*DpDCepr(k,l);
              }
            }
          }
        }
    
        static STensor43 CeprToF;
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                CeprToF(i,j,k,l) = 2.*Fepr(k,i)*invFp0(j,l);
              }
            }
          }
        }
    
        STensor3& DgammaDF = q1->_DgammaDF;
        static STensor43 DKcorDF;
    
        STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
        if (Gamma > 0){
          STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
          STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
        }
        else{
          STensorOperation::zero(DgammaDF);
          STensorOperation::zero(dFpdF);
        }
    
        static STensor43 DinvFpDF; //
        for (int i=0; i<3; i++)
          for (int s=0; s<3; s++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                DinvFpDF(i,s,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int j=0; j<3; j++)
                    DinvFpDF(i,s,k,l) -= Fpinv(i,m)*dFpdF(m,j,k,l)*Fpinv(j,s);
    
              }
    
        for (int m=0; m<3; m++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                dFedF(m,j,k,l) = _I(m,k)*Fpinv(l,j);
                for (int s=0; s<3; s++)
                  dFedF(m,j,k,l) += F(m,s)*DinvFpDF(s,j,k,l);
              }
    
        static STensor63 DlnDF;
        if (_order != 1){
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  for (int p=0; p<3; p++){
                    for (int q=0; q<3; q++){
                      DlnDF(i,j,k,l,p,q) = 0.;
                      for (int r=0; r<3; r++){
                        for (int s=0; s<3; s++){
                          for (int a=0; a<3; a++){
                            DlnDF(i,j,k,l,p,q) += DDlnDDCe(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
    
        }
        else{
          STensorOperation::zero(DlnDF);
        }
    
    
        static STensor43 dSdF;
        STensorOperation::zero(dSdF);
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++)
                for (int m=0; m<3; m++)
                  for (int n=0; n<3; n++){
                    dSdF(i,j,k,l) += DKcorDF(m,n,k,l)*DlnDCe(m,n,i,j);
                    dSdF(i,j,k,l) += KS(m,n)*DlnDF(m,n,i,j,k,l);
                  }
    
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                Tangent(i,j,k,l) = 0.;
                for (int m=0; m<3; m++){
                  for (int n=0; n<3; n++){
                    Tangent(i,j,k,l) +=dFedF(i,m,k,l)*S(m,n)*Fpinv(j,n);
                    Tangent(i,j,k,l) += Fe(i,m)*dSdF(m,n,k,l)*Fpinv(j,n);
                    Tangent(i,j,k,l) += Fe(i,m)*S(m,n)*DinvFpDF(j,n,k,l);
                  }
                }
              }
        dDeformationEnergydF(*q1,  P, q1->getRefToDElasticEnergyPartdF(), Fp1, dFedF);
        dViscousEnergydF(*q0,  *q1,  DlnDCe, Fe, q1->getRefToDViscousEnergyPartdF(), dFedF);
        dPlasticEnergydF(N, Gamma, Dgamma, dKcorDcepr, DGDCepr, DdevNDCepr, DtrNDCepr, KS, CeprToF,q1->getRefToDPlasticEnergyPartdF());
    
        STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
    
        if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
          DirrEnergDF=q1->getConstRefToDElasticEnergyPartdF();
        }
        else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
                 (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
          DirrEnergDF=q1->getConstRefToDPlasticEnergyPartdF();
        }
        else{
          STensorOperation::zero(DirrEnergDF);
        }
    
      };
    };
    
    void mlawPowerYieldHyper::dDeformationEnergydF(const IPHyperViscoElastoPlastic& q,  const STensor3 &P, STensor3 &dElasticEnergydF, const STensor3 &Fp, const STensor43 &dFedF) const
    {
       STensorOperation::zero(dElasticEnergydF);
       for(int i=0;i<3;i++){
         for(int j=0;j<3;j++){
           for(int k=0;k<3;k++){
             for(int l=0;l<3;l++){
               for (int m=0; m<3; m++){
                 dElasticEnergydF(i,j) += P(k,m)*Fp(l,m)*dFedF(k,l,i,j);
               }
             }
          }
        }
      }
    }
    
    void mlawPowerYieldHyper::dViscousEnergydF(const IPHyperViscoElastoPlastic& q0,  const IPHyperViscoElastoPlastic& q,  const STensor43 &dlnCdC, const STensor3 &Fe, STensor3 &dViscousEnergydF, const STensor43 &dFedF) const
    {
      STensorOperation::zero(dViscousEnergydF);
      static STensor3 dViscousEnergyPartdFe;
      STensorOperation::zero(dViscousEnergyPartdFe);
      mlawHyperViscoElastic::dViscousEnergydF(dynamic_cast<const IPHyperViscoElastic &>(q0),dynamic_cast<const IPHyperViscoElastic &>(q),dlnCdC, Fe, dViscousEnergyPartdFe);
      for(int i=0; i<3; i++)
      {
        for(int J=0; J<3; J++)
        {
          for(int k=0; k<3; k++)
          {
            for(int L=0; L<3; L++)
            {
              dViscousEnergydF(i,J)+=dViscousEnergyPartdFe(k,L)*dFedF(k,L,i,J);
            }
          }
        }
      }
    }
    void mlawPowerYieldHyper::dPlasticEnergydF(const STensor3 &N, double Gamma, double Dgamma, const STensor43 &dKcorDcepr, const STensor3 &DGDCepr, const STensor43 &DdevNDCepr, const STensor3 &DtrNDCepr, const STensor3 &KS, const STensor43 &CeprToF, STensor3 &dPlasticEnergydF) const
    {
      STensorOperation::zero(dPlasticEnergydF);
      if (Dgamma > 0){
        static STensor3 DirrEnergDCepr;
        double dotKSN = dot(KS,N);
        DirrEnergDCepr = DGDCepr;
        DirrEnergDCepr *= dotKSN;
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                DirrEnergDCepr(i,j) += (Gamma*dKcorDcepr(k,l,i,j)*N(k,l)+ Gamma*KS(k,l)*(DdevNDCepr(k,l,i,j)+_I(k,l)*DtrNDCepr(i,j)/3.));
              }
            }
          }
        }
        STensorOperation::multSTensor3STensor43(DirrEnergDCepr,CeprToF,dPlasticEnergydF);
      }
    }
    
    void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF) const{
      /* compute elastic predictor */
      STensor3& Fp1 = q1->_Fp;
      const STensor3& Fp0 = q0->_Fp;
    
      Fp1 = Fp0; // plastic deformation tensor
      q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
      q1->_epspCompression = q0->_epspCompression;
      q1->_epspTraction = q0->_epspTraction;
      q1->_epspShear = q0->_epspShear;
      q1->_backsig = q0->_backsig; // backstress tensor
      q1->_DgammaDt = 0.; // plastic rate --> failure
    
      static STensor3 Fpinv, Ce, Fepr;
      STensorOperation::inverseSTensor3(Fp1,Fpinv);
      STensorOperation::multSTensor3(F,Fpinv,Fepr);
      STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
    
      static STensor3 invFp0; // plastic predictor
      invFp0= Fpinv;
      STensor3& Fe = q1->_Fe;
      Fe = Fepr;
    
      static STensor43 DlnDCepr, DlnDCe;
      static STensor63 DDlnDDCe;
      static STensor43 dexpAdA; // estimation of dexpA/dA
    
      STensor3& Ee = q1->_Ee;
      STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCepr,&DDlnDDCe);
      Ee *= 0.5;
      DlnDCe = DlnDCepr;
    
      // update A, B
      double Ge, Ke;
      
      
      
      double A_v_pr=0., B_d_pr = 0., dApr=0.;
      static STensor3 dB_devpr;
      STensorOperation::zero(dB_devpr);
      viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
    
      static STensor3 devKpr; // dev corotational kirchoff stress predictor
      static double ppr; // pressure predictor
    
      STensorOperation::decomposeDevTr(q1->_kirchhoff,devKpr,ppr);
      ppr /= 3.;
      double keqpr = sqrt(1.5*devKpr.dotprod());
    
      static STensor3 devK;
      devK= devKpr;  // dev corotational kirchoff stress
      double p = ppr; // pressure
    
       // hardening
      this->hardening(q0,q1);
      fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
      this->getYieldCoefficients(q1,a);
    
      static STensor3 devN; // dev part of yield normal
      double trN = 0.; // trace part of yield normal
      static STensor3 N; // yield normal
    
      STensorOperation::zero(devN);
      STensorOperation::zero(N);
    
      double sigVM = keqpr;
      double Dgamma = 0.;
      double Gamma = 0.;
    
      double g0 = 0.;
      double z = a(1)/a(2)*pow(sigVM,1.-_n);
      double A = sqrt(1.5+1.*z*z/(3.*_n*_n));
    
      static fullMatrix<double> invJ(2,2); // inversed jacobian
      invJ.setAll(0.);
    
      double dAdz =0.;
      double dzdDgamma = 0.;
      double dzdsigVM = 0.;
      double dg0dsigVM = -1./(3.*Ge);
    
      if (q1->dissipationIsBlocked()){
        q1->getRefToDissipationActive() = false;
      }
      else{
        double f = a(2)*pow(sigVM,_n) - (a(1)*ppr+a(0));
        if (f>_tol){
          q1->getRefToDissipationActive() = true;
          // plasticity
          int ite = 0;
          int maxite = 1000; // maximal number of iters
          double val  = 1.5*a(2)*_n*pow(3.*fabs(p),_n-2.);
          q1->_nup = (val*p+ a(1)/3.)/(val*2.*p - a(1)/3.);
          double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
          A*=kk;
    
          fullVector<double> R(2); // residual of eqs to estimate Dgamma, Dpression and Gamma
    
          R(0) = Dgamma - g0*A;
          R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
    
          double res = R.norm();
    
          while(res>_tol or ite<1){
    
            this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
    
            static fullMatrix<double> J(2,2);
            J.setAll(0.);
    
            dAdz = kk*z/(3.*_n*_n*A);
            dzdDgamma = (Da(1)*a(2) - a(1)*Da(2))/(a(2)*a(2))*pow(sigVM,1.-_n);
            dzdsigVM = a(1)/a(2)*(1.-_n)*pow(sigVM,-_n);
    
    
    
            J(0,0) = 1. - g0*dAdz*dzdDgamma;
            J(0,1) = -dg0dsigVM*A-g0*dAdz*dzdsigVM;
    
            J(1,0) =  Da(2)*pow(sigVM,_n) - Da(1)*(ppr+Ke*z*g0/_n) - Da(0) - a(1)*Ke*dzdDgamma*g0/_n;
            J(1,1) = a(2)*_n*pow(sigVM,_n-1) - a(1)*Ke/_n*(dzdsigVM*g0+z*dg0dsigVM);
    
    
    
            double detJ = J(0,0)*J(1,1) - J(1,0)*J(0,1);
            if (detJ == 0.) Msg::Error("the corrected system can not be solved: mlawPowerYieldHyper::predictorCorrector_associatedFlow");
    
            invJ(0,0) = J(1,1)/detJ;
            invJ(0,1) = -J(0,1)/detJ;
            invJ(1,0) = -J(1,0)/detJ;
            invJ(1,1) = J(0,0)/detJ;
    
            double solDgamma = -invJ(0,0)*R(0)-invJ(0,1)*R(1);
            double solsigVM = -invJ(1,0)*R(0)-invJ(1,1)*R(1);
    
           // bool isSolved = J.luSolve(residual,sol);
           // if (!isSolved) Msg::Error("the corrected system can not be solved: mlawPowerYieldHyper::predictorCorrector");
            //sol.print("sol");
    
            if (solDgamma >0.1)
              Dgamma+= 0.1;
            else if (Dgamma+ solDgamma< 0.)
              Dgamma /= 2.;
            else
              Dgamma += solDgamma;
    
            if (sigVM + solsigVM< 0.)
              sigVM /= 2.;
            else
              sigVM += solsigVM;
    
    
            updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
            this->hardening(q0,q1);
            this->getYieldCoefficients(q1,a);
    
            g0 = (keqpr- sigVM)/(3.*Ge);
            z = a(1)/a(2)*pow(sigVM,1.-_n);
            A = kk*sqrt(1.5+1.*z*z/(3.*_n*_n));
    
            R(0) = Dgamma - g0*A;
            R(1) = a(2)*pow(sigVM,_n) - a(1)*(ppr+Ke*z*g0/_n) - a(0);
    
            res = R.norm();
    
            ite++;
    
            //Msg::Info("iter = %d, res = %e,  Dgamma = %e, sigVM/keqpr = %e",ite,res, Dgamma, sigVM/keqpr);
    
            if(ite > maxite){
              Msg::Error("No convergence for plastic correction in mlawPowerYieldHyper iter = %d, res = %e!!",ite,res);
              P(0,0) = P(1,1) = P(2,2) = sqrt(-1.);
              break;
            }
          }
    
          Gamma = z*g0/(_n*a(1));
    
          // update
          p += Ke*Gamma*a(1);
    
          double ff =a(2)*_n*pow(sigVM,_n-2.);
          devK*= 1./(1.+3.*Ge*Gamma*ff);
    
          // estimate yield normal
          devN = devK;
          devN *=  1.5*ff;
          trN = -a(1);
    
          N = devN;
          N(0,0) += trN/3.;
          N(1,1) += trN/3.;
          N(2,2) += trN/3.;
    
          // estimate exp(GammaN)
          static STensor3 expGN;
          static STensor3 GammaN;
          GammaN = N;
          GammaN *= Gamma;
          STensorOperation::expSTensor3(GammaN,_order,expGN,&dexpAdA);
    
          // update plastic deformation tensor
          STensorOperation::multSTensor3(expGN,Fp0,Fp1);
          // update IP
          updateEqPlasticDeformation(q1,q0,q1->_nup,Dgamma);
    
          // update elastic deformation tensor, corotational stress
          STensorOperation::inverseSTensor3(Fp1,Fpinv);
          STensorOperation::multSTensor3(F,Fpinv,Fe);
          STensorOperation::multSTensor3FirstTranspose(Fe,Fe,Ce);
          STensorOperation::logSTensor3(Ce,_order,Ee,&DlnDCe,&DDlnDDCe);
          Ee *= 0.5;
    
          updateViscoElasticFlow(q0,q1,Ke,Ge);
        }
        else{
          q1->getRefToDissipationActive() = false;
        }
      }
    
      // corotational Kirchhoff stress tenor
    
      STensor3& KS = q1->_kirchhoff;
      KS = devK;
      KS(0,0) += p;
      KS(1,1) += p;
      KS(2,2) += p;
    
      // first Piola Kirchhoff stress
      static STensor3 S;
      S*=(0.);
      for(int i=0; i<3; i++)
        for(int j=0; j<3; j++)
          for(int k=0; k<3; k++)
            for(int l=0; l<3; l++)
              S(i,j)+= KS(k,l)*DlnDCe(k,l,i,j);
    
      STensorOperation::zero(P);
      for(int i=0; i<3; i++)
        for(int j=0; j<3; j++)
          for(int k=0; k<3; k++)
            for(int l=0; l<3; l++)
              P(i,j) += Fe(i,k)*S(k,l)*Fpinv(j,l);
    
      // defo energy
      q1->_elasticEnergy=deformationEnergy(*q1);
      q1->getRefToPlasticEnergy() = q0->plasticEnergy();
      if (Gamma > 0){
        double dotKSN = dot(KS,N);
        q1->getRefToPlasticEnergy() += Gamma*dotKSN;
      }
    
      if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
        q1->getRefToIrreversibleEnergy() = q1->defoEnergy();
      }
      else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
               (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
        q1->getRefToIrreversibleEnergy() = q1->plasticEnergy();
      }
      else{
        q1->getRefToIrreversibleEnergy() = 0.;
      }
    
      if (stiff){
        static STensor3 DpprDCepr;
        static STensor43 DdevKprDCepr;
        STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpprDCepr);
        DpprDCepr*= (0.5*Ke);
        STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKprDCepr);
        DdevKprDCepr*= Ge;
    
        static STensor3 DpDCepr;
        static STensor43 DdevKDCepr;
        DpDCepr = DpprDCepr;
        DdevKDCepr = DdevKprDCepr;
    
        static STensor43 dFpDCepr;
        static STensor3 DgamaDCepr;
        static STensor3 DGDCepr;
        static STensor3 DtrNDCepr;
        static STensor43 DdevNDCepr;
        if (Dgamma > 0.){
          static STensor3 dg0dCepr;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++){
              dg0dCepr(i,j) = 0.;
              for (int k=0; k<3; k++)
                for (int l=0; l<3; l++)
                  dg0dCepr(i,j) += devKpr(k,l)*DdevKprDCepr(k,l,i,j)/(2.*Ge*keqpr);
    
            }
    
          static STensor3 DPhigammaDCepr;
          DPhigammaDCepr = dg0dCepr;
          DPhigammaDCepr *= -A;
    
    
          static STensor3 DPhiSigVMDCepr;
          DPhiSigVMDCepr = (dg0dCepr);
          DPhiSigVMDCepr *= Ke*z/_n;
          DPhiSigVMDCepr += DpprDCepr;
          DPhiSigVMDCepr *= -a(1);
    
          static STensor3 DsigVMDCepr;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++){
              DgamaDCepr(i,j) = -invJ(0,0)*DPhigammaDCepr(i,j) - invJ(0,1)*DPhiSigVMDCepr(i,j);
              DsigVMDCepr(i,j) = -invJ(1,0)*DPhigammaDCepr(i,j) - invJ(1,1)*DPhiSigVMDCepr(i,j);
            }
    
    
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++){
              DGDCepr(i,j) = (g0/(_n*a(1)))*(dzdDgamma*DgamaDCepr(i,j)+dzdsigVM*DsigVMDCepr(i,j)) +
                        z/(_n*a(1))*(dg0dsigVM*DsigVMDCepr(i,j)+ dg0dCepr(i,j)) - z*g0/(_n*a(1)*a(1))*Da(1)*DgamaDCepr(i,j);
            }
    
    
    
          DtrNDCepr = (DgamaDCepr);
          DtrNDCepr  *= -Da(1);
    
    
          DdevNDCepr = (DdevKprDCepr);
          double B = a(2)*pow(sigVM,_n-2.);
          double u = 1.+3.*Ge*Gamma*_n*B;
          double u2 = u*u;
          static STensor3 DdevNDB;
          DdevNDB = (devKpr);
          DdevNDB *= (1.5*_n/u2);
          static STensor3 DdevNDGamma;
          DdevNDGamma = (devKpr);
          DdevNDGamma *= (-1.5*_n*B*3.*Ge*_n*B/u2);
    
          DdevNDCepr *= (1.5*_n*B/u);
          for (int i=0; i<3; i++)
              for (int j=0; j<3; j++)
                for (int k=0; k<3; k++)
                  for (int l=0; l<3; l++){
                    DdevNDCepr(i,j,k,l) +=DdevNDB(i,j)*(Da(2)*pow(sigVM,_n-2.)*DgamaDCepr(k,l)+
                                                        a(2)*(_n-2.)*pow(sigVM,_n-3.)*DsigVMDCepr(k,l));
                    DdevNDCepr(i,j,k,l) += DdevNDGamma(i,j)*DGDCepr(k,l);
                  }
    
    
    
          // compute dFpdF
          static STensor43 temp1;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++)
              for (int k=0; k<3; k++)
                for (int l=0; l<3; l++)
                  temp1(i,j,k,l) = N(i,j)*DGDCepr(k,l)+ Gamma*DdevNDCepr(i,j,k,l)+ Gamma/3.*_I(i,j)*DtrNDCepr(k,l);
    
          static STensor43 EprFp0;
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++)
              for (int k=0; k<3; k++)
                for (int l=0; l<3; l++){
                  EprFp0(i,j,k,l) = 0.;
                  for (int s=0; s<3; s++){
                    EprFp0(i,j,k,l) += dexpAdA(i,s,k,l)*Fp0(s,j);
                  }
                }
    
          STensorOperation::multSTensor43(EprFp0,temp1,dFpDCepr);
    
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              DpDCepr(i,j) -= Ke*(DGDCepr(i,j)*trN+Gamma*DtrNDCepr(i,j));
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  DdevKDCepr(i,j,k,l) -=  2.*Ge*(DGDCepr(k,l)*devN(i,j)+Gamma*DdevNDCepr(i,j,k,l));
                }
              }
            }
          }
        }
        else{
          STensorOperation::zero(DgamaDCepr);
          STensorOperation::zero(dFpDCepr);
          STensorOperation::zero(DGDCepr);
          STensorOperation::zero(DdevNDCepr);
          STensorOperation::zero(DtrNDCepr);
        }
    
    
        static STensor43 dKcorDcepr;
        dKcorDcepr = DdevKDCepr;
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                dKcorDcepr(i,j,k,l) += _I(i,j)*DpDCepr(k,l);
              }
            }
          }
        }
    
        static STensor43 CeprToF;
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            for (int k=0; k<3; k++){
              for (int l=0; l<3; l++){
                CeprToF(i,j,k,l) = 2.*Fepr(k,i)*invFp0(j,l);
              }
            }
          }
        }
    
        STensor3& DgammaDF = q1->_DgammaDF;
        static STensor43 DKcorDF;
    
        STensorOperation::multSTensor3STensor43(DgamaDCepr,CeprToF,DgammaDF);
        STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
        STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
    
        static STensor43 DinvFpDF; //
        for (int i=0; i<3; i++)
          for (int s=0; s<3; s++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                DinvFpDF(i,s,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int j=0; j<3; j++)
                    DinvFpDF(i,s,k,l) -= Fpinv(i,m)*dFpdF(m,j,k,l)*Fpinv(j,s);
    
              }
    
        for (int m=0; m<3; m++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                dFedF(m,j,k,l) = _I(m,k)*Fpinv(l,j);
                for (int s=0; s<3; s++)
                  dFedF(m,j,k,l) += F(m,s)*DinvFpDF(s,j,k,l);
              }
    
        static STensor63 DlnDF;
        if (_order != 1){
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  for (int p=0; p<3; p++){
                    for (int q=0; q<3; q++){
                      DlnDF(i,j,k,l,p,q) = 0.;
                      for (int r=0; r<3; r++){
                        for (int s=0; s<3; s++){
                          for (int a=0; a<3; a++){
                            DlnDF(i,j,k,l,p,q) += DDlnDDCe(i,j,k,l,r,s)*2.*Fe(a,r)*dFedF(a,s,p,q);
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
    
        }
        else{
          STensorOperation::zero(DlnDF);
        }
    
    
        static STensor43 dSdF;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                dSdF(i,j,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int n=0; n<3; n++){
                    dSdF(i,j,k,l) += DKcorDF(m,n,k,l)*DlnDCe(m,n,i,j);
                    if (_order != 1){
                      dSdF(i,j,k,l) += KS(m,n)*DlnDF(m,n,i,j,k,l);
                    }
                  }
    
              }
    
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                Tangent(i,j,k,l) = 0.;
                for (int m=0; m<3; m++){
                  for (int n=0; n<3; n++){
                    Tangent(i,j,k,l) += dFedF(i,m,k,l)*S(m,n)*Fpinv(j,n);
                    Tangent(i,j,k,l) += Fe(i,m)*dSdF(m,n,k,l)*Fpinv(j,n);
                    Tangent(i,j,k,l) += Fe(i,m)*S(m,n)*DinvFpDF(j,n,k,l);
                  }
                }
              }
    
    
        STensor3& DirrEnergDF = q1->_DirreversibleEnergyDF;
        if (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DEFO_ENERGY){
          for(int i=0;i<3;i++){
            for(int j=0;j<3;j++){
              DirrEnergDF(i,j) = 0.;
              for(int k=0;k<3;k++){
                for(int l=0;l<3;l++){
                  for (int m=0; m<3; m++){
                    DirrEnergDF(i,j) += P(k,m)*Fp1(l,m)*dFedF(k,l,i,j);
                  }
                }
              }
            }
          }
        }
        else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
                 (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)){
          if (Dgamma > 0){
            static STensor3 DirrEnergDCepr;
            double dotKSN = dot(KS,N);
            DirrEnergDCepr = DGDCepr;
            DirrEnergDCepr *= dotKSN;
            for (int i=0; i<3; i++){
              for (int j=0; j<3; j++){
                for (int k=0; k<3; k++){
                  for (int l=0; l<3; l++){
                    DirrEnergDCepr(i,j) += (Gamma*dKcorDcepr(k,l,i,j)*N(k,l)+ Gamma*KS(k,l)*(DdevNDCepr(k,l,i,j)+_I(k,l)*DtrNDCepr(i,j)/3.));
                  }
                }
              }
            }
            STensorOperation::multSTensor3STensor43(DirrEnergDCepr,CeprToF,DirrEnergDF);
          }
          else{
            STensorOperation::zero(DirrEnergDF);
          }
        }
        else{
          STensorOperation::zero(DirrEnergDF);
        }
      }
    };
    
    void mlawPowerYieldHyper::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent,STensor43& dFedF, STensor43& dFpdF,
                                STensor43* elasticTangent) const{
      if (_tangentByPerturbation){
        if (_nonAssociatedFlow){
          this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF);
        }
        else{
          this->predictorCorrector_associatedFlow(F,q0,q1,P,false,Tangent,dFedF,dFpdF);
        }
    
        if (stiff){
          tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1);
        }
    
      }
      else{
        if (_nonAssociatedFlow){
          this->predictorCorrector_nonAssociatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF);
        }
        else{
          this->predictorCorrector_associatedFlow(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF);
        }
      }
      // compute mechanical tengent
      if (elasticTangent!= NULL)
      {
        STensor3 Ptmp(0.);
        static IPHyperViscoElastoPlastic q1tmp(*q1);
        q1tmp.operator=(*(dynamic_cast<const IPVariable*> (q1)));
        mlawHyperViscoElastic::predictorCorrector_ViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), Ptmp, q0, &q1tmp,
                                      true, *elasticTangent);
      }
    
    };
    
    
    
    mlawPowerYieldHyper::mlawPowerYieldHyper(const int num,const double E,const double nu, const double rho,
                            const double tol,
                            const bool matrixbyPerturbation, const double pert):
                            mlawHyperViscoElastic(num,E,nu,rho,matrixbyPerturbation,pert),_tol(tol),
                            _n(1.5),_nonAssociatedFlow(true),_b(0.3),
                            _viscosity(NULL),_p(1.),_compression(NULL),_traction(NULL),_kinematic(NULL){
    };
    
    
    mlawPowerYieldHyper::mlawPowerYieldHyper(const mlawPowerYieldHyper& src):mlawHyperViscoElastic(src),
                            _tol(src._tol),_n(src._n),_nonAssociatedFlow(src._nonAssociatedFlow),_b(src._b), _p(src._p){
      _viscosity = NULL;
      if (src._viscosity) _viscosity = src._viscosity->clone();
    
      _compression = NULL;
      if (src._compression) _compression = src._compression->clone();
    
      _traction = NULL;
      if (src._traction) _traction = src._traction->clone();
    
      _kinematic= NULL;
      if (src._kinematic) _kinematic = src._kinematic->clone();
    };
    
    mlawPowerYieldHyper& mlawPowerYieldHyper::operator=(const materialLaw& source){
    
      mlawHyperViscoElastic::operator=(source);
      const mlawPowerYieldHyper* src =dynamic_cast<const mlawPowerYieldHyper*>(&source);
      if(src != NULL){
    
        _tol = src->_tol;
        _n = src->_n;
        _nonAssociatedFlow = src->_nonAssociatedFlow;
        _b = src->_b;
        _p = src->_p;
    
        if(_viscosity != NULL) delete _viscosity;
        if (src->_viscosity != NULL){ _viscosity = src->_viscosity->clone();}
    
        if(_compression != NULL) delete _compression;
        if (src->_compression!=NULL){ _compression = src->_compression->clone();}
    
        if(_traction != NULL) delete _traction;
        if (src->_traction!=NULL){ _traction = src->_traction->clone();}
    
        if(_kinematic != NULL) delete _kinematic;
        if (src->_kinematic!=NULL){ _kinematic = src->_kinematic->clone();}
    
      }
      return *this;
    }
    
    mlawPowerYieldHyper::~mlawPowerYieldHyper(){
      if (_compression) delete _compression; _compression = NULL;
      if (_traction) delete _traction; _traction = NULL;
      if (_kinematic) delete _kinematic; _kinematic = NULL;
      if (_viscosity) delete _viscosity; _viscosity = NULL;
    };
    
    
    void mlawPowerYieldHyper::setPowerFactor(const double n) {
      _n = n;
    };
    void mlawPowerYieldHyper::nonAssociatedFlowRuleFactor(const double b){
      _b = b;
    };
    void mlawPowerYieldHyper::setPlasticPoissonRatio(const double nup){
      _b = 4.5*(1.-2.*nup)/(1.+nup);
    };
    void mlawPowerYieldHyper::setNonAssociatedFlow(const bool flag) {
      _nonAssociatedFlow = flag;
      if (_nonAssociatedFlow){
        Msg::Info("non associated flow is used");
      }
      else{
        Msg::Info("associated flow is used, kinematic hardening is not considered");
      }
    };
    
    void mlawPowerYieldHyper::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
    {
      IPVariable* ipvi = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
      IPVariable* ipv1 = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
      IPVariable* ipv2 = new IPHyperViscoElastoPlastic(_compression,_traction,_kinematic,_N);
      if(ips != NULL) delete ips;
      ips = new IP3State(state_,ipvi,ipv1,ipv2);
    };
    
    
    mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const int num,const double E,const double nu, const double rho,
                        const double tol, const bool matrixbyPerturbation, const double pert):
                        mlawPowerYieldHyper(num,E,nu,rho,tol,matrixbyPerturbation,pert){
      _failureCriterion = new mlawHyperelasticFailureTrivialCriterion(num);
    };
    
    
    
    mlawPowerYieldHyperWithFailure::mlawPowerYieldHyperWithFailure(const mlawPowerYieldHyperWithFailure& src):
        mlawPowerYieldHyper(src){
      _failureCriterion = NULL;
      if (src._failureCriterion){
        _failureCriterion = src._failureCriterion->clone();
      }
    };
    mlawPowerYieldHyperWithFailure::~mlawPowerYieldHyperWithFailure(){
      if (_failureCriterion != NULL) {
        delete _failureCriterion;
        _failureCriterion = NULL;
      }
    };
    
    void mlawPowerYieldHyperWithFailure::checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{
      IPHyperViscoElastoPlastic *q1 = dynamic_cast<IPHyperViscoElastoPlastic*>(ipv);
      const IPHyperViscoElastoPlastic *q0 = dynamic_cast<const IPHyperViscoElastoPlastic*>(ipvprev);
      if (q1 != NULL){
        double& r = q1->getRefToFailureOnset();
        r = q0->getFailureOnset();
        double failCr = _failureCriterion->getFailureCriterion(q0,q1);
        if (failCr > r){
          r = failCr;
        }
        // check failure
        if (q0->inPostFailureStage()){
          q1->getRefToInPostFailureStage() = true;
        }
        else{
          if (r > 0.){
            q1->getRefToInPostFailureStage() = true;
          }
          else{
            q1->getRefToInPostFailureStage() = false;
          }
        }
      }
    };
    
    void mlawPowerYieldHyperWithFailure::predictorCorrector(const STensor3& F, const IPHyperViscoElastoPlastic *q0, IPHyperViscoElastoPlastic *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF,
                                STensor43* elasticTangent) const{
      mlawPowerYieldHyper::predictorCorrector(F,q0,q1,P,stiff,Tangent,dFedF,dFpdF, elasticTangent);
    
      // compute failure r and DrDF
      double& gF = q1->getRefToFailurePlasticity();
    	gF = q0->getFailurePlasticity();
    
      if (q1->dissipationIsBlocked()){
        if (stiff){
          STensorOperation::zero(q1->getRefToDFailurePlasticityDF());
        }
      }
      else{
        if (q1->inPostFailureStage()){
          gF += q1->_epspbarre - q0->_epspbarre;
          if (stiff){
            q1->getRefToDFailurePlasticityDF() = q1->_DgammaDF;
          }
        }
        else{
          if (stiff){
            STensorOperation::zero(q1->getRefToDFailurePlasticityDF());
          }
        }
    
      }
    };
    
    
    void mlawPowerYieldHyperWithFailure::setFailureCriterion(const mlawHyperelasticFailureCriterionBase& fCr){
      if (_failureCriterion) delete _failureCriterion;
      _failureCriterion = fCr.clone();
    };