diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
index e5f8c9693c49ad7b143e7422e3a781f6758fcc74..9a2db9172fee0a70a0281008f9a03b75261b46f8 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVM.cpp
@@ -19,30 +19,26 @@ IPNonLinearTVM::IPNonLinearTVM(const J2IsotropicHardening* comp,
                      _dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.), 
                      _trCorKirinf_TVE(0.), _devCorKirinf_TVE(0.), _corKirExtra(0.), 
                      _psiMax(0.), _mullinsDamage(0.), _DpsiDT(0.), _DpsiDE(0.), _DmullinsDamage_Dpsi_cap(0.), _DmullinsDamage_DT(0.), _P_cap(0.),
-                     _pressure(0.), _Av_TVE(1.), _Bd_TVE(1.), _dAv_dTrEe(0.), _dAv_dTrEe_TVE(0.), _intAv_TVE(0.), _intBd_TVE(0.), 
+                     _pressure(0.), _Av_TVE(1.), _Bd_TVE(1.), _dAv_dTrEe(0.), _dBd_dTrEe(0.), _dAv_dTrEe_TVE(0.), _intAv_TVE(0.), _intBd_TVE(0.), 
                      _intAv_1_TVE(0.), _intBd_1_TVE(0.), _intAv_2_TVE(0.), _intBd_2_TVE(0.), _dBd_dDevEe(0.), _dBd_dDevEe_TVE(0.),
                      _viscousDissipatedEnergy(0.), _mullinsDissipatedEnergy(0.){
         
-    _DA_DT.clear();
-    _DDA_DTT.clear();
-    _DA_DdevE.clear();
-    _DDA_DTDdevE.clear();
-    _DB_DT.clear();
-    _DDB_DTT.clear();
-    _DB_DtrE.clear();
-    _DDB_DTDtrE.clear();
-    _Av_TVE_vector.clear();
-    _Bd_TVE_vector.clear();
-    _dAv_dTrEe_TVE_vector.clear();
-    _intAv_TVE_vector.clear();
-    _intBd_TVE_vector.clear();
-    _dBd_dDevEe_TVE_vector.clear();
+    _DA_DT.clear(); _DDA_DTT.clear(); 
+    _DA_DtrE.clear(); _DDA_DTDtrE.clear();
+    _DA_DdevE.clear(); _DDA_DTDdevE.clear();
+    _DB_DT.clear(); _DDB_DTT.clear();
+    _DB_DtrE.clear(); _DDB_DTDtrE.clear();
+    _Av_TVE_vector.clear(); _dAv_dTrEe_TVE_vector.clear();
+    _Bd_TVE_vector.clear(); _dBd_dTrEe_TVE_vector.clear(); _dBd_dDevEe_TVE_vector.clear();
+    _intAv_TVE_vector.clear(); _intBd_TVE_vector.clear();
         
     for (int i=0; i < N; i++){
         STensor3 el(0.);
         STensor43 el43(0.);
         _DA_DT.push_back(el);
         _DDA_DTT.push_back(el);
+        _DA_DtrE.push_back(el);
+        _DDA_DTDtrE.push_back(el);
         _dBd_dDevEe_TVE_vector.push_back(el);
         _DB_DT.push_back(0.);
         _DDB_DTT.push_back(0.);
@@ -51,6 +47,7 @@ IPNonLinearTVM::IPNonLinearTVM(const J2IsotropicHardening* comp,
         _Av_TVE_vector.push_back(0.);
         _Bd_TVE_vector.push_back(0.);
         _dAv_dTrEe_TVE_vector.push_back(0.);
+        _dBd_dTrEe_TVE_vector.push_back(0.);
         _intAv_TVE_vector.push_back(0.);
         _intBd_TVE_vector.push_back(0.);
         _DA_DdevE.push_back(el43);
@@ -63,7 +60,7 @@ IPNonLinearTVM::IPNonLinearTVM(const J2IsotropicHardening* comp,
 };
 
 IPNonLinearTVM::IPNonLinearTVM(const IPNonLinearTVM& src): IPHyperViscoElastoPlastic(src),
-            _DA_DT(src._DA_DT), _DDA_DTT(src._DDA_DTT), _DA_DdevE(src._DA_DdevE), _DDA_DTDdevE(src._DDA_DTDdevE),
+            _DA_DT(src._DA_DT), _DDA_DTT(src._DDA_DTT), _DA_DdevE(src._DA_DdevE), _DDA_DTDdevE(src._DDA_DTDdevE), _DA_DtrE(src._DA_DT), _DDA_DTDtrE(src._DDA_DTT),
             _DB_DT(src._DB_DT), _DDB_DTT(src._DDB_DTT), _DB_DtrE(src._DB_DtrE), _DDB_DTDtrE(src._DDB_DTDtrE),
             _T(src._T),_devDE(src._devDE), _trDE(src._trDE), _DE(src._DE), _DcorKirDT(src._DcorKirDT), _DDcorKirDTT(src._DDcorKirDTT), _DDcorKirDTDE(src._DDcorKirDTDE),
             _mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdE(src._DmechSrcTVEdE),
@@ -73,11 +70,11 @@ IPNonLinearTVM::IPNonLinearTVM(const IPNonLinearTVM& src): IPHyperViscoElastoPla
             _mullinsDamage(src._mullinsDamage), _DpsiDT(src._DpsiDT), _DpsiDE(src._DpsiDE),
             _DmullinsDamage_Dpsi_cap(src._DmullinsDamage_Dpsi_cap), _DmullinsDamage_DT(src._DmullinsDamage_DT), _P_cap(src._P_cap),
             _pressure(src._pressure), _Av_TVE(src._Av_TVE), _Bd_TVE(src._Bd_TVE), 
-            _dAv_dTrEe(src._dAv_dTrEe), _dAv_dTrEe_TVE(src._dAv_dTrEe_TVE), 
+            _dAv_dTrEe(src._dAv_dTrEe), _dBd_dTrEe(src._dBd_dTrEe), _dAv_dTrEe_TVE(src._dAv_dTrEe_TVE), 
             _intAv_TVE(src._intAv_TVE), _intBd_TVE(src._intBd_TVE),
             _intAv_1_TVE(src._intAv_1_TVE), _intBd_1_TVE(src._intBd_1_TVE), _intAv_2_TVE(src._intAv_2_TVE), _intBd_2_TVE(src._intBd_2_TVE),
             _dBd_dDevEe(src._dBd_dDevEe), _dBd_dDevEe_TVE(src._dBd_dDevEe_TVE),
-            _Av_TVE_vector(src._Av_TVE_vector), _Bd_TVE_vector(src._Bd_TVE_vector), _dAv_dTrEe_TVE_vector(src._dAv_dTrEe_TVE_vector),
+            _Av_TVE_vector(src._Av_TVE_vector), _Bd_TVE_vector(src._Bd_TVE_vector), _dAv_dTrEe_TVE_vector(src._dAv_dTrEe_TVE_vector), _dBd_dTrEe_TVE_vector(src._dBd_dTrEe_TVE_vector),
             _intAv_TVE_vector(src._intAv_TVE_vector), _intBd_TVE_vector(src._intBd_TVE_vector), _dBd_dDevEe_TVE_vector(src._dBd_dDevEe_TVE_vector),
             _viscousDissipatedEnergy(src._viscousDissipatedEnergy), _mullinsDissipatedEnergy(src._mullinsDissipatedEnergy){
                 
@@ -95,6 +92,8 @@ IPNonLinearTVM& IPNonLinearTVM::operator=(const IPVariable& src)
   {
     _DA_DT = psrc->_DA_DT; 
     _DDA_DTT = psrc->_DDA_DTT; 
+    _DA_DtrE = psrc->_DA_DtrE; 
+    _DDA_DTDtrE = psrc->_DDA_DTDtrE; 
     _DA_DdevE = psrc->_DA_DdevE; 
     _DDA_DTDdevE = psrc->_DDA_DTDdevE; 
     _DB_DT = psrc->_DB_DT; 
@@ -128,6 +127,7 @@ IPNonLinearTVM& IPNonLinearTVM::operator=(const IPVariable& src)
     _Av_TVE = psrc->_Av_TVE;
     _Bd_TVE = psrc->_Bd_TVE;
     _dAv_dTrEe = psrc->_dAv_dTrEe;
+    _dBd_dTrEe = psrc->_dBd_dTrEe;
     _dAv_dTrEe_TVE = psrc->_dAv_dTrEe_TVE;
     _intAv_TVE = psrc->_intAv_TVE;
     _intBd_TVE = psrc->_intBd_TVE;
@@ -140,6 +140,7 @@ IPNonLinearTVM& IPNonLinearTVM::operator=(const IPVariable& src)
     _Av_TVE_vector = psrc->_Av_TVE_vector;
     _Bd_TVE_vector = psrc->_Bd_TVE_vector;
     _dAv_dTrEe_TVE_vector = psrc->_dAv_dTrEe_TVE_vector;
+    _dBd_dTrEe_TVE_vector = psrc->_dBd_dTrEe_TVE_vector;
     _intAv_TVE_vector = psrc->_intAv_TVE_vector;
     _intBd_TVE_vector = psrc->_intBd_TVE_vector;
     _dBd_dDevEe_TVE_vector = psrc->_dBd_dDevEe_TVE_vector;
@@ -173,6 +174,8 @@ void IPNonLinearTVM::restart(){
   IPHyperViscoElastoPlastic::restart();
   restartManager::restart(_DA_DT);
   restartManager::restart(_DDA_DTT);
+  restartManager::restart(_DA_DtrE);
+  restartManager::restart(_DDA_DTDtrE);
   restartManager::restart(_DA_DdevE);
   restartManager::restart(_DDA_DTDdevE);
   restartManager::restart(_DB_DT);
@@ -207,6 +210,7 @@ void IPNonLinearTVM::restart(){
   restartManager::restart(_Av_TVE);
   restartManager::restart(_Bd_TVE);
   restartManager::restart(_dAv_dTrEe);
+  restartManager::restart(_dBd_dTrEe);
   restartManager::restart(_dAv_dTrEe_TVE);
   restartManager::restart(_intAv_TVE);
   restartManager::restart(_intBd_TVE);
@@ -219,6 +223,7 @@ void IPNonLinearTVM::restart(){
   restartManager::restart(_Av_TVE_vector);
   restartManager::restart(_Bd_TVE_vector);
   restartManager::restart(_dAv_dTrEe_TVE_vector);
+  restartManager::restart(_dBd_dTrEe_TVE_vector);
   restartManager::restart(_intAv_TVE_vector);
   restartManager::restart(_intBd_TVE_vector);
   restartManager::restart(_dBd_dDevEe_TVE_vector);
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVM.h b/NonLinearSolver/internalPoints/ipNonLinearTVM.h
index c784ceca1b5cf6c4512c8052bb49b2bf99f8c42d..feaafe8519e6b5765b624d7f7767e7d8e3fcc844 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVM.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVM.h
@@ -35,13 +35,13 @@ class IPNonLinearTVM : public IPHyperViscoElastoPlastic{
         STensor3 _devCorKirinf_TVE;
         STensor3 _corKirExtra;
         double _pressure;
-        double _dAv_dTrEe;
+        double _dAv_dTrEe, _dBd_dTrEe;
         double _Av_TVE, _Bd_TVE, _dAv_dTrEe_TVE;
         double _intAv_TVE, _intBd_TVE;
         double _intAv_1_TVE, _intAv_2_TVE;
         STensor3 _dBd_dDevEe;
         STensor3 _dBd_dDevEe_TVE, _intBd_1_TVE, _intBd_2_TVE;
-        std::vector<double> _Av_TVE_vector, _Bd_TVE_vector, _dAv_dTrEe_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector;
+        std::vector<double> _Av_TVE_vector, _Bd_TVE_vector, _dAv_dTrEe_TVE_vector, _dBd_dTrEe_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector;
         std::vector<STensor3> _dBd_dDevEe_TVE_vector; 
         
         // DcorKirDT
@@ -50,7 +50,7 @@ class IPNonLinearTVM : public IPHyperViscoElastoPlastic{
         STensor43 _DDcorKirDTDE;
         
         // Derivatives of A and B
-        std::vector<STensor3> _DA_DT, _DDA_DTT;     
+        std::vector<STensor3> _DA_DT, _DDA_DTT, _DA_DtrE, _DDA_DTDtrE;     
         std::vector<STensor43> _DA_DdevE, _DDA_DTDdevE;     
         std::vector<double> _DB_DT, _DDB_DTT, _DB_DtrE, _DDB_DTDtrE;                        
         
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index cc120251f1ebb889146863a13dff0e18bc287ae0..7b2aed65b54050d24c305eea06c3a23186d5ce78 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -10,6 +10,87 @@
 #include "STensorOperations.h"
 #include "nonLinearMechSolver.h"
 
+mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
+                          const bool matrixbyPerturbation, const double pert):
+    materialLaw(num,true), _rho(rho),_tangentByPerturbation(matrixbyPerturbation),_perturbationfactor(pert),
+    _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.), 
+    _volCorrection(0.), _xivolCorrection(1.), _zetavolCorrection(1.), _devCorrection(0.), _thetadevCorrection(1.), _pidevCorrection(1.), 
+    _compCorrection(1.), _compCorrection_2(0.), _compCorrection_3(0.), _tensionCompressionRegularisation(5.), 
+    _extraBranchType(Bilogarithmic), _extraBranchNLType(tanhType){
+        
+   _Idev = _I4;
+   STensor3 mIon3(-1./3);
+   STensor43 mIIon3;
+   tensprod(_I,mIon3, mIIon3);
+   _Idev += mIIon3;
+
+  double _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
+
+  _mu = 0.5*E/(1.+nu);
+  _K = E/3./(1.-2.*nu);
+
+  Cel=0.;
+  Cel(0,0,0,0) = _lambda + 2.*_mu;
+  Cel(1,1,0,0) = _lambda;
+  Cel(2,2,0,0) = _lambda;
+  Cel(0,0,1,1) = _lambda;
+  Cel(1,1,1,1) = _lambda + 2.*_mu;
+  Cel(2,2,1,1) = _lambda;
+  Cel(0,0,2,2) = _lambda;
+  Cel(1,1,2,2) = _lambda;
+  Cel(2,2,2,2) = _lambda + 2.*_mu;
+
+  Cel(1,0,1,0) = _mu;
+  Cel(2,0,2,0) = _mu;
+  Cel(0,1,0,1) = _mu;
+  Cel(2,1,2,1) = _mu;
+  Cel(0,2,0,2) = _mu;
+  Cel(1,2,1,2) = _mu;
+
+  Cel(0,1,1,0) = _mu;
+  Cel(0,2,2,0) = _mu;
+  Cel(1,0,0,1) = _mu;
+  Cel(1,2,2,1) = _mu;
+  Cel(2,0,0,2) = _mu;
+  Cel(2,1,1,2) = _mu;
+};
+mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): materialLaw(src),
+      _rho(src._rho),
+      _mu(src._mu),_K(src._K),Cel(src.Cel),
+      _order(src._order),
+      _tangentByPerturbation(src._tangentByPerturbation),_perturbationfactor(src._perturbationfactor),
+      _N(src._N),_viscoMethod(src._viscoMethod),_Ki(src._Ki),_ki(src._ki),_Gi(src._Gi),_gi(src._gi),
+      _I4(src._I4), _I(src._I), _Idev(src._Idev), _volCorrection(src._volCorrection), _xivolCorrection(src._xivolCorrection), 
+      _zetavolCorrection(src._zetavolCorrection), _devCorrection(src._devCorrection), _thetadevCorrection(src._thetadevCorrection), 
+      _pidevCorrection(src._pidevCorrection), _compCorrection(src._compCorrection), _compCorrection_2(src._compCorrection_2),  _compCorrection_3(src._compCorrection_3),  
+      _tensionCompressionRegularisation(src._tensionCompressionRegularisation), _extraBranchType(src._extraBranchType), _extraBranchNLType(src._extraBranchNLType)
+{
+
+};
+
+mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &source)
+{
+  materialLaw::operator=(source);
+  const mlawHyperViscoElastic* src =dynamic_cast<const mlawHyperViscoElastic*>(&source);
+  if(src != NULL)
+  {
+      _rho = src->_rho,
+      _mu = src->_mu; _K = src->_K; Cel = src->Cel;
+      _order = src->_order;
+      _tangentByPerturbation = src->_tangentByPerturbation; _perturbationfactor = src->_perturbationfactor;
+      _N = src->_N; _viscoMethod = src->_viscoMethod; _Ki = src->_Ki; _ki = src->_ki; _Gi = src->_Gi; _gi = src->_gi;
+      _zetavolCorrection=src->_zetavolCorrection; _devCorrection=src->_devCorrection;
+       _thetadevCorrection=src->_thetadevCorrection; _pidevCorrection=src->_pidevCorrection;
+       _extraBranchType = src->_extraBranchType;
+       _compCorrection = src->_compCorrection;
+       _compCorrection_2 = src->_compCorrection_2;
+       _compCorrection_3 = src->_compCorrection_3;
+       _extraBranchNLType = src->_extraBranchNLType;
+       _tensionCompressionRegularisation = src->_tensionCompressionRegularisation;
+  }
+  return *this;
+}
+
 void mlawHyperViscoElastic::setStrainOrder(const int order){
   _order = order;
   Msg::Info("order = %d is used to approximate log and exp");
@@ -108,8 +189,8 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
   fclose(file);
 };
 
-void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &intA, double &B_d, STensor3 &dB_vddev, double &intB) const
-{
+void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &intA, double &B_d, STensor3 &dB_vddev, double &intB, 
+                                                    double* dB_dTrEe) const{
   if(_extraBranchNLType == expType)
   {
     A_v = getVolumeCorrection()*pow(exp(getXiVolumeCorrection()/3.*tr*tr)-1.,getZetaVolumeCorrection());
@@ -197,7 +278,8 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
         double x = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
         A_v = -1. + getVolumeCorrection()* 1/sqrt(x) ;
         if (tr < 0.){
-            A_v = -1. + getVolumeCorrection()* 1/sqrt(x) * (1 + _compCorrection/getVolumeCorrection()) ;
+             A_v = -1. + getVolumeCorrection()* 1/sqrt(x) * (1. + _compCorrection/getVolumeCorrection()) ;
+            // A_v += (_compCorrection*tanh(_compCorrection_2*tr*tr+_compCorrection_3)) ;
         }
         if(A_v > -1.){    
             dA_vdE = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5);
@@ -205,9 +287,9 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
             intA -= getVolumeCorrection()*sqrt(getZetaVolumeCorrection())/getXiVolumeCorrection(); // value at trEe = 0.
             
             if (tr < 0.){
-                dA_vdE = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5) * (1 + _compCorrection/getVolumeCorrection());
-                intA = - tr*tr/2. + getVolumeCorrection()*sqrt(x)/getXiVolumeCorrection() * (1 + _compCorrection/getVolumeCorrection()); // integral of A_v * trEe
-                intA -= getVolumeCorrection()*sqrt(getZetaVolumeCorrection())/getXiVolumeCorrection() * (1 + _compCorrection/getVolumeCorrection()); // value at trEe = 0.
+                dA_vdE = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5) * (1. + _compCorrection/getVolumeCorrection());
+                intA = - tr*tr/2. + getVolumeCorrection()*sqrt(x)/getXiVolumeCorrection() * (1. + _compCorrection/getVolumeCorrection()); // integral of A_v * trEe
+                intA -= getVolumeCorrection()*sqrt(getZetaVolumeCorrection())/getXiVolumeCorrection() * (1. + _compCorrection/getVolumeCorrection()); // value at trEe = 0.
             }
         }
         else{
@@ -218,6 +300,9 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
         
         double y = getThetaDevCorrection()*dev.dotprod() + getPiDevCorrection();
         B_d = -1. + getDevCorrection()*1/sqrt(y);
+        if (tr < 0.){
+            // B_d += (_compCorrection*tanh(_compCorrection_2*tr*tr+_compCorrection_3)) ;
+        }
         
         if(B_d > -1.){ 
             STensorOperation::zero(dB_vddev);
@@ -232,9 +317,50 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
             intB = 0.;
         }
 
-        // Msg::Error(" Inside evaluatePhiPCorrection, A_v = %e, B_d = %e !!", A_v, B_d);
+    }
+    else if(_extraBranchNLType == TensionCompressionRegularisedType){
+        
+        // Regularising function
+        double m = _tensionCompressionRegularisation;
+        double expmtr(0.);
+        if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
+        double sigmoid = 1/(1.+expmtr);
+        
+        double x = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
+        A_v = -1. + sigmoid*1./sqrt(x) + (1.-sigmoid)*(_compCorrection/sqrt(x));
+  
+        dA_vdE = - sigmoid*getXiVolumeCorrection()*tr/pow(x,1.5) - (1.-sigmoid)*(getXiVolumeCorrection()*_compCorrection*tr/pow(x,1.5)) 
+                   + (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x) - (_compCorrection*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x);
+            
+        if(getXiVolumeCorrection()>0.){
+            double integrand = 1./getXiVolumeCorrection()*sqrt(x); // integral of A_v * trEe
+            integrand -= ( 1./getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ); // value at trEe = 0.
+            intA = - tr*tr/2. + sigmoid*integrand + (1.-sigmoid)*integrand;
+        }
+        else{
+            intA = - tr*tr/2.;
+        }
+            
+        double y = getThetaDevCorrection()*dev.dotprod() + getPiDevCorrection();
+        B_d = -1. + sigmoid*1./sqrt(y) + (1.-sigmoid)*(_compCorrection/sqrt(y)); 
 
-    
+        STensorOperation::zero(dB_vddev);
+        dB_vddev = dev;
+        dB_vddev *= (-sigmoid*getThetaDevCorrection()/pow(y,1.5) - (1.-sigmoid)*(getThetaDevCorrection()*_compCorrection/pow(y,1.5)));
+        
+        if(dB_dTrEe!=NULL){
+            *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y) - (_compCorrection*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y);
+        }    
+        if(getThetaDevCorrection()>0.){
+            double integrand = 1./getThetaDevCorrection()*sqrt(y);  // integral of B_d * devEe
+            integrand -= (1./getThetaDevCorrection()*sqrt(getPiDevCorrection()) );  // value at devEe = 0.   
+            intB = -dev.dotprod()/2. + sigmoid*integrand + (1.-sigmoid)*integrand;
+        }
+        else{
+            intB = -dev.dotprod()/2.;
+        } 
+        
+        // Msg::Error(" Inside evaluatePhiPCorrection, A_v = %e, B_d = %e !!", A_v, B_d);
     }
     else{
       Msg::Error("_extraBranchNLType %d has not been defined",_extraBranchNLType);
@@ -379,83 +505,6 @@ void mlawHyperViscoElastic::dViscousEnergydF(const IPHyperViscoElastic& q0, cons
   }
 }
 
-
-mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const double nu, const double rho,
-                          const bool matrixbyPerturbation, const double pert):
-    materialLaw(num,true), _rho(rho),_tangentByPerturbation(matrixbyPerturbation),_perturbationfactor(pert),
-    _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.), 
-    _volCorrection(0.), _xivolCorrection(1.), _zetavolCorrection(1.), _devCorrection(0.), _thetadevCorrection(1.), _pidevCorrection(1.), _compCorrection(0.), 
-    _extraBranchType(Bilogarithmic), _extraBranchNLType(tanhType)
-{
-   _Idev = _I4;
-   STensor3 mIon3(-1./3);
-   STensor43 mIIon3;
-   tensprod(_I,mIon3, mIIon3);
-   _Idev += mIIon3;
-
-  double _lambda = (E*nu)/(1.+nu)/(1.-2.*nu);
-
-  _mu = 0.5*E/(1.+nu);
-  _K = E/3./(1.-2.*nu);
-
-  Cel=0.;
-  Cel(0,0,0,0) = _lambda + 2.*_mu;
-  Cel(1,1,0,0) = _lambda;
-  Cel(2,2,0,0) = _lambda;
-  Cel(0,0,1,1) = _lambda;
-  Cel(1,1,1,1) = _lambda + 2.*_mu;
-  Cel(2,2,1,1) = _lambda;
-  Cel(0,0,2,2) = _lambda;
-  Cel(1,1,2,2) = _lambda;
-  Cel(2,2,2,2) = _lambda + 2.*_mu;
-
-  Cel(1,0,1,0) = _mu;
-  Cel(2,0,2,0) = _mu;
-  Cel(0,1,0,1) = _mu;
-  Cel(2,1,2,1) = _mu;
-  Cel(0,2,0,2) = _mu;
-  Cel(1,2,1,2) = _mu;
-
-  Cel(0,1,1,0) = _mu;
-  Cel(0,2,2,0) = _mu;
-  Cel(1,0,0,1) = _mu;
-  Cel(1,2,2,1) = _mu;
-  Cel(2,0,0,2) = _mu;
-  Cel(2,1,1,2) = _mu;
-};
-mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src): materialLaw(src),
-      _rho(src._rho),
-      _mu(src._mu),_K(src._K),Cel(src.Cel),
-      _order(src._order),
-      _tangentByPerturbation(src._tangentByPerturbation),_perturbationfactor(src._perturbationfactor),
-      _N(src._N),_viscoMethod(src._viscoMethod),_Ki(src._Ki),_ki(src._ki),_Gi(src._Gi),_gi(src._gi),
-      _I4(src._I4), _I(src._I), _Idev(src._Idev), _volCorrection(src._volCorrection), _xivolCorrection(src._xivolCorrection), 
-      _zetavolCorrection(src._zetavolCorrection), _devCorrection(src._devCorrection), _thetadevCorrection(src._thetadevCorrection), 
-      _pidevCorrection(src._pidevCorrection), _compCorrection(src._compCorrection), _extraBranchType(src._extraBranchType), _extraBranchNLType(src._extraBranchNLType)
-{
-
-};
-
-mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &source)
-{
-  materialLaw::operator=(source);
-  const mlawHyperViscoElastic* src =dynamic_cast<const mlawHyperViscoElastic*>(&source);
-  if(src != NULL)
-  {
-      _rho = src->_rho,
-      _mu = src->_mu; _K = src->_K; Cel = src->Cel;
-      _order = src->_order;
-      _tangentByPerturbation = src->_tangentByPerturbation; _perturbationfactor = src->_perturbationfactor;
-      _N = src->_N; _viscoMethod = src->_viscoMethod; _Ki = src->_Ki; _ki = src->_ki; _Gi = src->_Gi; _gi = src->_gi;
-      _zetavolCorrection=src->_zetavolCorrection; _devCorrection=src->_devCorrection;
-       _thetadevCorrection=src->_thetadevCorrection; _pidevCorrection=src->_pidevCorrection;
-       _extraBranchType = src->_extraBranchType;
-       _compCorrection = src->_compCorrection;
-       _extraBranchNLType = src->_extraBranchNLType;
-  }
-  return *this;
-}
-
 double mlawHyperViscoElastic::soundSpeed() const
 {
   double _E=9*_K*_mu/(3*_K+_mu);
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index a7d0aaa2d6db2915abd60f3a683a31b0ea2d2970..c0fe1a8a337dd8a3c7b7ea74188c7ba3cd48eee0 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.h
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.h
@@ -22,7 +22,7 @@ class mlawHyperViscoElastic : public materialLaw{
   public:
     enum viscoelasticType{Maxwell=0,KelvinVoight=1};
     enum extraBranchType{Bilogarithmic=0};
-    enum extraBranchNLType{expType=0,tanhType=1,sigmoidType=2,hyperType=3};
+    enum extraBranchNLType{expType=0,tanhType=1,sigmoidType=2,hyperType=3,TensionCompressionRegularisedType=4};
   #ifndef SWIG
   protected:
     STensor43 _I4; // symetric 4th order unit tensor
@@ -31,14 +31,6 @@ class mlawHyperViscoElastic : public materialLaw{
     double _rho; // density
     double _mu; // 2nd lame parameter (=G)
     double _K; // bulk modulus
-    double _volCorrection; //Vc correction on volume potential
-    double _xivolCorrection; //Xi correction on volume potential
-    double _zetavolCorrection; //Zeta correction on volume potential
-    
-    double _devCorrection; //Dc correction on volume potential
-    double _thetadevCorrection; //Theta correction on volume potential
-    double _pidevCorrection; //pi correction on volume potential
-    double _compCorrection; // compression correction on volume potential
 
     STensor43 Cel; // elastic hook tensor
     int _order ; // to approximated log and exp of a tensor
@@ -54,6 +46,18 @@ class mlawHyperViscoElastic : public materialLaw{
     std::vector<double> _Gi;
     std::vector<double> _gi;
     
+    // extraBranch parameters
+    double _volCorrection; //Vc correction on volume potential
+    double _xivolCorrection; //Xi correction on volume potential
+    double _zetavolCorrection; //Zeta correction on volume potential
+    
+    double _devCorrection; //Dc correction on volume potential
+    double _thetadevCorrection; //Theta correction on volume potential
+    double _pidevCorrection; //pi correction on volume potential
+    double _compCorrection; // compression correction 
+    double _compCorrection_2; // compression correction 
+    double _compCorrection_3; // compression correction 
+    double _tensionCompressionRegularisation; // regularisation parameter for tension-compression asymmetry in extrabranch
     extraBranchType _extraBranchType;
     extraBranchNLType _extraBranchNLType;
 
@@ -75,7 +79,8 @@ class mlawHyperViscoElastic : public materialLaw{
     void predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0_, IPHyperViscoElastic *q1,
                                   const bool stiff, STensor43& Tangent) const;
 
-    void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr,  double &A_v, double &dAprdEepr, double &intA, double &B_d, STensor3 &dB_vddevEe, double &intB) const;
+    void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr,  double &A_v, double &dAprdEepr, double &intA, double &B_d, STensor3 &dB_vddevEe, double &intB, 
+                                double* dB_dTrEe = NULL) const;
     //void evaluateA_dA(double trEe, double &A_v, double &dAprdEepr) const;
     
   #endif // SWIG
@@ -99,7 +104,9 @@ class mlawHyperViscoElastic : public materialLaw{
     virtual double getDevCorrection() const {return _devCorrection;};
     virtual double getThetaDevCorrection() const {return _thetadevCorrection;};
     virtual double getPiDevCorrection() const {return _pidevCorrection;};
-    virtual void setExtraBranch_CompressionParameter(double compCorrection) {_compCorrection = compCorrection;};
+    virtual void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3) {
+                                    _compCorrection = compCorrection; _compCorrection_2 = compCorrection_2; _compCorrection_3 = compCorrection_3;};
+    virtual void setTensionCompressionRegularisation(double tensionCompressionRegularisation) {_tensionCompressionRegularisation = tensionCompressionRegularisation;};
     virtual void setExtraBranchType(const int type);
     virtual void setExtraBranchNLType(const int type);
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
index fd7328de3c08f18217f41729367fc3b2529c3c91..b94b227aba7ab73a8ff2bedd5cfc2ad9712c3ce0 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -381,10 +381,18 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   
   // update A, B -> ExtraBranch calculations are within
   double Ke(0.), Ge(0.), Ke_pr(0.), Ge_pr(0.);
-  static STensor43 Ge_Tensor_pr, Ge_Tensor, Bd_stiffnessTerm;
-  STensorOperation::zero(Ge_Tensor_pr);
-  STensorOperation::zero(Ge_Tensor);
-  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
+  static STensor43 Ge_Tensor_pr, Ge_Tensor, Bd_stiffnessTerm, Ge_TrEeTensor_pr, Ge_TrEeTensor;
+  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){
+    mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,stiff,Bd_stiffnessTerm);
+  }
+  else{ // 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);
+    Ge_TrEeTensor_pr = Ge_TrEeTensor;
+  }
+  
   Ge_Tensor = _I4;
   Ge_Tensor *= Ge*2; // Because the function does not do it
   Ge_Tensor += Bd_stiffnessTerm;
@@ -465,7 +473,12 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   devPhi = devPhipr;
   
   // Initialise the rest
-  getDho(Gamma,Cepr,Ceinvpr,Kepr,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
+  if (_extraBranchNLType != TensionCompressionRegularisedType && _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);
+  }
+  else{ // 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);
+  }
 
   double PhiEqpr2 = 1.5*devPhi.dotprod();
   double PhiEqpr = sqrt(PhiEqpr2);      // -> second invariant
@@ -524,9 +537,9 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                         DdevPhidGamma_RHS(i,j) += ( -0.5*2*_b*Ke*(-ptilde-Gamma*dPhiPdGamma)*_I(i,j) - 3*pow(kk,1)*Hb/Cxdev*devPhi(i,j) ); // pow(kk,2) DEBUG
                         for (int k=0; k<3; k++)
                             for (int l=0; l<3; l++){
-                                DdevPhidDgamma_RHS(i,j) += 0.5* 2.*_b/3.*Gamma*Dho(i,j,k,l)*_I(k,l)*Hp; 
-                                DdevPhidGamma_RHS(i,j) += 0.5* ( -3.*Ge_Tensor(i,j,k,l)*devPhi(k,l) + 2.*_b/3.*Gamma*Dho(i,j,k,l)*_I(k,l)*dPhiPdGamma ); 
-                                DdevPhidGamma_RHS(i,j) += 0.5*Dho(i,j,k,l)*N(k,l); // DEBUG
+                                DdevPhidDgamma_RHS(i,j) += 0.5* 2.*_b/3.*Gamma*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*_I(k,l)*Hp; 
+                                DdevPhidGamma_RHS(i,j) += 0.5* ( -3.*Ge_Tensor(i,j,k,l)*devPhi(k,l) + 2.*_b/3.*Gamma*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*_I(k,l)*dPhiPdGamma ); 
+                                DdevPhidGamma_RHS(i,j) += 0.5*(-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*N(k,l); // DEBUG
                             }
                     }
                 for (int i=0; i<3; i++)
@@ -610,8 +623,15 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                 }
                 
                 // Update Phi
-                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);
+                if (_extraBranchNLType != TensionCompressionRegularisedType && _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);
+                }
+                else{ // 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);
+                }
+                                
                 PhiEq = sqrt(1.5*devPhi.dotprod());
                 
                 // Update A
@@ -655,8 +675,14 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
             }
 
             // Correct Phi
-            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);
+            if (_extraBranchNLType != TensionCompressionRegularisedType && _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);
+            }
+            else{ // 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);
+            }
             PhiEq = sqrt(1.5*devPhi.dotprod());
     
             // Correct Normal, H = expGN
@@ -827,6 +853,8 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                                 DcorKirDEe(i,j,k,l) += Ge_Tensor(i,j,p,q)*_Idev(p,q,k,l);
                             }
                     }
+        DcorKirDEe_pr += Ge_TrEeTensor_pr; // TC Assymmetry
+        DcorKirDEe += Ge_TrEeTensor; // TC Assymmetry
         
         static STensor43 DcorKir_pr_DCepr, DcorKirDCepr;
         STensorOperation::multSTensor43(DcorKirDEe_pr,DlnDCepr,DcorKir_pr_DCepr);
@@ -940,7 +968,8 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                     for (int p=0; p<3; p++)
                       for (int q=0; q<3; q++){
                           
-                        DdevphiDCepr_RHS(i,j,k,l) +=  0.5* ( Dho(i,j,p,q)* ( 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) );
+                        DdevphiDCepr_RHS(i,j,k,l) +=  0.5* ( Ge_TrEeTensor(i,j,p,q)*0.5*DlnDCepr(p,q,k,l) + 
+                                        (-Ge_TrEeTensor(i,j,p,q) + Dho(i,j,p,q))* ( 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) );
                         
                         for (int r=0; r<3; r++)
                           for (int s=0; s<3; s++){
@@ -1015,8 +1044,9 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                     for (int p=0; p<3; p++)
                       for (int q=0; q<3; q++){
                           
-                        DdevphiDCepr_RHS(i,j,k,l) +=  (0.5* ( - 3.*Ge_Tensor(i,j,p,q)*devPhi(p,q)*DGDCepr(k,l) +
-                                                    Dho(i,j,p,q)* ( N(p,q)*DGDCepr(k,l) + 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) ));
+                        DdevphiDCepr_RHS(i,j,k,l) +=  (0.5* ( - 3.*Ge_Tensor(i,j,p,q)*devPhi(p,q)*DGDCepr(k,l) 
+                                                    + (-Ge_TrEeTensor(i,j,p,q) + Dho(i,j,p,q)) * ( N(p,q)*DGDCepr(k,l) + 2.*_b*Gamma/3.*_I(p,q)*DphiPDCepr(k,l) ) 
+                                                    +  Ge_TrEeTensor(i,j,p,q) *0.5*DlnDCepr(p,q,k,l) ) );
                         for (int r=0; r<3; r++)
                           for (int s=0; s<3; s++){
                               
@@ -1332,7 +1362,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                 for (int k=0; k<3; k++)
                   for (int l=0; l<3; l++){
                     DdevphiDT_RHS(i,j) += 0.5*( -3*dGammaDT*Ge_Tensor(i,j,k,l)*devPhi(k,l) 
-                                                + Dho(i,j,k,l)*( dGammaDT*N(k,l) + 2.*_b/3.*dPhiPdT*_I(k,l) ) 
+                                                + (-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*( dGammaDT*N(k,l) + 2.*_b/3.*dPhiPdT*_I(k,l) ) 
                                                 + Ce(i,k)*DcorKirDT(k,l)*Ceinv(l,j) );  
                   }
               }
@@ -1388,7 +1418,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                 for (int k=0; k<3; k++)
                   for (int l=0; l<3; l++){
                     DdevphiDT_RHS(i,j) += 0.5*( -3*dGammaDT*Ge_Tensor(i,j,k,l)*devPhi(k,l) 
-                                                + Dho(i,j,k,l)*( dGammaDT*N(k,l) + 2.*_b/3.*Gamma*dPhiPdT*_I(k,l) ) 
+                                                + (-Ge_TrEeTensor(i,j,k,l) + Dho(i,j,k,l))*( dGammaDT*N(k,l) + 2.*_b/3.*Gamma*dPhiPdT*_I(k,l) ) 
                                                 + Ce(i,k)*DcorKirDT(k,l)*Ceinv(l,j) );  
                   }
               }
@@ -1768,7 +1798,10 @@ void mlawNonLinearTVENonLinearTVP::get_G1_Tensor(const STensor3& Cepr, const STe
 void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS,
                                   const double& Ke, const STensor43& Ge_Tensor, 
                                   const double& kk, const double& Hb, const double& Cxtr, const double& Cxdev,
-                                  const STensor3& expGN, const STensor43& dexpAdA, STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv) const{
+                                  const STensor3& expGN, const STensor43& dexpAdA, 
+                                  STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, STensor43* Ge_TrEeTensor) const{
+
+  // Ge_TrEeTensor -> Dependence of devCorKir on TrEe. DO NOT CHANGE ITS VALUE!
 
   // Get H = expGN is a symmetric tensor, HT, Hinv
   static STensor3 HT, Hinv, HTinv; 
@@ -1861,6 +1894,9 @@ void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& C
                 for (int q=0; q<3; q++)
                     DcorKirDEe(i,j,k,l) += Ge_Tensor(i,j,p,q)*_Idev(p,q,k,l);
         }
+ if (Ge_TrEeTensor!= NULL){
+    DcorKirDEe += (*Ge_TrEeTensor);
+ }
  
   // Dho -> 2nd term
   for (int i=0; i<3; i++)
@@ -1924,12 +1960,19 @@ void mlawNonLinearTVENonLinearTVP::getDho(const double& Gamma, const STensor3& C
     for (int j=0; j<3; j++)
       for (int k=0; k<3; k++)
         for (int l=0; l<3; l++){
-          Dho2_u(i,j,k,l) += 3*Gamma*pow(kk,1)*Hb/Cxdev * _I4(i,j,k,l); // DEBUG
+          Dho2_u(i,j,k,l) += 3*Gamma*pow(kk,1)*Hb/Cxdev * _I4(i,j,k,l); // DEBUG pow(kk,2)
           for (int p=0; p<3; p++)
             for (int q=0; q<3; q++){
                 Dho2_u(i,j,k,l) += 1.5*Gamma*(Ge_Tensor(i,j,p,q)*_I4(p,q,k,l) - Dho(i,j,p,q)*_I4(p,q,k,l));
         }
       }
+  if (Ge_TrEeTensor!= NULL){
+    static STensor43 Ge_TrEeTensor_temp;
+    STensorOperation::zero(Ge_TrEeTensor_temp);
+    Ge_TrEeTensor_temp = (*Ge_TrEeTensor);
+    Ge_TrEeTensor_temp *= (1.5*Gamma);
+    Dho2_u += Ge_TrEeTensor_temp;
+  }
   STensorOperation::inverseSTensor43(Dho2_u,Dho2_u_inv);
   
   // Get Dho2_v, Dho2_v_inv
@@ -1944,7 +1987,8 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
                                             double& Ke, double& Ge, STensor43& Ge_Tensor, 
                                             double& ptilde, STensor3& devPhi,
                                             STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
-                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv) const{
+                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv,
+                                            STensor43* Ge_TrEeTensor) const{
 
   // This function calculates Phi iteratively. 
     // At every iteration it will calculate Ee, corKir and N
@@ -1978,7 +2022,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
   // Initialise corKir
   static STensor3 KS;
   static STensor43 Bd_stiffnessTerm;
-  mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
+  if (Ge_TrEeTensor == NULL){
+    mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
+  }
+  else{ // TC asymmetry
+    mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
+  }
   Ge_Tensor = _I4;
   Ge_Tensor *= Ge*2; // *2 because the function doesnt do it
   Ge_Tensor += Bd_stiffnessTerm;
@@ -2023,8 +2072,13 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
   }
 
   // Initialise Dho, Dho2, Dho2inv
-  getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
-
+  if (Ge_TrEeTensor == NULL){
+    getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
+  }
+  else{ // TC asymmetry
+    getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,Ge_TrEeTensor);
+  }
+  
   // Initialise DPhi
   static STensor3 DPhi;
 
@@ -2091,7 +2145,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
           Ee(i,j) = Eepr(i,j) - Gamma*N(i,j);
   
       // update corKir
-      mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,false,Bd_stiffnessTerm);
+      if (Ge_TrEeTensor == NULL){
+        mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm);
+      }
+      else{ // TC asymmetry
+        mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(Ee,q0->_Ee,q0,q1,Ke,Ge,T0,T,true,Bd_stiffnessTerm,Ge_TrEeTensor);
+      }  
       Ge_Tensor = _I4;
       Ge_Tensor *= Ge*2; // *2 because the function doesnt do it
       Ge_Tensor += Bd_stiffnessTerm;
@@ -2118,7 +2177,12 @@ void mlawNonLinearTVENonLinearTVP::getIterated_DPhi(const double& T0, const doub
           J(i,j) += ( Phi(i,j) - MS(i,j) + q1->_backsig(i,j) );
 
       // update Dho2, Dho2inv
-      getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
+      if (Ge_TrEeTensor == NULL){
+        getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv);
+      }
+      else{ // TC asymmetry
+        getDho(Gamma,Cepr,Ceinvpr,KS,Ke,Ge_Tensor,kk,Hb,Cxtr,Cxdev,expGN,dexpAdA,Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,Ge_TrEeTensor);
+      }
 
       // update J_tol
       /*
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
index 4b1e5980d58d323e6f49f35dee4a4db70dd869f3..0b2fc4f1b25d887cb9f0beb8e7232f9b44c36473 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.h
@@ -27,12 +27,14 @@ class mlawNonLinearTVENonLinearTVP : public mlawNonLinearTVP{
                                             double& Ke, double& Ge, STensor43& Ge_Tensor, 
                                             double& ptilde, STensor3& devPhi,
                                             STensor3& Phi, STensor3& N, STensor3& expGN, STensor43& dexpAdA,
-                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv) const;
+                                            STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv , double& Dho2_v_inv,
+                                            STensor43* Ge_TrEeTensor = NULL) const;
                                             
         virtual void getDho(const double& Gamma, const STensor3& Cepr, const STensor3& Ceinvpr, const STensor3& KS,
                                   const double& Ke, const STensor43& Ge_Tensor, 
                                   const double& kk, const double& Hb, const double& Cxtr, const double& Cxdev,
-                                  const STensor3& expGN, const STensor43& dexpAdA, STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv) const;
+                                  const STensor3& expGN, const STensor43& dexpAdA, 
+                                  STensor43& Dho, STensor43& Dho2inv, STensor43& Dho2_u_inv, double& Dho2_v_inv, STensor43* Ge_TrEeTensor = NULL) const;
 
         virtual void get_G1_Tensor(const STensor3& Cepr, const STensor3& expGN, const STensor43& DCeinvprDCepr, const STensor3& KS, STensor43& G1) const;
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index 7356f3b1b312482a11a8be03f5f9ae6d46ac24ef..800287d6854a190e786bc69cdc71b2a18b336bfe 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -406,8 +406,9 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                                                                 + STensorOperation::doubledot(q0._A[i],q0._A[i])/pow(1/q0._Bd_TVE,2));
                 
                     // 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){ 
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){ 
                     Psy_mech += _Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i];
                 }
             }
@@ -427,9 +428,29 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                 
         // extraBranch terms
         if (_useExtraBranch){
-            *DpsiDT += KT*q._intA*_temFunc_elasticCorrection_Bulk->getDiff(T) + dKdT*q._intA + 
-                        2.*GT*q._intB*_temFunc_elasticCorrection_Shear->getDiff(T) + 2.*dGdT*q._intB;
-            (*DpsiDE) += q._corKirExtra;
+            if(_extraBranchNLType != TensionCompressionRegularisedType){
+                *DpsiDT += KT*q._intA*_temFunc_elasticCorrection_Bulk->getDiff(T) + dKdT*q._intA + 
+                            2.*GT*q._intB*_temFunc_elasticCorrection_Shear->getDiff(T) + 2.*dGdT*q._intB;
+                (*DpsiDE) += q._corKirExtra;
+            }
+            else{
+                // Regularising function
+                double m = _tensionCompressionRegularisation;
+                double expmtr(0.);
+                if (exp(-m*trEe)<1.e+10){ expmtr = exp(-m*trEe);}
+                double sigmoid = 1/(1.+expmtr);
+                double CTasymm = sigmoid + _compCorrection*(1.-sigmoid);
+                double dCTasymmDtrE = m*expmtr/pow((1.+expmtr),2.) - _compCorrection*m*expmtr/pow((1.+expmtr),2.);
+                
+                *DpsiDT += (KT*q._intA*_temFunc_elasticCorrection_Bulk->getDiff(T) + dKdT*q._intA + 
+                            2.*GT*q._intB*_temFunc_elasticCorrection_Shear->getDiff(T) + 2.*dGdT*q._intB)*CTasymm;
+                (*DpsiDE) += q._corKirExtra;
+                
+                double psy_inf = (KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe) + KT*q._intA + 2.*GT*q._intB);
+                for (int j=0; j<3; j++)
+                    for (int k=0; k<3; k++)
+                        (*DpsiDE)(j,k) += psy_inf*dCTasymmDtrE*_I(j,k);
+            }
         }
         
         // viscoelastic terms
@@ -498,9 +519,10 @@ 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){ 
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){ 
                     
                     double dTrEei_DT(0.); // branchElasticStrain trace
                     static STensor3 dDevEei_DT;
@@ -508,16 +530,43 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                     for (int k=0; k<3; k++)
                         for (int l=0; l<3; l++)
                             dDevEei_DT(k,l) = q0._A[i](k,l) * Dexp_rec_g_DT + Dexp_mid_g_DT*q._devDE(k,l);
+                    
+                    if (_ExtraBranch_TVE_option == 2){
+                        *DpsiDT += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*dTrEei_DT + _Gi[i]*2.*q._Bd_TVE_vector[i]*STensorOperation::doubledot(q._A[i],dDevEei_DT);
+                        // STensorOperation::multSTensor43(q._DA_DdevE[i],_Idev,temp);
+                        for (int j=0; j<3; j++)
+                            for (int k=0; k<3; k++){
+                                (*DpsiDE)(j,k) += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*exp_mid_k*_I(j,k);
+                                for (int p=0; p<3; p++)
+                                    for (int r=0; r<3; r++)
+                                        (*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
                             
-                    *DpsiDT += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*dTrEei_DT + _Gi[i]*2.*q._Bd_TVE_vector[i]*STensorOperation::doubledot(q._A[i],dDevEei_DT);
-                    STensorOperation::multSTensor43(q._DA_DdevE[i],_Idev,temp);
-                    for (int j=0; j<3; j++)
-                        for (int k=0; k<3; k++){
-                            (*DpsiDE)(j,k) += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*exp_mid_k*_I(j,k);
-                            for (int p=0; p<3; p++)
-                                for (int r=0; r<3; r++)
-                                    (*DpsiDE)(j,k) += 2*_Gi[i]*q._Bd_TVE_vector[i]*q._A[i](p,r)*exp_mid_g*_Idev(p,r,j,k) ; 
+                        // Regularising function
+                        double m = _tensionCompressionRegularisation;
+                        double expmtr(0.);
+                        if (exp(-m*q._B[i])<1.e+10){ expmtr = exp(-m*q._B[i]);}
+                        double sigmoid = 1/(1.+expmtr);
+                        double CTasymm = sigmoid + _compCorrection*(1.-sigmoid);
+                        double dCTasymmDtrE = m*expmtr/pow((1.+expmtr),2.) - _compCorrection*m*expmtr/pow((1.+expmtr),2.);
+                        
+                        (*DpsiDT) += (_Ki[i]*q._Av_TVE_vector[i]*q._B[i]*dTrEei_DT + 
+                                    _Gi[i]*2.*q._Bd_TVE_vector[i]*STensorOperation::doubledot(q._A[i],dDevEei_DT)) * CTasymm;
+                                    
+                        (*DpsiDT) += (_Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i]) * dCTasymmDtrE * dTrEei_DT;
+                                    
+                        for (int j=0; j<3; j++)
+                            for (int k=0; k<3; k++){
+                                (*DpsiDE)(j,k) += _Ki[i]*q._Av_TVE_vector[i]*q._B[i]*exp_mid_k*_I(j,k) * CTasymm;
+                                (*DpsiDE)(j,k) += (_Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i]) * dCTasymmDtrE * exp_mid_k*_I(j,k);
+                                for (int p=0; p<3; p++)
+                                    for (int r=0; r<3; r++)
+                                        (*DpsiDE)(j,k) += 2*_Gi[i]*q._Bd_TVE_vector[i]*q._A[i](p,r)*exp_mid_g*_Idev(p,r,j,k) * CTasymm; 
                         }
+                        
+                    }
                 }
               }
 
@@ -567,9 +616,19 @@ void mlawNonLinearTVM::corKirInfinity(const IPNonLinearTVM *q1, const STensor3&
 void mlawNonLinearTVM::evaluateElasticCorrection(const double trE, const STensor3 &devE, const double& T, 
                                                  double &A_v, double &dA_vdE, double &intA, double &dAdT,
                                                  double &B_d, STensor3 &dB_vddev, double &intB, double &dBdT,
-                                                 double *ddAdTT, double *ddBdTT, double *ddA_vdTdE, STensor3 *ddB_vdTdevE) const{
+                                                 double *ddAdTT, double *ddBdTT, double *ddA_vdTdE, STensor3 *ddB_vdTdevE,
+                                                 double *dB_dTrEe, double *dB_dTdTrEe) const{
     
-    mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB);
+    if (_extraBranchNLType != TensionCompressionRegularisedType){ 
+        mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB);
+    }
+    else{
+        mlawHyperViscoElastic::evaluatePhiPCorrection(trE,devE,A_v,dA_vdE,intA,B_d,dB_vddev,intB,dB_dTrEe);
+        *dB_dTrEe *= _temFunc_elasticCorrection_Shear->getVal(T);
+        if(dB_dTdTrEe!=NULL){
+            *dB_dTdTrEe = (*dB_dTrEe)*_temFunc_elasticCorrection_Bulk->getDiff(T);
+        }
+    }
     
     A_v *= _temFunc_elasticCorrection_Bulk->getVal(T);
     dA_vdE *= _temFunc_elasticCorrection_Bulk->getVal(T);
@@ -601,7 +660,8 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
                         double* DsigV_dTrEe, STensor43* DsigD_dDevEe,
                         double* DsigV_dT, STensor3* DsigD_dT,
                         double* DDsigV_dTT, STensor3* DDsigD_dTT,
-                        double* DDsigV_dTdTrEe, STensor43* DDsigD_dTdDevEe) const
+                        double* DDsigV_dTdTrEe, STensor43* DDsigD_dTdDevEe,
+                        STensor3* DsigD_dTrEe, STensor3* DsigD_dTdTrEe) const
 {
     // DsigDEe = DcorKirDE
     
@@ -612,10 +672,16 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
     STensorOperation::decomposeDevTr(Ee,devEe,trEe);
     
     double A(0.), B(0.);
-    double dA_dTrEe(0.), dA_dTdTrEe(0.);
+    double dA_dTrEe(0.), dA_dTdTrEe(0.), dB_dTrEe(0.), dB_dTdTrEe(0.);
     STensor3 dB_dDevEe, ddB_dTdDevEe;
     double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.), ddA_dTT(0.), ddB_dTT(0.);
-    evaluateElasticCorrection(trEe, devEe, T, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, &ddA_dTT, &ddB_dTT, &dA_dTdTrEe, &ddB_dTdDevEe);
+    
+    if (_extraBranchNLType != TensionCompressionRegularisedType){
+        evaluateElasticCorrection(trEe, devEe, T, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, &ddA_dTT, &ddB_dTT, &dA_dTdTrEe, &ddB_dTdDevEe);
+    }
+    else {
+        evaluateElasticCorrection(trEe, devEe, T, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, &ddA_dTT, &ddB_dTT, &dA_dTdTrEe, &ddB_dTdDevEe, &dB_dTrEe, &dB_dTdTrEe);
+    }
     
     if (A <= -1. + 1.e-5){ // saturated
         double trEe_0 = q0->_Ee.trace();
@@ -662,6 +728,10 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
             *DsigDEe *= (2.*Gextra*B);
             STensorOperation::prodAdd(devEe, dB_dDevEe, 2.*Gextra,  *DsigDEe);            
             STensorOperation::prodAdd(_I,_I,Kextra*(A+trEe*dA_dTrEe),*DsigDEe);
+            
+            if (_extraBranchNLType == TensionCompressionRegularisedType){
+                STensorOperation::prodAdd(devEe, _I, 2.*Gextra*dB_dTrEe, *DsigDEe);
+            }
         }
         
         // DsigDT
@@ -707,6 +777,16 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
             STensorOperation::prodAdd(devEe, dB_dDevEe, 2.*dGex_dT, *DDsigD_dTdDevEe);
             STensorOperation::prodAdd(devEe, ddB_dTdDevEe, 2.*Gextra, *DDsigD_dTdDevEe);
         }
+        if(DsigD_dTrEe!=NULL){
+            STensorOperation::zero(*DsigD_dTrEe);
+            *DsigD_dTrEe = devEe;
+            *DsigD_dTrEe *= (2.*Gextra*dB_dTrEe);
+        }
+        if(DsigD_dTdTrEe!=NULL){
+            STensorOperation::zero(*DsigD_dTdTrEe);
+            *DsigD_dTdTrEe = devEe;
+            *DsigD_dTdTrEe *= (2.*dGex_dT*dB_dTrEe + 2*Gextra*dB_dTdTrEe);
+        } 
     };
   }
   else
@@ -717,7 +797,8 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
 
 void mlawNonLinearTVM::extraBranch_nonLinearTVE(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) const{
+                        double& Bd, STensor3& dB_dDevEe, double &intB, STensor3 &intB1, STensor3 &intB2, double &dB_dT, double &ddB_dTT, STensor3& ddB_dTdDevEe,
+                        double* dB_dTrEe) const{
 
     // Similar code to previous implementations of extraBranch except we dont do -1 + Av, have to make this a unique implementation.
     
@@ -798,7 +879,8 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const doubl
         double x = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
         Av = getVolumeCorrection()* 1/sqrt(x) ;
         if (tr < 0.){
-            Av = getVolumeCorrection()* 1/sqrt(x) * (1 + _compCorrection/getVolumeCorrection());
+            Av = getVolumeCorrection()* 1/sqrt(x) * (1. + _compCorrection/getVolumeCorrection());
+            // Av += (_compCorrection*tanh(_compCorrection_2*tr*tr+_compCorrection_3)) ;
         }
         
         dA_dTrEe = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5);
@@ -810,10 +892,10 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const doubl
             intA = 1.;
         }
         if (tr < 0.){
-            dA_dTrEe = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5) * (1 + _compCorrection/getVolumeCorrection());
+            dA_dTrEe = - getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x,1.5) * (1. + _compCorrection/getVolumeCorrection());
             if(getXiVolumeCorrection()>0.){
-                intA = getVolumeCorrection()/getXiVolumeCorrection()*sqrt(x) * (1 + _compCorrection/getVolumeCorrection()); 
-                intA -= ( getVolumeCorrection()/getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ) * (1 + _compCorrection/getVolumeCorrection()); 
+                intA = getVolumeCorrection()/getXiVolumeCorrection()*sqrt(x) * (1. + _compCorrection/getVolumeCorrection()); 
+                intA -= ( getVolumeCorrection()/getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ) * (1. + _compCorrection/getVolumeCorrection()); 
             }
             else{
                 intA = 1.;
@@ -833,6 +915,53 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const doubl
         else{
             intB = 1.;
         }   
+        
+    }
+    else if (_ExtraBranch_TVE_option == 3){
+        
+        // Tension-Compression Asymmetry 
+        // Standard form => A_v = a/(sqrt(b*tr*tr+c))
+        
+        // Regularising function
+        double m = _tensionCompressionRegularisation;
+        double expmtr(0.);
+        if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
+        double sigmoid = 1/(1.+expmtr);
+        
+        // Av
+        double x = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
+        Av = sigmoid*1./sqrt(x) + (1.-sigmoid)*(_compCorrection/sqrt(x));
+        
+        dA_dTrEe = - sigmoid*getXiVolumeCorrection()*tr/pow(x,1.5) - (1.-sigmoid)*(getXiVolumeCorrection()*_compCorrection*tr/pow(x,1.5)) 
+                   + (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x) - (_compCorrection*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(x);
+        if(getXiVolumeCorrection()>0.){
+            double integrand = 1./getXiVolumeCorrection()*sqrt(x); // integral of A_v * trEe
+            integrand -= ( 1./getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ); // value at trEe = 0.
+            intA = sigmoid*integrand + (1.-sigmoid)*integrand;
+        }
+        else{
+            intA = 1.;
+        }
+        
+        // Bd    
+        double y = getThetaDevCorrection()*dev.dotprod() + getPiDevCorrection();
+        Bd = sigmoid*1./sqrt(y) + (1.-sigmoid)*(_compCorrection/sqrt(y)); 
+
+        STensorOperation::zero(dB_dDevEe);
+        dB_dDevEe = dev;
+        dB_dDevEe *= (-sigmoid*getThetaDevCorrection()/pow(y,1.5) - (1.-sigmoid)*(getThetaDevCorrection()*_compCorrection/pow(y,1.5)));
+        
+        if(dB_dTrEe!=NULL){
+            *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y) - (_compCorrection*m*expmtr/pow((1.+expmtr),2.))*1./sqrt(y);
+        }    
+        if(getThetaDevCorrection()>0.){
+            double integrand = 1./getThetaDevCorrection()*sqrt(y);  // integral of B_d * devEe
+            integrand -= (1./getThetaDevCorrection()*sqrt(getPiDevCorrection()) );  // value at devEe = 0.   
+            intB = sigmoid*integrand + (1.-sigmoid)*integrand;
+        }
+        else{
+            intB = 1.;
+        }   
     }
     else{ // Test
         // Av = 1.;
@@ -898,7 +1027,7 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const STensor3& Ee,  const doubl
 void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
           const IPNonLinearTVM *q0, IPNonLinearTVM *q1,
           double& Ke, double& Ge, double& DKDTsum, double& DGDTsum,
-          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm) const{
+          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2) const{
 
     // Set Temperature to calculate temperature updated properties
     double T_set = setTemp(T0,T1,1);  // 1->T1, 2->T_mid, 3->T_midshift
@@ -952,27 +1081,15 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
         STensorOperation::zero(ddB_dTdDevEe);
         
         if (_useExtraBranch_TVE){
-            STensorOperation::zero(Bd_stiffnessTerm);
+            STensorOperation::zero(Bd_stiffnessTerm); //This term is for the additional term in Ge that is a tensor
             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);
-                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
-            
-            // Evaluate Bd_stiffnessTerm
-            STensorOperation::zero(temp100);
-            for (int k=0; k<3; k++)
-                for (int l=0; l<3; l++)
-                    for (int p=0; p<3; p++)
-                        for (int q=0; q<3; q++)
-                            temp100(k,l) = q0->_dBd_dDevEe_TVE(p,q)*_Idev(p,q,k,l);
+                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
+            }
+            if(_ExtraBranch_TVE_option == 3){
+                STensorOperation::zero(*Bd_stiffnessTerm2); //This term is for the additional term in Ge because Bd is a function of trEe
             }
         }
         
@@ -1016,27 +1133,13 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
             else{ // all extraBranches
                 
                 if(_ExtraBranch_TVE_option == 1){
-                    
-                    /* Test
-                    double Temp_a(0.), Temp_b(0.);
-                    static STensor3 Term2_test_a, Term2_test_b; 
-                    STensorOperation::doubleContractionSTensor3(q0->_dBd_dDevEe_TVE,devDE,Temp_a);
-                    Term2_test_a = q0->_A[i];
-                    Term2_test_a *= Temp_a;
-                    STensorOperation::doubleContractionSTensor3(q0->_A[i],q0->_dBd_dDevEe_TVE,Temp_b);
-                    Term2_test_b = devDE;
-                    Term2_test_b *= Temp_b;
-                    // Apparently, Term2_test_a and Term2_test_b are not equal
-                    */
                 
                     // Evaluate Bd_stiffnessTerm
                     for (int k=0; k<3; k++)
                         for (int l=0; l<3; l++)
                             for (int p=0; p<3; p++)
-                                for (int q=0; q<3; q++){
-                                    Bd_stiffnessTerm(k,l,p,q) +=  2*_Gi[i]*q0->_A[i](k,l)*q0->_dBd_dDevEe_TVE(p,q) * 1/q0->_Bd_TVE * exp_mid_g; // 1st try
-                                    // Bd_stiffnessTerm(k,l,p,q) += 2*_Gi[i]* q0->_A[i](k,l)*temp100(p,q) * 1/q0->_Bd_TVE * exp_mid_g; // 2nd try
-                                }
+                                for (int q=0; q<3; q++)
+                                    Bd_stiffnessTerm(k,l,p,q) +=  2*_Gi[i]*q0->_A[i](k,l)*q0->_dBd_dDevEe_TVE(p,q) * 1/q0->_Bd_TVE * exp_mid_g;
                 
                     // Ge for DcorKirDE
                     Ge += ( exp_mid_g *_Gi[i] * q0->_Bd_TVE ); // 2* (..) - Dont multiply by 2 if using lambda and mu in Hooks tensor function
@@ -1065,7 +1168,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     p += ( _Ki[i]*q1->_B[i] );
                     
                 }
-                else if (_ExtraBranch_TVE_option == 2){
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
                     
                     // Single Convolution - get branchElasticStrain
                     for (int k=0; k<3; k++)
@@ -1079,8 +1182,16 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     branchElasticStrain(0,0) += 1./3.*q1->_B[i];
                     branchElasticStrain(1,1) += 1./3.*q1->_B[i];
                     branchElasticStrain(2,2) += 1./3.*q1->_B[i];
-                    mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                    if (_ExtraBranch_TVE_option == 2){
+                        mlawNonLinearTVM::extraBranch_nonLinearTVE(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,
+                                                                Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
+                        q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
+                    }
                     
                     q1->_Av_TVE_vector[i] = Av;
                     q1->_Bd_TVE_vector[i] = Bd;
@@ -1089,18 +1200,24 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     q1->_intAv_TVE_vector[i] = intA;
                     q1->_intBd_TVE_vector[i] = intB;
                     
-                    // Evaluate Bd_stiffnessTerm
+                    // Evaluate Bd_stiffnessTerm and/or Bd_stiffnessTerm2
                     for (int k=0; k<3; k++)
                         for (int l=0; l<3; l++)
                             for (int p=0; p<3; p++)
-                                for (int q=0; q<3; q++)
+                                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){
+                                        (*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)
+                                    }
+                                }
                                 
                     // Ge for DcorKirDE
                     Ge += ( exp_mid_g *_Gi[i] * Bd); // 2* (..) - Dont multiply by 2 if using lambda and mu in Hooks tensor function
+                            // exp_mid_g = d(dev branchElasticStrain)_d(devEe)
                     
                     // Ke for DcorKirDE
-                    Ke += ( exp_mid_k *_Ki[i] * (dA_dTrEe*q1->_B[i] + Av) ); 
+                    Ke += ( exp_mid_k *_Ki[i] * (dA_dTrEe*q1->_B[i] + Av) );  // exp_mid_k = d(tr branchElasticStrain)_d(trEe)
                     
                     // Add to deviatoric stress
                     devK += ( 2.*_Gi[i]*Bd*q1->_A[i] );
@@ -1120,6 +1237,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
             std::vector<double>& DDB_DTDtrE = q1->_DDB_DTDtrE;
             std::vector<STensor3>& DA_DT = q1->_DA_DT;
             std::vector<STensor3>& DDA_DTT = q1->_DDA_DTT;
+            std::vector<STensor3>& DA_DtrE = q1->_DA_DtrE;
+            std::vector<STensor3>& DDA_DTDtrE = q1->_DDA_DTDtrE;
             std::vector<STensor43>& DA_DdevE = q1->_DA_DdevE;
             std::vector<STensor43>& DDA_DTDdevE = q1->_DDA_DTDdevE;
                     
@@ -1196,11 +1315,11 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
                     // 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
-                    static STensor3 dDevEei_DT, ddDevEei_DTT;
+                    static STensor3 dDevEei_DT, ddDevEei_DTT; // branchElasticStrain deviator
                     dTrEei_DT = q0->_B[i] * Dexp_rec_k_DT + Dexp_mid_k_DT*trDE;
                     ddTrEei_DTT = q0->_B[i] * DDexp_rec_k_DTT + DDexp_mid_k_DTT*trDE;
                     for (int k=0; k<3; k++)
@@ -1219,6 +1338,17 @@ 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); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
+                            
+                            if (_ExtraBranch_TVE_option == 3){
+                                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);  
+                                
+                                // _Idev is not required !! _I is not pre-multiplied !!
+                                DA_DtrE[i](k,l) = q1->_dBd_dTrEe_TVE_vector[i] * exp_mid_k * q1->_A[i](k,l); // *_I(r,s);
+                                DDA_DTDtrE[i](k,l) = q1->_dBd_dTrEe_TVE_vector[i] * (exp_mid_k*dDevEei_DT(k,l) + Dexp_mid_k_DT*q1->_A[i](k,l)); // *_I(r,s);
+                            }
+                                                                
                             for (int r=0; r<3; r++)
                                 for (int s=0; s<3; s++){
                                     DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dDevEe_TVE_vector[i](r,s)*dDevEei_DT(r,s); 
@@ -1283,7 +1413,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){      
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){      
                     // 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
@@ -1297,6 +1427,11 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     for (int k=0; k<3; k++)
                         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){
+                                DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
+                            }
+                            
                             for (int r=0; r<3; r++)
                                 for (int s=0; s<3; s++){
                                     DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dDevEe_TVE_vector[i](r,s)*dDevEei_DT(r,s); 
@@ -1574,6 +1709,12 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
                                     DQiDEe[i](k,l,r,s) = _Ki[i]*q1->_DB_DtrE[i]*_I(k,l)*_I(r,s);
                                     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){
+                                        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);
+                                    }
+                                    
                                     for (int p=0; p<3; p++){
                                         for (int q=0; q<3; q++){
                                             DQiDEe[i](k,l,r,s) += 2*_Gi[i]*q1->_DA_DdevE[i](k,l,p,q)*_Idev(p,q,r,s); // _Idev is not pre-multiplied in this case
@@ -1668,6 +1809,9 @@ 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){
+                                    DDcorKirDTDEe(j,k,r,s) += 2*_Gi[i]*q1->_DDA_DTDtrE[i](j,k)*_I(r,s);
+                                }
                                 for (int p=0; p<3; p++)
                                     for (int q=0; q<3; q++){
                                         DDcorKirDTDEe(j,k,r,s) += 2*_Gi[i]*q1->_DDA_DTDdevE[i](j,k,p,q)*_Idev(p,q,r,s);
@@ -1888,13 +2032,20 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
     
     // Stresses and Effective Moduli
     double Ke(0.), Ge(0.), DKDTsum(0.), DGDTsum(0.);
-    static STensor43 Bd_stiffnessTerm;
+    static STensor43 Bd_stiffnessTerm, Bd_stiffnessTerm2;
     STensorOperation::zero(Bd_stiffnessTerm);
-    ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm); // Updates the values of moduli to the current temperatures and calculate stresses.
+    STensorOperation::zero(Bd_stiffnessTerm2);
+    if(_ExtraBranch_TVE_option != 3){
+        ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm); 
+                // Updates the values of moduli to the current temperatures and calculate stresses.
+    }
+    else{
+        ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm,&Bd_stiffnessTerm2); 
+    }
     
-    // additional branch
+    // additional hyperelastic branch
     static double DsigV_dT(0.), DDsigV_dTT(0.), DDsigV_dTdTrEe(0.);
-    static STensor3 sigExtra, DsigExtraDT, DsigD_dT, DDsigD_dTT;
+    static STensor3 sigExtra, DsigExtraDT, DsigD_dT, DDsigD_dTT, DDsigD_dTdTrEe;
     static STensor43 DsigExtraDE, DDsigD_dTdDevEe;
     STensorOperation::zero(sigExtra);
     STensorOperation::zero(DsigExtraDE);
@@ -1902,9 +2053,16 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
     STensorOperation::zero(DsigD_dT);
     STensorOperation::zero(DDsigD_dTT);
     STensorOperation::zero(DDsigD_dTdDevEe);
+    STensorOperation::zero(DDsigD_dTdTrEe);
     if (_useExtraBranch){
-        mlawNonLinearTVM::extraBranchLaw(E,T,q0,q1,sigExtra,stiff,&DsigExtraDE,&DsigExtraDT,NULL,NULL,
+        if (_extraBranchNLType != TensionCompressionRegularisedType){
+            mlawNonLinearTVM::extraBranchLaw(E,T,q0,q1,sigExtra,stiff,&DsigExtraDE,&DsigExtraDT,NULL,NULL,
                                             &DsigV_dT, &DsigD_dT, &DDsigV_dTT, &DDsigD_dTT, &DDsigV_dTdTrEe, &DDsigD_dTdDevEe);
+        }
+        else{
+            mlawNonLinearTVM::extraBranchLaw(E,T,q0,q1,sigExtra,stiff,&DsigExtraDE,&DsigExtraDT,NULL,NULL,
+                                            &DsigV_dT, &DsigD_dT, &DDsigV_dTT, &DDsigD_dTT, &DDsigV_dTdTrEe, &DDsigD_dTdDevEe, NULL, &DDsigD_dTdTrEe);
+        }
         q1->_corKirExtra = sigExtra;
     }
     q1->_kirchhoff += sigExtra;
@@ -1980,8 +2138,12 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
         for (int i=0; i<3; i++)
             for (int j=0; j<3; j++)
                 for (int k=0; k<3; k++)
-                    for (int l=0; l<3; l++)
+                    for (int l=0; l<3; l++){
                         DDcorKirDTDEe(i,j,k,l) += DDsigV_dTdTrEe*_I(i,j)*_I(k,l) + DDsigD_dTdDevEe(i,j,k,l);
+                        if (_extraBranchNLType == TensionCompressionRegularisedType){
+                            DDcorKirDTDEe(i,j,k,l) += DDsigD_dTdTrEe(i,j)*_I(k,l);
+                        }
+                    }
     }
 
     // Add corKirinf terms in IP
@@ -2057,6 +2219,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
                                 for (int q=0; q<3; q++){
                                     DcorKDE(i,j,k,l) += Bd_stiffnessTerm(i,j,p,q)*_Idev(p,q,k,l);
                                 }
+           DcorKDE += Bd_stiffnessTerm2; // for _ExtraBranch_TVE_option == 3 
         }
         
         
@@ -2278,10 +2441,13 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
 void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, const STensor3& Ee0,
           const IPNonLinearTVM *q0, IPNonLinearTVM *q1,
           double& Ke, double& Ge,
-          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm) const{
+          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2) const{
 
     //
     STensorOperation::zero(Bd_stiffnessTerm);
+    if(_extraBranchNLType == TensionCompressionRegularisedType && Bd_stiffnessTerm2!= NULL){
+        STensorOperation::zero(*Bd_stiffnessTerm2);
+    }
     
     // Set Temperature to calculate temperature updated properties
     double T_set = setTemp(T0,T1,1);  // 1->T1, 2->T_mid, 3->T_midshift
@@ -2318,7 +2484,14 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
         double A(0.), B(0.), dA_dTrEe(0.);
         STensor3 dB_dDevEe;
         double intA(0.), intB(0.), dA_dT(0.), dB_dT(0.);
-        evaluateElasticCorrection(trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT);
+        if(_extraBranchNLType != TensionCompressionRegularisedType){
+            evaluateElasticCorrection(trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT);
+        }
+        else{
+            double dB_dTrEe(0.);
+            evaluateElasticCorrection(trEe, devEe, T1, A, dA_dTrEe, intA, dA_dT, B, dB_dDevEe, intB, dB_dT, NULL ,NULL, NULL, NULL, &dB_dTrEe);
+            q1->_dBd_dTrEe = dB_dTrEe;
+        }
     
         if (A <= -1. + 1.e-5){ // saturated
             double trEe_0 = q0->_Ee.trace();
@@ -2358,8 +2531,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
         for (int k=0; k<3; k++)
             for (int l=0; l<3; l++)
                 for (int p=0; p<3; p++)
-                    for (int q=0; q<3; q++)
+                    for (int q=0; q<3; q++){
                         Bd_stiffnessTerm(k,l,p,q) +=  2*GT*devEe(k,l)*dB_dDevEe(p,q);
+                        if(_extraBranchNLType == TensionCompressionRegularisedType){
+                            (*Bd_stiffnessTerm2)(k,l,p,q) +=  2*GT*devEe(k,l)*_I(p,q)*q1->_dBd_dTrEe;
+                        }
+                    }
                         
     }
     
@@ -2386,16 +2563,10 @@ 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);
-                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
+                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
             }
         }
         
@@ -2468,7 +2639,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){
+             else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
                     
                 // Single Convolution - get branchElasticStrain
                 for (int k=0; k<3; k++)
@@ -2482,8 +2653,16 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 branchElasticStrain(0,0) += 1./3.*q1->_B[i];
                 branchElasticStrain(1,1) += 1./3.*q1->_B[i];
                 branchElasticStrain(2,2) += 1./3.*q1->_B[i];
-                mlawNonLinearTVM::extraBranch_nonLinearTVE(branchElasticStrain,T1,Av,dA_dTrEe,intA,intA1,intA2,dA_dT,ddA_dTT,ddA_dTdTrEe,
+                if (_ExtraBranch_TVE_option == 2){
+                    mlawNonLinearTVM::extraBranch_nonLinearTVE(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,
+                                                                Bd,dB_dDevEe,intB,intB1,intB2,dB_dT,ddB_dTT,ddB_dTdDevEe,&dB_dTrEe);
+                    q1->_dBd_dTrEe_TVE_vector[i] = dB_dTrEe;
+                }
                     
                 q1->_Av_TVE_vector[i] = Av;
                 q1->_Bd_TVE_vector[i] = Bd;
@@ -2496,8 +2675,12 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 for (int k=0; k<3; k++)
                     for (int l=0; l<3; l++)
                         for (int p=0; p<3; p++)
-                            for (int q=0; q<3; q++)
-                                Bd_stiffnessTerm(k,l,p,q) +=  2*_Gi[i]*q1->_A[i](k,l)*dB_dDevEe(p,q) * exp_mid_g; 
+                            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){
+                                    (*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;
+                                } 
+                            }
                                 
                 // Ge for DcorKirDE
                 Ge += ( exp_mid_g *_Gi[i] * Bd); // 2* (..) - Dont multiply by 2 if using lambda and mu in Hooks tensor function
@@ -2523,6 +2706,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
             std::vector<double>& DDB_DTDtrE = q1->_DDB_DTDtrE;
             std::vector<STensor3>& DA_DT = q1->_DA_DT;
             std::vector<STensor3>& DDA_DTT = q1->_DDA_DTT;
+            std::vector<STensor3>& DA_DtrE = q1->_DA_DtrE;
+            std::vector<STensor3>& DDA_DTDtrE = q1->_DDA_DTDtrE;
             std::vector<STensor43>& DA_DdevE = q1->_DA_DdevE;
             std::vector<STensor43>& DDA_DTDdevE = q1->_DDA_DTDdevE;
                     
@@ -2597,7 +2782,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){
                     // 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
@@ -2620,6 +2805,17 @@ 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); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
+                            
+                            if (_ExtraBranch_TVE_option == 3){
+                                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);  
+                                
+                                // _Idev is not required !! _I is not pre-multiplied !!
+                                DA_DtrE[i](k,l) = q1->_dBd_dTrEe_TVE_vector[i] * exp_mid_k * q1->_A[i](k,l); // *_I(r,s);
+                                DDA_DTDtrE[i](k,l) = q1->_dBd_dTrEe_TVE_vector[i] * (exp_mid_k*dDevEei_DT(k,l) + Dexp_mid_k_DT*q1->_A[i](k,l)); // *_I(r,s);
+                            }                                   
+                                    
                             for (int r=0; r<3; r++)
                                 for (int s=0; s<3; s++){
                                     DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dDevEe_TVE_vector[i](r,s)*dDevEei_DT(r,s); 
@@ -2684,7 +2880,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){      
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3){      
                     // 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
@@ -2698,6 +2894,11 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                     for (int k=0; k<3; k++)
                         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){
+                                DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
+                            }
+                            
                             for (int r=0; r<3; r++)
                                 for (int s=0; s<3; s++){
                                     DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dDevEe_TVE_vector[i](r,s)*dDevEei_DT(r,s); 
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
index 0c3fc0177e7f499a30a06ee2d01ffa3918529faa..55e6bc3d7e5b3b7d93b52a3c11e4d2472ad0ff5c 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
@@ -83,30 +83,33 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
                                                  double &A_v, double &dAdE, double &intA, double &dAdT,
                                                  double &B_d, STensor3 &dB_vddev, double &intB, double &dBdT,
                                                  double *ddAdTT = NULL, double *ddBdTT = NULL,
-                                                 double *ddA_vdTdE = NULL, STensor3 *ddB_vdTdevE = NULL) const;
+                                                 double *ddA_vdTdE = NULL, STensor3 *ddB_vdTdevE = NULL,
+                                                 double *dB_dTrEe = NULL, double *dB_dTdTrEe = NULL) const;
     
     virtual void extraBranchLaw(const STensor3& Ee,  const double& T, const IPNonLinearTVM *q0, IPNonLinearTVM *q1, 
                         STensor3& sig, const bool stiff = false, STensor43* DsigDEe = NULL, STensor3* DsigDT = NULL, 
                         double* DsigV_dTrEe = NULL, STensor43* DsigD_dDevEe = NULL,
                         double* DsigV_dT = NULL, STensor3* DsigD_dT = NULL,
                         double* DDsigV_dTT = NULL, STensor3* DDsigD_dTT = NULL, 
-                        double* DDsigV_dTdTrEe = NULL, STensor43* DDsigD_dTdDevEe = NULL) const;
+                        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, 
                         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) const;
+                        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;
                         
     void corKirInfinity(const IPNonLinearTVM *q1, const STensor3& devEe, const double& trEe, const double T, STensor3& CorKirDevInf, double& CorKirTrInf) const;
     
     void ThermoViscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
           const IPNonLinearTVM *q0, IPNonLinearTVM *q1,
           double& Ke, double& Ge, double& DKDTsum, double& DGDTsum,
-          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm) const;
+          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL) const;
           
     void ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, const STensor3& Ee0,
           const IPNonLinearTVM *q0, IPNonLinearTVM *q1,
           double& Ke, double& Ge,
-          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm) const;
+          const double T0, const double T1, const bool stiff, STensor43& Bd_stiffnessTerm, STensor43* Bd_stiffnessTerm2 = NULL) const;
               
     void isotropicHookTensor(const double K, const double G, STensor43& L) const;
                                   
diff --git a/dG3D/benchmarks/nonLinearTVE_allExtraBranches_uniaxial/runTVE.py b/dG3D/benchmarks/nonLinearTVE_allExtraBranches_uniaxial/runTVE.py
index 968cbfd8c8ba36e616e9dadfaf21f028316d6cdd..32dcd320548bad715ffa44a676166bdd85f635a3 100644
--- a/dG3D/benchmarks/nonLinearTVE_allExtraBranches_uniaxial/runTVE.py
+++ b/dG3D/benchmarks/nonLinearTVE_allExtraBranches_uniaxial/runTVE.py
@@ -60,7 +60,7 @@ law1.setViscoelasticMethod(0)
 law1.setShiftFactorConstantsWLF(C1,C2)
 law1.setReferenceTemperature(Tref)
 
-law1.setExtraBranchNLType(2)
+law1.setExtraBranchNLType(3)
 law1.useExtraBranchBool(True)
 law1.useExtraBranchBool_TVE(True,2)
 law1.setVolumeCorrection(1.48385807e-02, 4.86512658e-01, 2.03662889e-04, 5.47347064e-02, 6.57278720e+00, 2.77066402e-03)
diff --git a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
index 698bf3d4ec2619ffc6313f554fd450fb198e1e25..caf8ca522facba0151c80366d46a80b4b754d112 100644
--- a/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
+++ b/dG3D/benchmarks/nonLinearTVP_allExtraBranches_cube/cubeTVP.py
@@ -76,7 +76,6 @@ hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 2., 10.)
 # hardent.setTemperatureFunction_hexp(ShiftFactor_tempfunc)
 
 law1 = NonLinearTVENonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp,False,1e-8)
-law1.setExtraBranchNLType(2)
 
 law1.setCompressionHardening(hardenc)
 law1.setTractionHardening(hardent)
@@ -90,6 +89,7 @@ law1.setViscoelasticMethod(0)
 law1.setShiftFactorConstantsWLF(C1,C2)
 law1.setReferenceTemperature(Tref)
 
+law1.setExtraBranchNLType(3)
 law1.useExtraBranchBool(True)
 law1.useExtraBranchBool_TVE(True,2)
 law1.setVolumeCorrection(1.48385807e-02, 4.86512658e+00, 2.03662889e-04, 5.47347064e-02, 6.57278720e+00, 2.77066402e-03)
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index a3004ada9f70cea288653a9e36fd6f9eb2cdfa2f..dcdbf9b1ae8108935b7657fb44943e85edc6aee1 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -247,7 +247,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
 
   _stressLaw->stress(ipvcur,ipvprev,false,checkfrac,dTangent);
   
-#if 0
+#if 1
   _stressLaw->stress(ipvcur,ipvprev,true,checkfrac,dTangent);
   ipvcur->getConstRefToDeformationGradient().print("F Analytique"); // FLE
   ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Analytique"); // FLE
@@ -288,6 +288,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
         Msg::Info("d Mechanical source dp Analytique: %e", ipvcur->getRefTodMechanicalSourcedNonLocalVariable()(j,i));
         Msg::Info("dEMFieldSource dp Analytique: %e",ipvcur->getConstRefTodEMFieldSourcedNonLocalVariable()(j,i));
     }
+  }
   for (int i=0; i< numExtraDof; i++)
   {
     Msg::Info("Field Analytique: %e",ipvcur->getConstRefToField(i)); // FLE
@@ -323,7 +324,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
         ipvcur->getConstRefTodMechanicalSourcedVectorCurl()[i][k].print("dMechanicalSource dMagneticVectorCurl Analytique");
         ipvcur->getRefTodEMFieldSourcedVectorPotential()[i][k].print("dEMFieldSource dMagneticVectorPotential Analytique");
         ipvcur->getRefTodEMFieldSourcedVectorCurl()[i][k].print("dEMFieldSource dMagneticVectorCurl Analytique");
-    }*/
+    }
     for (unsigned int k = 0; k < numNonLocalVars; ++k)
     {
         Msg::Info("dp dExtraDofDiffusionField Analytique: %e",ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(k,i));
@@ -614,7 +615,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
         }
 
     }
-#if 0
+#if 1
     ipvcur->getConstRefToDeformationGradient().print("F Numerique"); // FLE
     ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Numerique"); // FLE
     ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");
@@ -687,7 +688,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
           ipvcur->getConstRefTodMechanicalSourcedVectorCurl()[i][k].print("dMechanicalSource dMagneticVectorCurl Numerique");
           ipvcur->getRefTodEMFieldSourcedVectorPotential()[i][k].print("dEMFieldSource dMagneticVectorPotential Numerique");
           ipvcur->getRefTodEMFieldSourcedVectorCurl()[i][k].print("dEMFieldSource dMagneticVectorCurl Numerique");
-      }*/
+      }
         for (unsigned int k = 0; k < numNonLocalVars; ++k)
         {
             Msg::Info("dp dExtraDofDiffusionField Numerique: %e",ipvcur->getRefTodLocalVariableDExtraDofDiffusionField()(k,i));
@@ -3742,8 +3743,11 @@ void HyperViscoElasticDG3DMaterialLaw::setViscoElasticData(const std::string fil
 void HyperViscoElasticDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc ,double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
 };
-void HyperViscoElasticDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void HyperViscoElasticDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
+};
+void HyperViscoElasticDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
 
 void HyperViscoElasticDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
@@ -3859,9 +3863,16 @@ void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setViscoElasticData(const
 void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc ,double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc, _xivc, _zetavc, _dc, _thetadc, _pidc);
 };
-void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
+};
+void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setExtraBranchNLType(const int method){
+  _viscoLaw.setExtraBranchNLType(method);
+};
+void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
+
 void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setCompressionHardening(const J2IsotropicHardening& comp){
   _viscoLaw.setCompressionHardening(comp);
 };
@@ -10114,8 +10125,14 @@ void NonLinearTVMDG3DMaterialLaw::useExtraBranchBool_TVE(const bool useExtraBran
 void NonLinearTVMDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVMDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void NonLinearTVMDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
+};
+void NonLinearTVMDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
+};
+void NonLinearTVMDG3DMaterialLaw::setExtraBranchNLType(const int method){
+  _viscoLaw.setExtraBranchNLType(method);
 };
 void NonLinearTVMDG3DMaterialLaw::setMullinsEffect(const mullinsEffect& mullins){
   _viscoLaw.setMullinsEffect(mullins);
@@ -10350,9 +10367,16 @@ void NonLinearTVEDG3DMaterialLaw::useExtraBranchBool(const bool useExtraBranch){
 void NonLinearTVEDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVEDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void NonLinearTVEDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
+  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
+};
+void NonLinearTVEDG3DMaterialLaw::setExtraBranchNLType(const int method){
+  _viscoLaw.setExtraBranchNLType(method);
+};
+void NonLinearTVEDG3DMaterialLaw::setTensionCompressionRegularisation(double tensionCompressionRegularisation){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
+
 // Added extra argument in the NonLinearTVEDG3DIPVariable contructor - Tinitial   (added from LinearThermoMech)
 void NonLinearTVEDG3DMaterialLaw::createIPState(IPStateBase* &ips, bool hasBodyForce, const bool* state_,const MElement *ele, const int nbFF_, const IntPt *GP, const int gpt) const
 {
@@ -10601,7 +10625,6 @@ void NonLinearTVPDG3DMaterialLaw::setViscosityEffect(const viscosityLaw& v, cons
 void NonLinearTVPDG3DMaterialLaw::setMullinsEffect(const mullinsEffect& mullins){
   _viscoLaw.setMullinsEffect(mullins);
 };
-
 void NonLinearTVPDG3DMaterialLaw::setYieldPowerFactor(const double n){
   _viscoLaw.setPowerFactor(n);
 };
@@ -10636,10 +10659,15 @@ void NonLinearTVPDG3DMaterialLaw::useExtraBranchBool(const bool useExtraBranch){
 void NonLinearTVPDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void NonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, 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){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
-
 // 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
 {
@@ -10910,10 +10938,15 @@ void NonLinearTVENonLinearTVPDG3DMaterialLaw::useExtraBranchBool_TVE(const bool
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setVolumeCorrection(double _vc, double _xivc, double _zetavc, double _dc, double _thetadc, double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
-void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection){
-  _viscoLaw.setExtraBranch_CompressionParameter(compCorrection);
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, 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){
+  _viscoLaw.setTensionCompressionRegularisation(tensionCompressionRegularisation); 
 };
-
 // 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 bedc4a60f28f1b95c220d6146ed65ad911ea23f3..a0efb8fee0d4c88b6baf66ee1b009f15c8fc8558 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -893,8 +893,9 @@ class HyperViscoElasticDG3DMaterialLaw : public dG3DMaterialLaw{
     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);
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
     void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
     #ifndef SWIG
      HyperViscoElasticDG3DMaterialLaw(const HyperViscoElasticDG3DMaterialLaw &source);
     virtual ~HyperViscoElasticDG3DMaterialLaw(){}
@@ -946,8 +947,10 @@ class HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw  : public dG3DMaterialLaw
     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);
-        
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
+    void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    
     void setCompressionHardening(const J2IsotropicHardening& comp);
     void setTractionHardening(const J2IsotropicHardening& trac);
     void setKinematicHardening(const kinematicHardening& kin);
@@ -3240,7 +3243,9 @@ class NonLinearTVMDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     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);
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
+    void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
     
     void setReferenceTemperature(const double Tref);
     void setShiftFactorConstantsWLF(const double C1, const double C2);
@@ -3319,7 +3324,9 @@ class NonLinearTVEDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     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);
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
+    void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
     
     void setReferenceTemperature(const double Tref);
     void setTempFuncOption(const double TemFuncOpt);
@@ -3433,7 +3440,9 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     // 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);
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
+    void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
     
     void setCompressionHardening(const J2IsotropicHardening& comp);
     void setTractionHardening(const J2IsotropicHardening& trac);
@@ -3538,8 +3547,10 @@ class NonLinearTVENonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // pub
     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);
-
+    void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
+    void setExtraBranchNLType(int type);
+    void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
+    
     void setCompressionHardening(const J2IsotropicHardening& comp);
     void setTractionHardening(const J2IsotropicHardening& trac);
     void setKinematicHardening(const kinematicHardening& kin);