diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
index a351ebed3a258ec9016fcbc9c7af2c059597b92f..ac55da68ee35a585b609c1062a336fe140445be9 100644
--- a/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
+++ b/NonLinearSolver/materialLaw/mlawHyperelastic.cpp
@@ -16,6 +16,7 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const int num,const double E,const
     _viscoMethod(Maxwell),_N(0.),_order(1),_Ki(0),_ki(0),_Gi(0),_gi(0), _I4(1.,1.), _I(1.), 
     _volCorrection(1.), _xivolCorrection(1.), _zetavolCorrection(1.), _devCorrection(1.), _thetadevCorrection(1.), _pidevCorrection(1.), 
     _compCorrection(1.), _compCorrection_2(0.), _compCorrection_3(0.), _tensionCompressionRegularisation(5.), 
+	_volCorrection2(1.), _xivolCorrection2(1.), _zetavolCorrection2(1.), _devCorrection2(1.), _thetadevCorrection2(1.), _pidevCorrection2(1.),
     _extraBranchType(Bilogarithmic), _extraBranchNLType(tanhType){
         
    _Idev = _I4;
@@ -63,7 +64,9 @@ mlawHyperViscoElastic::mlawHyperViscoElastic(const mlawHyperViscoElastic& src):
       _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)
+      _tensionCompressionRegularisation(src._tensionCompressionRegularisation), _extraBranchType(src._extraBranchType), _extraBranchNLType(src._extraBranchNLType),
+	  _volCorrection2(src._volCorrection2), _xivolCorrection2(src._xivolCorrection2),
+      _zetavolCorrection2(src._zetavolCorrection2), _devCorrection2(src._devCorrection2), _thetadevCorrection2(src._thetadevCorrection2)
 {
 
 };
@@ -79,14 +82,14 @@ mlawHyperViscoElastic& mlawHyperViscoElastic::operator=(const materialLaw &sourc
       _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;
+      _volCorrection=src->_volCorrection; _xivolCorrection=src->_xivolCorrection;
+      _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;
+       _compCorrection = src->_compCorrection; _compCorrection_2 = src->_compCorrection_2; _compCorrection_3 = src->_compCorrection_3;
        _extraBranchNLType = src->_extraBranchNLType;
        _tensionCompressionRegularisation = src->_tensionCompressionRegularisation;
+       _volCorrection2=src->_volCorrection2; _xivolCorrection2=src->_xivolCorrection2;
+       _zetavolCorrection2=src->_zetavolCorrection2; _devCorrection2=src->_devCorrection2; _thetadevCorrection2=src->_thetadevCorrection2; _pidevCorrection2=src->_pidevCorrection2;
   }
   return *this;
 }
@@ -366,13 +369,66 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
         else{
             intB = -dev.dotprod()/2.;
         } 
+    }
+    else if(_extraBranchNLType == hyper_exp_TCasymm_Type){
         
-        /*
-        if (psiInfCorrector != NULL){
-            double temp1000 = getVolumeCorrection()/sqrt(getXiVolumeCorrection()) * asinh(sqrt(getXiVolumeCorrection()/getZetaVolumeCorrection())*tr);
-            *psiInfCorrector = (sigmoid*temp1000 + (1.-sigmoid)*temp1000); // approximation
-        }*/
-        // Msg::Error(" Inside evaluatePhiPCorrection, A_v = %e, B_d = %e !!", A_v, B_d);
+        // Regularising function
+        double m = _tensionCompressionRegularisation;
+        double expmtr = exp(-m*tr);
+        // if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
+        double sigmoid = 1/(1.+expmtr);
+
+        double x1 = getXiVolumeCorrection()*tr*tr + getZetaVolumeCorrection();
+        double x2 = _xivolCorrection2*tr*tr + _zetavolCorrection2;
+
+        A_v = -1. + sigmoid*(getVolumeCorrection()/sqrt(x1)) + (1-sigmoid)*_compCorrection*(getVolumeCorrection()/sqrt(x1) + _volCorrection2*exp(x2));
+
+        dA_vdE = sigmoid*(-getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x1,1.5))
+        		   + (1.-sigmoid)*_compCorrection*( -getVolumeCorrection()*getXiVolumeCorrection()*tr/pow(x1,1.5) + 2*_volCorrection2*_xivolCorrection2*tr*exp(x2))
+                   + (m*expmtr/pow((1.+expmtr),2.))*getVolumeCorrection()*1./sqrt(x1) - (m*expmtr/pow((1.+expmtr),2.))*_compCorrection*(getVolumeCorrection()/sqrt(x1) + _volCorrection2*exp(x2)) ;
+
+        if(getXiVolumeCorrection()>0. && _xivolCorrection2>0.){
+            double integrand1 = 1./getXiVolumeCorrection()*sqrt(x1); // integral of A_v * trEe
+            integrand1 -= ( 1./getXiVolumeCorrection()*sqrt(getZetaVolumeCorrection()) ); // value at trEe = 0.
+            integrand1 *= getVolumeCorrection();
+
+            double integrand2 = exp(x2)/(2*_xivolCorrection2);
+            integrand2 -= exp(_zetavolCorrection2)/(2*_xivolCorrection2);
+			integrand2 *= _volCorrection2;
+
+            intA = - tr*tr/2. + sigmoid*integrand1 + (1.-sigmoid)*_compCorrection*(integrand1 + integrand2);
+        }
+        else{
+            intA = - tr*tr/2.;
+        }
+
+        double y1 = getThetaDevCorrection()*dev.dotprod() + getPiDevCorrection();
+        double y2 = _thetadevCorrection2*dev.dotprod() + _pidevCorrection2;
+        B_d = -1. + sigmoid*(getDevCorrection()/sqrt(y1)) + (1-sigmoid)*_compCorrection*(getDevCorrection()/sqrt(y1) + _devCorrection2*exp(y2));
+
+        STensorOperation::zero(dB_vddev);
+        dB_vddev = dev;
+        dB_vddev *= ( sigmoid*(-getDevCorrection()*getThetaDevCorrection()/pow(y1,1.5))
+        				+ (1.-sigmoid)*_compCorrection*( -getDevCorrection()*getThetaDevCorrection()/pow(y1,1.5) + 2*_devCorrection2*_thetadevCorrection2*exp(y2)) );
+
+        if(dB_dTrEe!=NULL){
+            *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*getDevCorrection()*1./sqrt(y1) - (m*expmtr/pow((1.+expmtr),2.))*_compCorrection*(getDevCorrection()/sqrt(y1) + _devCorrection2*exp(y2));
+        }
+        if(getThetaDevCorrection()>0. && _thetadevCorrection2>0.){
+
+            double integrand1 = 1./getThetaDevCorrection()*sqrt(y1);  // integral of B_d * devEe
+            integrand1 -= (1./getThetaDevCorrection()*sqrt(getPiDevCorrection()) );  // value at devEe = 0.
+            integrand1 *= getDevCorrection();
+
+            double integrand2 = exp(y2)/(2*_thetadevCorrection2);
+            integrand2 -= exp(_pidevCorrection2)/(2*_thetadevCorrection2);
+			integrand2 *= _devCorrection2;
+
+			intB = -dev.dotprod()/2. + sigmoid*integrand1 + (1.-sigmoid)*_compCorrection*(integrand1 + integrand2);
+        }
+        else{
+            intB = -dev.dotprod()/2.;
+        }
     }
     else{
       Msg::Error("_extraBranchNLType %d has not been defined",_extraBranchNLType);
diff --git a/NonLinearSolver/materialLaw/mlawHyperelastic.h b/NonLinearSolver/materialLaw/mlawHyperelastic.h
index 7df27b06289819de19cd6256287bbe37d1c51a3a..46ab9a2e5b810ee42c512e37eb3c646a862b9b34 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,TensionCompressionRegularisedType=4};
+    enum extraBranchNLType{expType=0,tanhType=1,sigmoidType=2,hyperType=3,TensionCompressionRegularisedType=4,hyper_exp_TCasymm_Type=5};
   #ifndef SWIG
   protected:
     STensor43 _I4; // symetric 4th order unit tensor
@@ -61,6 +61,13 @@ class mlawHyperViscoElastic : public materialLaw{
     extraBranchType _extraBranchType;
     extraBranchNLType _extraBranchNLType;
 
+    double _volCorrection2; //Vc correction on volume potential
+    double _xivolCorrection2; //Xi correction on volume potential
+    double _zetavolCorrection2; //Zeta correction on volume potential
+    double _devCorrection2; //Dc correction on volume potential
+    double _thetadevCorrection2; //Theta correction on volume potential
+    double _pidevCorrection2; //pi correction on volume potential
+
   protected:
     virtual double deformationEnergy(const IPHyperViscoElastic& q) const ;
     virtual void dDeformationEnergydF(const IPHyperViscoElastic& q,  const STensor3 &P, STensor3 &dElasticEnergydF) const ;
@@ -106,6 +113,8 @@ class mlawHyperViscoElastic : public materialLaw{
     virtual double getPiDevCorrection() const {return _pidevCorrection;};
     virtual void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3) {
                                     _compCorrection = compCorrection; _compCorrection_2 = compCorrection_2; _compCorrection_3 = compCorrection_3;};
+    virtual void setAdditionalVolumeCorrections(double V3, double V4, double V5, double D3, double D4, double D5) {
+    								_volCorrection2 = V3; _xivolCorrection2 = V4; _zetavolCorrection2 = V5; _devCorrection2 = D3; _thetadevCorrection2 = D4; _pidevCorrection2 = D5;};
     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 7b81fbf99829b11a98692428008afbb6367d54c6..26c5b923fe07694fab93b34acc3b229e7f5b844d 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVENonLinearTVP.cpp
@@ -385,7 +385,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   STensorOperation::zero(Ge_Tensor_pr); STensorOperation::zero(Ge_Tensor);
   STensorOperation::zero(Ge_TrEeTensor_pr); STensorOperation::zero(Ge_TrEeTensor);
   
-  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+  if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
       // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
 	  mlawNonLinearTVM::ThermoViscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm,&Ge_TrEeTensor);
         Ge_TrEeTensor_pr = Ge_TrEeTensor;
@@ -473,7 +473,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
   devPhi = devPhipr;
   
   // Initialise the rest
-  if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+  if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
     // 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);
   }
@@ -626,7 +626,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
                 }
                 
                 // Update Phi
-                if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     // 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);
@@ -679,7 +679,7 @@ void mlawNonLinearTVENonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow_nonL
             }
 
             // Correct Phi
-            if (_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                 // TC asymmetry -> for either case of TensionCompressionRegularisedType and _ExtraBranch_TVE_option == 3
                 getIterated_DPhi(T0,T,q0,q1,Gamma,Cxtr,Cxdev,Cepr,Eepr,trXn,devXn,Ke,Ge,Ge_Tensor,ptilde,devPhi,Phi,N,expGN,dexpAdA,
                     Dho,Dho2inv,Dho2_u_inv,Dho2_v_inv,&Ge_TrEeTensor);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
index 101e7fe590953c0cb51e69facc2d85c89b116ead..b7d3305e5fdd2a87c4865691ba0341ee560e73aa 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.cpp
@@ -21,7 +21,8 @@ mlawNonLinearTVM::mlawNonLinearTVM(const int num,const double E,const double nu,
      mlawPowerYieldHyper(num, E, nu, rho, tol, matrixbyPerturbation, pert),
      _Tinitial(Tinitial),_scalarAlpha(Alpha),_scalarK(KThCon),_Cp(Cp),
      _thermalEstimationPreviousConfig(thermalEstimationPreviousConfig), _mullinsEffect(NULL), _elasticPotential(NULL), 
-     _useExtraBranch(false), _useExtraBranch_TVE(false), _ExtraBranch_TVE_option(2), _V1(0.),_V2(0.),_D1(0.),_D2(0.),_Ci(0.){
+     _useExtraBranch(false), _useExtraBranch_TVE(false), _ExtraBranch_TVE_option(2), _V1(0.),_V2(0.),_D1(0.),_D2(0.),_Ci(0.),
+	 _V3(0.),_V4(0.),_V5(0.),_D3(0.),_D4(0.),_D5(0.){
 
   _G = _mu;
 
@@ -62,6 +63,7 @@ mlawNonLinearTVM::mlawNonLinearTVM(const mlawNonLinearTVM& src): mlawPowerYieldH
   _useExtraBranch_TVE = src._useExtraBranch_TVE;
   _ExtraBranch_TVE_option = src._ExtraBranch_TVE_option;
   _V1 = src._V1; _V2 = src._V2; _D1 = src._D1; _D2 = src._D2; _Ci = src._Ci;
+  _V3 = src._V3; _V4 = src._V4; _V5 = src._V5; _D3 = src._D3; _D4 = src._D4; _D5 = src._D5;
 
   _temFunc_K = NULL;                                                                    // bulk modulus
   if (src._temFunc_K != NULL){ _temFunc_K = src._temFunc_K->clone();}
@@ -112,6 +114,7 @@ mlawNonLinearTVM& mlawNonLinearTVM::operator=(const materialLaw &source){
     _useExtraBranch_TVE = src->_useExtraBranch_TVE;
     _ExtraBranch_TVE_option = src->_ExtraBranch_TVE_option;
     _V1 = src->_V1; _V2 = src->_V2; _D1 = src->_D1; _D2 = src->_D2; _Ci = src->_Ci;
+    _V3 = src->_V3; _V4 = src->_V4; _V5 = src->_V5; _D3 = src->_D3; _D4 = src->_D4; _D5 = src->_D5;
     
     if(_temFunc_K != NULL) delete _temFunc_K;                                                   // bulk modulus
     if (src->_temFunc_K != NULL){ _temFunc_K = src->_temFunc_K->clone();}
@@ -193,6 +196,15 @@ void mlawNonLinearTVM::setCorrectionsAllBranchesTVE(const int i, const double V1
     Msg::Info("setting: element=%d, V1 = %e, V2 = %e, D1 = %e, D2 = %e",i-1,V1,V2,D1,D2);
   }
 };
+void mlawNonLinearTVM::setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  if (i> _N or i<1)
+    Msg::Error("This setting is invalid %d > %d",i,_N);
+  else{
+    _V3[i-1] = V3; _V4[i-1] = V4; _V5[i-1] = V5; _D3[i-1] = D3; _D4[i-1] = D4; _D5[i-1] = D5;
+
+    Msg::Info("setting: element=%d, V3 = %e, V4 = %e, V5 = %e, D3 = %e, D4 = %e, D5 = %e",i-1,V3,V4,V5,D3,D4,D5);
+  }
+};
 
 void mlawNonLinearTVM::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
   if (i> _N or i<1)
@@ -432,7 +444,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                     // Psy_mech += q._psi_branch[i];*/
                     Msg::Info("setting: _ExtraBranch_TVE_option=%d integrates the viscoelastic stress and currently has no definition for deformation energy. Don't use it with     mullin's effect",_ExtraBranch_TVE_option);
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     Psy_mech += _Ki[i]*q._intAv_TVE_vector[i] + 2*_Gi[i]*q._intBd_TVE_vector[i];
                 }
             }
@@ -551,7 +563,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                     Msg::Info("setting: _ExtraBranch_TVE_option=%d integrates the viscoelastic stress and currently has no definition for deformation energy. Don't use it with     mullin's effect",_ExtraBranch_TVE_option);
                     
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     
                     double dTrEei_DT(0.); // branchElasticStrain trace
                     static STensor3 dDevEei_DT;
@@ -570,7 +582,7 @@ double mlawNonLinearTVM::freeEnergyMechanical(const IPNonLinearTVM& q0, IPNonLin
                                         (*DpsiDE)(j,k) += 2*_Gi[i]*q._Bd_TVE_vector[i]*q._A[i](p,r)*exp_mid_g*_Idev(p,r,j,k); 
                         }
                     }
-                    else { // _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4
+                    else { // _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5
                     
                         // NOTE : _Av_TVE_vector, _Bd_TVE_vector, _intAv_TVE_vector, _intBd_TVE_vector have CTasymm already
                         //          Therefore, to derive CTasymm, first divide to remove it and then derive.
@@ -779,7 +791,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
             STensorOperation::prodAdd(devEe, dB_dDevEe, 2.*Gextra,  *DsigDEe);            
             STensorOperation::prodAdd(_I,_I,Kextra*(A+trEe*dA_dTrEe),*DsigDEe);
             
-            if (_extraBranchNLType == TensionCompressionRegularisedType){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                 STensorOperation::prodAdd(devEe, _I, 2.*Gextra*dB_dTrEe, *DsigDEe);
             }
         }
@@ -807,7 +819,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
         if(DsigD_dT!=NULL){
             STensorOperation::zero(*DsigD_dT);
             double temp = 2*dGex_dT*B;
-            if (_extraBranchNLType == TensionCompressionRegularisedType){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                 temp += 2*Gextra*dB_dTrEe*DtrEeDT;
             }
             *DsigD_dT = devEe;
@@ -822,7 +834,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
         if(DDsigD_dTT!=NULL){
             STensorOperation::zero(*DDsigD_dTT);
             double temp = 2*ddGex_dTT*B;
-            if (_extraBranchNLType == TensionCompressionRegularisedType){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                 temp += 2*dGex_dT*(2*dB_dTrEe*DtrEeDT) + 2*Gextra*dB_dTrEe*DDtrEeDTT; // Assuming ddB_ddTrEe, ddB_dTdTrEe = 0.
             }
             *DDsigD_dTT = devEe;
@@ -836,7 +848,7 @@ void mlawNonLinearTVM::extraBranchLaw(const STensor3& Ee, const double& T, const
             STensorOperation::zero(*DDsigD_dTdDevEe);
             *DDsigD_dTdDevEe = _I4;
             double temp = (2.*dGex_dT*B);
-            if (_extraBranchNLType == TensionCompressionRegularisedType){
+            if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                 temp += 2*Gextra*dB_dTrEe*DtrEeDT;
             }
             *DDsigD_dTdDevEe *= temp;
@@ -1090,6 +1102,67 @@ void mlawNonLinearTVM::extraBranch_nonLinearTVE(const int i, const STensor3& Ee,
                 intB = 1.;
             }
         }
+    else if (_ExtraBranch_TVE_option == 5){
+
+                // Sigmoid + Exponential + Tension-Compression Asymmetry
+                // Standard form => A_v = a/(sqrt(b*tr*tr+c)) + d*exp(f*tr*tr+g)
+
+                // Regularising function
+                double m = _tensionCompressionRegularisation;
+                double expmtr = exp(-m*tr);
+                // if (exp(-m*tr)<1.e+10){ expmtr = exp(-m*tr);}
+                double sigmoid = 1/(1.+expmtr);
+
+                // Av
+                double x1 = _V1[i]*tr*tr + _V2[i];
+                double x2 = _V4[i]*tr*tr + _V5[i];
+                Av = sigmoid*(getVolumeCorrection()/sqrt(x1)) + (1-sigmoid)*_Ci[i]*(getVolumeCorrection()/sqrt(x1) + _V3[i]*exp(x2));
+
+                dA_dTrEe = sigmoid*(-getVolumeCorrection()*_V1[i]*tr/pow(x1,1.5)) + (1.-sigmoid)*_Ci[i]*( -getVolumeCorrection()*_V1[i]*tr/pow(x1,1.5) + 2*_V3[i]*_V4[i]*tr*exp(x2))
+                           + (m*expmtr/pow((1.+expmtr),2.))*getVolumeCorrection()*1./sqrt(x1) - (m*expmtr/pow((1.+expmtr),2.))*_Ci[i]*(getVolumeCorrection()/sqrt(x1) + _V3[i]*exp(x2)) ;
+
+                if(_V1[i]>0. && _V4[i]>0.){
+                    double integrand1 = 1./_V1[i]*sqrt(x1); // integral of A_v * trEe
+                    integrand1 -= ( 1./_V1[i]*sqrt(_V2[i]) ); // value at trEe = 0.
+                    integrand1 *= getVolumeCorrection();
+
+                    double integrand2 = exp(x2)/(2*_V4[i]);
+                    integrand2 -= exp(_V5[i])/(2*_V4[i]);
+					integrand2 *= _V3[i];
+
+                    intA = sigmoid*integrand1 + (1.-sigmoid)*_Ci[i]*(integrand1 + integrand2);
+                }
+                else{
+                    intA = 1.;
+                }
+
+                // Bd
+                double y1 = _D1[i]*dev.dotprod() + _D2[i];
+                double y2 = _D4[i]*dev.dotprod() + _D5[i];
+                Bd = sigmoid*(getDevCorrection()/sqrt(y1)) + (1-sigmoid)*_Ci[i]*(getDevCorrection()/sqrt(y1) + _D3[i]*exp(y2));
+
+                STensorOperation::zero(dB_dDevEe);
+                dB_dDevEe = dev;
+                dB_dDevEe *= ( sigmoid*(-getDevCorrection()*_D1[i]/pow(y1,1.5)) + (1.-sigmoid)*_Ci[i]*( -getDevCorrection()*_D1[i]/pow(y1,1.5) + 2*_D3[i]*_D4[i]*exp(y2)) );
+
+                if(dB_dTrEe!=NULL){
+                    *dB_dTrEe = (m*expmtr/pow((1.+expmtr),2.))*getDevCorrection()*1./sqrt(y1) - (m*expmtr/pow((1.+expmtr),2.))*_Ci[i]*(getDevCorrection()/sqrt(y1) + _D3[i]*exp(y2));
+                }
+                if(_D1[i]>0. && _D4[i]>0.){
+                    double integrand1 = 1./_D1[i]*sqrt(y1);  // integral of B_d * devEe
+                    integrand1 -= (1./_D1[i]*sqrt(_D2[i]) );  // value at devEe = 0.
+                    integrand1 *= getDevCorrection();
+
+                    double integrand2 = exp(y2)/(2*_D4[i]);
+                    integrand2 -= exp(_D5[i])/(2*_D4[i]);
+					integrand2 *= _D3[i];
+
+					intB = sigmoid*integrand1 + (1.-sigmoid)*_Ci[i]*(integrand1 + integrand2);
+                }
+                else{
+                    intB = 1.;
+                }
+            }
     else{ // Test
         // Av = 1.;
         // intA1 = tr;
@@ -1159,7 +1232,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
 
     // NEW
     STensorOperation::zero(Bd_stiffnessTerm);
-    if(_extraBranchNLType == TensionCompressionRegularisedType && Bd_stiffnessTerm2!= NULL){
+    // if(_extraBranchNLType == TensionCompressionRegularisedType  || _extraBranchNLType == hyper_exp_TCasymm_Type && Bd_stiffnessTerm2!= NULL){
+    if(Bd_stiffnessTerm2!= NULL){
         STensorOperation::zero(*Bd_stiffnessTerm2);
     }
     
@@ -1243,7 +1317,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                 for (int p=0; p<3; p++)
                     for (int q=0; q<3; q++){
                         Bd_stiffnessTerm(k,l,p,q) +=  2*GT*devEe(k,l)*dB_dDevEe(p,q);
-                        if(_extraBranchNLType == TensionCompressionRegularisedType){
+                        if(_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                             (*Bd_stiffnessTerm2)(k,l,p,q) +=  2*GT*devEe(k,l)*_I(p,q)*q1->_dBd_dTrEe;
                         }
                     }
@@ -1365,7 +1439,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                     p += ( _Ki[i]*q1->_B[i] );
                     
                 }
-                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     
                     // Single Convolution - get branchElasticStrain
                     for (int k=0; k<3; k++)
@@ -1403,7 +1477,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                             for (int p=0; p<3; p++)
                                 for (int q=0; q<3; q++){
                                     Bd_stiffnessTerm(k,l,p,q) +=  2*_Gi[i]*q1->_A[i](k,l)*dB_dDevEe(p,q) * exp_mid_g; 
-                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                         (*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)
                                     }
@@ -1517,7 +1591,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     // Here the derivatives are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                     // 2nd order derivatives of _Av_TVE_vector and _Bd_TVE_vector are ignored
                     
@@ -1542,7 +1616,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                 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);  
@@ -1624,7 +1698,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                                                                                 q0->_A[i](k,l) * 1/q0->_Bd_TVE * temp2); 
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     // Here the derivative are made as: B[i] = Av[i]*B[i]; A[i] = Bd[i]*A[i]; Above, BranchElasticStrain = A[i] + 1/3*B[i]*_I
                     // 2nd order derivatives of _Av_TVE_vector and _Bd_TVE_vector are ignored
                     
@@ -1640,7 +1714,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
                         for (int l=0; l<3; l++){
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                             }
                             
@@ -1930,7 +2004,7 @@ void mlawNonLinearTVM::getMechSource_nonLinearTVE_term(const IPNonLinearTVM *q0,
                                     DDQiDTDEe[i](k,l,r,s) = _Ki[i]*q1->_DDB_DTDtrE[i]*_I(k,l)*_I(r,s);
                                     dot_DGammaiDEe[i](k,l,r,s) = 1/3*(1.- exp_mid_k)*_I(k,l)*_I(r,s)  + _Idev(k,l,r,s);
                                     
-                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                                    if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                         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);
                                     }
@@ -2045,7 +2119,7 @@ void mlawNonLinearTVM::getTVEdCorKirDT(const IPNonLinearTVM *q0, IPNonLinearTVM
                     for (int l=0; l<3; l++){
                         DDcorKirDTDEe(i,j,k,l) += DDsigV_dTdTrEe*_I(i,j)*_I(k,l);
                         
-                        if (_extraBranchNLType == TensionCompressionRegularisedType){
+                        if (_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                             DDcorKirDTDEe(i,j,k,l) += DDsigD_dTdTrEe(i,j)*_I(k,l);
                         }
                         for (int r=0; r<3; r++)
@@ -2070,7 +2144,7 @@ void mlawNonLinearTVM::getTVEdCorKirDT(const IPNonLinearTVM *q0, IPNonLinearTVM
                             else{ // all extraBranches, because _Idev was not multiplied to _DDA_DTDdevE in the predictorCorrector function
                                 DDcorKirDTDEe(j,k,r,s) += _Ki[i]*q1->_DDB_DTDtrE[i]*_I(j,k)*_I(r,s);
                                 
-                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                     DDcorKirDTDEe(j,k,r,s) += 2*_Gi[i]*q1->_DDA_DTDtrE[i](j,k)*_I(r,s);
                                 }
                                 
@@ -2298,7 +2372,7 @@ void mlawNonLinearTVM::predictorCorrector_ThermoViscoElastic(
     static STensor43 Bd_stiffnessTerm, Bd_stiffnessTerm2;
     STensorOperation::zero(Bd_stiffnessTerm);
     STensorOperation::zero(Bd_stiffnessTerm2);
-    if(_extraBranchNLType == TensionCompressionRegularisedType || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+    if(_extraBranchNLType == TensionCompressionRegularisedType  || _extraBranchNLType == hyper_exp_TCasymm_Type || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
         ThermoViscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,DKDTsum,DGDTsum,T0,T,stiff,Bd_stiffnessTerm,&Bd_stiffnessTerm2); 
                 // Updates the values of moduli to the current temperatures and calculate stresses.
     }
@@ -2681,7 +2755,8 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
 
     //
     STensorOperation::zero(Bd_stiffnessTerm);
-    if(_extraBranchNLType == TensionCompressionRegularisedType && Bd_stiffnessTerm2!= NULL){
+    // if(_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type && Bd_stiffnessTerm2!= NULL){
+    if(Bd_stiffnessTerm2!= NULL){
         STensorOperation::zero(*Bd_stiffnessTerm2);
     }
     
@@ -2763,7 +2838,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 for (int p=0; p<3; p++)
                     for (int q=0; q<3; q++){
                         Bd_stiffnessTerm(k,l,p,q) +=  2*GT*devEe(k,l)*dB_dDevEe(p,q);
-                        if(_extraBranchNLType == TensionCompressionRegularisedType){
+                        if(_extraBranchNLType == TensionCompressionRegularisedType || _extraBranchNLType == hyper_exp_TCasymm_Type){
                             (*Bd_stiffnessTerm2)(k,l,p,q) +=  2*GT*devEe(k,l)*_I(p,q)*q1->_dBd_dTrEe;
                         }
                     }
@@ -2869,7 +2944,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                 // Add to volumetric stress
                 p += ( _Ki[i]*q1->_B[i] );
              }
-             else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+             else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     
                 // Single Convolution - get branchElasticStrain
                 for (int k=0; k<3; k++)
@@ -2907,7 +2982,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                         for (int p=0; p<3; p++)
                             for (int q=0; q<3; q++){
                                 Bd_stiffnessTerm(k,l,p,q) += 2*_Gi[i]*q1->_A[i](k,l)*dB_dDevEe(p,q) * exp_mid_g;
-                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                                if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                     (*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;
                                 } 
                             }
@@ -3017,7 +3092,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                                 }
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     // 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
@@ -3041,7 +3116,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             DDA_DTT[i](k,l) = q1->_Bd_TVE_vector[i]*ddDevEei_DTT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                 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);  
@@ -3123,7 +3198,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                                                                                 q0->_A[i](k,l) * 1/q0->_Bd_TVE * temp2); 
                         }
                  }
-                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                 else if (_ExtraBranch_TVE_option == 2 || _ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                     // 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
@@ -3138,7 +3213,7 @@ void mlawNonLinearTVM::ThermoViscoElasticPredictor_forTVP(const STensor3& Ee, co
                         for (int l=0; l<3; l++){
                             DA_DT[i](k,l) = q1->_Bd_TVE_vector[i]*dDevEei_DT(k,l); 
                             
-                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4){
+                            if (_ExtraBranch_TVE_option == 3 || _ExtraBranch_TVE_option == 4 || _ExtraBranch_TVE_option == 5){
                                 DA_DT[i](k,l) += q1->_A[i](k,l)*q1->_dBd_dTrEe_TVE_vector[i]*dTrEei_DT;
                             }
                             
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
index b34b05e624a54a0a4b8b01ba676b6e1f892bae51..472cf6687936f3678bdb86c6e9928ef3711e5f12 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVM.h
@@ -69,8 +69,14 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
   // Extrabranch Corrector parameters
     std::vector<double> _V1;
     std::vector<double> _V2;
+    std::vector<double> _V3;
+    std::vector<double> _V4;
+    std::vector<double> _V5;
     std::vector<double> _D1;
     std::vector<double> _D2;
+    std::vector<double> _D3;
+    std::vector<double> _D4;
+    std::vector<double> _D5;
     std::vector<double> _Ci;
 
   #ifndef SWIG  
@@ -194,6 +200,7 @@ class mlawNonLinearTVM : public mlawPowerYieldHyper{
     virtual void setViscoElasticData_Bulk(const int i, const double Ki, const double ki);
     virtual void setViscoElasticData_Shear(const int i, const double Gi, const double gi);
     virtual void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    virtual void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     virtual void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     virtual void setViscoElasticData(const std::string filename);
     
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 4e93b930656dae7e1bffe429c1f5ecfa5a061e01..8a35f90ae4fd723266d7817935c7a311288fdcbf 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -3743,6 +3743,9 @@ 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::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void HyperViscoElasticDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -3863,6 +3866,9 @@ 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::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void HyperViscoElastoPlasticPowerYieldDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -10273,6 +10279,9 @@ void NonLinearTVMDG3DMaterialLaw::useExtraBranchBool_TVE(const bool useExtraBran
 void NonLinearTVMDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
+void NonLinearTVMDG3DMaterialLaw::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVMDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -10285,6 +10294,9 @@ void NonLinearTVMDG3DMaterialLaw::setExtraBranchNLType(const int method){
 void NonLinearTVMDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
   _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
 };
+void NonLinearTVMDG3DMaterialLaw::setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalCorrectionsAllBranchesTVE(i,V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVMDG3DMaterialLaw::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
   _viscoLaw.setCompressionCorrectionsAllBranchesTVE(i,Ci);
 };
@@ -10521,6 +10533,9 @@ 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::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVEDG3DMaterialLaw::setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -10813,6 +10828,9 @@ void NonLinearTVPDG3DMaterialLaw::useExtraBranchBool(const bool useExtraBranch){
 void NonLinearTVPDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
+void NonLinearTVPDG3DMaterialLaw::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -11092,6 +11110,9 @@ void NonLinearTVENonLinearTVPDG3DMaterialLaw::useExtraBranchBool_TVE(const bool
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc){
   _viscoLaw.setVolumeCorrection(_vc,_xivc,_zetavc,_dc,_thetadc,_pidc);
 };
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalVolumeCorrections(V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3){
   _viscoLaw.setExtraBranch_CompressionParameter(compCorrection,compCorrection_2,compCorrection_3);
 };
@@ -11104,6 +11125,9 @@ void NonLinearTVENonLinearTVPDG3DMaterialLaw::setTensionCompressionRegularisatio
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2){
   _viscoLaw.setCorrectionsAllBranchesTVE(i,V1,V2,D1,D2);
 };
+void NonLinearTVENonLinearTVPDG3DMaterialLaw::setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5){
+  _viscoLaw.setAdditionalCorrectionsAllBranchesTVE(i,V3,V4,V5,D3,D4,D5);
+};
 void NonLinearTVENonLinearTVPDG3DMaterialLaw::setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci){
   _viscoLaw.setCompressionCorrectionsAllBranchesTVE(i,Ci);
 };
diff --git a/dG3D/src/dG3DMaterialLaw.h b/dG3D/src/dG3DMaterialLaw.h
index 0c2e5cba07b02f64be553e750057036396e7e0c4..1ef36a8e9ac105dd16855aaccbe406f711ada5ea 100644
--- a/dG3D/src/dG3DMaterialLaw.h
+++ b/dG3D/src/dG3DMaterialLaw.h
@@ -894,6 +894,7 @@ 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 setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
     void setExtraBranchNLType(int type);
     void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
@@ -948,6 +949,7 @@ 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 setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(double compCorrection, double compCorrection_2, double compCorrection_3);
     void setExtraBranchNLType(int type);
     void setTensionCompressionRegularisation(double tensionCompressionRegularisation);
@@ -3280,10 +3282,12 @@ 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(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
     void setExtraBranchNLType(const int type);
     void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
     void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     
     void setReferenceTemperature(const double Tref);
@@ -3363,6 +3367,7 @@ class NonLinearTVEDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     void setViscoElasticData_Alpha(const int i, const double Alphai);
     void setViscoElasticData(const std::string filename);
     void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
     void setExtraBranchNLType(const int type);
     void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
@@ -3479,6 +3484,7 @@ class NonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // public Material
     // void setTemperatureFunction_BranchThermalDilationCoefficient(const scalarFunction& Tfunc,int i);
     void useExtraBranchBool(const bool useExtraBranch);
     void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
     void setExtraBranchNLType(const int type);
     void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
@@ -3586,10 +3592,12 @@ class NonLinearTVENonLinearTVPDG3DMaterialLaw : public dG3DMaterialLaw{   // pub
     void useExtraBranchBool(const bool useExtraBranch);
     void useExtraBranchBool_TVE(const bool useExtraBranch_TVE, const int ExtraBranch_TVE_option);
     void setVolumeCorrection(const double _vc, const double _xivc, const double _zetavc, const double _dc, const double _thetadc, const double _pidc);
+    void setAdditionalVolumeCorrections(const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setExtraBranch_CompressionParameter(const double compCorrection, const double compCorrection_2, const double compCorrection_3);
     void setExtraBranchNLType(const int type);
     void setTensionCompressionRegularisation(const double tensionCompressionRegularisation);
     void setCorrectionsAllBranchesTVE(const int i, const double V1, const double V2, const double D1, const double D2);
+    void setAdditionalCorrectionsAllBranchesTVE(const int i, const double V3, const double V4, const double V5, const double D3, const double D4, const double D5);
     void setCompressionCorrectionsAllBranchesTVE(const int i, const double Ci);
     
     void setCompressionHardening(const J2IsotropicHardening& comp);