diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 937e8e80d67c133732cd1555ed149c4a759bc17e..3206b2ac6331291a40ab66f8d8ca75bcd8328af7 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -23,7 +23,8 @@ mlawNonLinearTVP::mlawNonLinearTVP(const int num,const double E,const double nu,
       // 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)
@@ -40,6 +41,14 @@ mlawNonLinearTVP::mlawNonLinearTVP(const mlawNonLinearTVP& src): mlawNonLinearTV
   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){
@@ -54,9 +63,14 @@ mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
         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){
+        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;
 };
@@ -64,6 +78,8 @@ mlawNonLinearTVP& mlawNonLinearTVP::operator=(const materialLaw& source){
 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??
@@ -108,20 +124,20 @@ void mlawNonLinearTVP::getModifiedMandel(const STensor3& C, const IPNonLinearTVP
     M = 0.5*(corKir + temp2);
 };
 
-void mlawNonLinearTVP::checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const{
+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, temp5, temp6;
-    STensorOperation::multSTensor3(C, S, temp1);
-    STensorOperation::multSTensor3(S, C, temp2);
-    STensorOperation::multSTensor3(corKir, C, temp3);
-    STensorOperation::multSTensor3(C, corKir, temp4);
+    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 
-    temp5 = temp1;
-    temp5 -= temp2;
+    commuteCheck = temp1;
+    commuteCheck -= temp2;
     
     // Check corKir (just for fun)
     temp6 = temp3;
@@ -130,14 +146,13 @@ void mlawNonLinearTVP::checkCommutavity(const STensor3& C, const STensor3& S, co
     // 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) = 0; // 1; 
+    coeffs(0) = 0; // 1; //DEBUGGING
     coeffs(1) = 0;
     
     // add dependency on p later
@@ -170,13 +185,16 @@ void mlawNonLinearTVP::getYieldCoefficientTempDers(const IPNonLinearTVP *q1, con
   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) = DDa1Dmm*DmDT*DmDT + Da1Dm*DDmDTT + 3.*2.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc*sigc)*dsigcdT*dsigcdT - 3.*(pow(m,_n)-1.)/(m+1.)/(sigc*sigc)*ddsigcdTT;
+  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;
     
@@ -547,7 +565,7 @@ void mlawNonLinearTVP::getMechSourceTVP(const IPNonLinearTVP *q0, IPNonLinearTVP
     const std::vector<double>& dIsoHardForcedT =  q1->_dIsoHardForcedT;
     const std::vector<STensor3>& dIsoHardForcedF =  q1->_dIsoHardForcedF;
     
-    double R, dRdT, ddRdTT;
+    double R(0.), dRdT(0.), ddRdTT(0.);
     STensor3 dRdF, ddRdFdT; STensorOperation::zero(dRdF); STensorOperation::zero(ddRdFdT);
     
     for (int i=0; i<3; i++){
@@ -1004,10 +1022,11 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
                                   double &dmechanicalSourcedT,
                                   STensor3 &dmechanicalSourceF) const{
                                       
-  // compute elastic predictor 
+  // Compute elastic predictor first
   STensor3& Fp1 = q1->_Fp;
   const STensor3& Fp0 = q0->_Fp;
-
+  
+  // Initialise predictor values
   Fp1 = Fp0; // plastic deformation tensor
   q1->_epspbarre = q0->_epspbarre; // plastic equivalent strain
   q1->_epspCompression = q0->_epspCompression;
@@ -1017,23 +1036,38 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   q1->_DbackSigDT = q0->_DbackSigDT; // DbackSigDT
   q1->_DgammaDt = 0.;
 
+  // Get hardening stuff 
+  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.), dHb(0.), dHbdT(0.), ddHbddT(0.);
+  if (q1->_ipKinematic != NULL){
+    Hb = q1->_ipKinematic->getDR(); // kinematic hardening parameter
+    dHb = q1->_ipKinematic->getDDR(); // kinematic hardening parameter derivative (dHb/dgamma ??)
+    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");
+  
+  
+  // Get Cepr, Ceinvpr
   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;
 
+  // Predictor - TVE
   static STensor43 DlnDCepr, DlnDCe;
   static STensor63 DDlnDDCe;
   static STensor3 expGN;
@@ -1045,95 +1079,89 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   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); 
   getTVEdCorKirDT(q0,q1,T0,T);
   
-  // NEW
+  // Get Mepr
   static STensor3 Me, Mepr, Kepr;  // Ke is corKir
-  
   getModifiedMandel(Ce,q0,q1);
   Me = q1->_ModMandel;
   Mepr = Me;
   Kepr = q1-> _kirchhoff;
-  // const STensor3& dKeprDT = q1->getConstRefToDcorKirDT();
-  // NEW
 
-  // NEW
+  // Get Xn
   static STensor3 devXn;
   static double trXn;
   STensorOperation::decomposeDevTr(q0->_backsig,devXn,trXn); // needed for chaboche model
-  // NEW
   
+  // Initialise Phipr
   static STensor3 PhiPr;
-  PhiPr = q1->_ModMandel; // _kirchhoff; //NEW
+  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");
-
+  
+  // Initialise Normal
   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));
+  
+  // Get plastic poisson ratio
+  q1->_nup = (9.-2.*_b)/(18.+2.*_b);
+  double kk = 1./sqrt(1.+2.*q1->_nup*q1->_nup);
+  
+  // Initialise variables
+  // double Gamma = 0.;  
+  double& Gamma = q1->_Gamma;   // flow rule parameter 
+  double dDgammaDGamma = 0.;
+  double Dgamma = 0.; // eqplastic strain
+  static fullVector<double> m(2);
+  double Px(0.);
+  getChabocheCoeffs(m,0,q1);
+  getPlasticPotential(ptildepr, PhiEqpr, Px);  // CHANGE -> this is wrong, change to backstress instead of phi 
 
   double DfDGamma = 0.;
   double dfdDgamma = 0.;  
   double u = 1.;
   double v = 1.;
+  
+  // Initialise Cx, Gt, Kt
+  double Cx = m(0)*Dgamma + m(1)*Px - 1/Hb * dHbdT * (T-T0) + 1;   // CHECK 
+  double Gt = Ge;
+  double Kt = Ke;
+  
+  // Initialise ptilde and devPhi (phi)
+  ptilde = ptildepr + 1/3*trXn; // current effective pressure  -> first invariant
+  ptilde *= 1/v;
+  devPhi = devPhipr + devXn;
+  devPhi *= 1/u;
+  
+  if (abs(Cx) > 0.0){
+    Gt += pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
+    Kt += pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
+    
+    ptilde -= 1/3*trXn*1/Cx;
+    devPhi -= devXn;
+    devPhi *= 1/Cx;
+  }
+  
+  double PhiEqpr2 = 1.5*devPhi.dotprod();
+  double PhiEqpr = sqrt(PhiEqpr2);      // -> second invariant
+  double PhiEq = PhiEqpr;   // current effective deviator 
  
+  // Get f and A
+  double f = a(2)*pow(PhiEq,_n) - (a(1)*ptilde+a(0)); 
   double A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
 
-  double dDgammaDGamma = 0.;
-  double Dgamma = 0.; // eqplastic strain
-  
-// NEW
-  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
-// NEW  
-
 
+  // NR Loop
   if (q1->dissipationIsBlocked()){
     q1->getRefToDissipationActive() = false;
   }
@@ -1144,60 +1172,94 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
       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.);
+          
+        // Viscosity  
+        double eta(0.), Deta(0.), DetaDT(0.);
         if (_viscosity != NULL)
           _viscosity->get(q1->_epspbarre,T,eta,Deta,DetaDT);
         double etaOverDt = eta/this->getTimeStep();
-
-    // NEW
-        double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A);
-        // double dAdGamma = -(72.*Gt*PhiEq*PhiEq/u+ 16.*Kt*_b*_b*_b*ptilde*ptilde/(3.*v))/(2.*A); // -> the same
-    // NEW
         
+        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); 
-    // NEW
-        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; // CHANGE THIS DERIVATIVE -> tensor subtracting scalar
-        double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma;
+        
+        // Get He, Hp
+        double dCxdDgamma = m(0); // CHANGE to add derivatives of m(0) + dm0Ddgamma*Dgamma
+        
+        // double Stemp = STensorOperation::doubledot(devPhi,devXn);
+        // double He = (PhiEq*6*Gamma*kk*kk*Hb - Stemp*3.*1./PhiEq)/(2*Cx*Cx*u) * dCxdDgamma; // CHANGE THIS DERIVATIVE -> tensor subtracting scalar //CHECK
+        // double Hp = ptilde*(2*_b*Gamma*kk*kk*Hb - trXn)/(3*Cx*Cx*v) * dCxdDgamma; // CHECK
+        
+        static STensor3 DdevPhidDgamma, term2; 
+        STensorOperation::zero(DdevPhidDgamma);
+        if (abs(Cx) > 0.0){
+            // 1st term of He
+            DdevPhidDgamma = 3*Gamma*pow(kk,1)*Hb*devPhipr;  // pow(kk,2) CHANGE
+            DdevPhidDgamma -= devXn*(1 + 6*Gamma*(Ge + pow(kk,1)*Hb/2)); // pow(kk,2) CHANGE
+            DdevPhidDgamma *= 1/(Cx*Cx*u*u) * dCxdDgamma;
+            
+            // 2nd term of He
+            term2 = 3*Gamma*pow(kk,1)*Cx*devPhi; // pow(kk,2) CHANGE
+            term2 *= 1/(Cx*Cx*u) * dHb;
+            
+            DdevPhidDgamma -= term2;
+        }
+        double Stemp = STensorOperation::doubledot(devPhi,DdevPhidDgamma);
+        double He = 1.5*Stemp/PhiEq;
+        
+        double Hp = ( 2*_b*Gamma*pow(kk,1)*Hb*ptildepr - trXn*(1 + 2*_b*Gamma*(Ke + pow(kk,1)*Hb/3)) )*1/(3*Cx*Cx*v*v)*dCxdDgamma;
+        Hp -= (2*_b*Gamma*pow(kk,1)*Cx*ptilde)*1/(3*Cx*Cx*v) * dHb;
+        
         dfdDgamma += a(2)*_n*pow(PhiEq,(_n-1))*He - a(1)*Hp;
-    // NEW
+        // 
+        
         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.);
+          DfDGamma -= pow(etaOverDt,_p)*_p*pow(Gamma,(_p-1.);
 
-        double dGamma = -f/DfDGamma; // WHAT 
+        double dGamma = -f/DfDGamma;    // NR 
 
         if (Gamma + dGamma <=0.){
-            Gamma /= 2.;                // WHAT
+            Gamma /= 2.;                // NR
         }
         else
           Gamma += dGamma;
 
         //Msg::Error("Gamma = %e",Gamma);
-
+        
+        // Correct Phi, A, delgamma
+        Cx = m(0)*Dgamma + m(1)*Px - 1/Hb * dHbdT * (T-T0) + 1;
         u = 1.+6.*Gt*Gamma;
         v = 1.+2.*_b*Kt*Gamma; 
-    // NEW       
-        // PhiEq = PhiEqpr/u; // -> original
-        // ptilde = ptildepr/v; // -> original
+        
+        ptilde = ptildepr + 1/3*trXn; // current effective pressure  -> first invariant
+        ptilde *= 1/v;
+        devPhi = devPhipr + devXn;
+        devPhi *= 1/u;
+        
         devPhi = (devPhipr + devXn * (1-1/Cx));
         devPhi *= 1./u;  
         PhiEq = sqrt(1.5*devPhi.dotprod());
-        ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;
-     // NEW   
+        ptilde = (ptildepr + 1/3*trXn*(1-1/Cx))/v;  
+        
+        if (abs(Cx) > 0.0){
+            Gt += pow(kk,1) * Hb/(2*Cx); // pow(kk,2) CHANGE
+            Kt += pow(kk,1) * Hb/(3*Cx); // pow(kk,2) CHANGE
+    
+            ptilde -= 1/3*trXn*1/Cx;
+            devPhi -= devXn;
+            devPhi *= 1/Cx;
+        }
         
         A = sqrt(6.*PhiEq*PhiEq+4.*_b*_b/3.*ptilde*ptilde);
         Dgamma = kk*Gamma*A;
@@ -1312,15 +1374,16 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   // STensorOperation::multSTensor3STensor43(KS,DlnDCe,S);
   
   // NEW 
-  getModifiedMandel(Ce, q0, q1); // update Mandel // this step is probably wrong -> Mandel is already being updated above in the NR loop. Should corKir come from Mandel then?
-                                    // WHAT ???-
+  getModifiedMandel(Ce, q0, q1); // update Mandel
+  
   const STensor3& MS = q1->_ModMandel;
   // second Piola Kirchhoff stress
   static STensor3 S, Ceinv;
   STensorOperation::inverseSTensor3(Ce,Ceinv);
   STensorOperation::multSTensor3(Ceinv,q1->_ModMandel,S);
   // NEW
-
+   
+  // first PK 
   for(int i=0; i<3; i++)
     for(int j=0; j<3; j++){
       P(i,j) = 0.;
@@ -1362,55 +1425,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
    STensorOperation::multSTensor3SVector3(Cinv,gradT,fluxT);
    fluxT *= (-KThConT*J);
 
-   // thermalEnergy -> assuming stiff = false  
-   double CpT; 
-   if (stiff){
-       double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
-       getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
-       CpT = -T*DDpsiTVMdTT;
-    }
-   else{
-       getCp(CpT,T);
-    }
- 
-   q1->_thermalEnergy = CpT*T;  
-    
-   // thermalSource
-   if (this->getTimeStep() > 0.){
-    thermalSource = -CpT*(T-T0)/this->getTimeStep();
-    }
-   else
-    thermalSource = 0.;
-
-  // NEW - mechSourceTVP
-  double& Wm = q1->getRefToMechSrcTVP();
-  double& Wm_TVE = q1->getRefToMechSrcTVE();
-  
-  // Have to do this to update DcorKirDT in IP, and mechSourceTVE and derivatives
-  mechanicalSource = 0.;
-  mlawNonLinearTVE::getMechSourceTVE(q0,q1,T0,T,DKDTsum,DGDTsum,_I4,&Wm_TVE);
-  // STensor3 temp101 = q1->_DcorKirDT;
-  mechanicalSource += Wm_TVE;
- 
-  STensor3 P_TVE(0.);
-  SVector3 fluxT_TVE(0.);
-  double thermalSource_TVE(0.);
-  double mechanicalSource_TVE(0.);
- /* mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(q0->getConstRefToFe(), q1->getConstRefToFe(), P_TVE, q0, q1, Tangent, dFedF, dFedT,
-                                                      T0, T, gradT0, gradT, fluxT_TVE, dPdT, dfluxTdgradT, dfluxTdT, dfluxTdF, 
-                                                      thermalSource_TVE, dthermalSourcedT, dthermalSourcedF, 
-                                                      mechanicalSource, dmechanicalSourcedT, dmechanicalSourceF, 
-                                                      stiff);*/
-
-
-  getMechSourceTVP(q0,q1,T0,T,_I,_I4,_I,&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);
+   // 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){
@@ -1601,7 +1619,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
        temp8 *= (ptildepr + 1/3*trXn*(1-1/Cx));
        temp8 *= 1/pow(v,2);
        
-       DphiPDCepr -= temp7;
+       DphiPDCepr -= temp7;  // CHECK sign
        DphiPDCepr -= temp8;
        
        // DdevphiDCepr
@@ -1613,15 +1631,15 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
        STensorOperation::zero(temp13);
        STensorOperation::zero(temp14);
        
-       temp13 = devXn;
-       temp13 *= 1/pow(Cx,2);
+       temp13 = devXn;       // This must be temp12?
+       temp13 *= 1/pow(Cx,2);  // This must be temp12?
        
        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);
+           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
          }
        } 
@@ -1629,7 +1647,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
        STensorOperation::prod(temp13, temp14, 1., temp11);
        temp11 *= 1/pow(u,2);
        
-       DdevphiDCepr -= temp10;
+       DdevphiDCepr -= temp10;  // CHECK sign
        DdevphiDCepr -= temp11;
        
        //9. get DtrNDCepr DdevNdCepr
@@ -1794,7 +1812,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
          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);
+              DphiPDF(i,j) += DphiPDCepr(k,l)*CeprToF(k,l,i,j);
           }
       }
     
@@ -1808,7 +1826,8 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
     }
     else{
       STensorOperation::zero(DgammaDF);
-      STensorOperation::zero(dFpdF);
+      STensorOperation::zero(DGammaDF);
+      // STensorOperation::zero(dFpdF);
     }
 
     // 16. Everything else and Tangent 
@@ -1973,7 +1992,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
     static double temp18;
     static STensor3 temp19;
      
-    dPhiPdT = dPhipprDT/v - 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this
+    dPhiPdT = dPhipprDT/v - 1/3*trXn/pow(Cx,2)*dCxdT/v; // no correction to this  // CHECK the sign of the 2nd term
     temp18 = (ptildepr+1/3*trXn*(1-1/Cx))*(2*_b*Gamma*dKxdT)/pow(v,2); // Gamma derivatives will change this
     dPhiPdT -= temp18;
     
@@ -1981,7 +2000,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
     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);
+        dDevPhiDT(i,j) = (dDevPhiprDT(i,j) - devXn(i,j)/pow(Cx,2)*dCxdT)/u - temp19(i,j);  // CHECK the sign of the 2nd term
       }
     }
    
@@ -2010,7 +2029,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
         getYieldCoefficientTempDers(q1,T,DaDT,DDaDTT);
         
         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);
+        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));
         }
@@ -2043,12 +2062,12 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
         DdevNdT = dDevPhiDT;
         DdevNdT *= 3;
         
-        // 10. get dFpDT
+        // 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) + Gamma/3.*_I(i,j)*DtrNdT;
+            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++)
@@ -2063,10 +2082,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
         
         for (int i=0; i<3; i++){
           for (int j=0; j<3; j++){
-            dFpDT(i,j) = 0.;
+            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);
+                  dFpdT(i,j) += EprFp0(i,j,k,l)*temp20(k,l);
               }
             }
           }
@@ -2092,12 +2111,16 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
     else{
         // elastic
         STensorOperation::zero(DgammaDT);
-        STensorOperation::zero(dFpDT);
+        // 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 
     
@@ -2109,7 +2132,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
       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);
+            DinvFpdT(i,j) -= Fpinv(i,p)*dFpdT(p,q)*Fpinv(q,j);
           }
     
     STensorOperation::zero(dFedT);
@@ -2173,6 +2196,10 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
       }
   
   
+  }
+  
+  if(stiff){
+  
   // TVP - flux, thermalEnergy, thermalSource
     if (!_thermalEstimationPreviousConfig){
       
@@ -2207,22 +2234,65 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
         }
       }
     }
+  }
 
+  // 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; // DEBUGGING
+
+  getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,&Wm);
+  // mechanicalSource += Wm;    // DEBUGGING 
+    
+  // 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 energy
-    q1->_thermalEnergy = CpT*T; 
+    // 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 = -CpT/this->getTimeStep() - (CpT-CpT_0)/this->getTimeStep();
+     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();
@@ -2258,16 +2328,16 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
     
     dmechanicalSourcedT = 0.;
     STensorOperation::zero(dmechanicalSourceF);
-    dmechanicalSourceF += dWmdF_TVE;
-    dmechanicalSourcedT += dWmdT_TVE;
+    // dmechanicalSourceF += dWmdF_TVE; // DEBUGGING 
+    // dmechanicalSourcedT += dWmdT_TVE; // DEBUGGING
     
     double& dWmdT = q1->getRefTodMechSrcTVPdT();
     STensor3& dWmdF = q1->getRefTodMechSrcTVPdF();
 
-    getMechSourceTVP(q0,q1,T0,T,dFpDT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT); 
+    getMechSourceTVP(q0,q1,T0,T,dFpdT,dFpdF,DphiPDF,NULL,&dWmdF,&dWmdT); 
     
-    dmechanicalSourceF += dWmdF;
-    dmechanicalSourcedT += dWmdT;
+    // dmechanicalSourceF += dWmdF;   // DEBUGGING 
+    // dmechanicalSourcedT += dWmdT;  // DEBUGGING 
 }
 };
   
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
index eaef7b8ecd93f1edeed564c893dbd70d258b94d1..922286a716d4ae3e6cbc8f83a340f0379460d07f 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.h
@@ -22,15 +22,16 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
         
         scalarFunction* _temFunc_Sy0; // temperature dependence of initial yield stress
         scalarFunction* _temFunc_H; // temperature dependence of hardening stress    
-    
+        scalarFunction* _temFunc_visc; // temperature dependence of viscosity law
+        scalarFunction* _temFunc_Hb; // temperature dependence of kinematic hardening modulus  
+        
         std::vector<double> _HR; // isotropic hardening moduli for R
         
     protected:
     
-        // NEW
         virtual void hardening(const IPNonLinearTVP* q0, IPNonLinearTVP* q1, const double& T) const;
         virtual void getModifiedMandel(const STensor3& C, const IPNonLinearTVP *q0_, IPNonLinearTVP *q1) const;
-        virtual void checkCommutavity(const STensor3& C, const STensor3& S, const IPNonLinearTVP *q1) const;
+        virtual void checkCommutavity(STensor3& commuteCheck, const STensor3& Ce, const STensor3& S, const IPNonLinearTVP *q1) const;
         virtual void getChabocheCoeffs(fullVector<double>& coeffs, const double& opt, const IPNonLinearTVP *q1) const;
         virtual void getPlasticPotential(const double& phiP, const double& phiE, double& Px) const;
         virtual void getYieldCoefficientTempDers(const IPNonLinearTVP *q1, const double& T, fullVector<double>& DcoeffsDT, fullVector<double>& DDcoeffsDTT) const;
@@ -43,7 +44,6 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
                                       
         virtual void getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP *q1, const double& T0, const double& T,
                                       double *psiTVM = NULL, double *DpsiTVMdT = NULL, double *DDpsiTVMdTT = NULL) const;
-        // NEW
         
         virtual void 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, 
@@ -159,6 +159,19 @@ class mlawNonLinearTVP : public mlawNonLinearTVE{
         void setTemperatureFunction_Hardening(const scalarFunction& Tfunc){
             if (_temFunc_H != NULL) delete _temFunc_H;
             _temFunc_H = Tfunc.clone();
+            // _compression->setTemperatureFunction(Tfunc);
+            // _traction->setTemperatureFunction(Tfunc); // - this is difficult, because of multiple parameters
+        }
+        void setTemperatureFunction_KinematicHardening(const scalarFunction& Tfunc){
+            if (_temFunc_Hb != NULL) delete _temFunc_Hb;
+            _temFunc_Hb = Tfunc.clone();
+            // _kinematic->setTemperatureFunction_K(Tfunc); 
+                    // - this is not a pure virtual parent class, on the other hand the viscosity/kine harde type isnt set here but in the .py file.
+        }
+        void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc){
+            if (_temFunc_visc != NULL) delete _temFunc_visc;
+            _temFunc_visc = Tfunc.clone();
+            // _viscosity->setTemperatureFunction(Tfunc);
         }
         
         void setTaylorQuineyFactor(const double f){_TaylorQuineyFactor = f;};
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 9b814e9329d2d206a2b836707a240e03524f0f91..074c235685f1464342f4e05f1a75d5052b929888 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -245,7 +245,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
   int numCurlDof = ipvcur->getNumConstitutiveCurlVariable();
 
   _stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
-#if 1
+#if 0
   _stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
   ipvcur->getConstRefToTangentModuli().print("dPdE Analytique");
   for (int k=0; k< numNonLocalVars; k++)
@@ -607,7 +607,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
         }
 
     }
-#if 0
+#if 1
     ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");
     for (int k=0; k< numNonLocalVars; k++)
     {
@@ -8575,7 +8575,6 @@ void NonLinearTVPDG3DMaterialLaw::setTaylorQuineyFactor(const double f){
 void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_InitialYieldStress(const scalarFunction& Tfunc){
   _viscoLaw.setTemperatureFunction_InitialYieldStress(Tfunc);
 };
-/*
 void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_Hardening(const scalarFunction& Tfunc){
   _viscoLaw.setTemperatureFunction_Hardening(Tfunc);
 };
@@ -8585,7 +8584,6 @@ void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_KinematicHardening(cons
 void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_Viscosity(const scalarFunction& Tfunc){
   _viscoLaw.setTemperatureFunction_Viscosity(Tfunc);
 };
-*/
 
 // Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial   (added from LinearThermoMech)
 void NonLinearTVPDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index a70797517e2d65a6a40ca3246fd15ba738df44c6..234321640135c11b24053a1687bf9c35e0bb1a8a 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -3000,9 +3000,9 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
 
     void setTaylorQuineyFactor(const double f);
     void setTemperatureFunction_InitialYieldStress(const scalarFunction& Tfunc);
-    /*void setTemperatureFunction_Hardening(const scalarFunction& Tfunc);
+    void setTemperatureFunction_Hardening(const scalarFunction& Tfunc);
     void setTemperatureFunction_KinematicHardening(const scalarFunction& Tfunc);
-    void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc);*/
+    void setTemperatureFunction_Viscosity(const scalarFunction& Tfunc);
 
     #ifndef SWIG
      NonLinearTVPDG3DMaterialLaw(const NonLinearTVPDG3DMaterialLaw &source);