new 3D mesher is now efficient (in serial)

CMakeLists.txt

@@ -1312,7 +1312,7 @@ endif(NOWARN)
 # disable compile optimization on some known problematic files
 check_cxx_compiler_flag("-O0" NOOPT)
 if(NOOPT)
-  file(GLOB_RECURSE NOOPT_SRC Numeric/robustPredicates.cpp Mesh/BDS.cpp
+  file(GLOB_RECURSE NOOPT_SRC Mesh/BDS.cpp
                               Parser/Gmsh.tab.cpp contrib/Tetgen*/predicates.cxx)
   foreach(FILE ${NOOPT_SRC})
     get_source_file_property(PROP ${FILE} COMPILE_FLAGS)

Mesh/delaunay3d.cpp

@@ -346,7 +346,6 @@ static void delaunayCavity2 (Tet *t,
 }
 
 static void delaunayCavity (Tet *t,
-			    Tet *prev,
 			    Vertex *v,
 			    cavityContainer &cavity,
 			    connContainer &bnd,
@@ -361,7 +360,7 @@ static void delaunayCavity (Tet *t,
     }
     else if (neigh->inSphere(v,thread)){
       if (!(neigh->isSet(thread, iPnt))) {
-	delaunayCavity (neigh, t, v, cavity,bnd,thread, iPnt); 
+	delaunayCavity (neigh, v, cavity,bnd,thread, iPnt); 
       }
     }
     else {
@@ -529,7 +528,7 @@ void delaunayTrgl (const unsigned int numThreads,
   double totCavityGlob=0;
 #endif
 
-  double t1,t2=0,t3=0,t4=0;
+  //  double t1,t2=0,t3=0,t4=0;
 
   //  checkLocalDelaunayness(allocator, 0, "initial");
 
@@ -571,7 +570,8 @@ void delaunayTrgl (const unsigned int numThreads,
 #ifdef _HAVE_NUMA
       allocatedVerts [K] = (Vertex*)numa_alloc_local (locSizeK[K]*sizeof(Vertex));
 #else
-      allocatedVerts [K] = (Vertex*)calloc (locSizeK[K],sizeof(Vertex));
+      //      allocatedVerts [K] = (Vertex*)calloc (locSizeK[K],sizeof(Vertex));
+      allocatedVerts [K] = new Vertex [locSizeK[K]];
 #endif
       for (unsigned int iP=0 ; iP < locSizeK[K] ; iP++){
 	allocatedVerts[K][iP] = *(assignTo[K+myThread*NPTS_AT_ONCE][iP]);
@@ -591,7 +591,7 @@ void delaunayTrgl (const unsigned int numThreads,
       std::vector<Tet*> t(NPTS_AT_ONCE);
       //	  double c1 = Cpu();
       // FIND SEEDS
-      t1 = Cpu();
+      //      t1 = Cpu();
       for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
 	vToAdd[K] = iPGlob <  locSizeK[K] ? &allocatedVerts[K][iPGlob] : NULL;
 	if(vToAdd[K]){
@@ -606,10 +606,10 @@ void delaunayTrgl (const unsigned int numThreads,
 	  }
 	}
       }
-      t2+= Cpu() - t1;
+      //      t2+= Cpu() - t1;
       //      double c1 = Cpu();
       // BUILD CAVITIES
-      t1 = Cpu();
+      //      t1 = Cpu();
       for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
 	if(vToAdd[K]){
 	  cavityContainer &cavityK = cavity[K];
@@ -618,6 +618,7 @@ void delaunayTrgl (const unsigned int numThreads,
 	  for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
 	  cavityK.clear(); bndK.clear();
 	  delaunayCavity2(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K);
+	  //delaunayCavity(t[K], vToAdd[K], cavityK, bndK, myThread, K);
 	  // verify the cavity
 	  for (unsigned int i=0; i< bndK.size(); i++) {
 	    const double val =   robustPredicates::orient3d((double*)bndK[i].f.V[0],
@@ -633,7 +634,7 @@ void delaunayTrgl (const unsigned int numThreads,
 	}
       }
 
-      t3 += Cpu() - t1;
+      //      t3 += Cpu() - t1;
 
 #pragma omp barrier
 
@@ -641,7 +642,7 @@ void delaunayTrgl (const unsigned int numThreads,
 	if (!vToAdd[K])ok[K]=false;
 	else ok[K] = canWeProcessCavity (cavity[K], myThread, K);
       }
-      t1 = Cpu();
+      //      t1 = Cpu();
 
       for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
 	if (ok[K]){
@@ -676,7 +677,7 @@ void delaunayTrgl (const unsigned int numThreads,
 	  }
 	}
       }
-      t4 += Cpu() - t1;
+      //      t4 += Cpu() - t1;
     }
 #ifdef _VERBOSE
     #pragma omp critical
@@ -698,7 +699,7 @@ void delaunayTrgl (const unsigned int numThreads,
   //checkLocalDelaunayness(T,"final");
 
 
-  printf(" %12.5E %12.5E  %12.5E tot  %12.5E \n",t2,t3,t4,t2+t3+t4);
+  //  printf(" %12.5E %12.5E  %12.5E tot  %12.5E \n",t2,t3,t4,t2+t3+t4);
 
 #ifdef _VERBOSE
   printf("average searches per point  %12.5E\n",totSearchGlob/Npts);

Mesh/delaunay3d_private.h

@@ -10,6 +10,9 @@
 #include <math.h>
 #include "robustPredicates.h"
 #include <stdio.h>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
 
 #ifndef MAX_NUM_THREADS_
 #define MAX_NUM_THREADS_ 1
@@ -292,7 +295,16 @@ class tetContainer {
   tetContainer (int nbThreads, int preallocSizePerThread) {
     // FIXME !!!
     if (nbThreads != 1) throw;
-    _perThread.push_back (new aBunchOfStuff<Tet>    (preallocSizePerThread) );
+    _perThread.resize(nbThreads);
+#pragma omp parallel num_threads(nbThreads)
+    {
+#ifdef _OPENMP
+      int  myThread = omp_get_thread_num();
+#else
+      int  myThread = 0;
+#endif
+      _perThread [myThread] = new aBunchOfStuff<Tet>    (preallocSizePerThread) ;
+    }
     //    _perThreadV.push_back(new aBunchOfStuff<Vertex> (preallocSizePerThread) );
   }
   inline Tet * newTet(int thread) {

Mesh/delaunay_refinement.cpp

@@ -311,8 +311,7 @@ void edgeBasedRefinement (const int numThreads,
 
   // fill up old Datastructures
 
-  tetContainer allocator (numThreads,3000000);
-
+  tetContainer allocator (numThreads,1000000);
 
   std::vector<Vertex *> _vertices;
   edgeContainer ec;
@@ -328,13 +327,13 @@ void edgeBasedRefinement (const int numThreads,
     }
 
 
-    FILE *f = fopen ("pts_init.dat","w");
-    fprintf(f,"%d\n",all.size());
-    for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
-      MVertex *mv = *it;      
-      fprintf(f,"%12.5E %12.5E %12.5E\n",mv->x(),mv->y(),mv->z());
-    }
-    fclose(f);
+    //    FILE *f = fopen ("pts_init.dat","w");
+    //    fprintf(f,"%d\n",all.size());
+    //    for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
+    //      MVertex *mv = *it;      
+    //      fprintf(f,"%12.5E %12.5E %12.5E\n",mv->x(),mv->y(),mv->z());
+    //    }
+    //    fclose(f);
 
 
     _vertices.resize(all.size());
@@ -347,22 +346,25 @@ void edgeBasedRefinement (const int numThreads,
       _ma[v] = mv;
       counter++;
     }
-    connSet faceToTet;
-    // FIXME MULTITHREADING
-    for (unsigned int i=0;i<T.size();i++){
-      MTetrahedron  *tt = T[i];
-      int i0 = tt->getVertex(0)->getIndex();
-      int i1 = tt->getVertex(1)->getIndex();
-      int i2 = tt->getVertex(2)->getIndex();
-      int i3 = tt->getVertex(3)->getIndex();
-      Tet *t = allocator.newTet(0) ; t->setVertices (_vertices[i0],_vertices[i1],_vertices[i2],_vertices[i3]);
-      computeAdjacencies (t,0,faceToTet);
-      computeAdjacencies (t,1,faceToTet);
-      computeAdjacencies (t,2,faceToTet);
-      computeAdjacencies (t,3,faceToTet);
-      delete tt;
+
+    {
+      connSet faceToTet;
+      // FIXME MULTITHREADING
+      for (unsigned int i=0;i<T.size();i++){
+	MTetrahedron  *tt = T[i];
+	int i0 = tt->getVertex(0)->getIndex();
+	int i1 = tt->getVertex(1)->getIndex();
+	int i2 = tt->getVertex(2)->getIndex();
+	int i3 = tt->getVertex(3)->getIndex();
+	Tet *t = allocator.newTet(0) ; t->setVertices (_vertices[i0],_vertices[i1],_vertices[i2],_vertices[i3]);
+	computeAdjacencies (t,0,faceToTet);
+	computeAdjacencies (t,1,faceToTet);
+	computeAdjacencies (t,2,faceToTet);
+	computeAdjacencies (t,3,faceToTet);
+	delete tt;
+      }
+      T.clear();
     }
-    T.clear();
   }
 
   // do not allow to saturate boundary edges
@@ -436,7 +438,7 @@ void edgeBasedRefinement (const int numThreads,
 		(t5-t4),
 		(t5-__t__),
 		allocator.size(0));
-      printf("%d calls to inSphere\n",Tet::in_sphere_counter);
+      //      printf("%d calls to inSphere\n",Tet::in_sphere_counter);
       iter++;
     }
   }

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];