diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 6875522ae4eb3c55c7f7020d4e6be35e989f9e4e..7ba803bdd4e3135286578067a6f957a71f7829a6 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -809,7 +809,7 @@ void delaunayTriangulation (const int numThreads,
   S.clear();
   delaunayTrgl(1,1, assignTo0[0].size(),assignTo0,allocator);
   delaunayTrgl(numThreads,nptsatonce, N,&assignTo[0], allocator);
-  __print ("finalTetrahedrization.pos",0, allocator);
+  //  __print ("finalTetrahedrization.pos",0, allocator);
 }
 
 
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index faa304de63eb37ea42a8667dab2b63e29216fe69..6046d2e66c735a56768393ded3fa4b06cbb81045 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -57,9 +57,9 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
 				double lc1 , double lc2 ,
 				double (*f)(const SPoint3 &p, void *),
 				void *data, std::vector< IPT > & _result,
-				double &dl, double epsilon = 1.e-5)
+				double &dl, std::stack<IPT> &_stack, double epsilon = 1.e-5)
 {
-  std::stack< IPT > _stack;
+  //  _stack.clear();
   _result.clear();
   // local parameters on the edge
   double t1 = 0.0;
@@ -112,12 +112,12 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
 }
 
 
-void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &p, void *), void *data) {
+void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &p, void *), void *data, std::stack<IPT> &temp) {
   std::vector< IPT > _result;
   double dl;
   SPoint3 p1 = e.first->point();
   SPoint3 p2 = e.second->point();
-  const double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl);
+  const double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl, temp);
   const int N = (int) (dN+0.1);
   const double interval = dN/N;
   double L = 0.0;
@@ -160,6 +160,7 @@ void saturateEdges ( edgeContainer &ec,
 		     int nbThreads,
 		     std::vector<Vertex*> &S,
 		     double (*f)(const SPoint3 &p, void *), void *data) {
+  std::stack<IPT> temp;
   AVGSEARCH= 0;
   // FIXME
   const int N = T.size(0);
@@ -171,7 +172,7 @@ void saturateEdges ( edgeContainer &ec,
 	Edge ed = t->getEdge(iEdge);
 	bool isNew = ec.addNewEdge(ed);
 	if (isNew){
-	  saturateEdge (ed, S, f, data);
+	  saturateEdge (ed, S, f, data, temp);
 	}
       }
     }
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 040db18bfbda297690dd7de31a7ab708d35b0e5b..352edc18aac790f7096feeb1d40137bb2bd4c0f4 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -592,7 +592,7 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> &regions)
 }
 
 // uncomment this to test the new code
-//#define NEW_CODE
+#define NEW_CODE
 
 static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions) {
   GRegion *gr = regions[0];
diff --git a/benchmarks/3d/Cube-01.geo b/benchmarks/3d/Cube-01.geo
index ebb2f0b2d8e7c947ba0ba7d2fa524d5d3f3df10d..239b5af5e4ba4565d9e0d830c9b15000f54d6965 100644
--- a/benchmarks/3d/Cube-01.geo
+++ b/benchmarks/3d/Cube-01.geo
@@ -18,3 +18,4 @@ Line Loop(5) = {2,3,4,1};
 Plane Surface(6) = {5};         
 Extrude Surface { 6, {0,0.0,1} };         
 //Characteristic Length {10} = lc/100;
+Physical Point(1)={1 };