3D Delaunnay-ization of a set of points is now fully gmsh based

Mesh/meshGRegionDelaunayInsertion.cpp

@@ -18,6 +18,7 @@
 #include "MTriangle.h"
 #include "Numeric.h"
 #include "Context.h"
+#include "HilbertCurve.h"
 
 int MTet4::radiusNorm = 2;
 static double LIMIT_ = 1;
@@ -330,15 +331,6 @@ void nonrecurFindCavity(std::list<faceXtet> & shell,
                      MVertex *v ,
                      MTet4 *t)
 {
-  // Msg::Info("tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
-  //     t->tet()->getVertex(1)->getNum(),
-  //     t->tet()->getVertex(2)->getNum(),
-  //     t->tet()->getVertex(3)->getNum());
-
-  // invariant : this one has to be inserted in the cavity
-  // consider this tet deleted
-  // remove its reference to its neighbors
-
   std::stack<MTet4*> _stack;
   _stack.push(t);
   while(!_stack.empty()){
@@ -1612,3 +1604,152 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
   MTet4::radiusNorm = 2;
   LIMIT_ = 1;
 }
+
+
+///// do a 3D delaunay mesh assuming a set of vertices
+
+static void initialCube (std::vector<MVertex*> &v,
+			 MVertex *box[8],
+			 std::vector<MTet4*> &t){
+  SBoundingBox3d bbox ;
+  for (size_t i=0;i<v.size();i++){
+    MVertex *pv = v[i];
+    bbox += SPoint3(pv->x(),pv->y(),pv->z());    
+  }
+  bbox *= 1.3;
+  box[0] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.min().z());
+  box[1] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.min().z());
+  box[2] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.min().z());
+  box[3] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.min().z());
+  box[4] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.max().z());
+  box[5] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.max().z());
+  box[6] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.max().z());
+  box[7] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.max().z());
+  std::vector<MTetrahedron*> t_box;
+  MTetrahedron *t0 = new MTetrahedron (box[2],box[7],box[3],box[1]);
+  MTetrahedron *t1 = new MTetrahedron (box[0],box[7],box[1],box[3]);
+  MTetrahedron *t2 = new MTetrahedron (box[6],box[1],box[7],box[2]);
+  MTetrahedron *t3 = new MTetrahedron (box[0],box[1],box[7],box[4]);
+  MTetrahedron *t4 = new MTetrahedron (box[1],box[4],box[5],box[7]);
+  MTetrahedron *t5 = new MTetrahedron (box[1],box[7],box[5],box[6]);
+  t.push_back(new MTet4(t0,0.0));
+  t.push_back(new MTet4(t1,0.0));
+  t.push_back(new MTet4(t2,0.0));
+  t.push_back(new MTet4(t3,0.0));
+  t.push_back(new MTet4(t4,0.0));
+  t.push_back(new MTet4(t5,0.0));
+  connectTets(t);
+}
+
+int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
+                            double P[3], double N[3], double);
+
+static MTet4* search4Tet (MTet4 *t, MVertex *v, int _size) {
+  if (t->inCircumSphere(v)) return t;
+  int ITER = 0;
+  SPoint3 p2 (v->x(),v->y(),v->z());
+  while (1){
+    SPoint3 p1 = t->tet()->barycenter();
+    int found = -1;
+    for (int i = 0; i < 4; i++){
+      MTet4 *neigh = t->getNeigh(i);
+      if (neigh){
+	faceXtet fxt (t, i);
+	double X[3] = {fxt.v[0]->x(),fxt.v[1]->x(),fxt.v[2]->x()};
+	double Y[3] = {fxt.v[0]->y(),fxt.v[1]->y(),fxt.v[2]->y()};
+	double Z[3] = {fxt.v[0]->z(),fxt.v[1]->z(),fxt.v[2]->z()};
+	//	printf("%d %d\n",i,intersect_line_triangle(X,Y,Z,p1,p2-p1));
+	if (intersect_line_triangle(X,Y,Z,p1,p1-p2, 1.e-12) > 0) {
+	  found = i;
+	  break;
+	}
+      }
+    }
+    if (found < 0){
+      return 0;      
+    }
+    t = t->getNeigh(found);
+    if (t->inCircumSphere(v)) {
+      //      printf("FOUND\n");
+      return t;
+    }
+    if (ITER++ > _size) break;
+  }
+  return 0;
+}
+
+
+MTet4 * getTetToBreak (MVertex *v, std::vector<MTet4*> &t, int &NB_GLOBAL_SEARCH){
+  // last inserted is used as starting point
+  // we know it is not deleted
+  MTet4 *start = t[t.size() - 1];
+  start = search4Tet (start,v,(int)t.size());
+  if (start)return start;
+  NB_GLOBAL_SEARCH++;
+  for (size_t i = 0;i<t.size();i++){
+    if (!t[i]->isDeleted() && t[i]->inCircumSphere(v))return t[i];
+  }  
+  return 0;
+}
+
+bool tetOnBox (MTetrahedron *t, MVertex *box[8]){
+  for (size_t i = 0;i<4;i++)
+    for (size_t j = 0;j<8;j++)
+      if (t->getVertex(i) == box[j])return true;
+  return false;
+}
+
+void delaunayMeshIn3D(std::vector<MVertex*> &v, std::vector<MTetrahedron*> &result)
+{
+  std::vector<MTet4*> t;
+  MVertex *box[8];
+  initialCube (v,box,t);
+
+  int NB_GLOBAL_SEARCH = 0;
+  SortHilbert(v);
+  clock_t t1 = clock();
+
+  for (size_t i=0;i<v.size();i++){
+    MVertex *pv = v[i];    
+    MTet4 * found = getTetToBreak (pv,t,NB_GLOBAL_SEARCH);
+    std::list<faceXtet> shell;
+    std::list<MTet4*> cavity;
+    recurFindCavity(shell, cavity, pv, found);
+    std::vector<MTet4*> extended_cavity;
+    std::list<faceXtet>::iterator it = shell.begin();
+    while (it != shell.end()){
+      MTetrahedron *tr = new MTetrahedron(it->getVertex(0), 
+					  it->getVertex(1), 
+					  it->getVertex(2), pv);
+      MTet4 *t4 = new MTet4(tr, 0.0);
+      t.push_back(t4);
+      extended_cavity.push_back(t4);
+      MTet4 *otherSide = it->t1->getNeigh(it->i1);
+      if (otherSide)
+	extended_cavity.push_back(otherSide);
+      ++it;
+    }
+    //    printf("connecting %i tets\n",extended_cavity.size());
+    connectTets_vector(extended_cavity.begin(),extended_cavity.end());
+  }
+
+  clock_t t2 = clock();
+  printf("%d global searches among %d CPU = %g\n",NB_GLOBAL_SEARCH,v.size(),(double)(t2-t1)/CLOCKS_PER_SEC);
+
+
+  //  FILE *f = fopen ("tet.pos","w");
+  //  fprintf(f,"View \"\"{\n");
+  for (size_t i = 0;i<t.size();i++){
+    if (t[i]->isDeleted() || tetOnBox (t[i]->tet(),box)) delete t[i]->tet();
+    else {
+      result.push_back(t[i]->tet());
+      //      t[i]->tet()->writePOS (f, false,false,true,false,false,false);
+    }
+    delete t[i];
+  }
+  //  fprintf(f,"};\n");
+  //  fclose(f);
+}
+
+
+

Mesh/meshGRegionDelaunayInsertion.h

@@ -173,6 +173,8 @@ class MTet4
 
 void connectTets(std::list<MTet4*> &);
 void connectTets(std::vector<MTet4*> &);
+// IN --> Vertices ----  OUT --> Tets
+void delaunayMeshIn3D(std::vector<MVertex*> &, std::vector<MTetrahedron*> &);
 void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000, bool _classify = true);
 void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
 GRegion *getRegionFromBoundingFaces(GModel *model,

Numeric/CMakeLists.txt

@@ -7,6 +7,7 @@ set(SRC
   Numeric.cpp
     fullMatrix.cpp
   BasisFactory.cpp
+  discreteFrechetDistance.cpp
     nodalBasis.cpp
 	polynomialBasis.cpp
 	pyramidalBasis.cpp
@@ -27,6 +28,7 @@ set(SRC
     GaussQuadraturePyr.cpp
     GaussLegendreSimplex.cpp
   GaussJacobi1D.cpp
+  HilbertCurve.cpp
  robustPredicates.cpp
   mathEvaluator.cpp
   Iso.cpp

Numeric/HilbertCurve.cpp

+#include "SBoundingBox3d.h"
+#include "MVertex.h"
+
+
+struct HilbertSort
+{
+// The code for generating table transgc
+// from: http://graphics.stanford.edu/~seander/bithacks.html.    
+  int transgc[8][3][8];
+  int tsb1mod3[8];
+  int maxDepth;
+  int Limit;
+  SBoundingBox3d bbox ;
+  void ComputeGrayCode(int n);
+  int Split(MVertex** vertices,
+	    int arraysize,int GrayCode0,int GrayCode1,
+	    double BoundingBoxXmin, double BoundingBoxXmax, 
+	    double BoundingBoxYmin, double BoundingBoxYmax, 
+	    double BoundingBoxZmin, double BoundingBoxZmax);
+  void Sort(MVertex** vertices, int arraysize, int e, int d, 
+	   double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax, 
+	   double BoundingBoxZmin, double BoundingBoxZmax, int depth);
+  HilbertSort (int m = 0, int l=1) : maxDepth(m),Limit(l)
+  {
+    ComputeGrayCode(3);
+  }
+  void MultiscaleSortHilbert(MVertex** vertices, int arraysize, 
+			     int threshold, double ratio, int *depth)
+  {
+    int middle;
+    
+    middle = 0;
+    if (arraysize >= threshold) {
+      (*depth)++;
+      middle = arraysize * ratio;
+      MultiscaleSortHilbert(vertices, middle, threshold, ratio, depth);
+    }
+    Sort (&(vertices[middle]),arraysize - middle,0,0,
+	  bbox.min().x(),bbox.max().x(),
+	  bbox.min().y(),bbox.max().y(),
+	  bbox.min().z(),bbox.max().z(),0);
+  }
+  
+  void Apply (std::vector<MVertex*> &v)
+  {
+    for (size_t i=0;i<v.size();i++){
+      MVertex *pv = v[i];
+      bbox += SPoint3(pv->x(),pv->y(),pv->z());    
+    }
+    bbox *= 1.01;
+    MVertex**pv = &v[0];
+    int depth;
+    MultiscaleSortHilbert(pv, v.size(), 128, 0.125,&depth);
+  }
+};
+
+
+void HilbertSort::ComputeGrayCode(int n)
+{
+  int gc[8], N, mask, travel_bit;
+  int e, d, f, k, g;
+  int v, c;
+  int i;
+
+  N = (n == 2) ? 4 : 8;
+  mask = (n == 2) ? 3 : 7;
+
+  // Generate the Gray code sequence.
+  for (i = 0; i < N; i++) {
+    gc[i] = i ^ (i >> 1);
+  }
+
+  for (e = 0; e < N; e++) {
+    for (d = 0; d < n; d++) {
+      // Calculate the end point (f).
+      f = e ^ (1 << d);  // Toggle the d-th bit of 'e'.
+      // travel_bit = 2**p, the bit we want to travel. 
+      travel_bit = e ^ f;
+      for (i = 0; i < N; i++) {
+        // // Rotate gc[i] left by (p + 1) % n bits.
+        k = gc[i] * (travel_bit * 2);
+        g = ((k | (k / N)) & mask);
+        // Calculate the permuted Gray code by xor with the start point (e).
+        transgc[e][d][i] = (g ^ e);
+      }
+      //      assert(transgc[e][d][0] == e);
+      //      assert(transgc[e][d][N - 1] == f);
+    } // d
+  } // e
+
+  // Count the consecutive '1' bits (trailing) on the right.
+  tsb1mod3[0] = 0;
+  for (i = 1; i < N; i++) {
+    v = ~i; // Count the 0s.
+    v = (v ^ (v - 1)) >> 1; // Set v's trailing 0s to 1s and zero rest
+    for (c = 0; v; c++) {
+      v >>= 1;
+    }
+    tsb1mod3[i] = c % n;
+  }
+}
+
+
+int HilbertSort::Split(MVertex** vertices,
+		       int arraysize,int GrayCode0,int GrayCode1,
+		       double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax, 
+		       double BoundingBoxZmin, double BoundingBoxZmax)
+{
+  MVertex* swapvert;
+  int axis, d;
+  double split;
+  int i, j;
+
+
+  // Find the current splitting axis. 'axis' is a value 0, or 1, or 2, which 
+  //   correspoding to x-, or y- or z-axis.
+  axis = (GrayCode0 ^ GrayCode1) >> 1; 
+
+  // Calulate the split position along the axis.
+  if (axis == 0) {
+    split = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+  } else if (axis == 1) {
+    split = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+  } else { // == 2
+    split = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+  }
+
+  // Find the direction (+1 or -1) of the axis. If 'd' is +1, the direction
+  //   of the axis is to the positive of the axis, otherwise, it is -1.
+  d = ((GrayCode0 & (1<<axis)) == 0) ? 1 : -1;
+
+
+  // Partition the vertices into left- and right-arrays such that left points
+  //   have Hilbert indices lower than the right points.
+  i = 0;
+  j = arraysize - 1;
+
+  // Partition the vertices into left- and right-arrays.
+  if (d > 0) {
+    do {
+      for (; i < arraysize; i++) {      
+        if (vertices[i]->point()[axis] >= split) break;
+      }
+      for (; j >= 0; j--) {
+        if (vertices[j]->point()[axis] < split) break;
+      }
+      // Is the partition finished?
+      if (i == (j + 1)) break;
+      // Swap i-th and j-th vertices.
+      swapvert = vertices[i];
+      vertices[i] = vertices[j];
+      vertices[j] = swapvert;
+      // Continue patitioning the array;
+    } while (true);
+  } else {
+    do {
+      for (; i < arraysize; i++) {      
+        if (vertices[i]->point()[axis] <= split) break;
+      }
+      for (; j >= 0; j--) {
+        if (vertices[j]->point()[axis] > split) break;
+      }
+      // Is the partition finished?
+      if (i == (j + 1)) break;
+      // Swap i-th and j-th vertices.
+      swapvert = vertices[i];
+      vertices[i] = vertices[j];
+      vertices[j] = swapvert;
+      // Continue patitioning the array;
+    } while (true);
+  }
+
+  return i;
+}
+
+// The sorting code is inspired by Tetgen 1.5
+
+void HilbertSort::Sort(MVertex** vertices, int arraysize, int e, int d, 
+		       double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax, 
+		       double BoundingBoxZmin, double BoundingBoxZmax, int depth)
+{
+  double x1, x2, y1, y2, z1, z2;
+  int p[9], w, e_w, d_w, k, ei, di;
+  int n = 3, mask = 7;
+
+  p[0] = 0;
+  p[8] = arraysize;
+
+  p[4] = Split(vertices, p[8], transgc[e][d][3], transgc[e][d][4], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+  p[2] = Split(vertices, p[4], transgc[e][d][1], transgc[e][d][2], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+  p[1] = Split(vertices, p[2], transgc[e][d][0], transgc[e][d][1], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+  p[3] = Split(&(vertices[p[2]]), p[4] - p[2], 
+	       transgc[e][d][2], transgc[e][d][3], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[2];
+  p[6] = Split(&(vertices[p[4]]), p[8] - p[4], 
+	       transgc[e][d][5], transgc[e][d][6], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[4];
+  p[5] = Split(&(vertices[p[4]]), p[6] - p[4], 
+	       transgc[e][d][4], transgc[e][d][5], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[4];
+  p[7] = Split(&(vertices[p[6]]), p[8] - p[6], 
+	       transgc[e][d][6], transgc[e][d][7], 
+	       BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[6];
+
+  if (maxDepth > 0) {
+    if ((depth + 1) == maxDepth) {
+      return;
+    }
+  }
+
+  // Recursively sort the points in sub-boxes.
+  for (w = 0; w < 8; w++) {
+    if ((p[w+1] - p[w]) > Limit) {
+      if (w == 0) {
+        e_w = 0;
+      } else {
+        k = 2 * ((w - 1) / 2); 
+        e_w = k ^ (k >> 1); 
+      }
+      k = e_w;
+      e_w = ((k << (d+1)) & mask) | ((k >> (n-d-1)) & mask);
+      ei = e ^ e_w;
+      if (w == 0) {
+        d_w = 0;
+      } else {
+        d_w = ((w % 2) == 0) ? tsb1mod3[w - 1] : tsb1mod3[w];
+      }
+      di = (d + d_w + 1) % n;
+      if (transgc[e][d][w] & 1) {
+        x1 = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+        x2 = BoundingBoxXmax;
+      } else {
+        x1 = BoundingBoxXmin;
+        x2 = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+      }
+      if (transgc[e][d][w] & 2) { // y-axis
+        y1 = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+        y2 = BoundingBoxYmax;
+      } else {
+        y1 = BoundingBoxYmin;
+        y2 = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+      }
+      if (transgc[e][d][w] & 4) { // z-axis
+        z1 = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+        z2 = BoundingBoxZmax;
+      } else {
+        z1 = BoundingBoxZmin;
+        z2 = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+      }
+      Sort(&(vertices[p[w]]), p[w+1] - p[w], ei, di, 
+                    x1, x2, y1, y2, z1, z2, depth+1);
+    } 
+  } 
+}
+
+
+
+void SortHilbert (std::vector<MVertex*>& v){
+  HilbertSort h;
+  h.Apply(v);
+}

Numeric/HilbertCurve.h

+#ifndef _HILBERT_CURVE_
+#define _HILBERT_CURVE_
+void SortHilbert (std::vector<MVertex*>&);
+#endif