diff --git a/NonLinearSolver/materialLaw/CoalescenceLaw.cpp b/NonLinearSolver/materialLaw/CoalescenceLaw.cpp
index 25928db4ae497fbc90628aae9fa4ea5f14c5f526..97a061c04169f17186b99d82adba5d7589c851d7 100644
--- a/NonLinearSolver/materialLaw/CoalescenceLaw.cpp
+++ b/NonLinearSolver/materialLaw/CoalescenceLaw.cpp
@@ -183,7 +183,7 @@ ThomasonCoalescenceLaw::ThomasonCoalescenceLaw(const int num, const double lamda
   _acceleratedRate(1.),
   _checkWithNormal(false), _onlyInterface(false)
 {
-  _ChiRegularizedFunction = new twoVariableExponentialSaturationScalarFunction(0.98,0.99);
+  _ChiRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-3, 1.e-6, 0.99, 0.999);
 };
 
 
@@ -194,7 +194,7 @@ ThomasonCoalescenceLaw::ThomasonCoalescenceLaw(const int num, const double lamda
   _acceleratedRate(accRate),
   _checkWithNormal(withNormal), _onlyInterface(onlyInt)
 {
- _ChiRegularizedFunction = new twoVariableExponentialSaturationScalarFunction(0.98,0.99);
+ _ChiRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-3, 1.e-6, 0.99, 0.999);
 };
 
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
index 8c0021d7891ee86d70327347de6015e50c984e02..0d31be8346840398f5c12497de35164011f4d4c6 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalDamageGurson.cpp
@@ -27,10 +27,10 @@ mlawNonLocalDamageGurson::mlawNonLocalDamageGurson(const int num,const double E,
   //_correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 0., 0.9*_fV_failure,0.999*_fV_failure);
 
   if (_correctedRegularizedFunction != NULL)  delete _correctedRegularizedFunction;
-  _correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 0., 0.98*_fV_failure,0.999*_fV_failure);
+  _correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.98*_fV_failure,0.999*_fV_failure);
 
   if (_localRegularizedFunction != NULL)  delete _localRegularizedFunction;
-  _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 0., 0.9*_localfVFailure,0.999*_localfVFailure);
+  _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.9*_localfVFailure,0.999*_localfVFailure);
 
   _failedTol = 0.995;
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
index 98cceb5bf9f3a1a3cbec9ddfa1a489f402ca47fb..3ad0b48d8a64f418a56cf7aa26279819a3059373 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorous.cpp
@@ -76,8 +76,8 @@ mlawNonLocalPorosity::mlawNonLocalPorosity(const int num,const double E,const do
   _temFunc_alp =new constantScalarFunction(1.);
   _triFunc_kw = new constantScalarFunction(1.);
   //---------------------------end
-  _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 0., 0.9*_localfVFailure,0.999*_localfVFailure);
-  _correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 0., 0.9*_fV_failure,0.999*_fV_failure);
+  _localRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.9*_localfVFailure,0.999*_localfVFailure);
+  _correctedRegularizedFunction = new fourVariableExponentialSaturationScalarFunction(1.e-6, 1.e-8, 0.9*_fV_failure,0.999*_fV_failure);
 }
 
 mlawNonLocalPorosity::mlawNonLocalPorosity(const mlawNonLocalPorosity &source) :
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
index 21d6743968f7ac6cfb2116732ba21fd6e7e1b47d..6d040b87edaf4fdd69c02bfe9f0ca67eb493598b 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoalescence.cpp
@@ -129,9 +129,12 @@ double mlawNonLocalPorousThomasonLaw::yieldFunction(const double kcorEq, const d
   // Get values in ipvs
   double Cft = q1->getConstRefToIPCoalescence().getConcentrationFactor();
   double CftOffset = q0->getConstRefToIPCoalescence().getCrackOffsetOnCft();
-  if (_CfTOffsetMethod==2 or _CfTOffsetMethod==3){
+  if (_CfTOffsetMethod==2){
     if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffset = 1.;};
-  };
+  }
+  else if (_CfTOffsetMethod==3){
+    if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffset = 1.;};
+  }
 
   double R0 = _j2IH->getYield0();
   if (_roundedMethod == 1){
@@ -172,10 +175,13 @@ void mlawNonLocalPorousThomasonLaw::computeYieldDerivatives( fullVector<double>&
   const IPThomasonCoalescence* q0Thom = dynamic_cast<const IPThomasonCoalescence*>(&(q0->getConstRefToIPCoalescence()));
   double Cft = q1Thom->getConcentrationFactor();
   double CftOffset = q0Thom->getCrackOffsetOnCft();
-  if (_CfTOffsetMethod==2 or _CfTOffsetMethod==3){
+  if (_CfTOffsetMethod==2){
     if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffset = 1.;};
-  };
-
+  }
+  else if (_CfTOffsetMethod==3){
+    if (!q0->getConstRefToIPCoalescence().getCoalescenceOnsetFlag()){CftOffset = 1.;};
+  }
+  
   double DCftDChi = q1Thom->getDConcentrationFactorDLigamentRatio();
   double DCftDW = q1Thom->getDConcentrationFactorDAspectRatio();
 
diff --git a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
index 769e1743c995afc0c7e4008d94ed8ad825f35aca..e8d467314626e37ffe6add5e79d1b587dff8c1b9 100644
--- a/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
+++ b/NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp
@@ -499,9 +499,9 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
             q1Thom->getRefToAccelerateRate() = 1.0;
             
             if (_CfTOffsetMethod == 3){
-              q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
+              q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft() + yieldThomason*R0/(Cft*R);
               
-              if (q1Thom->getRefToCrackOffsetOnCft() > 2.){Msg::Error("High value of Cft Offset = %e",q1Thom->getRefToCrackOffsetOnCft());};
+              if (q1Thom->getRefToCrackOffsetOnCft() > 2.){Msg::Error("High predicted value of Cft Offset = %e",q1Thom->getRefToCrackOffsetOnCft());};
             }
             
           }
diff --git a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
index f04f9714378a2a20cd1979b8ea5c13e0ad6d3253..34b3357a6e3156559e787f1a28a97bf8d1c1deb9 100644
--- a/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
+++ b/dG3D/src/dG3DCohesiveBandMaterialLaw.cpp
@@ -915,32 +915,27 @@ void NonlocalPorosityCohesiveBand3DLaw::checkCohesiveInsertion(IPStateBase* ipsm
 	IPCoalescence& ipcolm = ipvm->getIPNonLocalPorosity()->getRefToIPCoalescence();
 	IPCoalescence& ipcolp = ipvp->getIPNonLocalPorosity()->getRefToIPCoalescence();
 
-  if ((ipcolm.getCoalescenceActiveFlag() or ipcolp.getCoalescenceActiveFlag()) or forcedInsert){
-		printf("rank %d: cohesive element is inserted \n",Msg::GetCommRank());
+  // Check if crack insertion is requested (because coalescence occurs or it is forced to occur)
+  if ((ipcolm.getCoalescenceOnsetFlag() or ipcolp.getCoalescenceOnsetFlag()) or forcedInsert){
+    printf("rank %d: cohesive element is inserted \n",Msg::GetCommRank());
     fMinusCur->broken();
     fPlusCur->broken();
 
-		if (!forcedInsert){
-			if (!ipcolm.getCoalescenceActiveFlag()) {
-				ipcolm.operator=(*dynamic_cast<const IPCoalescence*>(&ipcolp));
-			}
-			else{
-				ipcolp.operator=(*dynamic_cast<const IPCoalescence*>(&ipcolm));
-			}
-		}
-		else{
-			if (_coalescenceLaw != NULL){
-				_coalescenceLaw->forceCoalescence(ipvm->getIPNonLocalPorosity()->getYieldPorosity(),ipvm->getIPNonLocalPorosity(),&ipcolm);
-				_coalescenceLaw->forceCoalescence(ipvp->getIPNonLocalPorosity()->getYieldPorosity(),ipvp->getIPNonLocalPorosity(),&ipcolp);
-			}
-			else
-				Msg::Warning("problem with forcedInsert in PorosityCohesiveBand3DLaw::checkCohesiveInsertion");
-		}
+    // Force IP to be in a coalescence state if not already in that state
+    if (!ipcolm.getCoalescenceOnsetFlag()){
+      _coalescenceLaw->forceCoalescence(ipvm->getIPNonLocalPorosity()->getYieldPorosity(), ipvm->getIPNonLocalPorosity(),&ipcolm);
+    }
+    if (!ipcolp.getCoalescenceOnsetFlag()){
+      _coalescenceLaw->forceCoalescence(ipvp->getIPNonLocalPorosity()->getYieldPorosity(), ipvp->getIPNonLocalPorosity(),&ipcolp);
+    }
 
+    
+    
     // activate nonlocal to local transition
     fMinusCur->setNonLocalToLocal(true);
     fPlusCur->setNonLocalToLocal(true);
 
+
 		// broken is checked via bulk terial law, see function stress
     // initialize cohesive law in negative part
 	  static SVector3 zeroVec;
diff --git a/dG3D/src/dG3DTerms.cpp b/dG3D/src/dG3DTerms.cpp
index 2f48def3eca32ab7eb1bccdc6e96cd84f8c18ef0..0e70e272e662674a9b66836e4383d3668cd8a3fa 100644
--- a/dG3D/src/dG3DTerms.cpp
+++ b/dG3D/src/dG3DTerms.cpp
@@ -909,8 +909,8 @@ void dG3DForceInter::get(MElement *ele, int npts, IntPt *GP, fullVector<double>
             gradm = ipvm->getConstRefToGradNonLocalVariable()[nlk];
             gradp = ipvp->getConstRefToGradNonLocalVariable()[nlk];
             if (_incrementNonlocalBased){
-              gradm = ipvmprev->getConstRefToGradNonLocalVariable()[nlk];
-              gradp = ipvpprev->getConstRefToGradNonLocalVariable()[nlk];
+              gradm -= ipvmprev->getConstRefToGradNonLocalVariable()[nlk];
+              gradp -= ipvpprev->getConstRefToGradNonLocalVariable()[nlk];
             };
 
             SVector3 cgGradm;