diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index 7fa320ea097e5cb372309c9542c2d041ff250087..0948b41e1b748b3419c034902c1235f5721cbacb 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -333,6 +333,8 @@ class IPField : public elementsField {
                     EQUIVALENT_EIGENSTRESS,EQUIVALENT_EIGENSTRAIN, 
                     EIGENSTRESS_XX, EIGENSTRESS_YY, EIGENSTRESS_ZZ, EIGENSTRESS_XY, EIGENSTRESS_XZ, EIGENSTRESS_YZ,
                     EIGENSTRAIN_XX, EIGENSTRAIN_YY, EIGENSTRAIN_ZZ, EIGENSTRAIN_XY, EIGENSTRAIN_XZ, EIGENSTRAIN_YZ,
+                    // moduli scalers e.g. mullins effect
+                    MAX_DEFO_ENERGY, BULK_SCALE_FACTOR, SHEAR_SCALE_FACTOR,
                     UNDEFINED};
     enum  Operator { MEAN_VALUE=1, MIN_VALUE, MAX_VALUE, CRUDE_VALUE};      
     static std::string ToString(const int i);
diff --git a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
index 5ae3728b9a0566549831a3e3747000b1f8b4dc1b..f53d5384b6f9d950b75f2773f3f012ac35ffc4d7 100644
--- a/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
+++ b/NonLinearSolver/internalPoints/ipMullinsEffect.cpp
@@ -10,7 +10,7 @@
 
 #include "ipMullinsEffect.h"
 
-IPMullinsEffect::IPMullinsEffect(): eta(0), DetaDpsi(0), DDetaDpsipsi(0), DetaDT(0.), DDetaDTT(0.)
+IPMullinsEffect::IPMullinsEffect(): eta(1.0), DetaDpsi(0.), DDetaDpsipsi(0.), DetaDT(0.), DDetaDTT(0.)
 {
 
 }
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
index 9a047d211abf10f6914424f19a97e8aa9417e541..fce0b2dc47e83dd30caac8a1bf4ac09fe89529fa 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.cpp
@@ -18,7 +18,7 @@ IPNonLinearTVE::IPNonLinearTVE(const J2IsotropicHardening* comp,
                     const kinematicHardening* kin, const int N, const mullinsEffect* mullins): IPHyperViscoElastoPlastic(comp,trac,kin,N), // IPHyperViscoElastic(N), 
                                         _T(298.15),_devDE(0.),_trDE(0.),_DE(0.),_DcorKirDT(0.),_DDcorKirDTT(0.),
                                         _mechSrcTVE(0.),_DmechSrcTVEdT(0.),_DmechSrcTVEdF(0.),
-                                        _thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.){
+                                        _thermalEnergy(0.),_dtShift(0.),_DdtShiftDT(0.),_DDdtShiftDDT(0.), _psiMax(0.){
     _devOGammai.clear();
     _trOGammai.clear();
     _devDOGammaiDT.clear();
@@ -66,7 +66,7 @@ IPNonLinearTVE::IPNonLinearTVE(const IPNonLinearTVE& src): IPHyperViscoElastoPla
             _mechSrcTVE(src._mechSrcTVE),_DmechSrcTVEdT(src._DmechSrcTVEdT),_DmechSrcTVEdF(src._DmechSrcTVEdF),
             _C(src._C),_D(src._D), _E(src._E), _F(src._F), _G(src._G),
             _thermalEnergy(src._thermalEnergy),
-            _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT){
+            _dtShift(src._dtShift),_DdtShiftDT(src._DdtShiftDT),_DDdtShiftDDT(src._DDdtShiftDDT), _psiMax(src._psiMax){
                 
     if (src._ipMullinsEffect != NULL)
         _ipMullinsEffect = dynamic_cast<IPMullinsEffect*>(src._ipMullinsEffect->clone());
@@ -110,6 +110,7 @@ IPNonLinearTVE& IPNonLinearTVE::operator=(const IPVariable& src)
     _dtShift = psrc->_dtShift;
     _DdtShiftDT = psrc->_DdtShiftDT;
     _DDdtShiftDDT = psrc->_DDdtShiftDDT;
+    _psiMax = psrc->_psiMax;
     
     if ( psrc->_ipMullinsEffect != NULL) {
       if (_ipMullinsEffect == NULL){
@@ -173,6 +174,7 @@ void IPNonLinearTVE::restart(){
   restartManager::restart(_dtShift);
   restartManager::restart(_DdtShiftDT);
   restartManager::restart(_DDdtShiftDDT);
+  restartManager::restart(_psiMax);
   
   if (_ipMullinsEffect != NULL)
    restartManager::restart(_ipMullinsEffect);
diff --git a/NonLinearSolver/internalPoints/ipNonLinearTVE.h b/NonLinearSolver/internalPoints/ipNonLinearTVE.h
index 4b55144a0c7eb563627ab62a8c26095169d43690..82c653b90783dc57f1e0db17db9b8d34a6a8cb36 100644
--- a/NonLinearSolver/internalPoints/ipNonLinearTVE.h
+++ b/NonLinearSolver/internalPoints/ipNonLinearTVE.h
@@ -62,6 +62,7 @@ class IPNonLinearTVE : public IPHyperViscoElastoPlastic{
         
         // Mullins Effect IP
         IPMullinsEffect* _ipMullinsEffect;
+        double _psiMax; // maximum strain energy
         
     public:
         IPNonLinearTVE(const J2IsotropicHardening* comp,
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
index a814172f69a60454e4fd5a0d20b2e3abf379045f..706f2ff6142e341ac7a2db5ae2e16f27a7efcd41 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.cpp
@@ -534,11 +534,11 @@ double mlawNonLinearTVE::setTemp(const double T0, const double T1, const int par
     return T_set;
 }
 
-double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const{
+double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const double T0, const double T) const{
 
     // Update the Properties to the current temperature (see below which ones are being used)
     double KT, GT, AlphaT; 
-    getK(&q,KT,Tc); getG(&q,GT,Tc); getAlpha(AlphaT,Tc);
+    getK(&q,KT,T); getG(&q,GT,T); getAlpha(AlphaT,T);
 
     double Psy_mech = 0.;
 
@@ -548,7 +548,7 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
     STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
 
     // Get thermal strain
-    double Eth = 3.*AlphaT*(Tc-_Tinitial);
+    double Eth = 3.*AlphaT*(T-_Tinitial);
 
     // Get Equilibrium freeEnergy_mech
     Psy_mech = KT*0.5*(trEe-Eth)*(trEe-Eth) + GT*STensorOperation::doubledot(devEe,devEe);
@@ -558,45 +558,57 @@ double mlawNonLinearTVE::freeEnergyMechanical(const IPNonLinearTVE& q, const dou
 
         std::vector<double> KiT,GiT,AlphaiT; 
         KiT.resize(_N,0.); GiT.resize(_N,0.); AlphaiT.resize(_N,0.); 
-        getKi(&q,KiT,Tc); getGi(&q,GiT,Tc); getAlphai(AlphaiT,Tc);
+        getKi(&q,KiT,T); getGi(&q,GiT,T); getAlphai(AlphaiT,T);
 
         for (int i=0; i<_N; i++){
             static STensor3 devEebranch;
             devEebranch = devEe;
+            
+            // viscous overstrain = _OGammai
+            // viscous strain (_Gammai) = viscous overstrain - Ee
+            double trGamma = q._trOGammai[i]; 
+            trGamma -= trEe; 
+            static STensor3 devGamma;
+            devGamma = q._devOGammai[i]; 
+            devGamma -= devEe; 
+            
+            devEebranch -= devGamma; // temporary tensor
 
-            // devEebranch -= q._A[i];
-            // Psy_mech += KiT[i]*0.5*(trEe-q._B[i])*(trEe-q._B[i])+GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch);
-
-            double trGamma = q._trOGammai[i];
-            static STensor3 devGamma = q._devOGammai[i];
-            devEebranch -= devGamma;
 
             if (_modelFlag == 2){  // Additional volumetric terms
-                trGamma += 3*(AlphaT-AlphaiT[i])*(Tc-_Tinitial);
+                trGamma += 3*(AlphaT-AlphaiT[i])*(T-_Tinitial);
             }
 
             Psy_mech += KiT[i]*0.5*(trEe-trGamma)*(trEe-trGamma) + GiT[i]*STensorOperation::doubledot(devEebranch,devEebranch);
         }
     }
-
+    
+    // Get thermal Energy
+    // Psy_mech += _Cp*( (T-T0) - T*log(T/T0) ); // this probably does nothing since we are interested only in the elastic energy
+    
     return Psy_mech;
 }
 
 void mlawNonLinearTVE::calculateMullinsEffectScaler(const IPNonLinearTVE* q0, IPNonLinearTVE *q1, const double T) const{
     
+    // Initialise
+    double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
+    
+    // Msg::Error(" Inside calculateMullinsEffectScaler, before updating eta = %e!!",eta);
+    
     // update eta - Mullins Effect Scalar
     if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
-    _mullinsEffect->mullinsEffectScaling(q0->_elasticEnergy, *q0->_ipMullinsEffect, q1->_elasticEnergy,*q1->_ipMullinsEffect);
+    _mullinsEffect->mullinsEffectScaling(q1->_psiMax, *q0->_ipMullinsEffect, q0->_elasticEnergy, q1->_elasticEnergy, *q1->_ipMullinsEffect);
     }
-
+    
     // get eta - Mullins Effect Scalar (at unloading intitiation-> eta = 1.; unloading end-> eta = eta_min; reloading -> eta = eta_min;)
-    double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
     if (q1->_ipMullinsEffect != NULL){
         eta = q1->_ipMullinsEffect->getEta(); 
         DetaDpsi = q1->_ipMullinsEffect->getDetaDpsi();
         DetaDT = q1->_ipMullinsEffect->getDetaDT();
         DDetaDTT = q1->_ipMullinsEffect->getDDetaDTT();
     }
+    // Msg::Error(" Inside calculateMullinsEffectScaler, after getting and updating eta = %e!!",eta);
     
     // update Bulk/Shear PropertyScaleFactor
     q1->getRefToElasticBulkPropertyScaleFactor() = eta;
@@ -1029,7 +1041,7 @@ void mlawNonLinearTVE::ThermoViscoElasticPredictor(const STensor3& Ee, const STe
           double& DKe, double& DGe, double& DKde, double& DGde,
           double& KTsum, double& GTsum, double& DKDTsum, double& DGDTsum,
           const double T0, const double T1, const bool calculateStress) const{
-              
+          
     // Mullins Effect Scaler -> should be before any calls to get functions
     if (_mullinsEffect != NULL && q1->_ipMullinsEffect != NULL){
        mlawNonLinearTVE::calculateMullinsEffectScaler(q0, q1, T1);
@@ -1502,7 +1514,14 @@ void mlawNonLinearTVE::predictorCorrector_ThermoViscoElastic(
     STensorOperation::multSTensor3(F,secondPK,P);                  // 1st PK
 
     // elastic energy
-    q1->_elasticEnergy = freeEnergyMechanical(*q1,T);     // deformationEnergy(*q1,T);
+    q1->_elasticEnergy = freeEnergyMechanical(*q1,T0,T);
+    if (q1->_elasticEnergy > q0->_psiMax){
+        q1->_psiMax = q1->_elasticEnergy;
+    }
+    else{
+        q1->_psiMax = q0->_psiMax;
+    }
+    // Msg::Error(" Inside TVE predCorr, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
 
     // thermal energy
     q1->_thermalEnergy = CpT*T;
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
index 90ae75b93e749f22388e8b8fce874349673cd7a4..26073b30ddfc6cf5ae81ca572e23700c702517fe 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVE.h
@@ -77,7 +77,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
                              double *t_shift = NULL, double *Ddt_shiftDT = NULL, double *DDdt_shiftDT = NULL,
                              const int opt_time = 1, const int opt_order = 2) const;
     virtual double linearInterpTemp(const double t_interp, const double t0, const double t1, const double T0, const double T1) const;
-    virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double Tc) const ;
+    virtual double freeEnergyMechanical(const IPNonLinearTVE& q, const double T0, const double T) const ;
     virtual double setTemp(const double T0, const double T1, const int para) const ;
     
     virtual void calculateMullinsEffectScaler(const IPNonLinearTVE* q0, IPNonLinearTVE *q1, const double T) const;
@@ -229,6 +229,7 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
 
     virtual void getK(const IPNonLinearTVE *q1, double& KT, const double Tc, double *dK = NULL, double *ddK = NULL) const{
         double K = _K*q1->getConstRefToElasticBulkPropertyScaleFactor();
+        // Msg::Error(" Inside getK, after getting and updating eta; K = %e!!",K);
         KT = K*_temFunc_K->getVal(Tc);
         if(dK!=NULL){
             *dK = K*_temFunc_K->getDiff(Tc);  // Here, *dK can be imagined as *(dK+0) or dK[0] if dK were a vector. See dKi below.
@@ -274,13 +275,12 @@ class mlawNonLinearTVE : public mlawPowerYieldHyper{
 
     // The 3rd argument is a pointer to a vector.
     virtual void getKi(const IPNonLinearTVE *q1, std::vector<double>& Ki_T, const double Tc, std::vector<double>* dKi=NULL, std::vector<double>* ddKi=NULL) const{
-        
         std::vector<double> Ki;
         Ki.clear(); Ki.resize(_N,0.);
         for (int i=0; i<_Ki.size();i++){
           Ki[i] = _Ki[i]*q1->getConstRefToElasticBulkPropertyScaleFactor();  
         }
-        
+        // Msg::Error(" Inside getKi, after getting and updating eta; Ki[0] = %e!!",Ki[0]);
         Ki_T.resize(_N,0.);
         for (int i=0; i<_Ki.size();i++){
             Ki_T[i] = Ki[i]*_temFunc_Ki[i]->getVal(Tc);
diff --git a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
index 873300d9660130da66ba4632e51691ed142aaf97..a67d67427ed9fa9e3bd63fd3394868bc9a997da6 100644
--- a/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLinearTVP.cpp
@@ -364,7 +364,7 @@ void mlawNonLinearTVP::getFreeEnergyTVM(const IPNonLinearTVP *q0, IPNonLinearTVP
     const std::vector<STensor3>& devDDOiDTT = q1->getConstRefToDevDDOiDTT();
     const std::vector<double>& trDDOiDTT = q1->getConstRefToTrDDOiDTT();
         
-    psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);
+    psiVE = mlawNonLinearTVE::freeEnergyMechanical(*q1,T0,T);
     
     DpsiVEdT = 1/9 * (1/KT * trKeinf * DtrKeinfDT - 1/(2*pow(KT,2)) * dKdT * pow(trKeinf,2));
     DDpsiVEdTT = 1/9 * (1/KT * (pow(DtrKeinfDT,2)+trKeinf*DDtrKeinfDTT) - 1/pow(KT,2)*dKdT*trKeinf*DtrKeinfDT 
@@ -3249,7 +3249,7 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
  double CpT;
  if (stiff){
      double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
-     getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
+     // getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
      // CpT = -T*DDpsiTVMdTT;
      getCp(CpT,T);
   }
@@ -3281,14 +3281,21 @@ void mlawNonLinearTVP::predictorCorrector_TVP_nonAssociatedFlow(const STensor3&
   double& psiTVM = q1->getRefToFreeEnergyTVM();
   // getFreeEnergyTVM(q0,q1,T0,T,&psiTVM,NULL);
 
-  q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T);     // deformationEnergy(*q1,T);
+  q1->_elasticEnergy = mlawNonLinearTVE::freeEnergyMechanical(*q1,T0,T);
+  if (q1->_elasticEnergy > q0->_psiMax){
+    q1->_psiMax = q1->_elasticEnergy;
+  }
+  else{
+    q1->_psiMax = q0->_psiMax;
+  }
+  // Msg::Error(" Inside TVP after correction, after updating q1->_elasticEnergy = %e!!",q1->_elasticEnergy);
   
 
   if (stiff){
 
-    const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT();
-    double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
-    getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
+    // const double& DDpsiTVMdTT_0 = q0->getConstRefToDDFreeEnergyTVMdTT();
+    // double& DDpsiTVMdTT = q1->getRefToDDFreeEnergyTVMdTT();
+    // getFreeEnergyTVM(q0,q1,T0,T,NULL,NULL,&DDpsiTVMdTT);
     // double CpT = -T*DDpsiTVMdTT;
     // double CpT_0 = -T0*DDpsiTVMdTT_0;
 
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.cpp b/NonLinearSolver/materialLaw/mullinsEffect.cpp
index 47f32a8b9178899e13a5d6e790d7eb1cd15432a7..14bda0551a31f7bfbe0a14538662792b48734e1c 100644
--- a/NonLinearSolver/materialLaw/mullinsEffect.cpp
+++ b/NonLinearSolver/materialLaw/mullinsEffect.cpp
@@ -85,44 +85,46 @@ void negativeExponentialScaler::setTemperatureFunction_r(const scalarFunction& T
   _temFunc_r = Tfunc.clone();
 }
 
+void negativeExponentialScaler::setTemperatureFunction_m(const scalarFunction& Tfunc){
+  if (_temFunc_m != NULL) delete _temFunc_m;
+  _temFunc_m = Tfunc.clone();
+}
+
 void negativeExponentialScaler::createIPVariable(IPMullinsEffect* &ipv) const
 {
   if(ipv != NULL) delete ipv;
   ipv = new IPMullinsEffect();
 }
 
-void negativeExponentialScaler::mullinsEffectScaling(double psi0,  const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const{
+void negativeExponentialScaler::mullinsEffectScaling(double psi_max,  const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const{
     
   // at unloading intitiation-> eta = 1.; unloading end-> eta = eta_min; reloading -> eta = eta_min;  
     
   double eta0 = ipvprev.getEta();
   
-  double psi = _psi;
-  
   double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
   
-  if (psi<psi0){
-    eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi0-psi)) ) );
+  if (_psi<psi_max && _psi<_psi0){
+    eta = 1. - _r*( 1. - exp( - sqrt(_m*(psi_max-_psi)) ) );
   }
   else{
-    eta = eta0;
+    eta = 1.;
   }
    
   ipv.set(eta,DetaDpsi,DDetaDpsipsi,DetaDT,DDetaDTT);
 }
 
-void negativeExponentialScaler::mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const{
+void negativeExponentialScaler::mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const{
     
   double eta0 = ipvprev.getEta();
    
-  double psi = _psi;
   double r = _r*_temFunc_r->getVal(T);
   double m = _m*_temFunc_m->getVal(T);
   
-  double eta = 0., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
+  double eta = 1., DetaDpsi = 0., DDetaDpsipsi = 0., DetaDT = 0., DDetaDTT = 0.;
   
-  if (psi<psi0){
-    eta = 1. - r*( 1. - exp( - sqrt(m*(psi0-psi)) ) );
+  if (_psi<psi_max && _psi<_psi0){
+    eta = 1. - r*( 1. - exp( - sqrt(m*(psi_max-_psi)) ) );
   }
   else{
     eta = eta0;
diff --git a/NonLinearSolver/materialLaw/mullinsEffect.h b/NonLinearSolver/materialLaw/mullinsEffect.h
index 9e06ee651f89844e61f4e9ad7a13eed49e76e590..dd10cdd9f5d0700c561931d38df10cdc037f50d6 100644
--- a/NonLinearSolver/materialLaw/mullinsEffect.h
+++ b/NonLinearSolver/materialLaw/mullinsEffect.h
@@ -33,14 +33,14 @@ class mullinsEffect {
     virtual void createIPVariable(IPMullinsEffect* &ipv) const=0;
     virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw)=0;
     virtual const bool isInitialized() {return _initialized;}
-    virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const=0;
-    virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const=0;
+    virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const=0;
+    virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const=0;
     virtual mullinsEffect * clone() const=0;
   #endif
 
 };
 
-// eta = 1 - r* ( 1 - exp(-sqrt(m*(psi0_max-psi0))) ) -> Ricker et al, 2021 (modified Zuniga, 2005)
+// eta = 1 - r* ( 1 - exp(-sqrt(m*(psi_max-psi))) ) -> Ricker et al, 2021 (modified Zuniga, 2005)
 class negativeExponentialScaler : public mullinsEffect{
   #ifndef SWIG
   protected:
@@ -61,8 +61,8 @@ class negativeExponentialScaler : public mullinsEffect{
     virtual scalingFunction getType() const{return mullinsEffect::negativeExponentialScaler; };
     virtual void createIPVariable(IPMullinsEffect* &ipv) const;
     virtual void initLaws(const std::map<int,mullinsEffect*> &maplaw) {}; //nothing now, we will see if we use the mapping
-    virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, IPMullinsEffect &ipv) const;
-    virtual void mullinsEffectScaling(double psi0, const IPMullinsEffect &ipvprev, double _psi, double T, IPMullinsEffect &ipv) const;
+    virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, IPMullinsEffect &ipv) const;
+    virtual void mullinsEffectScaling(double psi_max, const IPMullinsEffect &ipvprev, double _psi0, double _psi, double T, IPMullinsEffect &ipv) const;
     virtual mullinsEffect * clone() const;
     #endif // SWIG
 };
diff --git a/dG3D/src/dG3DIPVariable.cpp b/dG3D/src/dG3DIPVariable.cpp
index 88095add61d8a9d8d6859662ac0db7bf1bac5dc8..e48cf6989daaedc28266df9154bd52dcb0036ddb 100644
--- a/dG3D/src/dG3DIPVariable.cpp
+++ b/dG3D/src/dG3DIPVariable.cpp
@@ -5251,6 +5251,18 @@ double NonLinearTVEDG3DIPVariable::get(const int i) const{
   else if (i == IPField::THERMALFLUX_Z){return getConstRefToFlux()[0](2);}
   else if (i == IPField::FIELD_SOURCE){return getConstRefToFieldSource()(0);}
   else if (i == IPField::MECHANICAL_SOURCE){return getConstRefToMechanicalSource()(0);}
+  else if (i == IPField::DEFO_ENERGY)
+  {
+    return _ipViscoElastic->_elasticEnergy;
+  }
+  else if (i == IPField::MAX_DEFO_ENERGY)
+  {
+    return _ipViscoElastic->_psiMax;
+  }
+  else if (i == IPField::BULK_SCALE_FACTOR)
+  {
+    return _ipViscoElastic->getConstRefToElasticBulkPropertyScaleFactor();
+  }
 
   else
   {
@@ -5453,7 +5465,18 @@ double NonLinearTVPDG3DIPVariable::get(const int i) const{
   {
     return _ipViscoElastoPlastic->getConstRefToKinematicHardening().getR();
   }
-
+  else if (i == IPField::DEFO_ENERGY)
+  {
+    return _ipViscoElastoPlastic->_elasticEnergy;
+  }
+  else if (i == IPField::MAX_DEFO_ENERGY)
+  {
+    return _ipViscoElastoPlastic->_psiMax;
+  }
+  else if (i == IPField::BULK_SCALE_FACTOR)
+  {
+    return _ipViscoElastoPlastic->getConstRefToElasticBulkPropertyScaleFactor();
+  }
   else
   {
     return dG3DIPVariable::get(i);
diff --git a/dG3D/src/dG3DMaterialLaw.cpp b/dG3D/src/dG3DMaterialLaw.cpp
index 1c42e2af34e06f851cf51071384597f2f506a206..f7f61561b0884906a0518b5d03286b2b40ca6f50 100644
--- a/dG3D/src/dG3DMaterialLaw.cpp
+++ b/dG3D/src/dG3DMaterialLaw.cpp
@@ -616,7 +616,7 @@ void dG3DMaterialLawWithTangentByPerturbation::stress(IPVariable* ipv, const IPV
         }
 
     }
-#if 1
+#if 0
     ipvcur->getConstRefToDeformationGradient().print("F Numerique"); // FLE
     ipvcur->getConstRefToFirstPiolaKirchhoffStress().print("P Numerique"); // FLE
     ipvcur->getConstRefToTangentModuli().print("dPdE Numerique");