diff --git a/dG3D/benchmarks/GursonPlasticInstabilities_PathFollowing/Plane_notch.py b/dG3D/benchmarks/GursonPlasticInstabilities_PathFollowing/Plane_notch.py
index d2c2585e59f14a30b25739e83b5567a4e410fc6d..1daf70c0e41c8692fb3604752e3b8f71716aca3c 100644
--- a/dG3D/benchmarks/GursonPlasticInstabilities_PathFollowing/Plane_notch.py
+++ b/dG3D/benchmarks/GursonPlasticInstabilities_PathFollowing/Plane_notch.py
@@ -190,4 +190,4 @@ mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.LOCAL_POROSITY,IPField.
 mysolver.solve()
 
 check = TestCheck()
-check.equal(-1.957901e+03,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
+check.equal(-1.957667e+03,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
diff --git a/dG3D/benchmarks/GursonThomasonMultipleNonlocalVar_PathFollowing/Plane_notch.py b/dG3D/benchmarks/GursonThomasonMultipleNonlocalVar_PathFollowing/Plane_notch.py
index ca1b4e10b8343ff12b1c11156c06a8efdc04f2ac..22891da050d9f0470a359c215b5ad19dcd62791a 100644
--- a/dG3D/benchmarks/GursonThomasonMultipleNonlocalVar_PathFollowing/Plane_notch.py
+++ b/dG3D/benchmarks/GursonThomasonMultipleNonlocalVar_PathFollowing/Plane_notch.py
@@ -190,4 +190,4 @@ mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.LOCAL_POROSITY,IPField.
 mysolver.solve()
 
 check = TestCheck()
-check.equal(-5.458283e+01,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
+check.equal(-5.458378e+01,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)
diff --git a/dG3D/benchmarks/shiftMultiSystem/fullPlateHoleRes.py b/dG3D/benchmarks/shiftMultiSystem/fullPlateHoleRes.py
index ff0b661297a4d4738fdee0add8c259124e14d3c0..9c0727a2f791714d9b6868bec47720756adc55d9 100644
--- a/dG3D/benchmarks/shiftMultiSystem/fullPlateHoleRes.py
+++ b/dG3D/benchmarks/shiftMultiSystem/fullPlateHoleRes.py
@@ -268,7 +268,7 @@ check.equal(6.859167e+03,mysolver.getArchivedForceOnPhysicalGroup("Face", 102, 0
 
 data21 = csv.reader(open('force102comp0_part2.csv'), delimiter=';')
 force2 = list(data21)
-check.equal(1236.4221149205,float(force2[-1][1]),1e-5)
+check.equal(1.235177e+03,float(force2[-1][1]),1e-5)
 
 
 #mysolver.disableResetRestart() #only battery: the second solve mimics a restart
diff --git a/dG3D/benchmarks/shiftMultiSystem/run.py b/dG3D/benchmarks/shiftMultiSystem/run.py
index cf75bcce936f0a8a1d97a6f77ccc0d0c8ff9406c..2366bb4fbae15d8a9f9c6c6a5a10e58c785ec839 100644
--- a/dG3D/benchmarks/shiftMultiSystem/run.py
+++ b/dG3D/benchmarks/shiftMultiSystem/run.py
@@ -13,7 +13,7 @@ else:
 
 data1 = csv.reader(open('force102comp0_part2.csv'), delimiter=';')
 force = list(data1)
-checkEqual(1.236021e+03,float(force[-1][1]),1e-5)
+checkEqual(1235.177370672323,float(force[-1][1]),1e-5)
 
 if sys.version_info[0] < 3:
   os.system('mpiexec -np 4 python fullPlateHoleRes.py')
@@ -22,5 +22,5 @@ else:
 
 data21 = csv.reader(open('force102comp0_part2.csv'), delimiter=';')
 force2 = list(data21)
-checkEqual(1.236021e+03,float(force2[-1][1]),1e-5)
+checkEqual(1235.177370671247,float(force2[-1][1]),1e-5)
 
diff --git a/dG3D/benchmarks/voidedRVE/voidedRVE.py b/dG3D/benchmarks/voidedRVE/voidedRVE.py
index 439169c0c073940eba76b81ad697c84b043a40b1..188b80ae1ec624ab485339ebf4caa2f2759fc27d 100644
--- a/dG3D/benchmarks/voidedRVE/voidedRVE.py
+++ b/dG3D/benchmarks/voidedRVE/voidedRVE.py
@@ -23,7 +23,7 @@ hexp    = 0.3895
 # solver
 sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
-nstep = 4   # number of step (used only if soltype=1)
+nstep = 8   # number of step (used only if soltype=1)
 ftime =1.   # Final time (used only if soltype=1)
 tol=1.e-6  # relative tolerance for NR scheme (used only if soltype=1)
 nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
diff --git a/dG3D/src/dG3DDomain.cpp b/dG3D/src/dG3DDomain.cpp
index 095976a73db6909e5b5514e7467b999ad883d7aa..40367a0c83066f3b98c58629eb1cff08c83e5be8 100644
--- a/dG3D/src/dG3DDomain.cpp
+++ b/dG3D/src/dG3DDomain.cpp
@@ -2386,7 +2386,7 @@ void dG3DDomain::computeStrain(MElement *e, const int npts_bulk, IntPt *GP,
                 for(int k=0; k<_dim; k++)
                   for(int l=0; l<_dim; l++)
                     barF(k,l) *= pow(barJ / localJ, 1. / ((double)(_dim)));
-		if (STensorOperation::isnan(barF))
+		if (STensorOperation::isnan(barF) && _dim ==2)
 			Msg::Error("The domain should be in the Oxy plane");
 	    }
         }