diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
index bd98287baae05ba83d36f08e861f582895aeef04..7b81fbf99829b11a98692428008afbb6367d54c6 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -385,14 +385,12 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   STensorOperation::zero(Ge_Tensor_pr); STensorOperation::zero(Ge_Tensor);
   STensorOperation::zero(Ge_TrEeTensor_pr); STensorOperation::zero(Ge_TrEeTensor);
   
-  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3){
+  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
       // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
-        // mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
 	  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
         Ge_TrEeTensor_pr = Ge_TrEeTensor;
   }
   else{ 
-        // mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
 	  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm);
   }
     
@@ -475,7 +473,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   devPhi = devPhipr;
   
   // Initialise the rest
-  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3){
+  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
     // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
     getDho(Gamma,Cepr,Ceinvpr,Kepr,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
   }
@@ -628,7 +626,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                 }
                 
                 // Update Phi
-                if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3){
+                if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
                     getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
                                     Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
@@ -681,7 +679,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
             }
 
             // Correct Phi
-            if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                 // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
                 getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
                     Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index a21489b30b807e4a881cf6f030e6eec4ca8f85c1..101e7fe590953c0cb51e69facc2d85c89b116ead 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -21,7 +21,7 @@ mlawNonLinearTVM::mlawNonLinearTVM(const int num,const double E,const double nu,
      mlawPowerYieldHyper(num, E, nu, rho, tol, matrixbyPerturbation, pert),
      _Tinitial(Tinitial),_scalarAlpha(Alpha),_scalarK(KThCon),_Cp(Cp),
      _thermalEstimationPreviousConfig(thermalEstimationPreviousConfig), _mullinsEffect(NULL), _elasticPotential(NULL), 
-     _useExtraBranch(false), _useExtraBranch_TVE(false), _ExtraBranch_TVE_option(2){
+     _useExtraBranch(false), _useExtraBranch_TVE(false), _ExtraBranch_TVE_option(2), _V1(0.),_V2(0.),_D1(0.),_D2(0.),_Ci(0.){
 
   _G = _mu;
 
@@ -61,7 +61,8 @@ mlawNonLinearTVM::mlawNonLinearTVM(const mlawNonLinearTVM& src): mlawPowerYieldH
   _useExtraBranch = src._useExtraBranch;
   _useExtraBranch_TVE = src._useExtraBranch_TVE;
   _ExtraBranch_TVE_option = src._ExtraBranch_TVE_option;
-    
+  _V1 = src._V1; _V2 = src._V2; _D1 = src._D1; _D2 = src._D2; _Ci = src._Ci;
+
   _temFunc_K = NULL;                                                                    // bulk modulus
   if (src._temFunc_K != NULL){ _temFunc_K = src._temFunc_K->clone();}
 
@@ -110,6 +111,7 @@ mlawNonLinearTVM& mlawNonLinearTVM::operator=(const materialLaw &source){
     _useExtraBranch = src->_useExtraBranch;
     _useExtraBranch_TVE = src->_useExtraBranch_TVE;
     _ExtraBranch_TVE_option = src->_ExtraBranch_TVE_option;
+    _V1 = src->_V1; _V2 = src->_V2; _D1 = src->_D1; _D2 = src->_D2; _Ci = src->_Ci;
     
     if(_temFunc_K != NULL) delete _temFunc_K;                                                   // bulk modulus
     if (src->_temFunc_K != NULL){ _temFunc_K = src->_temFunc_K->clone();}
@@ -159,10 +161,10 @@ void mlawNonLinearTVM::setViscoElasticNumberOfElement(const int N){
   Msg::Info("Numer of Spring/Dashpot for viscoelastic model: %d",_N);
   _Ki.clear(); _ki.clear();
   _Gi.clear(); _gi.clear();
-  _Ki.resize(_N,0.);
-  _ki.resize(_N,0.);
-  _Gi.resize(_N,0.);
-  _gi.resize(_N,0.);
+  _Ki.resize(_N,0.); _ki.resize(_N,0.);
+  _Gi.resize(_N,0.); _gi.resize(_N,0.);
+  _V1.clear(); _V2.clear(); _D1.clear(); _D2.clear(); _Ci.clear();
+  _V1.resize(_N,0.); _V2.resize(_N,0.); _D1.resize(_N,0.); _D2.resize(_N,0.); _Ci.resize(_N,0.);
 };
 
 void mlawNonLinearTVM::setViscoElasticData(const int i, const double Ei, const double taui){
@@ -182,6 +184,26 @@ void mlawNonLinearTVM::setViscoElasticData(const int i, const double Ei, const d
   }
 };
 
+void mlawNonLinearTVM::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
+  if (i> _N or i<1)
+    Msg::Error("This setting is invalid %d > %d",i,_N);
+  else{
+    _V1[i-1] = V1; _V2[i-1] = V2; _D1[i-1] = D1; _D2[i-1] = D2;
+
+    Msg::Info("setting: element=%d, V1 = %e, V2 = %e, D1 = %e, D2 = %e",i-1,V1,V2,D1,D2);
+  }
+};
+
+void mlawNonLinearTVM::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
+  if (i> _N or i<1)
+    Msg::Error("This setting is invalid %d > %d",i,_N);
+  else{
+    _Ci[i-1] = Ci;
+
+    Msg::Info("setting: element=%d, Ci = %e",i-1,Ci);
+  }
+};
+
 void mlawNonLinearTVM::setViscoElasticData_Bulk(const int i, const double Ki, const double ki){
   if (i> _N or i<1)
     Msg::Error("This setting is invalid %d > %d",i,_N);
@@ -410,7 +432,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                     // Psy_mech += q._psi_branch[i];*/
                     Msg::Info("setting: _ExtraBranch_TVE_option=%d integrates the viscoelastic stress and currently has no definition for deformation energy. Don't use it with     mullin's effect",_ExtraBranch_TVE_option);
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){ 
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     Psy_mech += _Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i];
                 }
             }
@@ -529,7 +551,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                     Msg::Info("setting: _ExtraBranch_TVE_option=%d integrates the viscoelastic stress and currently has no definition for deformation energy. Don't use it with     mullin's effect",_ExtraBranch_TVE_option);
                     
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){ 
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     
                     double dTrEei_DT(0.); // branchElasticStrain trace
                     static STensor3 dDevEei_DT;
@@ -548,7 +570,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                                         (*DpsiDE)(j,k) += 2*_Gi[i]*q._Bd_TVE_vector[i]*q._A[i](p,r)*exp_mid_g*_Idev(p,r,j,k); 
                         }
                     }
-                    else { // _ExtraBranch_TVE_option == 3
+                    else { // _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4
                     
                         // NOTE : _Av_TVE_vector, _Bd_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector have CTasymm already
                         //          Therefore, to derive CTasymm, first divide to remove it and then derive.
@@ -838,7 +860,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
   }
 };
 
-void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const double& T, 
+void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,  const double& T,
                         double& Av, double& dA_dTrEe, double& intA, double& intA1, double& intA2, double& dA_dT, double& ddA_dTT, double& ddA_dTdTrEe,
                         double& Bd, STensor3& dB_dDevEe, double &intB, STensor3 &intB1, STensor3 &intB2, double &dB_dT, double &ddB_dTT, STensor3& ddB_dTdDevEe,
                         double* dB_dTrEe) const{
@@ -1014,6 +1036,60 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const doubl
             intB = 1.;
         }   
     }
+    else if (_ExtraBranch_TVE_option == 4){
+
+            // Tension-Compression Asymmetry
+            // Standard form => A_v = a/(sqrt(b*tr*tr+c))
+
+            // Regularising function
+            double m = _tensionCompressionRegularisation;
+            double expmtr = exp(-m*tr);
+            // if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
+            double sigmoid = 1/(1.+expmtr);
+
+            // Av
+            double x = _V1[i]*tr*tr + _V2[i];
+            Av = sigmoid*1./sqrt(x) + (1.-sigmoid)*(_Ci[i]/sqrt(x));
+            Av *= getVolumeCorrection();
+
+            dA_dTrEe = - sigmoid*_V1[i]*tr/pow(x,1.5) - (1.-sigmoid)*(_V1[i]*_Ci[i]*tr/pow(x,1.5))
+                       + (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x);
+            dA_dTrEe *= getVolumeCorrection();
+
+            if(_V1[i]>0.){
+                double integrand = 1./_V1[i]*sqrt(x); // integral of A_v * trEe
+                integrand -= ( 1./_V1[i]*sqrt(_V2[i]) ); // value at trEe = 0.
+                intA = sigmoid*integrand + _Ci[i]*(1.-sigmoid)*integrand;
+                intA *= getVolumeCorrection();
+            }
+            else{
+                intA = 1.;
+            }
+
+            // Bd
+            double y = _D1[i]*dev.dotprod() + _D2[i];
+            Bd = sigmoid*1./sqrt(y) + (1.-sigmoid)*(_Ci[i]/sqrt(y));
+            Bd *= getDevCorrection();
+
+            STensorOperation::zero(dB_dDevEe);
+            dB_dDevEe = dev;
+            dB_dDevEe *= (-sigmoid*_D1[i]/pow(y,1.5) - (1.-sigmoid)*(_D1[i]*_Ci[i]/pow(y,1.5)));
+            dB_dDevEe *= getDevCorrection();
+
+            if(dB_dTrEe!=NULL){
+                *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y) - (_Ci[i]*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y);
+                *dB_dTrEe *= getDevCorrection();
+            }
+            if(_D1[i]>0.){
+                double integrand = 1./_D1[i]*sqrt(y);  // integral of B_d * devEe
+                integrand -= (1./_D1[i]*sqrt(_D2[i]) );  // value at devEe = 0.
+                intB = sigmoid*integrand + _Ci[i]*(1.-sigmoid)*integrand;
+                intB *= getDevCorrection();
+            }
+            else{
+                intB = 1.;
+            }
+        }
     else{ // Test
         // Av = 1.;
         // intA1 = tr;
@@ -1201,7 +1277,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
             // NEW
             
             if(_ExtraBranch_TVE_option == 1){
-                mlawNonLinearTVM::extraBranch_nonLinearTVE(Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+                mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
                 q1->_Av_TVE = Av; q1->_Bd_TVE = Bd; q1->_dAv_dTrEe_TVE = dA_dTrEe; q1->_dBd_dDevEe_TVE = dB_dDevEe;
                 q1->_intAv_TVE = intA; q1->_intBd_TVE = intB; q1->_intAv_1_TVE = intA1; q1->_intBd_1_TVE = intB1;  // this is zero
                 q1->_intAv_2_TVE = intA2; q1->_intBd_2_TVE = intB2;  // this is zero
@@ -1289,7 +1365,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     p += ( _Ki[i]*q1->_B[i] );
                     
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     
                     // Single Convolution - get branchElasticStrain
                     for (int k=0; k<3; k++)
@@ -1304,12 +1380,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     branchElasticStrain(1,1) += 1./3.*q1->_B[i];
                     branchElasticStrain(2,2) += 1./3.*q1->_B[i];
                     if (_ExtraBranch_TVE_option == 2){
-                        mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                        mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
                                                                 Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
                     }
                     else{
                         double dB_dTrEe(0.);
-                        mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                        mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
                                                                 Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
                         q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
                     }
@@ -1327,7 +1403,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                             for (int p=0; p<3; p++)
                                 for (int q=0; q<3; q++){
                                     Bd_stiffnessTerm(k,l,p,q) +=  2*_Gi[i]*q1->_A[i](k,l)*dB_dDevEe(p,q) * exp_mid_g; 
-                                    if (_ExtraBranch_TVE_option == 3){
+                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                         (*Bd_stiffnessTerm2)(k,l,p,q) += 2*_Gi[i]*q1->_dBd_dTrEe_TVE_vector[i]*q1->_A[i](k,l)*_I(p,q) * exp_mid_k;
                                             // exp_mid_k = d(tr branchElasticStrain)_d(trEe)
                                     }
@@ -1441,7 +1517,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     // Here the derivatives are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                     // 2nd order derivatives of _Av_TVE_vector and _Bd_TVE_vector are ignored
                     
@@ -1466,7 +1542,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                                 DDA_DTT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*ddTrEei_DTT +
                                                     2*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT*dDevEei_DT(k,l);  
@@ -1548,7 +1624,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                                                                                 q0->_A[i](k,l) * 1/q0->_Bd_TVE * temp2); 
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){      
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     // Here the derivative are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                     // 2nd order derivatives of _Av_TVE_vector and _Bd_TVE_vector are ignored
                     
@@ -1564,7 +1640,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                         for (int l=0; l<3; l++){
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                             }
                             
@@ -1768,7 +1844,7 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
     // Add visco terms
     if ((_Ki.size() > 0) or (_Gi.size() > 0)){
         
-        if (_ExtraBranch_TVE_option == 2){ 
+        if (_ExtraBranch_TVE_option != 1){
         
             // TVE Mechanical Source Term
             
@@ -1854,7 +1930,7 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
                                     DDQiDTDEe[i](k,l,r,s) = _Ki[i]*q1->_DDB_DTDtrE[i]*_I(k,l)*_I(r,s);
                                     dot_DGammaiDEe[i](k,l,r,s) = 1/3*(1.- exp_mid_k)*_I(k,l)*_I(r,s)  + _Idev(k,l,r,s);
                                     
-                                    if (_ExtraBranch_TVE_option == 3){
+                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                         DQiDEe[i](k,l,r,s) +=  2*_Gi[i]*q1->_DA_DtrE[i](k,l)*_I(r,s);
                                         DDQiDTDEe[i](k,l,r,s) +=  2*_Gi[i]*q1->_DDA_DTDtrE[i](k,l)*_I(r,s);
                                     }
@@ -1994,7 +2070,7 @@ void mlawNonLinearTVM::getTVEdCorKirDT(const IPNonLinearTVM *q0, IPNonLinearTVM
                             else{ // all extraBranches, because _Idev was not multiplied to _DDA_DTDdevE in the predictorCorrector function
                                 DDcorKirDTDEe(j,k,r,s) += _Ki[i]*q1->_DDB_DTDtrE[i]*_I(j,k)*_I(r,s);
                                 
-                                if (_ExtraBranch_TVE_option == 3){
+                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                     DDcorKirDTDEe(j,k,r,s) += 2*_Gi[i]*q1->_DDA_DTDtrE[i](j,k)*_I(r,s);
                                 }
                                 
@@ -2222,7 +2298,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
     static STensor43 Bd_stiffnessTerm, Bd_stiffnessTerm2;
     STensorOperation::zero(Bd_stiffnessTerm);
     STensorOperation::zero(Bd_stiffnessTerm2);
-    if(_ExtraBranch_TVE_option == 3 || _extraBranchNLType == TensionCompressionRegularisedType){
+    if(_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
         ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm,&Bd_stiffnessTerm2); 
                 // Updates the values of moduli to the current temperatures and calculate stresses.
     }
@@ -2716,7 +2792,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
         
         if (_useExtraBranch_TVE){                            
             if(_ExtraBranch_TVE_option == 1){
-                mlawNonLinearTVM::extraBranch_nonLinearTVE(Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
+                mlawNonLinearTVM::extraBranch_nonLinearTVE(0,Ee,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
                 q1->_Av_TVE = Av; q1->_Bd_TVE = Bd; q1->_dAv_dTrEe_TVE = dA_dTrEe; q1->_dBd_dDevEe_TVE = dB_dDevEe;
                 q1->_intAv_TVE = intA; q1->_intBd_TVE = intB;
                 q1->_intAv_1_TVE = intA1; q1->_intBd_1_TVE = intB1;  // this is zero
@@ -2793,7 +2869,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 // Add to volumetric stress
                 p += ( _Ki[i]*q1->_B[i] );
              }
-             else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
+             else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     
                 // Single Convolution - get branchElasticStrain
                 for (int k=0; k<3; k++)
@@ -2808,12 +2884,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 branchElasticStrain(1,1) += 1./3.*q1->_B[i];
                 branchElasticStrain(2,2) += 1./3.*q1->_B[i];
                 if (_ExtraBranch_TVE_option == 2){
-                    mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                    mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
                                                                 Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe);
                 }
                 else{
                     double dB_dTrEe(0.);
-                    mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                    mlawNonLinearTVM::extraBranch_nonLinearTVE(i,branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
                                                                 Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
                     q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
                 }
@@ -2831,7 +2907,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                         for (int p=0; p<3; p++)
                             for (int q=0; q<3; q++){
                                 Bd_stiffnessTerm(k,l,p,q) += 2*_Gi[i]*q1->_A[i](k,l)*dB_dDevEe(p,q) * exp_mid_g;
-                                if (_ExtraBranch_TVE_option == 3){
+                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                     (*Bd_stiffnessTerm2)(k,l,p,q) += 2*_Gi[i]*q1->_dBd_dTrEe_TVE_vector[i]*q1->_A[i](k,l)*_I(p,q) * exp_mid_k;
                                 } 
                             }
@@ -2941,7 +3017,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     // Here the derivative are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                     
                     double dTrEei_DT(0.), ddTrEei_DTT(0.); // branchElasticStrain trace
@@ -2965,7 +3041,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                                 DDA_DTT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*ddTrEei_DTT +
                                                     2*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT*dDevEei_DT(k,l);  
@@ -3047,7 +3123,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                                                                                 q0->_A[i](k,l) * 1/q0->_Bd_TVE * temp2); 
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){      
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                     // Here the derivative are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                      
                     double dTrEei_DT(0.); // branchElasticStrain trace
@@ -3062,7 +3138,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                         for (int l=0; l<3; l++){
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                             }
                             
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
index 64d2a130b01d3c35577e4e5442c4edf6f431ae19..b34b05e624a54a0a4b8b01ba676b6e1f892bae51 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
@@ -66,6 +66,13 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
     scalarFunction* _temFunc_elasticCorrection_Bulk;
     scalarFunction* _temFunc_elasticCorrection_Shear;
     
+  // Extrabranch Corrector parameters
+    std::vector<double> _V1;
+    std::vector<double> _V2;
+    std::vector<double> _D1;
+    std::vector<double> _D2;
+    std::vector<double> _Ci;
+
   #ifndef SWIG  
   protected:
   // Const Functions for constitutive function - To update any variable, pass it as an argument by reference.
@@ -94,7 +101,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
                         double* DDsigV_dTdTrEe = NULL, STensor43* DDsigD_dTdDevEe = NULL,
                         STensor3* DsigD_dTrEe = NULL, STensor3* DsigD_dTdTrEe = NULL) const;
     
-    virtual void extraBranch_nonLinearTVE(const STensor3& Ee,  const double& T, 
+    virtual void extraBranch_nonLinearTVE(const int i, const STensor3& Ee,  const double& T,
                         double& Av, double& dA_dTrEe, double& intA, double& intA1, double& intA2, double& dA_dT, double& ddA_dTT, double& ddA_dTdTrEe,
                         double& Bd, STensor3& dB_dDevEe, double &intB, STensor3 &intB1, STensor3 &intB2, double &dB_dT, double &ddB_dTT, STensor3& ddB_dTdDevEe,
                         double* dB_dTrEe = NULL) const;
@@ -186,6 +193,8 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
     virtual void setViscoElasticData(const int i, const double Ei, const double taui);
     virtual void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
     virtual void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
+    virtual void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    virtual void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     virtual void setViscoElasticData(const std::string filename);
     
   // Non-Const Functions for temp and temp_functions that are called from .py file.  ---- These WILL change the values of the variables defined in the class object.
@@ -339,4 +348,4 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
     virtual void checkInternalState(IPVariable* ipv, const IPVariable* ipvprev) const{}; // do nothing
     #endif // SWIG
 };
-#endif //mlawNonLinearTVM
\ No newline at end of file
+#endif //mlawNonLinearTVM
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index bb1f96e5a58056f33e28d35c213f9be63d353855..4e93b930656dae7e1bffe429c1f5ecfa5a061e01 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -10270,18 +10270,24 @@ void NonLinearTVMDG3DMaterialLaw::useExtraBranchBool(const bool useExtraBranch){
 void NonLinearTVMDG3DMaterialLaw::useExtraBranchBool_TVE(const bool useExtraBranch_TVE, const int ExtraBranch_TVE_option){
   _viscoLaw.useExtraBranchBool_TVE(useExtraBranch_TVE, ExtraBranch_TVE_option);
 };
-void NonLinearTVMDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+void NonLinearTVMDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVMDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+void NonLinearTVMDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
-void NonLinearTVMDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+void NonLinearTVMDG3DMaterialLaw::setTensionCompressionRegularisation(const double tensionCompressionRegularisation){
   _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
 void NonLinearTVMDG3DMaterialLaw::setExtraBranchNLType(const int method){
   _viscoLaw.setExtraBranchNLType(method);
 };
+void NonLinearTVMDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
+  _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
+};
+void NonLinearTVMDG3DMaterialLaw::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
+  _viscoLaw.setCompressionCorrectionsAllBranchesTVE(i,Ci);
+};
 void NonLinearTVMDG3DMaterialLaw::setMullinsEffect(const mullinsEffect& mullins){
   _viscoLaw.setMullinsEffect(mullins);
 };
@@ -10804,16 +10810,16 @@ void NonLinearTVPDG3DMaterialLaw::setTemperatureFunction_Viscosity(const scalarF
 void NonLinearTVPDG3DMaterialLaw::useExtraBranchBool(const bool useExtraBranch){
   _viscoLaw.useExtraBranchBool(useExtraBranch);
 };
-void NonLinearTVPDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+void NonLinearTVPDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+void NonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
 void NonLinearTVPDG3DMaterialLaw::setExtraBranchNLType(const int method){
   _viscoLaw.setExtraBranchNLType(method);
 };
-void NonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+void NonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisation(const double tensionCompressionRegularisation){
   _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
 // Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial   (added from LinearThermoMech)
@@ -11083,18 +11089,25 @@ void NonLinearTVENonLinearTVPDG3DMaterialLaw::useExtraBranchBool(const bool useE
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::useExtraBranchBool_TVE(const bool useExtraBranch_TVE, const int ExtraBranch_TVE_option){
   _viscoLaw.useExtraBranchBool_TVE(useExtraBranch_TVE, ExtraBranch_TVE_option);
 };
-void NonLinearTVENonLinearTVPDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranchNLType(const int method){
   _viscoLaw.setExtraBranchNLType(method);
 };
-void NonLinearTVENonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisation(const double tensionCompressionRegularisation){
   _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
+  _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
+};
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
+  _viscoLaw.setCompressionCorrectionsAllBranchesTVE(i,Ci);
+};
+
 // Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial   (added from LinearThermoMech)
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::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 20aaacf2d10764416e83b10c9744b6b24dd45424..0c2e5cba07b02f64be553e750057036396e7e0c4 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -3279,10 +3279,12 @@ class NonLinearTVMDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
     void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
     void setViscoElasticData(const std::string filename);
-    void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
-    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
-    void setExtraBranchNLType(int type);
-    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
+    void setExtraBranchNLType(const int type);
+    void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
+    void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     
     void setReferenceTemperature(const double Tref);
     void setShiftFactorConstantsWLF(const double C1, const double C2);
@@ -3360,10 +3362,10 @@ class NonLinearTVEDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
     void setViscoElasticData_Alpha(const int i, const double Alphai);
     void setViscoElasticData(const std::string filename);
-    void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
-    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
-    void setExtraBranchNLType(int type);
-    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
+    void setExtraBranchNLType(const int type);
+    void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
     
     void setReferenceTemperature(const double Tref);
     void setTempFuncOption(const double TemFuncOpt);
@@ -3476,10 +3478,10 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     // void setTemperatureFunction_BranchShearModuli(const scalarFunction& Tfunc,int i);
     // void setTemperatureFunction_BranchThermalDilationCoefficient(const scalarFunction& Tfunc,int i);
     void useExtraBranchBool(const bool useExtraBranch);
-    void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
-    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
-    void setExtraBranchNLType(int type);
-    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
+    void setExtraBranchNLType(const int type);
+    void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
     
     void setCompressionHardening(const J2IsotropicHardening& comp);
     void setTractionHardening(const J2IsotropicHardening& trac);
@@ -3583,10 +3585,12 @@ class NonLinearTVENonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // pub
     // void setTemperatureFunction_BranchThermalDilationCoefficient(const scalarFunction& Tfunc,int i);
     void useExtraBranchBool(const bool useExtraBranch);
     void useExtraBranchBool_TVE(const bool useExtraBranch_TVE, const int ExtraBranch_TVE_option);
-    void setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc);
-    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
-    void setExtraBranchNLType(int type);
-    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
+    void setExtraBranchNLType(const int type);
+    void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
+    void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     
     void setCompressionHardening(const J2IsotropicHardening& comp);
     void setTractionHardening(const J2IsotropicHardening& trac);