diff --git a/NonLinearSolver/internalPoints/ipField.cpp b/NonLinearSolver/internalPoints/ipField.cpp
index 922e0fdd50d563c98d63c07ee68fd0e32a3ea7f4..ef374c7d9c7a8a81eaa47b1a3447b447584b0c8d 100644
--- a/NonLinearSolver/internalPoints/ipField.cpp
+++ b/NonLinearSolver/internalPoints/ipField.cpp
@@ -1160,6 +1160,11 @@ std::string IPField::ToString(const int i){
   else if (i == DEPLAST_ZZ_DEPSILON_XX) return "DEPLAST_ZZ_DEPSILON_XX";
   else if (i == DEPLAST_XY_DEPSILON_XX) return "DEPLAST_XY_DEPSILON_XX";
   
+  //strain increments (added on 21/05/2025)
+  else if (i == DEPSILON_XX) return "DEPSILON_XX";
+  else if (i == DEPSILON_YY) return "DEPSILON_YY";
+  else if (i == DEPSILON_XY) return "DEPSILON_XY";
+  
   else{
     Msg::Warning("This IP field %d is not defined in IPField::ToString(%d)",i,i);
     return "UNDEFINED";
diff --git a/NonLinearSolver/internalPoints/ipField.h b/NonLinearSolver/internalPoints/ipField.h
index c19e5143813beb01aab955ba95b145b295ea7bb4..3b0817f82627df830d8e0dfb72111f5264984ca8 100644
--- a/NonLinearSolver/internalPoints/ipField.h
+++ b/NonLinearSolver/internalPoints/ipField.h
@@ -358,6 +358,8 @@ class IPField : public elementsField {
                     DSIGMAREC_XX_DEPSILON_XX, DSIGMAREC_YY_DEPSILON_XX,DSIGMAREC_ZZ_DEPSILON_XX,DSIGMAREC_XY_DEPSILON_XX,
                     //DEplastDE (we start in 1D so only 4 components for the plastic strain: Eplast_XX, Eplast_YY, Eplast_ZZ, Eplast_XY and 1 component for epsilon: epsilon_XX)
                     DEPLAST_XX_DEPSILON_XX, DEPLAST_YY_DEPSILON_XX, DEPLAST_ZZ_DEPSILON_XX, DEPLAST_XY_DEPSILON_XX,
+                    //strain increments
+                    DEPSILON_XX, DEPSILON_YY, DEPSILON_XY,
                     UNDEFINED};
     enum  Operator { MEAN_VALUE=1, MIN_VALUE, MAX_VALUE, CRUDE_VALUE};
     static std::string ToString(const int i);
diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.cpp b/dG3D/src/nonLocalDamageDG3DIPVariable.cpp
index 4082b2d30b85e193f015ed2c14de8e7117ef434b..e8a6b4c02baf2567b3e06e03db93efe8a3f40680 100644
--- a/dG3D/src/nonLocalDamageDG3DIPVariable.cpp
+++ b/dG3D/src/nonLocalDamageDG3DIPVariable.cpp
@@ -237,7 +237,7 @@ void nonLocalDamageDG3DIPVariable::restart()
 /* Beginning IPvariable for nonlocalDamageTorchANNDG3DIPVariable*/
 nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable(const std::vector< CLengthLaw *> cll, const int n_stress, const int n_elTang, const int n_NonDamTang, const double initial_h, const int numNonlocal, const bool createBodyForceHO, const bool oninter):
       dG3DIPVariable(createBodyForceHO, oninter,numNonlocal), _n_stress(n_stress), _n_elTang(n_elTang), _n_NonDamTang(n_NonDamTang), _numNonlocal(numNonlocal), _initial_h(initial_h), _nlD(1,6), _D(1,6), _nlD36comp(6,6), _D36comp(6,6),
-      _kinematicVariables(1,6),   
+      _kinematicVariables(1,6), _dExx(0.), _dEyy(0.), _dExy(0.),  
       _defoEnergy(0.), _plasticEnergy(0.), _damageEnergy(0.), _irreversibleEnergy(0.), _DirrEnergyDF(0.), _DirrEnergyDVar_nl0(0.), _DirrEnergyDVar_nl1(0.)
 
 {
@@ -262,7 +262,7 @@ nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::nonLocalDamageTorchANNDG3DIP
     _D(source._D), _nlD36comp(source._nlD36comp), 
     _D36comp(source._D36comp), _DTensor43(source._DTensor43), _nlDTensor43(source._nlDTensor43), 
      _internalVars_stress(source._internalVars_stress), _internalVars_elTang(source._internalVars_elTang), _internalVars_NonDamTang(source._internalVars_NonDamTang), 
-    _kinematicVariables(source._kinematicVariables),_irreversibleEnergy(source._irreversibleEnergy), _DirrEnergyDF(source._DirrEnergyDF),_DirrEnergyDVar_nl0(source._DirrEnergyDVar_nl0), _DirrEnergyDVar_nl1(source._DirrEnergyDVar_nl1), _plasticEnergy(source._plasticEnergy), 
+    _kinematicVariables(source._kinematicVariables),_irreversibleEnergy(source._irreversibleEnergy), _DirrEnergyDF(source._DirrEnergyDF),_DirrEnergyDVar_nl0(source._DirrEnergyDVar_nl0), _DirrEnergyDVar_nl1(source._DirrEnergyDVar_nl1), _plasticEnergy(source._plasticEnergy), _dExx(source._dExx), _dEyy(source._dEyy), _dExy(source._dExy),
      _defoEnergy(source._defoEnergy), _damageEnergy(source._damageEnergy), _ipMeca(NULL), _DirrEnergyDVar_nl(source._DirrEnergyDVar_nl)
 
 {
@@ -332,6 +332,9 @@ nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable& nonLocalDamageTorchANNDG3DIP
     _Cel = psrc->_Cel;   //elastic stiffness without damage
     _DS_nlDE = psrc->_DS_nlDE; //DStressreconstructedDStrain
     _DEplastDE = psrc->_DEplastDE;   //DEplastDE
+    _dExx = psrc->_dExx;   //strain increment following xx (added on 21/05/25)
+    _dEyy = psrc->_dEyy;   //strain increment following yy (added on 21/05/25)
+    _dExy = psrc->_dExy;   //strain increment following xy (added on 21/05/25)
     _defoEnergy = psrc -> _defoEnergy; //added for energy-based path Following
     _plasticEnergy = psrc -> _plasticEnergy; //added for energy-based path Following
     _damageEnergy = psrc -> _damageEnergy; //added for energy-based path Following    
@@ -730,6 +733,22 @@ double nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable::get(const int comp) c
   {
     return _plasticEnergy;
   }
+  else if (comp == IPField::PLASTIC_ENERGY)
+  {
+    return _plasticEnergy;
+  }
+  else if (comp == IPField::DEPSILON_XX)
+  {
+    return _dExx;
+  }
+  else if (comp == IPField::DEPSILON_YY)
+  {
+    return _dEyy;
+  }
+  else if (comp == IPField::DEPSILON_XY)
+  {
+    return _dExy;
+  }
   else
   {
     return dG3DIPVariable::get(comp);
diff --git a/dG3D/src/nonLocalDamageDG3DIPVariable.h b/dG3D/src/nonLocalDamageDG3DIPVariable.h
index 406266045aeb095dafb4659a4418b6339091c0fb..95f784df80b3404c978d13df990e3243ab09a6a0 100644
--- a/dG3D/src/nonLocalDamageDG3DIPVariable.h
+++ b/dG3D/src/nonLocalDamageDG3DIPVariable.h
@@ -1592,6 +1592,10 @@ class nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable : public dG3DIPVariable
     //DStressreconstructedDStrain (maybe this variable can be removed because defined elsewhere)
     STensor43 _DS_nlDE;
     
+    double _dExx; //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    double _dEyy; //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    double _dExy; //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    
     double _defoEnergy;
     double _plasticEnergy;        //added for energy-based path Following
     double _damageEnergy;
@@ -1738,6 +1742,13 @@ class nonLocalDamageTorchANNDG3DIPVariableDG3DIPVariable : public dG3DIPVariable
     virtual double & getRefToDamageEnergy() {return _damageEnergy;};   //added for energy-based path Following
     virtual const double & getConstRefToDamageEnergy() const {return _damageEnergy;};   //added for energy-based path Following
  //////////////////////////////////////////////////////////////////////////////////
+ 
+    virtual const double & getConstRefTodExx() const {return _dExx;};   //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    virtual double & getRefTodExx() {return _dExx;};   
+    virtual const double & getConstRefTodEyy() const {return _dEyy;};   //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    virtual double & getRefTodEyy() {return _dEyy;}; 
+    virtual const double & getConstRefTodExy() const {return _dExy;};   //added on 21/05/25 to check if the strain increment is not too big for the RNN
+    virtual double & getRefTodExy() {return _dExy;}; 
    
     virtual const std::vector<IPCLength*> &getConstRefToIPCLength() const
     {
diff --git a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
index cc869bf0448bd960517a25c89da5d3e941bcb8f2..4e11be29bf28b712d3048d48ab14dc61f6ac9add 100644
--- a/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
+++ b/dG3D/src/nonLocalDamageDG3DMaterialLaw.cpp
@@ -3016,6 +3016,25 @@ void NonlocalDamageTorchANNBasedDG3DMaterialLaw::stress(IPVariable* ipv, const I
   STensor3 E1Tensor;
   STensorOperation::fromFullVectorToSTensor3(E1Vector, E1Tensor);
   
+  //beginning step -1: we check if the strain increments are below the limit increment that the RNN can handle (below 4e-4 for each component for the specifics RNN we have trained so far, this value could be higher with others trained RNN), for the moment it is done for a 2D RNN with only epsilon_XX, epsilon_YY and epsilon_XY as inputs (added on 21/05/25)
+  double incr_limit=0.0004;
+  double& dExx=ipvcur->getRefTodExx();
+  double& dEyy=ipvcur->getRefTodEyy();
+  double& dExy=ipvcur->getRefTodExy();
+  dExx=E1(0,0)-E0(0,0); 
+  dEyy=E1(0,1)-E0(0,1);
+  dExy=E1(0,3)-E0(0,3); 
+  if (abs(dExx)>incr_limit){
+    std::cout << "Warning: beyond the limit increment for epsilon_XX" << "\n";
+  }
+  if (abs(dEyy)>incr_limit){
+    std::cout << "Warning: beyond the limit increment for epsilon_YY" << "\n";
+  }
+  if (abs(dExy)>incr_limit){
+    std::cout << "Warning: beyond the limit increment for epsilon_XY" << "\n";
+  }
+  //end step -1
+  
   /*if (E1Tensor(0,0)>=0.04){
     std::cout << "Exx = " << E1Tensor(0,0) << "\n";   //we check if Exx becomes out of the RNN training range
   }*/