From 13ffc6f8b5bf05faac1678d7e052197a9368ec45 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 15 Jan 2016 16:11:39 +0000
Subject: [PATCH] more optimizations

---
 Mesh/delaunay3d.cpp          |  2 +-
 Mesh/delaunay_refinement.cpp | 11 ++++++-----
 Mesh/meshGRegion.cpp         |  2 +-
 benchmarks/3d/Cube-01.geo    |  1 +
 4 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 6875522ae4..7ba803bdd4 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 faa304de63..6046d2e66c 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 040db18bfb..352edc18aa 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 ebb2f0b2d8..239b5af5e4 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 };
-- 
GitLab