Skip to content
Snippets Groups Projects
Select Git revision
  • 298d74b3c28a96c6c074008e97a5db10fffd7604
  • master default protected
  • 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
  • convert_fdivs
  • tmp_jcjc24
  • 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

GModelIO_GEOM.cpp

Blame
  • mlawNonLinearTVP.cpp 96.91 KiB
    //
    // C++ Interface: Material Law
    //
    // Description: Non-Linear Thermo-Visco-Mechanics (Thermo-ViscoElasto-ViscoPlasto Law) with Non-Local Damage Interface (soon.....)
    //
    // Author: <Ujwal Kishore J - FLE_Knight>, (C) 2022
    //
    // Copyright: See COPYING file that comes with this distribution
    //
    //
    
    #include "STensorOperations.h"
    #include "nonLinearMechSolver.h"
    #include "mlawNonLinearTVP.h"
    
    
    mlawNonLinearTVP::mlawNonLinearTVP(const int num,const double E,const double nu, const double rho, const double tol,
                                       const double Tinitial, const double Alpha, const double KThCon, const double Cp,
                                       const bool matrixbyPerturbation, const double pert, const bool thermalEstimationPreviousConfig):
                                       mlawNonLinearTVE(num, E, nu, rho, tol, Tinitial, Alpha, KThCon, Cp, matrixbyPerturbation, pert, thermalEstimationPreviousConfig),
                                       _TaylorQuineyFactor(0.9), _HR(1.){
    
          // by default, no temperature dependence
          _temFunc_Sy0 = new constantScalarFunction(1.);
          _temFunc_H = new constantScalarFunction(1.);
          _temFunc_Hb = new constantScalarFunction(1.);
          _temFunc_visc = new constantScalarFunction(1.);
             
    };
    
    // , mlawPowerYieldHyper(src)
    mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTVE(src) {
        
        _TaylorQuineyFactor = src._TaylorQuineyFactor;
        _HR = src._HR;
    
      _temFunc_Sy0 = NULL; // temperature dependence of initial yield stress
      if (src._temFunc_Sy0!= NULL)
        _temFunc_Sy0 = src._temFunc_Sy0->clone();
    
      _temFunc_H = NULL; // temperature dependence of hardening stress
      if (src._temFunc_H != NULL)
        _temFunc_H = src._temFunc_H->clone();
      
        _temFunc_Hb = NULL; // temperature dependence of kinematic hardening
      if (src._temFunc_Hb != NULL)
        _temFunc_Hb = src._temFunc_Hb->clone();
        
        _temFunc_visc = NULL; // temperature dependence of viscosity
      if (src._temFunc_visc != NULL)
        _temFunc_visc = src._temFunc_visc->clone();
      
    };
    
    mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
        
        mlawNonLinearTVE::operator=(source);
        // mlawPowerYieldHyper::operator=(source);
        const mlawNonLinearTVP* src =dynamic_cast<const mlawNonLinearTVP*>(&source);
        if(src != NULL){
            _TaylorQuineyFactor = src->_TaylorQuineyFactor;
            _HR = src->_HR;
            if(_temFunc_Sy0 != NULL) delete _temFunc_Sy0; // temperature dependence of initial yield stress
            if (src->_temFunc_Sy0!= NULL)
                _temFunc_Sy0 = src->_temFunc_Sy0->clone();
            if(_temFunc_H != NULL) delete _temFunc_H; // temperature dependence of hardening stress
            if (src->_temFunc_H != NULL)
                _temFunc_H = src->_temFunc_H->clone();
            if(_temFunc_Hb != NULL) delete _temFunc_Hb; // temperature dependence of kinematic hardening
            if (src->_temFunc_Hb != NULL)
                _temFunc_Hb = src->_temFunc_Hb->clone();
            if(_temFunc_visc != NULL) delete _temFunc_visc; // temperature dependence of viscosity
            if (src->_temFunc_visc != NULL)
                _temFunc_visc = src->_temFunc_visc->clone();
             
        }
        return *this;
    };
    
    mlawNonLinearTVP::~mlawNonLinearTVP(){
        if (_temFunc_Sy0 != NULL) delete _temFunc_Sy0;  _temFunc_Sy0 = NULL;
        if (_temFunc_H != NULL) delete _temFunc_H; _temFunc_H= NULL;  
        if (_temFunc_Hb != NULL) delete _temFunc_Hb; _temFunc_Hb= NULL;  
        if (_temFunc_visc != NULL) delete _temFunc_visc; _temFunc_visc= NULL;  
    };
    
    // CHECK IF IPSTATE NEEDS the same definition as in mlawPowerYieldHyper --- WHAT??
    
    // NEW
    
    void mlawNonLinearTVP::setIsotropicHardeningCoefficients(const double HR1, const double HR2, const double HR3){
    
        _HR.clear();
        _HR.resize(3,0.);
        _HR[0] = HR1;
        _HR[1] = HR2;
        _HR[2] = HR3;
    };
    
    void mlawNonLinearTVP::hardening(const IPNonLinearTVP* q0, IPNonLinearTVP* q1, const double& T) const{
      //Msg::Error("epspCompression = %e, epspTRaction = %e, epspShear = %e",q->_epspCompression,q->_epspTraction,q->_epspShear);
      if (_compression != NULL && q1->_ipCompression != NULL){
        _compression->hardening(q0->_epspCompression, *q0->_ipCompression, q1->_epspCompression,T,*q1->_ipCompression);
      }
    
      if (_traction!= NULL && q1->_ipTraction != NULL){
        _traction->hardening(q0->_epspTraction,*q0->_ipTraction, q1->_epspTraction,T,*q1->_ipTraction);
      }
    
      if (_kinematic!= NULL && q1->_ipKinematic != NULL)
        _kinematic->hardening(q0->_epspbarre,*q0->_ipKinematic,q1->_epspbarre,T,*q1->_ipKinematic);
    };
    
    void mlawNonLinearTVP::getModifiedMandel(const STensor3& C, const IPNonLinearTVP *q0, IPNonLinearTVP *q1) const{
        
        STensor3& M = q1->getRefToModifiedMandel();
        const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
        
        STensor3 Cinv;
        STensorOperation::inverseSTensor3(C, Cinv);
        
        STensor3 temp1, temp2;
        STensorOperation::multSTensor3(C, corKir, temp1);
        STensorOperation::multSTensor3(temp1, Cinv, temp2);
        
        M = 0.5*(corKir + temp2);
    };
    
    void mlawNonLinearTVP::checkCommutavity(STensor3& commuteCheck,const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const{
        
        const STensor3& M = q1->getConstRefToModifiedMandel();
        const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
        
        STensor3 temp1, temp2, temp3, temp4, temp6;
        STensorOperation::multSTensor3(Ce, S, temp1);
        STensorOperation::multSTensor3(S, Ce, temp2);
        STensorOperation::multSTensor3(corKir, Ce, temp3);
        STensorOperation::multSTensor3(Ce, corKir, temp4);
        
        // Check S 
        commuteCheck = temp1;
        commuteCheck -= temp2;
        
        // Check corKir (just for fun)
        temp6 = temp3;
        temp6 -= temp4;
        
        // RANK OF A ZERO TENSOR is ZERO - maybe use this?
    };
    
    
    void mlawNonLinearTVP::getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const{
    
        double p = q1->_epspbarre; // eq.plastic strain
        
        coeffs.resize(2);
        
        coeffs(0) = 1; //DEBUGGING
        coeffs(1) = 0;
        
        // add dependency on p later
    };
    
    void mlawNonLinearTVP::getPlasticPotential(const double& phiP, const double& phiE, double& Px) const{
        
        Px = pow(phiE,2) + _b*pow(phiP,2);
    };
    
    void mlawNonLinearTVP::getYieldCoefficientTempDers(const IPNonLinearTVP *q1, const double& T, fullVector<double>& DcoeffsDT, fullVector<double>& DDcoeffsDTT) const{
      double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
      sigc = q1->_ipCompression->getR();
      Hc = q1->_ipCompression->getDR();
      dsigcdT = q1->_ipCompression->getDRDT();
      ddsigcdTT = q1->_ipCompression->getDDRDTT();
      
    
      double sigt(0.), Ht(0.), dsigtdT(0.), ddsigtdTT(0.);
      sigt = q1->_ipTraction->getR();
      Ht = q1->_ipTraction->getDR();
      dsigtdT = q1->_ipTraction->getDRDT();
      ddsigtdTT = q1->_ipTraction->getDDRDTT();
    
    
      DcoeffsDT.resize(3);
      DDcoeffsDTT.resize(3);
      double m = sigt/sigc;
      double DmDT = (dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc); 
      double DDmDTT = (ddsigtdTT*sigc - sigt*ddsigcdTT)/(sigc*sigc) - 2*(dsigtdT*sigc- sigt*dsigcdT)/(sigc*sigc*sigc)*dsigcdT;
      double Da1Dm = 3./sigc*(_n*pow(m,_n-1.)/(m+1.) - (pow(m,_n)-1.)/(m+1.)/(m+1.));
      double DDa1Dmm = 3./sigc*(_n*(_n-1)*pow(m,_n-2.)/(m+1.) - _n*(pow(m,_n)-1.)/(m+1.)/(m+1.) - _n*pow(m,_n-1)/(m+1.)/(m+1.) + 2.*(pow(m,_n)-1.)/(m+1.)/(m+1.)/(m+1.));
      double term = _n*pow(m,_n) + _n*pow(m,_n-1.) - pow(m,_n);
      double dterm = pow(_n,2)*pow(m,_n-1.) + _n*(_n-1)*pow(m,_n-2.) - _n*pow(m,_n-1);
    
      DcoeffsDT(2) = -_n*pow(sigc,-_n-1.)*dsigcdT;
      DcoeffsDT(1) = Da1Dm*DmDT - 3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*dsigcdT;
      DcoeffsDT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DmDT;
      
      DDcoeffsDTT(2) = -_n*pow(sigc,-_n-1.)*ddsigcdTT -_n*(-_n-1)*pow(sigc,-_n-2.)*dsigcdT*dsigcdT;
      DDcoeffsDTT(1) = -1/sigc*dsigcdT*DcoeffsDT(1) + 3/sigc*( (dterm/(m+1.)/(m+1.) -2*term/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT + term/(m+1.)/(m+1.)*DDmDTT
                       + 1/sigc*( (1/sigc*pow(dsigtdT,2) - ddsigcdTT) * (pow(m,_n)-1.)/(m+1.)  + dsigtdT * term/(m+1.)/(m+1.) ) );
      DDcoeffsDTT(0) = ((_n*pow(m,_n-1)+1.)/(m+1) - (pow(m,_n)+m)/(m+1.)/(m+1.))*DDmDTT 
                        + ((_n*(_n-1)*pow(m,_n-2))/(m+1) + 2*(pow(m,_n)+m)/(m+1.)/(m+1.)/(m+1.))*DmDT*DmDT;
        
    };
    
    void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T,
                                          double *psiTVM, double *DpsiTVMdT, double *DDpsiTVMdTT) const{
    
        // get some Things
        const double& Gamma = q1->_Gamma;
        const double& DgammaDT = q1->_DgammaDT;
        const double& dGammaDT = q1->_DGammaDT;
        const STensor3& Me = q1->_ModMandel;
        const STensor3& X = q1->_backsig;
        const STensor3& dXdT  = q1->_DbackSigDT;
        
        STensor3 devMe, devX, dDevXdT, dDevMeDT;
        double trMe, trX, dTrXdT, dTrMeDT;
        STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe); 
        STensorOperation::decomposeDevTr(q1->_backsig,devX,trX); 
        STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT); 
        STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
        
        // Start
        double psiVE, psiVP, psiT = 0.;
        double DpsiVEdT, DpsiVPdT, DpsiTdT = 0.;
        double DDpsiVEdTT, DDpsiVPdTT, DpsiTdTT = 0.;
        
        double CpT, DCpDT, DDCpDTT;
        mlawNonLinearTVE::getCp(CpT,T,&DCpDT,&DDCpDTT);
        
        // psiT 
        psiT = CpT*((T-_Tinitial)-T*log(T/_Tinitial));
        DpsiTdT = -CpT*log(T/_Tinitial) - psiT/CpT*DCpDT;
        DpsiTdTT = -CpT/T - psiT/CpT*DDCpDTT;
        
        //
        // psiVP
        
        // get R, dRdT
        fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
        getYieldCoefficients(q1,a);
        getYieldCoefficientDerivatives(q1,q1->_nup,Da);
        getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
        double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
        sigc = q1->_ipCompression->getR();
        Hc = q1->_ipCompression->getDR();
        dsigcdT = q1->_ipCompression->getDRDT();
        ddsigcdTT = q1->_ipCompression->getDDRDTT();
        
        double R1(0.), R2(0.), R3(0.), dR1dT(0.), dR2dT(0.), dR3dT(0.), ddR1dTT(0.), ddR2dTT(0.), ddR3dTT(0.);
        
        const std::vector<double>& IsoHardForce =  q1->_IsoHardForce;
        const std::vector<double>& dIsoHardForcedT =  q1->_dIsoHardForcedT;
        const std::vector<double>& ddIsoHardForcedTT =  q1->_ddIsoHardForcedTT;
        
        R1 = IsoHardForce[0]; R2 = IsoHardForce[1]; R3 = IsoHardForce[2];
        dR1dT = dIsoHardForcedT[0]; dR2dT = dIsoHardForcedT[1]; dR3dT = dIsoHardForcedT[2];
        ddR1dTT = ddIsoHardForcedTT[0]; ddR2dTT = ddIsoHardForcedTT[1]; ddR3dTT = ddIsoHardForcedTT[2];
        
         // get Kinematic Hardening stuff
        double Hb(0.), dHbdT(0.), ddHbddT(0.);
        if (q1->_ipKinematic != NULL){
            Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
            dHbdT = q1->_ipKinematic->getDRDT();
            ddHbddT = q1->_ipKinematic->getDDRDDT();
        }
    
        double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
        static STensor3 alpha;
        alpha = X;
        alpha *= 1/(pow(kk,2)*Hb);
        
        // ddXdTT 
        static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
        STensorOperation::zero(ddXdTT);
        double delT = T-T0;
        if (delT>1e-10){
         for (int i=0; i<3; i++)
           for (int j=0; j<3; j++){
             ddXdTT(i,j) += 2*dXdT(i,j)/(T-T0) - 2*(X(i,j)-q0->_backsig(i,j))/pow((T-T0),2);
           }
        }
        
        static STensor3 dAlphadT, ddAlphadTT;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) - X(i,j)*1/(pow(kk*Hb,2))*dHbdT;
            ddAlphadTT(i,j) = ddXdTT(i,j)*1/(pow(kk,2)*Hb) - 2*dXdT(i,j)*1/(pow(kk*Hb,2))*dHbdT 
                            - X(i,j)*( -2/(pow(kk*Hb,3))*dHbdT*dHbdT +  1/(pow(kk*Hb,2))*ddHbddT);
          }
    
        // finally
        STensorOperation::doubleContractionSTensor3(X,alpha,psiVP);
        psiVP *= 0.5;
        psiVP += 0.5*( 1/_HR[0]*R1*a(1)*sigc + 1/_HR[1]*pow(R2,2) + 1/_HR[2]*pow(R3,2) );
        
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            DpsiVPdT += 0.5*(dXdT(i,j)*alpha(i,j)+X(i,j)*dAlphadT(i,j));
            DDpsiVPdTT += 0.5*(ddXdTT(i,j)*alpha(i,j) + 2*dXdT(i,j)*dAlphadT(i,j) + X(i,j)*ddAlphadTT(i,j));
          }
        DpsiVPdT += 0.5*( 1/_HR[0]*(dR1dT*a(1)*sigc + R1*(DaDT(1)*sigc + a(1)*dsigcdT)) + 2/_HR[1]*R2*dR2dT + 2/_HR[2]*R3*dR3dT );
        
        DDpsiVPdTT += 0.5*( 1/_HR[0]*(ddR1dTT*a(1)*sigc + 2*dR1dT*(DaDT(1)*sigc + a(1)*dsigcdT) + R1*(DDaDTT(1)*sigc + 2*DaDT(1)*dsigcdT + a(1)*ddsigcdTT))
                    + 2/_HR[1]*R2*ddR2dTT + 2/_HR[2]*R3*ddR3dTT  + 2/_HR[1]*dR2dT*dR2dT + 2/_HR[2]*dR3dT*dR3dT); // CHANGE!!  2nd derivative of R and X are difficult
        
        //
        // psiVE
        
        double KT, GT, cteT, dKdT, ddKdTT, dGdT, ddGdTT, DcteDT, DDcteDTT; getK(KT,T,&dKdT,&ddKdTT); getG(GT,T,&dGdT,&ddGdTT); getAlpha(cteT,T,&DcteDT,&DDcteDTT);
        
        const STensor3& E = q1->getConstRefToElasticStrain();
        static STensor3 devKeinf, DdevKeinfDT, DDdevKeinfDTT, devEe;
        static double trKeinf, DtrKeinfDT, DDtrKeinfDTT, trEe;
        STensorOperation::decomposeDevTr(E,devEe,trEe);
        mlawNonLinearTVE::corKirInfinity(devEe,trEe,T,devKeinf,trKeinf);
        
        DtrKeinfDT = trKeinf*dKdT/KT - 3*KT*(DcteDT*(T-_Tinitial) + cteT);
        DDtrKeinfDTT = trKeinf*ddKdTT/KT - 6*dKdT*(DcteDT*(T-_Tinitial) + cteT)- 3*KT*(DDcteDTT*(T-_Tinitial) + 2*DcteDT);
        DdevKeinfDT = devKeinf;
        DdevKeinfDT *= dGdT/GT;
        DDdevKeinfDTT = devKeinf;
        DDdevKeinfDTT *= ddGdTT/GT;
    
        const std::vector<STensor3>& devOi = q1->getConstRefToDevViscoElasticOverStress();
        const std::vector<double>& trOi = q1->getConstRefToTrViscoElasticOverStress();
        const std::vector<STensor3>& devDOiDT = q1->getConstRefToDevDOiDT();
        const std::vector<double>& trDOiDT = q1->getConstRefToTrDOiDT();
        const std::vector<STensor3>& devDDOiDTT = q1->getConstRefToDevDDOiDTT();
        const std::vector<double>& trDDOiDTT = q1->getConstRefToTrDDOiDTT();
            
        psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);
        
        DpsiVEdT = 1/9 * (1/KT * trKeinf * DtrKeinfDT - 1/(2*pow(KT,2)) * dKdT * pow(trKeinf,2));
        DDpsiVEdTT = 1/9 * (1/KT * (pow(DtrKeinfDT,2)+trKeinf*DDtrKeinfDTT) - 1/pow(KT,2)*dKdT*trKeinf*DtrKeinfDT 
                            - 1/(2*pow(KT,2)) *(ddKdTT*pow(trKeinf,2) + 2*dKdT*trKeinf*DtrKeinfDT) 
                            + 1/pow(KT,3)*pow(dKdT,2)*pow(trKeinf,2)); 
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            DpsiVEdT += 2/GT * devKeinf(i,j)*DdevKeinfDT(i,j) - 1/pow(GT,2) * dGdT * devKeinf(i,j)*devKeinf(i,j);
            DDpsiVEdTT += 2/GT*(DdevKeinfDT(i,j)*DdevKeinfDT(i,j) + devKeinf(i,j)*DDdevKeinfDTT(i,j) - dGdT*devKeinf(i,j)*DdevKeinfDT(i,j)) 
                          - 1/pow(GT,2) * (ddGdTT*devKeinf(i,j)*devKeinf(i,j) + 2*dGdT*devKeinf(i,j)*DdevKeinfDT(i,j));
          }
        if ((_Ki.size() > 0) or (_Gi.size() > 0)){
          for (int k=0; k<_Gi.size(); k++){
            DpsiVEdT += 1/9 * 1/_Ki[k] * trOi[k] * trDOiDT[k];
            DDpsiVEdTT += 1/9 * 1/_Ki[k] * (trOi[k]*trDDOiDTT[k] + pow(trDOiDT[k],2));
            for (int i=0; i<3; i++)
              for (int j=0; j<3; j++){
                DpsiVEdT += 2/_Gi[k] * devOi[k](i,j) * devDOiDT[k](i,j);
                DDpsiVEdTT += 2/_Gi[k] * (devOi[k](i,j) * devDDOiDTT[k](i,j) + pow(devDOiDT[k](i,j),2));
              }    
            }
        }
        
        // FINALLY
        if(psiTVM!=NULL){
            *psiTVM = psiT + psiVE + psiVP;
        }
        if(DpsiTVMdT!=NULL){
            *DpsiTVMdT = DpsiTdT + DpsiVEdT + DpsiVPdT;
        }
        if(DDpsiTVMdTT!=NULL){
            *DDpsiTVMdTT = DpsiTdTT + DDpsiVEdTT + DDpsiVPdTT;
        }
    
    };
    
    void mlawNonLinearTVP::getIsotropicHardeningForce(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
                                          const STensor3& DphiPDF, std::vector<double>* ddRdTT, std::vector<STensor3>* ddRdFdT) const{
         
        const double& Gamma_0 = q0->_Gamma;
        const double& Gamma = q1->_Gamma;
        const double& DgammaDT = q1->_DgammaDT;
        const double& dGammaDT_0 = q0->_DGammaDT;
        const double& dGammaDT = q1->_DGammaDT; // = 0. for debugging //DEBUGGING
        const STensor3& DgammaDF = q1->_DgammaDF;
        const STensor3& DGammaDF = q1->_DGammaDF;
        const double& gamma0 = q0->_epspbarre;
        const double& gamma1 = q1->_epspbarre;
        
        static STensor3 devMe, devX1, dDevXdT, dDevMeDT;
        static double trMe, trX1, dTrXdT, dTrMeDT;
        static double trMe_0, trX1_0, dTrXdT_0, dTrMeDT_0;
        
        STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe); 
        STensorOperation::decomposeDevTr(q0->_ModMandel,devMe,trMe_0); 
        STensorOperation::decomposeDevTr(q1->_backsig,devX1,trX1); 
        STensorOperation::decomposeDevTr(q0->_backsig,devX1,trX1_0); 
        STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT); 
        STensorOperation::decomposeDevTr(q0->_DbackSigDT,dDevXdT,dTrXdT_0); 
        STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
        STensorOperation::decomposeDevTr(q0->_DModMandelDT,dDevMeDT,dTrMeDT_0);
        
        trMe /= 3; trX1 /= 3; dTrMeDT /= 3; dTrXdT /= 3;  
        trMe_0 /= 3; trX1_0 /= 3; dTrMeDT_0 /= 3; dTrXdT_0 /= 3;  
         
        // dTrMeDT = 0.; dTrXdT = 0.; // for debugging only //DEBUGGING
        
        // get Dgamma
        double Dgamma = gamma1 - gamma0;
        
        // get R
        double eta(0.),Deta(0.),DetaDT(0.),DDetaDTT(0.);
        fullVector<double> a(3), Da(3), DaDT(3), DDaDTT(3);
        getYieldCoefficients(q1,a);
        getYieldCoefficientDerivatives(q1,q1->_nup,Da);
        getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
        
        if (_viscosity != NULL)
            _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT,DDetaDTT);
        double etaOverDt = eta/this->getTimeStep();
        double viscoTerm = etaOverDt*Gamma;
        
        const std::vector<double>& R0 = q0->_IsoHardForce;
        std::vector<double>& R = q1->_IsoHardForce; 
        const std::vector<double>& dRdT_0 = q0->_dIsoHardForcedT; 
        std::vector<double>& dRdT = q1->_dIsoHardForcedT;
        R.resize(3,0.);
        dRdT.resize(3,0.);
        
        double sigc(0.), Hc(0.), dsigcdT(0.), ddsigcdTT(0.);
        sigc = q1->_ipCompression->getR();
        Hc = q1->_ipCompression->getDR();
        dsigcdT = q1->_ipCompression->getDRDT();
        ddsigcdTT = q1->_ipCompression->getDDRDTT();
        
        R[0] = - a(1) * (trMe-trX1) * sigc; 
        dRdT[0] = - DaDT(1)*(trMe-trX1)*sigc - a(1)*( (dTrMeDT-dTrXdT)*sigc + (trMe-trX1) * dsigcdT);
        R[1] = - a(0) * sigc; 
        dRdT[1] = - (DaDT(0) * sigc + a(0) * dsigcdT);
        if (Gamma>0 and etaOverDt>0){
           R[2] = - pow(eta*Gamma/this->getTimeStep(),_p) * sigc;
           dRdT[2] = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep() * sigc - pow(eta*Gamma/this->getTimeStep(),_p) * dsigcdT; 
        }
         
        // dRdF
        // static STensor3 dRdF;
        const std::vector<STensor3>& dRdF_0 = q0->_dIsoHardForcedF;
        std::vector<STensor3>& dRdF = q1->_dIsoHardForcedF;
        dRdF.resize(3,0.);
        
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dRdF[0](i,j) = -Da(1)*(trMe-trX1)*DgammaDF(i,j)*sigc - a(1)*(DphiPDF(i,j)*sigc + (trMe-trX1)*Hc*DgammaDF(i,j));  
            dRdF[1](i,j) = -Da(0)*DgammaDF(i,j)*sigc - a(0)*Hc*DgammaDF(i,j); 
            }
              
        if (Gamma>0 and etaOverDt>0){
          for (int i=0; i<3; i++)
            for (int j=0; j<3; j++){
              dRdF[2](i,j) = - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DGammaDF(i,j)/this->getTimeStep())*sigc - pow(eta*Gamma/this->getTimeStep(),_p) * Hc*DgammaDF(i,j); 
            }
        }
        
        // ddRdTT
        double delT = T-T0;
        if(ddRdTT!=NULL){
           
           ddRdTT->at(0) = - DDaDTT(1)*(trMe-trX1)*sigc - a(1)*( 2*(dTrMeDT-dTrXdT)*dsigcdT + (trMe-trX1) * ddsigcdTT );
           ddRdTT->at(1) = -( DDaDTT(0)*sigc + 2*DaDT(0)*dsigcdT + a(0)*ddsigcdTT);
           
           if (Gamma>0 and etaOverDt>0){
            ddRdTT->at(2) = - _p*(_p-1)*pow(eta*Gamma/this->getTimeStep(),_p-2)*pow(((DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep()),2) * sigc
                            - _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*((DDetaDTT*Gamma + 2*DetaDT*dGammaDT)/this->getTimeStep() * sigc + 
                            2*(DetaDT*Gamma + eta*dGammaDT)/this->getTimeStep()* dsigcdT)
                            - pow(eta*Gamma/this->getTimeStep(),_p) * ddsigcdTT;
           }
           
           if (delT>1e-10){
            double DphiPDT = dTrMeDT- dTrXdT; double DphiPDT_0 = dTrMeDT_0 - dTrXdT_0;
            double DDphiPDTT = (DphiPDT-DphiPDT_0)/delT; // (DphiPDT - (trMe-trMe_0+trX1_0-trX1)/delT)/delT; // 2*DphiPDT/delT - 2*(trMe-trMe_0+trX1_0-trX1)/pow(delT,2); // DEBUGGING
            ddRdTT->at(0) -= a(1) * DDphiPDTT * sigc;
            
            if (Gamma>0 and etaOverDt>0){
             double DDGammaDTT = (dGammaDT-dGammaDT_0)/delT; // 2*dGammaDT/delT - 2*(Gamma-Gamma_0)/pow(delT,2);
             ddRdTT->at(2) -= _p*pow(eta*Gamma/this->getTimeStep(),_p-1)*(eta*DDGammaDTT)/this->getTimeStep() * sigc;
            } 
           }
           q1->_ddIsoHardForcedTT = *ddRdTT;
        }
        // ddRdTT = 2*dRdT/(T-T0) - 2*(R-R0)/pow((T-T0),2); // ddGammaDTT is very difficult to estimate
        
        // ddRdTF
        if(ddRdFdT!=NULL){
     
          for (int k=0; k<3; k++){
            STensorOperation::zero(ddRdFdT->at(k));
            if(delT>1e-10){
              ddRdFdT->at(k) = dRdF[k];
              ddRdFdT->at(k) -= dRdF_0[k];
              ddRdFdT->at(k) *= 1/delT;
            }
          }  
        }
    };
                                          
    void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double T0, const double& T,
                                            const STensor3& dFpdT, const STensor43& dFpdF, const STensor3& DphiPDF,
                                            double *Wm, STensor3 *dWmdF, double *dWmdT) const{
        
        const double& Gamma = q1->_Gamma;                                        
        const double& DgammaDT = q1->_DgammaDT;
        const double& dGammaDT = q1->_DGammaDT;
        const STensor3& DgammaDF = q1->_DgammaDF;
        const STensor3& DGammaDF = q1->_DGammaDF;
        const double& gamma0 = q0->_epspbarre;
        const double& gamma1 = q1->_epspbarre;
        const STensor3& Fp0 = q0->_Fp;
        const STensor3& Fp1 = q1->_Fp;
        const STensor3& Me = q1->_ModMandel;
        const STensor3& X1 = q1->_backsig;
        const STensor3& X0 = q0->_backsig;
        const STensor3& dXdT = q1->_DbackSigDT; // Just for debugging = 0. //DEBUGGING
        const STensor3& dXdT_0  = q0->_DbackSigDT;
        const STensor43& dXdF_0  = q0->_DbackSigDF;
        const STensor43& dXdF  = q1->_DbackSigDF;
        const STensor3& dMeDT  = q1->_DModMandelDT;
        const STensor3& dMeDT_0  = q0->_DModMandelDT;
        const STensor43& dMedF  = q1->_DModMandelDF;
        const STensor43& dMedF_0  = q0->_DModMandelDF;
        
        static STensor3 devMe, devX1, dDevXdT, dDevMeDT;
        static double trMe, trX1, dTrXdT, dTrMeDT;
        
        STensorOperation::decomposeDevTr(q1->_ModMandel,devMe,trMe); 
        STensorOperation::decomposeDevTr(q1->_backsig,devX1,trX1); 
        STensorOperation::decomposeDevTr(q1->_DbackSigDT,dDevXdT,dTrXdT); 
        STensorOperation::decomposeDevTr(q1->_DModMandelDT,dDevMeDT,dTrMeDT);
        
        // get Dgamma
        double Dgamma = gamma1 - gamma0;
        
        // get Hb
        double Hb(0.), dHb(0), dHbdT(0.), ddHbddT(0.);
        if (q1->_ipKinematic != NULL){
            Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
            dHb = q1->_ipKinematic->getDDR();
            dHbdT = q1->_ipKinematic->getDRDT();
            ddHbddT = q1->_ipKinematic->getDDRDDT();
        }
    
        // get Dp
        static STensor3 Fpinv, Fp_dot, Dp;
        STensorOperation::inverseSTensor3(Fp1,Fpinv);
        Fp_dot = Fp1;
        Fp_dot -= Fp0;
        Fp_dot *= 1/this->getTimeStep();
        STensorOperation::multSTensor3(Fp_dot,Fpinv,Dp);
       
        // get alpha 
        double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
        static STensor3 alpha0, alpha1;
        alpha0 = X0; alpha1 = X1;
        alpha0 *= 1/(pow(kk,2)*Hb); alpha1 *= 1/(pow(kk,2)*Hb); // how to get alpha0-current or previous Hb??? WHAT?  
        alpha1 -= alpha0;
        alpha1 *= 1/this->getTimeStep();
        
        // get TempX - convenient backStress tensor
        STensor3 tempX;
        tempX = X1;
        tempX -= T*dXdT; 
        
        // get R, dRdT, dRdF
        std::vector<double> ddIsoHardForcedTT; ddIsoHardForcedTT.resize(3,0.);
        std::vector<STensor3> ddIsoHardForcedFdT; ddIsoHardForcedFdT.resize(3,0.);
        
        getIsotropicHardeningForce(q0,q1,T0,T,DphiPDF,&ddIsoHardForcedTT,&ddIsoHardForcedFdT);
        const std::vector<double>& IsoHardForce =  q1->_IsoHardForce;
        const std::vector<double>& dIsoHardForcedT =  q1->_dIsoHardForcedT;
        const std::vector<STensor3>& dIsoHardForcedF =  q1->_dIsoHardForcedF;
        
        double R(0.), dRdT(0.), ddRdTT(0.);
        STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
        
        for (int i=0; i<3; i++){
          R += IsoHardForce[i];
          dRdT += dIsoHardForcedT[i];
          ddRdTT += ddIsoHardForcedTT[i];
          dRdF += dIsoHardForcedF[i];
          ddRdFdT += ddIsoHardForcedFdT[i]; //DEBUGGING
        }
        
        R -= T*dRdT;
        
        // mechSourceTVP
        if(Wm!=NULL){
            double mechSourceTVP = 0.;
            if (this->getTimeStep() > 0){
              mechSourceTVP += STensorOperation::doubledot(Me,Dp);
              mechSourceTVP -= STensorOperation::doubledot(tempX,alpha1);
              mechSourceTVP -= (R*Dgamma)/this->getTimeStep();
            }
            *Wm = mechSourceTVP;
            // Second Term TBD - need dXdT, dRdT -> added above
        }
        
        // Some Generic Conversion Tensors
        static STensor43 dDpdFp, dFpinvdFp;
        static STensor3 dFpinvdT;
        
        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++){
                 dFpinvdFp(i,j,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int s=0; s<3; s++)
                    dFpinvdFp(i,j,k,l) -= Fpinv(i,m)*_I4(m,s,k,l)*Fpinv(s,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++){
                dDpdFp(i,j,k,l) = 0.;
                for (int m=0; m<3; m++)
                  dDpdFp(i,j,k,l) += 1/this->getTimeStep()*_I4(i,m,k,l)*Fpinv(m,j) + Fp_dot(i,m)*dFpinvdFp(m,j,k,l);
              }
        
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dFpinvdT(i,j) = 0.;
            for (int m=0; m<3; m++)
              for (int s=0; s<3; s++)
                dFpinvdT(i,j) -= Fpinv(i,m)*dFpdT(m,s)*Fpinv(s,j);
              }
              
        // dDpdF and dDpdT
        static STensor43 dDpdF;   
        // STensorOperation::zero(dDpdF);
        STensorOperation::multSTensor43(dDpdFp, dFpdF, dDpdF);
              
        static STensor3 dDpdT;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dDpdT(i,j) = 0.;
            for (int m=0; m<3; m++){
              dDpdT(i,j) += dFpdT(i,m)*Fpinv(m,j)*1/this->getTimeStep() + Fp_dot(i,m)*dFpinvdT(m,j); // dDpdFp(i,j,k,l)*dFpdT(k,l);  // maybe use Dp = Fp_dot Fpinv
          }}
          
       // dXdF is available 
       // const STensor43 dXdF = q1->_DbackSigDF;
       
       // ddXdTF -> make it -- CHANGE!!
       static STensor43 ddXdTF; // ddGammaDTF and ddQDTF are very difficult to estimate
       static STensor3 ddXdTT; // ddGammaDTT and ddQDTT are very difficult to estimate
       STensorOperation::zero(ddXdTF);
       STensorOperation::zero(ddXdTT);
       double delT = T-T0;
       if (delT>1e-10){
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            ddXdTT(i,j) += (dXdT(i,j)-dXdT_0(i,j))/delT; // 2*dXdT(i,j)/(T-T0) - 2*(X1(i,j)-X0(i,j))/pow((T-T0),2); // = 0. in DEBUGGING 
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
               ddXdTF(i,j,k,l) += (dXdF(i,j,k,l)-dXdF_0(i,j,k,l))/delT; // = 0. in DEBUGGING
              }
            }
          }
       
       // dAlphadF - from dXdF
       static STensor43 dAlphadF;
       STensorOperation::zero(dAlphadF);
        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++){
                dAlphadF(i,j,k,l) += dXdF(i,j,k,l)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDF(k,l)*1/(pow(kk*Hb,2))*dHb;
              }
              
       // dAlphadT - from dXdT   
       static STensor3 dAlphadT;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dAlphadT(i,j) = dXdT(i,j)*1/(pow(kk,2)*Hb) + X1(i,j)*DgammaDT*1/(pow(kk*Hb,2))*dHbdT;
            }
       
       
       // dMedF is available
       // const STensor43 dMedF = q1->_DModMandelDF;
    
        if(dWmdF!=NULL){
            STensor3 dmechSourceTVPdF(0.);
            if (this->getTimeStep() > 0){
              for (int i=0; i<3; i++)
                for (int j=0; j<3; j++){
                  dmechSourceTVPdF(i,j) = - (dRdF(i,j) - T*ddRdFdT(i,j))*Dgamma/this->getTimeStep() - R*DgammaDF(i,j)/this->getTimeStep();
                  for (int k=0; k<3; k++)
                    for (int l=0; l<3; l++){
                      dmechSourceTVPdF(i,j) += dMedF(i,j,k,l)*Dp(i,j) + Me(i,j)*dDpdF(i,j,k,l) 
                                            - (dXdF(i,j,k,l) - T*ddXdTF(i,j,k,l))*alpha1(i,j) - tempX(i,j)*dAlphadF(i,j,k,l)/this->getTimeStep();
                    }
                }
            }  
            *dWmdF = dmechSourceTVPdF;
        }
        
        if(dWmdT!=NULL){
            double dmechSourceTVPdT(0.);
            if (this->getTimeStep() > 0){
              for (int i=0; i<3; i++)
                for (int j=0; j<3; j++){
                  dmechSourceTVPdT += dMeDT(i,j)*Dp(i,j) + Me(i,j)*dDpdT(i,j) 
                                            + T*ddXdTT(i,j)*alpha1(i,j) - tempX(i,j)*dAlphadT(i,j)/this->getTimeStep()
                                            + T*ddRdTT*Dgamma/this->getTimeStep() - R*DgammaDT/this->getTimeStep() ;
                }
            }
            *dWmdT = dmechSourceTVPdT;
        }
    };
    
    // NEW
    
    void mlawNonLinearTVP::constitutive( // This function is just a placeholder, defined in the pure virtual class -> does nothing, its never called.
                                const STensor3& F0,         // previous deformation gradient (input @ time n)
                                const STensor3& F1,         // current deformation gradient (input @ time n+1)
                                      STensor3 &P1,                // current 1st Piola-Kirchhoff stress tensor (output)
                                const IPVariable *q0i,       // array of previous internal variables
                                      IPVariable *q1i,             // current array of internal variable (in ipvcur on output),
                                      STensor43 &Tangent,         // tangents (output)
                                const bool stiff,            // if true compute the tangents
                                      STensor43* elasticTangent, 
                                const bool dTangent,
                                      STensor63* dCalgdeps) const{
        static SVector3 gradT, temp2;
        static STensor3 temp3;
        static STensor33 temp33;    
        static double tempVal;
        static STensor43 dFpdF, dFedF;
        static STensor3 dFpdT, dFedT;
    
        const IPNonLinearTVP *q0 = dynamic_cast<const IPNonLinearTVP *> (q0i);
              IPNonLinearTVP *q1 = dynamic_cast<IPNonLinearTVP *>(q1i);
              
        if(elasticTangent != NULL) Msg::Error("mlawNonLinearTVP elasticTangent not defined"); 
        predictorCorrector_ThermoViscoPlastic(F0,F1,P1,q0,q1,Tangent,dFpdF, dFpdT, dFedF, dFedT,
                                              _Tref,_Tref,gradT,gradT,temp2,temp3,temp3,temp2,temp33,
                                              tempVal,tempVal,temp3,
                                              tempVal,tempVal,temp3,
                                              stiff, NULL);
    };
                
    void mlawNonLinearTVP::constitutive(
                                const STensor3& F0,                     // initial deformation gradient (input @ time n)
                                const STensor3& F1,                     // 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 mechanical tangents (output)
                                const double& T0,                       // previous temperature
                                const double& T,                        // temperature
                                const SVector3 &gradT0,                 // previous temperature gradient
                                const SVector3 &gradT,                  // temperature gradient
                                      SVector3 &fluxT,                  // temperature flux
                                      STensor3 &dPdT,                   // mechanical-thermal coupling
                                      STensor3 &dfluxTdgradT,           // thermal tengent
                                      SVector3 &dfluxTdT,
                                      STensor33 &dfluxTdF,              // thermal-mechanical coupling
                                      double &thermalSource,            // - Cp*dTdt
                                      double &dthermalSourcedT,         // thermal source
                                      STensor3 &dthermalSourcedF,
                                      double &mechanicalSource,         // mechanical source--> convert to heat
                                      double &dmechanicalSourcedT,
                                      STensor3 &dmechanicalSourceF,
                                const bool stiff,
                                      STensor43* elasticTangent) const{
    
        const IPNonLinearTVP *q0 = dynamic_cast<const IPNonLinearTVP *>(q0i);
              IPNonLinearTVP *q1 = dynamic_cast<IPNonLinearTVP *>(q1i);
    
        // Declaring here, coz I dont know better.
        static STensor43 dFedF, dFpdF;
        static STensor3 dFpdT, dFedT;
            
        if (!_tangentByPerturbation){
    
            this->predictorCorrector_ThermoViscoPlastic(F0, F1, P, q0, q1, Tangent, dFpdF, dFpdT, dFedF, dFedT,
                                                        T0, T, gradT0, gradT, fluxT, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
                                                        thermalSource, dthermalSourcedT, dthermalSourcedF,
                                                        mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF,
                                                        stiff, elasticTangent);
        }
            
        else{
    
            this->predictorCorrector_ThermoViscoPlastic(F0, F1, P, q0, q1, T0, T, gradT0, gradT, fluxT, thermalSource, mechanicalSource, elasticTangent);
            if (stiff)
                this->tangent_full_perturbation(F0,F1,P,q0,q1,T0, T, gradT0, gradT, fluxT, thermalSource, mechanicalSource,
                                                Tangent, dFpdF, dFpdT, dFedF, dFedT,
                                                dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF,
                                                dthermalSourcedT, dthermalSourcedF,
                                                dmechanicalSourcedT, dmechanicalSourceF);            
        }
    };
    
    
    void mlawNonLinearTVP::predictorCorrector_ThermoViscoPlastic(
                                                        const STensor3& F0,         // initial deformation gradient (input @ time n)
                                                        const STensor3& F1,         // updated deformation gradient (input @ time n+1)
                                                              STensor3 &P,                // updated 1st Piola-Kirchhoff stress tensor (output)
                                                        const IPNonLinearTVP *q0,       // array of initial internal variable
                                                              IPNonLinearTVP *q1,             // updated array of internal variable (in ipvcur on output),
                                                        const double& T0, // previous temperature
                                                        const double& T, // temperature
                                                        const SVector3 &gradT0, // previous temeprature gradient
                                                        const SVector3 &gradT, // temeprature gradient
                                                              SVector3 &fluxT, // temperature flux)
                                                              double &thermalSource,
                                                              double& mechanicalSource,
                                                              STensor43* elasticTangent) const{
      // temp variables
      static STensor43 Tangent, dFpdF, dFedF;
      static STensor3 dPdT, dfluxTdgradT, dthermalSourcedF, dmechanicalSourceF, dFpdT, dFedT;
      static STensor33 dfluxTdF;
      static SVector3 dfluxTdT;
      static double dthermalSourcedT, dmechanicalSourcedT;
      predictorCorrector_ThermoViscoPlastic(F0,F1,P,q0,q1,Tangent,dFpdF,dFpdT,dFedF,dFedT,
                                            T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,
                                            thermalSource,dthermalSourcedT,dthermalSourcedF,
                                            mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF,false,elasticTangent);
    };
    
    
    void mlawNonLinearTVP::tangent_full_perturbation(
                                const STensor3& F0,
                                const STensor3& F1,         // updated deformation gradient (input @ time n+1)
                                      STensor3& P,                // updated 1st Piola-Kirchhoff stress tensor (output)
                                const IPNonLinearTVP *q0,       // array of initial internal variable
                                      IPNonLinearTVP *q1,             // updated array of internal variable (in ipvcur on output),
                                const double& T0, // previous temperature
                                const double& T, // temperature
                                const SVector3 &gradT0, // previous temeprature gradient
                                const SVector3 &gradT, // temeprature gradient
                                      SVector3 &fluxT, // temperature flux)
                                      double &thermalSource,
                                      double &mechanicalSource,
                                      STensor43 &Tangent,         // mechanical tangents (output)
                                      STensor43 &dFpdF, // plastic tangent
                                      STensor3 &dFpdT, // plastic tangent
                                      STensor43 &dFedF, // elastic tangent
                                      STensor3 &dFedT, // elastic tangent
                                      STensor3 &dPdT, // mechanical-thermal coupling
                                      STensor3 &dfluxTdgradT, // thermal tengent
                                      SVector3 &dfluxTdT,
                                      STensor33 &dfluxTdF, // thermal-mechanical coupling
                                      double &dthermalSourcedT, // thermal source
                                      STensor3 &dthermalSourcedF,
                                      double &dmechanicalSourcedT,
                                      STensor3 &dmechanicalSourceF) const{
    
        static STensor3 Fplus, Pplus;
        static SVector3 fluxTPlus, gradTplus;
        static double thermalSourcePlus;
        static double mechanicalSourcePlus;
        static IPNonLinearTVP qPlus(*q0);
    
        // perturb F
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            Fplus = F1;
            Fplus(i,j) += _perturbationfactor;
            
            predictorCorrector_ThermoViscoPlastic(F0,Fplus,Pplus,q0,&qPlus,T0,T,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
            
            q1->_DgammaDF(i,j) = (qPlus._epspbarre - q1->_epspbarre)/_perturbationfactor;
            q1->_DirreversibleEnergyDF(i,j) = (qPlus._irreversibleEnergy - q1->_irreversibleEnergy)/_perturbationfactor;
            
            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);
                dFpdF(k,l,i,j) = (qPlus._Fp(k,l)-q1->_Fp(k,l))/_perturbationfactor;
                dFedF(k,l,i,j) = (qPlus._Fe(k,l)-q1->_Fe(k,l))/_perturbationfactor;
              }
              dfluxTdF(k,i,j) = (fluxTPlus(k) - fluxT(k))/_perturbationfactor;
            }
            
            q1->getRefToDElasticEnergyPartdF()(i,j)=(qPlus.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/_perturbationfactor;
            q1->getRefToDViscousEnergyPartdF()(i,j)=(qPlus.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/_perturbationfactor;
            q1->getRefToDPlasticEnergyPartdF()(i,j)=(qPlus.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/_perturbationfactor;
          
            dthermalSourcedF(i,j) = (thermalSourcePlus - thermalSource)/_perturbationfactor;
            dmechanicalSourceF(i,j) = (mechanicalSourcePlus - mechanicalSource)/_perturbationfactor;
          }
        }
        
        // perturb gradT
        double gradTpert = _perturbationfactor*T0/1e-3;
        for (int i=0; i<3; i++){
          gradTplus = gradT;
          gradTplus(i) += gradTpert;
          predictorCorrector_ThermoViscoPlastic(F0,F1,Pplus,q0,&qPlus,T0,T,gradT0,gradTplus,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
          for (int k=0; k<3; k++){
            dfluxTdgradT(k,i) = (fluxTPlus(k) - fluxT(k))/gradTpert;
          }
        }
    
        // perturb T
        double Tplus = T+ _perturbationfactor*T0;
        predictorCorrector_ThermoViscoPlastic(F0,F1,Pplus,q0,&qPlus,T0,Tplus,gradT0,gradT,fluxTPlus,thermalSourcePlus,mechanicalSourcePlus, NULL);
        
        q1->_DgammaDT = (qPlus._epspbarre - q1->_epspbarre)/(_perturbationfactor*T0);
        q1->_DirreversibleEnergyDT = (qPlus._irreversibleEnergy - q1->_irreversibleEnergy)/(_perturbationfactor*T0);
            
        for (int k=0; k<3; k++){
          for (int l=0; l<3; l++){
            dFpdT(k,l) = (qPlus._Fp(k,l) - q1->_Fp(k,l))/(_perturbationfactor*T0);
            dFedT(k,l) = (qPlus._Fe(k,l) - q1->_Fe(k,l))/(_perturbationfactor*T0);
            dPdT(k,l) = (Pplus(k,l) - P(k,l))/(_perturbationfactor*T0);
          }
          dfluxTdT(k) = (fluxTPlus(k) - fluxT(k))/(_perturbationfactor*T0);
        }
        
        q1->getRefToDElasticEnergyPartdT() = (qPlus.getConstRefToElasticEnergyPart()-q1->getConstRefToElasticEnergyPart())/(_perturbationfactor*T0);
        q1->getRefToDViscousEnergyPartdT() = (qPlus.getConstRefToViscousEnergyPart()-q1->getConstRefToViscousEnergyPart())/(_perturbationfactor*T0);
        q1->getRefToDPlasticEnergyPartdT() = (qPlus.getConstRefToPlasticEnergyPart()-q1->getConstRefToPlasticEnergyPart())/(_perturbationfactor*T0);
          
        dthermalSourcedT = (thermalSourcePlus - thermalSource)/(_perturbationfactor*T0);
        dmechanicalSourcedT = (mechanicalSourcePlus - mechanicalSource)/(_perturbationfactor*T0);
                                          
    };
    
    void mlawNonLinearTVP::predictorCorrector_ThermoViscoPlastic(
                                                        const STensor3& F0,         // initial deformation gradient (input @ time n)
                                                        const STensor3& F1,         // updated deformation gradient (input @ time n+1)
                                                              STensor3& P,                // updated 1st Piola-Kirchhoff stress tensor (output)
                                                        const IPNonLinearTVP *q0,       // array of initial internal variable
                                                              IPNonLinearTVP *q1,             // updated array of internal variable (in ipvcur on output),
                                                              STensor43 &Tangent,         // mechanical tangents (output)
                                                              STensor43 &dFpdF, // plastic tangent
                                                              STensor3 &dFpdT, // plastic tangent
                                                              STensor43 &dFedF, // elastic tangent
                                                              STensor3 &dFedT, // elastic tangent
                                                        const double& T0, // previous temperature
                                                        const double& T, // temperature
                                                        const SVector3 &gradT0, // previoustemeprature gradient
                                                        const SVector3 &gradT, // temeprature gradient
                                                              SVector3 &fluxT, // temperature flux
                                                              STensor3 &dPdT, // mechanical-thermal coupling
                                                              STensor3 &dfluxTdgradT, // thermal tengent
                                                              SVector3 &dfluxTdT,
                                                              STensor33 &dfluxTdF, // thermal-mechanical coupling
                                                              double &thermalSource,   // - cp*dTdt
                                                              double &dthermalSourcedT, // thermal source
                                                              STensor3 &dthermalSourcedF,
                                                              double &mechanicalSource, // mechanical source--> convert to heat
                                                              double &dmechanicalSourcedT,
                                                              STensor3 &dmechanicalSourceF,
                                                        const bool stiff,
                                                              STensor43* elasticTangent) const{
    
      if (_tangentByPerturbation){
        if (_nonAssociatedFlow){
          this->predictorCorrector_TVP_nonAssociatedFlow(F0,F1,q0,q1,P,false,Tangent,dFedF,dFpdF,dFedT,dFpdT, 
                                                            T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,      
                                                            thermalSource,dthermalSourcedT,dthermalSourcedF,
                                                            mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF);
        }
        else{
          this->predictorCorrector_TVP_AssociatedFlow(F1,q0,q1,P,false,Tangent,dFedF,dFpdF,T0,T);
        }
    
    //    if (stiff){                                                                      // WHAT?????  WHY TWICE? (The same line happens above as STIFF is always TRUE)
    //      tangent_full_perturbation(Tangent,dFedF,dFpdF,P,F,q0,q1);
    //    }
    
      }
      else{
        if (_nonAssociatedFlow){
          this->predictorCorrector_TVP_nonAssociatedFlow(F0,F1,q0,q1,P,stiff,Tangent,dFedF,dFpdF,dFedT,dFpdT, 
                                                            T0,T,gradT0,gradT,fluxT,dPdT,dfluxTdgradT,dfluxTdT,dfluxTdF,      
                                                            thermalSource,dthermalSourcedT,dthermalSourcedF,
                                                            mechanicalSource,dmechanicalSourcedT,dmechanicalSourceF);
        }
        else{
          this->predictorCorrector_TVP_AssociatedFlow(F1,q0,q1,P,stiff,Tangent,dFedF,dFpdF,T0,T);
        }
      }
      
      /*
      // compute mechanical tengent
      if (elasticTangent!= NULL)
      {
        STensor3 Ptmp(0.);
        static IPNonLinearTVP q1tmp(*q1);
        q1tmp.operator=(*(dynamic_cast<const IPVariable*> (q1)));
        mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), Ptmp, q0, &q1tmp, *elasticTangent, dFedF, dFedT,
                                                          T0, T, gradT0, gradT, fluxT, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF, 
                                                          thermalSource, dthermalSourcedT, dthermalSourcedF, 
                                                          mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF, 
                                                          true);
                                                          // mechSource requires taylorQuinney Term - CHANGE 
        // UPDATE FOR TVP
        // thermalSource derivatives -DEFINE LATER
        
        double CpT; 
        if (!stiff){ // if stiff, we dont need this technically
           getCp(CpT,T);
        }
        else{
            // thermalSource TVP  
            CpT = -T*q1->getConstRefToDDFreeEnergyTVMdTT();
        }
        
        // thermalSource
        if (this->getTimeStep() > 0.){
            thermalSource = -CpT*(T-T0)/this->getTimeStep();
        }
        else{thermalSource = 0.;}
      
        // mechSource - TVP
        mechanicalSource += q1->getConstRefToMechSrcTVP();
        dmechanicalSourcedT += q1->getConstRefTodMechSrcTVPdT();
        dmechanicalSourceF += q1->getConstRefTodMechSrcTVPdF();
    
      }*/
    
    };
    
    void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3& F0, const STensor3& F, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, STensor3& dFedT, STensor3& dFpdT, 
                                const double T0, 
                                const double T,                     
                                const SVector3 &gradT0,                 // previous temperature gradient
                                const SVector3 &gradT,                  // temperature gradient
                                      SVector3 &fluxT,                  // temperature flux
                                      STensor3 &dPdT,                   // mechanical-thermal coupling
                                      STensor3 &dfluxTdgradT,           // thermal tengent
                                      SVector3 &dfluxTdT,
                                      STensor33 &dfluxTdF,              // thermal-mechanical coupling
                                      double &thermalSource,            // - Cp*dTdt
                                      double &dthermalSourcedT,         // thermal source
                                      STensor3 &dthermalSourcedF,
                                      double &mechanicalSource,         // mechanical source--> convert to heat
                                      double &dmechanicalSourcedT,
                                      STensor3 &dmechanicalSourceF) 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->_DbackSigDT = q0->_DbackSigDT; // DbackSigDT
      q1->_DgammaDt = 0.;
    
    
      static STensor3 Fpinv, Ce, Fepr;
      STensorOperation::inverseSTensor3(Fp1,Fpinv);
      STensorOperation::multSTensor3(F,Fpinv,Fepr);
      STensorOperation::multSTensor3FirstTranspose(Fepr,Fepr,Ce);
      
      // NEW
      static STensor3 Cepr, Ceinvpr;
      Cepr = Ce;
      STensorOperation::inverseSTensor3(Cepr,Ceinvpr);
      // NEW
    
      static STensor3 invFp0; // plastic predictor
      invFp0= Fpinv;
      STensor3& Fe = q1->_Fe;
      Fe = Fepr;
    
      static STensor43 DlnDCepr, DlnDCe;
      static STensor63 DDlnDDCe;
      static STensor3 expGN;
      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 Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
      double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
      ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); 
      getTVEdCorKirDT(q0,q1,T0,T);
      
      static STensor3 Me, Mepr, Kepr;  // Ke is corKir
      
      getModifiedMandel(Ce,q0,q1);
      Me = q1->_ModMandel;
      Mepr = Me;
      Kepr = q1-> _kirchhoff;
    
      static STensor3 devXn;
      static double trXn;
      STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
      
      static STensor3 PhiPr;
      PhiPr = q1->_ModMandel; 
      PhiPr -= q1->_backsig;
      
      static STensor3 devPhipr,devPhi; // effective dev stress predictor
      static double ptildepr,ptilde;
      STensorOperation::decomposeDevTr(PhiPr,devPhipr,ptildepr);
      ptildepr/= 3.;
    
      ptilde = ptildepr; // current effective pressure  -> first invariant
      devPhi = devPhipr;
    
      double PhiEqpr2 = 1.5*devPhipr.dotprod();
      double PhiEqpr = sqrt(PhiEqpr2);      // -> second invariant
    
      // 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& Gamma = q1->_Gamma;   // flow rule parameter
      double PhiEq = PhiEqpr;   // current effective deviator  
    
    
       // hardening
      this->hardening(q0,q1,T);
      static fullVector<double> a(3), Da(3); // yield coefficients and derivatives in respect to plastic deformation
      this->getYieldCoefficients(q1,a);
    
      double Hb(0.), dHbdT(0.), ddHbddT(0.);
      if (q1->_ipKinematic != NULL){
        Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
        dHbdT = Hb * q1->_ipKinematic->getDRDT();  // make sure to normalise DRDT while defining in python file //NEW
        ddHbddT = Hb * q1->_ipKinematic->getDDRDDT();
      }
    
      //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
      
      static fullVector<double> m(2);
      getChabocheCoeffs(m,0,q1);
      double Cx(0.), Px(0.);
      getPlasticPotential(ptilde, PhiEq, Px); 
      Cx = m(0)*(q0->_epspbarre + Dgamma) + m(1)*Px - 1/Hb * dHbdT + 1;  
      double Gt = Ge + pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
      double Kt = Ke + pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
    
      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.),DetaDT(0.);
            if (_viscosity != NULL)
              _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
            double etaOverDt = eta/this->getTimeStep();
    
            double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
            
            dDgammaDGamma = kk*(A+Gamma*dAdGamma);  // mistake in the paper (VD 2016)
    
            this->getYieldCoefficientDerivatives(q1,q1->_nup,Da);
            dfdDgamma = Da(2)*pow(PhiEq,_n) - Da(1)*ptilde -Da(0); 
            
            double dCxdDgamma = m(0);
            double Stemp = STensorOperation::doubledot(devPhi,devXn);
            double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; 
            double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma;
            dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
    
            if (Gamma>0 and etaOverDt>0)
              dfdDgamma -= _p*pow(etaOverDt,_p-1.)*Deta/this->getTimeStep()*pow(Gamma,_p);  // THIS term is absent in the paper!! WHAT??
    
            DfDGamma = dfdDgamma*dDgammaDGamma - (_n*a(2)*6.*Gt)*pow(PhiEq,_n)/u + a(1)*ptilde*2.*_b*Kt/v;
            if (Gamma>0 and etaOverDt>0)
              DfDGamma -=pow(etaOverDt,_p)*_p*pow(Gamma,_p-1.);
    
            double dGamma = -f/DfDGamma;    // WHAT 
    
            if (Gamma + dGamma <=0.){
                Gamma /= 2.;                // WHAT
            }
            else
              Gamma += dGamma;
    
            //Msg::Error("Gamma = %e",Gamma);
    
            u = 1.+6.*Gt*Gamma;
            v = 1.+2.*_b*Kt*Gamma; 
    
            devPhi = (devPhipr + devXn * (1-1/Cx));
            devPhi *= 1./u;  
            PhiEq = sqrt(1.5*devPhi.dotprod());
            ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/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,T);
            getYieldCoefficients(q1,a);
            //a.print("a update");
    
            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);
    
            if (fabs(f) <_tol) break;
    
            if(ite > maxite){
              Msg::Error("No convergence for plastic correction in mlawNonLinearTVP 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 + devXn * (1-1/Cx));
          devPhi *= 1./u;  
          ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v; 
    
          // update normal
          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); // Why again???? WHAT??
    
          // 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; // Here We have assumed that, Ce commutes with Q (i.e. with M) 
          
          // update A, B
          ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); //,false); 
          getTVEdCorKirDT(q0,q1,T0,T); // update dCorKirdT
            // DKDTsum and DGDTsum = sum of bulk/shear moduli derivatives wrt T 
    
          // backstress
          static STensor3 DB; // increment
          DB = (N); // increment
          DB *= (pow(kk,1)*Hb*Gamma); // pow(kk,2) CHANGE
          DB += q0->_backsig;
          DB *= 1/Cx;
          DB -= q0->_backsig;
          q1->_backsig += DB; // update
          
          q1->_ModMandel = devPhi;
          q1->_ModMandel += q1->_backsig;
    
          q1->_ModMandel(0,0) += (ptilde);
          q1->_ModMandel(1,1) += (ptilde);
          q1->_ModMandel(2,2) += (ptilde);
    
        }
        else{
          q1->getRefToDissipationActive() = false;
        }
      }
      
      // update Stress
      const STensor3& KS = q1->_kirchhoff;
      getModifiedMandel(Ce, q0, q1); // update Mandel 
      const STensor3& MS = q1->_ModMandel;
      static STensor3 S, Ceinv; // get 2nd PK
      STensorOperation::inverseSTensor3(Ce,Ceinv);
      STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S); 
      
      // first PK
      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);
        }
    
      // check Commutavity
      static STensor3 Check;
      checkCommutavity(Check,Ce,S,q1);
      
      // defo energy
      // q1->getRefToElasticEnergy()=deformationEnergy(*q1);
      q1->getRefToViscousEnergyPart()=viscousEnergy(*q0,*q1)+q0->getConstRefToViscousEnergyPart();
      q1->getRefToPlasticEnergy() = q0->plasticEnergy();
      
      if (Gamma > 0){
        double dotMSN = dot(MS,N);
        q1->getRefToPlasticEnergy() += Gamma*dotMSN; // = Dp:Me
      }
    
      
       // fluxT
       double KThConT, DKThConDT; mlawNonLinearTVE::getKTHCon(KThConT,T,&DKThConDT);
       double J  = 1.;
       STensor3 Finv(0.);
       if (_thermalEstimationPreviousConfig){                                            // ADD  _thermalEstimationPreviousConfig
        STensorOperation::inverseSTensor3(F0,Finv);
        J = STensorOperation::determinantSTensor3(F0);
       }
       else{
        STensorOperation::inverseSTensor3(F,Finv);
        J = STensorOperation::determinantSTensor3(F);
       }
       
       static STensor3 Cinv;
       STensorOperation::multSTensor3SecondTranspose(Finv,Finv,Cinv);
       STensorOperation::multSTensor3SVector3(Cinv,gradT,fluxT);
       fluxT *= (-KThConT*J);
       
    
      // ThermSrc and MechSrc after dPdF and dPdT - But need the following tensors for the mechSrcTVP function call
      static STensor3 DphiPDF;
      STensorOperation::zero(dFpdF);
      STensorOperation::zero(dFpdT);
      
                                          
      // I didnt make this
      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){  
          
          // dP/dF
          
          //1. get DtrKeprDCepr and DdevKeprDCepr, also DKeprDCepr
          static STensor3 DpDCepr;
          static STensor43 DdevKDCepr;
          static STensor43 DdevKeprDCepr, DdevMeprDCepr, DKeprDCepr;  
          static STensor3 DpKeprDCepr, DpMeprDCepr;      // K = corKir; pK -> pressure of corKir
          STensorOperation::multSTensor3STensor43(_I,DlnDCepr,DpKeprDCepr);
          DpKeprDCepr*= (0.5*Ke);
          STensorOperation::multSTensor43(_Idev,DlnDCepr,DdevKeprDCepr);
          DdevKeprDCepr*= Ge;
          
          // DKeprDCepr
          DKeprDCepr = DdevKeprDCepr;
          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++){
                  DKeprDCepr(i,j,k,l) += _I(i,j)*DpKeprDCepr(k,l);
                }
    
          // initiate here - corrected corKir derivatives wrt Cepr 
          DpDCepr = DpKeprDCepr;
          DdevKDCepr = DdevKeprDCepr;
        
          //2. get DCeprDCepr and DCeinvprDCepr
          static STensor43 DCeprDCepr, DCeinvprDCepr;
          DCeprDCepr = _I4;
    
          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++){
                  DCeinvprDCepr(i,s,k,l) = 0.;
                  for (int m=0; m<3; m++)
                    for (int j=0; j<3; j++)
                      DCeinvprDCepr(i,s,k,l) -= Ceinvpr(i,m)*_I4(m,j,k,l)*Ceinvpr(j,s);
                }
           
          //3. get DtrMeprDCepr and DdevMeprDCepr
          static STensor43 DMeprDCepr;
          
          static STensor3 temp1;
          static STensor43 temp2, temp3;
          STensorOperation::multSTensor3(Kepr,Ceinvpr,temp1);
          
          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++){
                  temp2(i,s,k,l) = 0.;
                  temp3(i,s,k,l) = 0.;
                  for (int m=0; m<3; m++){
                      temp2(i,s,k,l) += _I4(i,m,k,l)*temp1(m,s);
                      temp3(i,s,k,l) += DKeprDCepr(i,m,k,l)*Ceinvpr(m,s) + Kepr(i,m)*DCeinvprDCepr(m,s,k,l);
                  }
              }
          
          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++){
                  DMeprDCepr(i,s,k,l) = 0.;
                  for (int m=0; m<3; m++)
                    DMeprDCepr(i,s,k,l) += Cepr(i,m)*temp3(m,s,k,l);
                }
           
          DMeprDCepr += temp2;
          DMeprDCepr += DKeprDCepr;
          DMeprDCepr *= 0.5;
          
          // Because, TrMepr = TrKepr 
          DpMeprDCepr = DpKeprDCepr;   
          
          // 
          STensorOperation::zero(DdevMeprDCepr);
          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++){
                  DdevMeprDCepr(i,s,k,l) += DMeprDCepr(i,s,k,l) - _I(i,s)*DpMeprDCepr(k,l);
                }
          
          //4. get DdevphiDCepr and DphiPprDCepr
          static STensor3 DphiPprDCepr, DphiPDCepr;
          static STensor43 DdevphiprDCepr, DdevphiDCepr;
          
          STensorOperation::zero(DphiPprDCepr);
          STensorOperation::zero(DphiPDCepr);
          STensorOperation::zero(DdevphiprDCepr);
          STensorOperation::zero(DdevphiDCepr);
          
          DphiPprDCepr = DpMeprDCepr;
          DdevphiprDCepr = DdevMeprDCepr;
          
          DdevphiDCepr = DdevphiprDCepr;
          DdevphiDCepr *= 1./u;
          DphiPDCepr = DphiPprDCepr;
          DphiPDCepr *= 1./v;
          
         //5. get DphiEprDCepr from DdevphiDCepr
          static STensor3 DphiEDCepr, DphiEDdevPhi;
          
          STensorOperation::zero(DphiEDdevPhi);
          DphiEDdevPhi = devPhi;
          DphiEDdevPhi *= 1.5/(PhiEq);
          
          STensorOperation::multSTensor3STensor43(DphiEDdevPhi,DdevphiDCepr,DphiEDCepr);
         
         
         // 6. to 11. (inside if loop-> Gamma > 0.)
          static STensor3 dAdCepr, dfDCepr;
          static STensor3 DGDCepr;
          static STensor3 DgammaDCepr;
          static STensor3 DtrNDCepr;
          static STensor43 DdevNDCepr;
          static STensor43 dFpDCepr;
          static STensor43 dCedCepr, dCeinvdCepr;
          static STensor43 dXdCepr;
          static STensor3 dTrXdCepr;
          
          if (Gamma >0){
              // plastic
              
          //6. get dAdCepr from DdevphiDCepr 
            double fact = a(2)*_n*pow(PhiEq,_n-1.);
            for (int i=0; i<3; i++){
              for (int j=0; j<3; j++){
                dAdCepr(i,j) = (4.*_b*_b*ptilde/(A*3.))*DphiPDCepr(i,j);
                dfDCepr(i,j) = -a(1)*DphiPDCepr(i,j);
                dAdCepr(i,j) += (6./A)*PhiEq*DphiEDCepr(i,j);
                dfDCepr(i,j) += fact*DphiEDCepr(i,j);
              }
            }
            
          //7. get DGDCepr, DgammaDCepr
            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++){
                DgammaDCepr(i,j) = kk*Gamma*dAdCepr(i,j)+ kk*dDgammaDGamma*DGDCepr(i,j);
              }
            }
          
         //8. update DdevphiDCepr and DphiEDCepr
           static STensor3 dCxdCepr, dPxdCepr;
           STensorOperation::zero(dPxdCepr);
           STensorOperation::zero(dCxdCepr);
           dPxdCepr *= m(1);
           dCxdCepr = DgammaDCepr;
           dCxdCepr *= m(0);
           dCxdCepr += dPxdCepr;
           
           // DphiEDCepr
           static STensor3 temp7, temp8, temp9;
           STensorOperation::zero(temp7);
           STensorOperation::zero(temp8);
           STensorOperation::zero(temp9);
           
           temp7 = dCxdCepr;
           temp7 *= (1/3*1/(pow(Cx,2))*trXn/v);
           
           temp9 = dCxdCepr;
           temp9 *= (2/3*_b*Gamma*pow(kk,1)*Hb/pow(Cx,2)); // pow(kk,2) CHANGE
           
           temp8 = DGDCepr;
           temp8 *= (2*_b*(Ke+pow(kk,1)*Hb/(3*Cx))); // pow(kk,2) CHANGE
           temp8 -= temp9;
           temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
           temp8 *= 1/pow(v,2);
           
           DphiPDCepr += temp7;
           DphiPDCepr -= temp8;
           
           // DdevphiDCepr
           static STensor43 temp10, temp11;
           static STensor3 temp12, temp13, temp14;
           // STensorOperation::zero(temp10);
           // STensorOperation::zero(temp11);
           STensorOperation::zero(temp12);
           STensorOperation::zero(temp13);
           STensorOperation::zero(temp14);
           
           temp12 = devXn;
           temp12 *= 1/pow(Cx,2);
           
           STensorOperation::prod(temp12, dCxdCepr, 1., temp10);
           temp10 *= 1/u;
           
           for (int i=0; i<3; i++){
             for (int j=0; j<3; j++){
               temp13(i,j) = devPhipr(i,j) + (1-1/Cx)*devXn(i,j);
               temp14(i,j) = 6*DGDCepr(i,j)*( Ge + pow(kk,1)*Hb/(2*Cx) ) - 6*Gamma*(pow(kk,1)*Hb/(2*pow(Cx,2))*dCxdCepr(i,j)); // pow(kk,2) CHANGE
             }
           } 
           
           STensorOperation::prod(temp13, temp14, 1., temp11);
           temp11 *= 1/pow(u,2);
           
           DdevphiDCepr += temp10;
           DdevphiDCepr -= temp11;
           
           //9. get DtrNDCepr DdevNdCepr
           DtrNDCepr = DphiPDCepr;
           DtrNDCepr *= (2*_b);
           
           DdevNDCepr = DdevphiDCepr;
           DdevNDCepr *= 3;
           
           // 10. get dFpDCepr and dCeinvdCepr
           
           // dFpDCepr
           static STensor43 temp15;
           static STensor3 temp15_tr;
           for (int i=0; i<3; i++)
             for (int j=0; j<3; j++){
                temp15_tr(i,j) = Gamma/3.*DtrNDCepr(i,j);
               for (int k=0; k<3; k++)
                 for (int l=0; l<3; l++)
                   temp15(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,temp15,dFpDCepr);
           
           
           // dCeinvdCepr; Here, A = GN
           static STensor43 dexpAdCepr, dexpAinvdCepr;
           STensorOperation::multSTensor43(dexpAdA,temp15,dexpAdCepr);
           
           static STensor3 expGNinv;
           STensorOperation::inverseSTensor3(expGN,expGNinv);
           
           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++){
                   dexpAinvdCepr(i,s,k,l) = 0.;
                   for (int m=0; m<3; m++)
                     for (int j=0; j<3; j++)
                       dexpAinvdCepr(i,s,k,l) -= expGNinv(i,m)*dexpAdCepr(m,j,k,l)*expGNinv(j,s);
                 }  
           
           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++){
                   dCedCepr(i,s,k,l) = 0.;
                   for (int m=0; m<3; m++)
                     for (int j=0; j<3; j++)
                       dCedCepr(i,s,k,l) +=  expGNinv(m,i)*_I4(m,j,k,l)*expGNinv(j,s) - expGNinv(m,i)*Cepr(m,j)*dexpAinvdCepr(j,s,k,l)
                                              - dexpAinvdCepr(m,i,k,l)*Cepr(m,j)*expGNinv(j,s); // assuming GN is symmetric
                 }    // This derivative is wrong somehow. Use d()/dF below instead.
           
                 
           // 10.1 get dXdCepr
           const STensor3 backSig = q1->_backsig;
           static double trX;
           static STensor3 devX;
           STensorOperation::decomposeDevTr(backSig,devX,trX);
           
           STensorOperation::zero(dXdCepr);
           STensorOperation::zero(dTrXdCepr);
           for (int i=0; i<3; i++)
             for (int j=0; j<3; j++){
               dTrXdCepr(i,j) = pow(kk,1)*Hb*temp15_tr(i,j)/Cx - trX*dCxdCepr(i,j)/Cx; // pow(kk,2) CHANGE
               for (int k=0; k<3; k++)
                 for (int l=0; l<3; l++)
                   dXdCepr(i,j,k,l) = pow(kk,1)*Hb*temp15(i,j,k,l)/Cx -  q1->_backsig(i,j)*dCxdCepr(k,l)/Cx; // CHECK, MIGHT HAVE TO DIFFERENTIATE Hb wrt gamma
                   // pow(kk,2) CHANGE
           }
            
           // 11. update DpDCepr, DdevKDCepr
           // (DpKeprDCepr DdevKeprDCepr)
           // 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));
               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{
          // elastic
          STensorOperation::zero(DgammaDCepr);
          STensorOperation::zero(dFpDCepr);
          STensorOperation::zero(DtrNDCepr);
          STensorOperation::zero(DdevNDCepr);
          STensorOperation::zero(DGDCepr);
          STensorOperation::zero(dXdCepr);
          STensorOperation::zero(dTrXdCepr);
          
          // STensorOperation::zero(dCedCepr);
          dCedCepr = _I4;                                           // DOUBT _ CHECK
        }
        
        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++){
                dCeinvdCepr(i,s,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int j=0; j<3; j++)
                    dCeinvdCepr(i,s,k,l) -= Ceinv(i,m)*dCedCepr(m,j,k,l)*Ceinv(j,s);
                 }
        
        // 12. get dKcorDcepr
        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);
              }
            }
          }
        }
        
        //13. get  CeprToF conversion tensor
        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);
              }
            }
          }
        }
        
        static STensor43 DKcorDF;
        STensorOperation::multSTensor43(dKcorDcepr,CeprToF,DKcorDF);
        
        // 14. get dSdCepr and dSdF; also get dMedCepr and dMedF
        // Done below -> need dCedF
        
        STensor43& dXdF = q1->getRefToDbackStressdF();
        static STensor3 dTrXdF;
        
        STensorOperation::multSTensor43(dXdCepr,CeprToF,dXdF); 
    
        // DphiPDF - use this for the pressure term in R
        // for the pressure term in R
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
             DphiPDF(i,j) = 0.;
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                  DphiPDF(i,j) += DphiPDCepr(k,l)*CeprToF(k,l,i,j);
              }
          }
        
        // 15. get DgammaDF
        STensor3& DgammaDF = q1->_DgammaDF;
        STensor3& DGammaDF = q1->_DGammaDF;
        if (Gamma > 0){
          STensorOperation::multSTensor3STensor43(DgammaDCepr,CeprToF,DgammaDF);
          STensorOperation::multSTensor3STensor43(DGDCepr,CeprToF,DGammaDF);
          STensorOperation::multSTensor43(dFpDCepr,CeprToF,dFpdF);
        }
        else{
          STensorOperation::zero(DgammaDF);
          STensorOperation::zero(DGammaDF);
          // STensorOperation::zero(dFpdF);
        }
    
        // 16. Everything else and Tangent 
        static STensor43 DinvFpDF, dCedF, dCeinvDF; //
        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);
              }
        
        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++){
                dCedF(i,j,k,l) = 0.;
                for (int p=0; p<3; p++){
                    dCedF(i,j,k,l) += Fe(p,i)*dFedF(p,j,k,l) + dFedF(p,i,k,l)*Fe(p,j); // This Works!
                }
              }
        
        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++){
                dCeinvDF(i,s,k,l) = 0.;
                for (int m=0; m<3; m++)
                  for (int j=0; j<3; j++)
                    dCeinvDF(i,s,k,l) -= Ceinv(i,m)*dCedF(m,j,k,l)*Ceinv(j,s);
              }   
        
        // 17. Tangent - dPdF      
        static STensor43 dSdCepr, dSdF, dMedCepr;
        STensor43& dMedF = q1->getRefToDModMandelDF();
        
        /* Something wrong with dCedCepr - good luck debugging
         * 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++){
                dMedCepr(i,j,k,l) = dKcorDcepr(i,j,k,l); // CHANGED
                for (int m=0; m<3; m++)
                  for (int s=0; s<3; s++)
                    dMedCepr(i,j,k,l) += dCedCepr(i,m,k,l)*KS(m,s)*Ceinv(s,j) + Ce(i,m)*(dKcorDcepr(m,s,k,l)*Ceinv(s,j) + KS(m,s)*dCeinvdCepr(s,j,k,l));
                                        // CHECK!! the last terms
              }
        dMedCepr *= 0.5;
        STensorOperation::multSTensor43(dMedCepr,CeprToF,dMedF);
        
        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++){
                dSdCepr(i,j,k,l) = 0.;
                for (int m=0; m<3; m++){
                    // dSdCepr(i,j,k,l) += dKcorDcepr(i,m,k,l)*Ceinv(m,j) + KS(i,m)*dCeinvdCepr(m,j,k,l); 
                    // dSdCepr(i,j,k,l) += dCeinvdCepr(i,m,k,l)*KS(m,j) + Ceinv(i,m)*dKcorDcepr(m,j,k,l);
                    dSdCepr(i,j,k,l) += dCeinvdCepr(i,m,k,l)*MS(m,j) + Ceinv(i,m)*dMedCepr(m,j,k,l);
                                     // KS is corKir
                }}
        // dSdCepr *= 0.5;
        STensorOperation::multSTensor43(dSdCepr,CeprToF,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++){
                dMedF(i,j,k,l) = DKcorDF(i,j,k,l); // CHANGED
                for (int m=0; m<3; m++)
                  for (int s=0; s<3; s++)
                    dMedF(i,j,k,l) += dCedF(i,m,k,l)*KS(m,s)*Ceinv(s,j) + Ce(i,m)*(DKcorDF(m,s,k,l)*Ceinv(s,j) + KS(m,s)*dCeinvDF(s,j,k,l));
                                        // CHECK!! the last terms
              }
        dMedF *= 0.5;
             
        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++){
                    // dSdF(i,j,k,l) += DKcorDF(i,m,k,l)*Ceinv(m,j) + KS(i,m)*dCeinvDF(m,j,k,l); 
                    // dSdF(i,j,k,l) += dCeinvDF(i,m,k,l)*KS(m,j) + Ceinv(i,m)*DKcorDF(m,j,k,l);
                    dSdF(i,j,k,l) += dCeinvDF(i,m,k,l)*MS(m,j) + Ceinv(i,m)*dMedF(m,j,k,l);
                                     // KS is corKir
                }}
    
        // dSdF *= 0.5;
        
        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);
                  }
                }
              }
              
              
        // dP/dT
        
        // 1. get dMeprDT from dKeprDT  -> DcorKirprDT from TVE
        static STensor3 dMeprDT, dDevMeprDT, DcorKirDT;
        double dpMeprDT(0.);
        const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
        STensorOperation::zero(dMeprDT);
        
        static STensor3 temp16, temp17;
        STensorOperation::multSTensor3(dKeprDT, Ceinvpr, temp16);
        STensorOperation::multSTensor3(Cepr, temp16, temp17);
    
        dMeprDT = 0.5*(dKeprDT + temp17);
        
        STensorOperation::decomposeDevTr(dMeprDT,dDevMeprDT,dpMeprDT);
        // dpMeprDT *= 1/3; WHAT IS THE PROBLEM HERE? WHY DOES IT GIVE ZERO?
        dpMeprDT = dMeprDT.trace()/3;
        
        DcorKirDT = dKeprDT; // update later
        
        // 2. get dPhipprDT and dDevPhiprDT
        static double dPhipprDT;
        static STensor3 dDevPhiprDT;
        
        dPhipprDT = dpMeprDT;
        dDevPhiprDT = dDevMeprDT;
        
        // 3. get dCxdT, dGxdT, dKxdT, dXdT
        static double dCxdT, dGxdT, dKxdT;
        
        dCxdT = 1/(pow(Hb,2))*pow(dHbdT,2)*(T-T0) - 1/Hb*ddHbddT*(T-T0) - 1/Hb*dHbdT; 
        dGxdT = DGDTsum + pow(kk,1)/(2*Cx)*dHbdT - (pow(kk,1)*Hb)/(2*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
        dKxdT = DKDTsum + pow(kk,1)/(3*Cx)*dHbdT - (pow(kk,1)*Hb)/(3*pow(Cx,2))*dCxdT; // pow(kk,2) CHANGE
        
        // need this for DmechsourceDT -> initialise here (X is backStress)
        STensor3& dXdT = q1->getRefToDbackStressdT();  
        STensorOperation::zero(dXdT);
        dXdT -= (q1->_backsig)*(dCxdT);
        dXdT *= 1/Cx;
        
        // 4. get dPhiPdT and dDevPhiDT
        static double dPhiPdT;
        static STensor3 dDevPhiDT;
        STensorOperation::zero(dDevPhiDT);
        
        static double temp18;
        static STensor3 temp19;
         
        dPhiPdT = dPhipprDT/v + 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this
        temp18 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*Gamma*dKxdT)/pow(v,2); // Gamma derivatives will change this
        dPhiPdT -= temp18;
        
        temp19 = (devPhipr + devXn*(1-1/Cx))*(6*Gamma*dGxdT);
        temp19 *= 1/pow(u,2);
        for (int i=0; i<3; i++){
          for (int j=0; j<3; j++){
            dDevPhiDT(i,j) = (dDevPhiprDT(i,j) + devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);
          }
        }
       
        // 5. get dPhiEdT
        static double dPhiEdT;
        STensorOperation::doubleContractionSTensor3(DphiEDdevPhi,dDevPhiDT,dPhiEdT);
        
        // 6. to 11. (inside if loop-> Gamma > 0.)
        static double dAdT, dfdT;
        static double dGammaDT;
        static double DgammaDT;
        static double DtrNdT;
        static STensor3 DdevNdT;
        
        double eta(0.),Deta(0.),DetaDT(0.);
        if (_viscosity != NULL)
            _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
        double etaOverDt = eta/this->getTimeStep();    
    
        if (Gamma >0){
            // plastic
            
            // 6. get DaDT, dAdT and dfdT
            fullVector<double> DaDT(3), DDaDTT(3);
            getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
            
            double a1, a2, a3;
            double a4, a5, a6;
            a1 = a(0); a2 = a(1); a3 = a(2);
            a4 = DaDT(0); a5 = DaDT(1); a6 = DaDT(2);
            
            dAdT = (4.*_b*_b*ptilde/(A*3.))*dPhiPdT + (6./A)*PhiEq*dPhiEdT;
            dfdT = DaDT(2)*pow(PhiEq,_n) + a(2)*_n*pow(PhiEq,_n-1.)*dPhiEdT - DaDT(1)*ptilde - a(1)*dPhiPdT - DaDT(0);
            if (this->getTimeStep()>0){
                dfdT -= (DetaDT*Gamma/this->getTimeStep())*_p*pow((eta*Gamma/this->getTimeStep()),(_p-1));
            }
            
            // 7. get dGammaDT, DgammaDT
            dGammaDT = (-dfdT-dfdDgamma*kk*Gamma*dAdT)/DfDGamma;
            DgammaDT = kk*Gamma*dAdT + kk*dDgammaDGamma*dGammaDT;
            
            // 8. update DdevphiDT and DphiPDT
            static double dPxdT, temp18_2;
            static STensor3 temp19_2;
            
            dPxdT = 0.*m(1);
            dCxdT += m(0)*DgammaDT;
            dCxdT += dPxdT; 
            
            dPhiPdT += temp18;
            temp18_2 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*(dGammaDT*(Ke+pow(kk,1)*Hb/(3*Cx)) + Gamma*dKxdT))/pow(v,2); // pow(kk,2) CHANGE
            dPhiPdT -= temp18_2;
            
            dDevPhiDT += temp19;
            temp19_2 = (devPhipr+devXn*(1-1/Cx))*(6*(dGammaDT*(Ge+pow(kk,1)*Hb/(2*Cx)) + Gamma*dGxdT)); // pow(kk,2) CHANGE
            temp19_2 *= 1/pow(u,2);
            dDevPhiDT -= temp19_2;
            
            // 9. get DtrNdT DdevNdT
            DtrNdT = dPhiPdT;
            DtrNdT *= (2*_b);
    
            DdevNdT = dDevPhiDT;
            DdevNdT *= 3;
            
            // 10. get dFpdT
            static STensor3 temp20;
            STensorOperation::zero(temp20);
            for (int i=0; i<3; i++)
              for (int j=0; j<3; j++)
                temp20(i,j) += N(i,j)*dGammaDT + Gamma*DdevNdT(i,j) + 1/3*Gamma*_I(i,j)*DtrNdT;
                
            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);
                    }
                  }
            
            for (int i=0; i<3; i++){
              for (int j=0; j<3; j++){
                dFpdT(i,j) = 0.;
                for (int k=0; k<3; k++){
                  for (int l=0; l<3; l++){
                      dFpdT(i,j) += EprFp0(i,j,k,l)*temp20(k,l);
                  }
                }
              }
            }
            
            // 10.1 dXdT - backstress temperature derivative -> correction
            
            for (int i=0; i<3; i++){
              for (int j=0; j<3; j++){
                dXdT(i,j) += 1/Cx*((pow(kk,1)*(dHbdT*Gamma*N(i,j) + Hb*temp20(i,j)))); // pow(kk,2) CHANGE
              }
            }
            
            // 11. DcorKirDT
            for (int i=0; i<3; i++){
              for (int j=0; j<3; j++){
                DcorKirDT(i,j) -= ( DKDTsum*Gamma*trN + Ke*(dGammaDT*trN+Gamma*DtrNdT))*_I(i,j);
                DcorKirDT(i,j) -=  2.*( DGDTsum*Gamma*devN(i,j) + Ge*(dGammaDT*devN(i,j)+Gamma*DdevNdT(i,j)));
              }
            }
            
        }
        else{
            // elastic
            STensorOperation::zero(DgammaDT);
            // STensorOperation::zero(dFpdT);
            STensorOperation::zero(DtrNdT);
            STensorOperation::zero(DdevNdT);
            STensorOperation::zero(dGammaDT);
            }
            
        // 11.1
        q1->_DGammaDT = dGammaDT;
        q1->_DgammaDT = DgammaDT;
        
        // 12. get dKcorDT
        // done above 
        
        // 13. get DinvFpdT, dFedT, dCedT, dCeinvDT
        static STensor3 DinvFpdT, dFedT, dCedT, dCeinvdT;
        
        STensorOperation::zero(DinvFpdT);
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int p=0; p<3; p++)
              for (int q=0; q<3; q++){
                DinvFpdT(i,j) -= Fpinv(i,p)*dFpdT(p,q)*Fpinv(q,j);
              }
        
        STensorOperation::zero(dFedT);
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++)
            for (int k=0; k<3; k++){
              dFedT(i,j) += F(i,k)*DinvFpdT(k,j);
            }
        
        STensorOperation::zero(dCedT);
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
              for (int p=0; p<3; p++){
                dCedT(i,j) += Fe(p,i)*dFedT(p,j) + dFedT(p,i)*Fe(p,j); // This Works!
              }
            }
    
        STensorOperation::zero(dCeinvdT);
        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++){
                dCeinvdT(i,s) -= Ceinv(i,k)*dCedT(k,l)*Ceinv(l,s);
              }
        
        // 13.1 get dMeDT -> needed for DmechSourceDT
        STensor3& dMeDT = q1->getRefToDModMandelDT();
        
        // STensorOperation::zero(dMeDT);
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dMeDT(i,j) = DcorKirDT(i,j);
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                dMeDT(i,j) += dCedT(i,k)*KS(k,l)*Ceinv(l,j) + Ce(i,k)*(DcorKirDT(k,l)*Ceinv(l,j) + KS(k,l)*dCeinvdT(l,j));
              }
          }
        dMeDT *= 0.5; 
        
        // 14. get dSdT
        static STensor3 dSdT;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dSdT(i,j) = 0.;
            for (int m=0; m<3; m++){
                  // dSdT(i,j) += DcorKirDT(i,m)*Ceinv(m,j) + KS(i,m)*dCeinvdT(m,j) + dCeinvdT(i,m)*KS(m,j) + Ceinv(i,m)*DcorKirDT(m,j);
                  dSdT(i,j) += Ceinv(i,m)*dMeDT(m,j) + dCeinvdT(i,m)*MS(m,j);
            }
          }
        // dSdT *= 0.5;
        
        // 15. Finally, get dPdT
        // static STensor3 dPdT;
        for (int i=0; i<3; i++)
          for (int j=0; j<3; j++){
            dPdT(i,j) = 0.;  
            for (int k=0; k<3; k++)
              for (int l=0; l<3; l++){
                dPdT(i,j) += (dFedT(i,k)*S(k,l)*Fpinv(j,l) + Fe(i,k)*dSdT(k,l)*Fpinv(j,l) + Fe(i,k)*S(k,l)*DinvFpdT(j,l));
              }
          }
      }
      
      if(stiff){
      
      // TVP - flux, thermalEnergy, thermalSource
        if (!_thermalEstimationPreviousConfig){
          
          static STensor3 DJDF;
          static STensor43 DCinvDF;
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              DJDF(i,j) = J*Finv(j,i);
              for (int k=0; k<3; k++){
                for (int l=0; l<3; l++){
                  DCinvDF(i,j,k,l) = 0.;
                  for (int p=0; p<3; p++){
                    for (int a=0; a<3; a++){
                      for (int b=0; b<3; b++){
                        DCinvDF(i,j,k,l) -= 2.*F(k,p)*Cinv(i,a)*Cinv(j,b)*_I4(a,b,p,l);
                      }
                    }
                  }
                }
              }
             }
           }
    
          for (int i=0; i<3; i++){
            for (int j=0; j<3; j++){
              for (int k=0; k<3; k++){
                dfluxTdF(i,j,k) += (DJDF(j,k)*fluxT(i)/J);
                for (int l=0; l<3; l++){
                  dfluxTdF(i,j,k) -= (J*DCinvDF(i,l,j,k)*gradT(l)*KThConT);
                }
              }
            }
          }
        }
      }
    
      // thermSrc and MechSrc
      
       // thermalEnergy 
       double CpT; 
       if (stiff){
           double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
           getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
           // CpT = -T*DDpsiTVMdTT;
           getCp(CpT,T);
        }
       else{
           getCp(CpT,T);
        }
     
       q1->_thermalEnergy = CpT*T;  
        
       // thermalSource
       if (this->getTimeStep() > 0.){
        thermalSource = -CpT*(T-T0)/this->getTimeStep();
        }
       else
        thermalSource = 0.;
    
      // mechSourceTVP
      double& Wm = q1->getRefToMechSrcTVP();
      double& Wm_TVE = q1->getRefToMechSrcTVE();
      
      mechanicalSource = 0.;
      mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE); // mechSourceTVE
      mechanicalSource += Wm_TVE;
    
      getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,&Wm);
      mechanicalSource += Wm;
        
      // freeEnergy and elastic energy
      double& psiTVM = q1->getRefToFreeEnergyTVM();
      getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
      
      q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);     // deformationEnergy(*q1,T);
      
    if (stiff){
        
        const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT();
        double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
        getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
        // double CpT = -T*DDpsiTVMdTT;
        // double CpT_0 = -T0*DDpsiTVMdTT_0;
    
        // thermal source derivatives
        double DCpDT(0.);
        mlawNonLinearTVE::getCp(CpT,T,&DCpDT);
        static STensor3 DCpDF;
        STensorOperation::zero(DCpDF);                                                       // CHANGE for DCpDF
        
        if (this->getTimeStep() > 0){
         dthermalSourcedT = -DCpDT*(T-T0)/this->getTimeStep() - CpT/this->getTimeStep();
         // dthermalSourcedT = -CpT/this->getTimeStep() - (CpT-CpT_0)/this->getTimeStep(); // thermalSource = -CpT*(T-T0)/this->getTimeStep();
         for(int i = 0; i< 3; i++){
           for(int j = 0; j< 3; j++){
             dthermalSourcedF(i,j) = -DCpDF(i,j)*(T-T0)/this->getTimeStep();
           }
         }
        }
        else{
         dthermalSourcedT = 0.;
         STensorOperation::zero(dthermalSourcedF);
        }
        
     
        // mechSourceTVP derivatives
        
        // conversion tensor for mechSourceTVE Terms
        static STensor43 DCeDFe, DEeDFe;
        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++){
                DCeDFe(i,j,k,l) = 0.;
                for (int p=0; p<3; p++){
                  DCeDFe(i,j,k,l) += 2*F(p,i)*dFedF(p,j,k,l);
                }
              }
              
        STensorOperation::multSTensor43(DlnDCe,DCeDFe,DEeDFe);
        DEeDFe *= 0.5;
        
        double& dWmdT_TVE = q1->getRefTodMechSrcTVEdT();
        STensor3& dWmdF_TVE = q1->getRefTodMechSrcTVEdF();
        mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,DEeDFe,NULL,&dWmdF_TVE,&dWmdT_TVE);
        
        dmechanicalSourcedT = 0.;
        STensorOperation::zero(dmechanicalSourceF);
        dmechanicalSourceF += dWmdF_TVE;
        dmechanicalSourcedT += dWmdT_TVE;
        
        double& dWmdT = q1->getRefTodMechSrcTVPdT();
        STensor3& dWmdF = q1->getRefTodMechSrcTVPdF();
    
        getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT); 
        
        dmechanicalSourceF += dWmdF;
        dmechanicalSourcedT += dWmdT;
    }
    };
      
    
    void mlawNonLinearTVP::predictorCorrector_TVP_AssociatedFlow(const STensor3& F, const IPNonLinearTVP *q0, IPNonLinearTVP *q1,
                                STensor3&P, const bool stiff, STensor43& Tangent, STensor43& dFedF, STensor43& dFpdF, const double T0, const double T) 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 Ke, Ge = 0.; double Kde, Gde = 0.; double KTsum, GTsum = 0.;
      double DKe, DGe = 0.; double DKde, DGde = 0.; double DKDTsum, DGDTsum = 0.;
      ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T); 
    
      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,T);
      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: mlawNonLinearTVP::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,T);
            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;
    
          // 
          ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,Kde,Gde,DKe,DGe,DKde,DGde,KTsum,GTsum,DKDTsum,DGDTsum,T0,T,false); 
        }
        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){
         // FILL THIS
          }
    };