switch plugin to in-house Delaunay (still slower, but we'll improve)

Plugin/Tetrahedralize.cpp

@@ -9,8 +9,8 @@
 #include "GmshMessage.h"
 #include "Tetrahedralize.h"
 
-#if defined(HAVE_TETGEN)
-#include "tetgen.h"
+#if defined(HAVE_MESH)
+#include "meshGRegionDelaunayInsertion.h"
 #endif
 
 StringXNumber TetrahedralizeOptions_Number[] = {
@@ -45,13 +45,12 @@ StringXNumber *GMSH_TetrahedralizePlugin::getOption(int iopt)
 
 class PointData {
  public:
-  std::vector<double> v;
+  MVertex *v;
+  std::vector<double> val;
   PointData(double x, double y, double z, int numVal)
   {
-    v.resize(3 + numVal);
-    v[0] = x;
-    v[1] = y;
-    v[2] = z;
+    v = new MVertex(x, y, z);
+    val.resize(numVal);
   }
 };
 
@@ -81,7 +80,7 @@ PView *GMSH_TetrahedralizePlugin::execute(PView *v)
       PointData p(x, y, z, numComp * numSteps);
       for(int step = 0; step < numSteps; step++)
         for(int comp = 0; comp < numComp; comp++)
-          data1->getValue(step, ent, ele, 0, comp, p.v[3 + numComp * step + comp]);
+          data1->getValue(step, ent, ele, 0, comp, p.val[numComp * step + comp]);
       points.push_back(p);
     }
   }
@@ -91,76 +90,80 @@ PView *GMSH_TetrahedralizePlugin::execute(PView *v)
     return v1;
   }
 
-#if !defined(HAVE_TETGEN)
-  Msg::Error("Gmsh has to be compiled with Tetgen support to run "
+#if !defined(HAVE_MESH)
+  Msg::Error("Gmsh has to be compiled with Mesh support to run "
              "Plugin(Tetrahedralize)");
   return v1;
 #else
-  // fill tetgen structure
-  tetgenio in, out;
-  in.mesh_dim = 3;
-  in.numberofpoints = points.size();
-  in.pointlist = new REAL[in.numberofpoints * 3];
+  std::vector<MVertex*> vertices(points.size());
   for(unsigned int i = 0; i < points.size(); i++){
-    in.pointlist[i * 3 + 0] = points[i].v[0];
-    in.pointlist[i * 3 + 1] = points[i].v[1];
-    in.pointlist[i * 3 + 2] = points[i].v[2];
-  }
-  try{
-    tetrahedralize((char*)"Q", &in, &out);
-  }
-  catch (int error){
-    Msg::Error("Could not tetrahedralize points");
-    return v1;
-  }
-
-  if(out.numberofpoints != (int)points.size()){
-    Msg::Error("Tetrahedralization added points (%d -> %d): aborting",
-               points.size(), out.numberofpoints);
-    return v1;
+    MVertex *v = points[i].v;
+    v->setIndex(i);
+    vertices[i] = v;
   }
+  std::vector<MTetrahedron*> tets;
+  delaunayMeshIn3D(vertices, tets);
 
   // create output
   PView *v2 = new PView();
   PViewDataList *data2 = getDataList(v2);
-  for(int i = 0; i < out.numberoftetrahedra; i++){
-    PointData *p[4];
-    p[0] = &points[out.tetrahedronlist[i * 4 + 0]];
-    p[1] = &points[out.tetrahedronlist[i * 4 + 1]];
-    p[2] = &points[out.tetrahedronlist[i * 4 + 2]];
-    p[3] = &points[out.tetrahedronlist[i * 4 + 3]];
+  for(unsigned int i = 0; i < tets.size(); i++){
+    MTetrahedron *t = tets[i];
+    int i0 = t->getVertex(0)->getIndex();
+    int i1 = t->getVertex(1)->getIndex();
+    int i2 = t->getVertex(2)->getIndex();
+    int i3 = t->getVertex(3)->getIndex();
+    int n = points.size();
+    if(i0 < 0 || i0 >= n || i1 < 0 || i1 >= n ||
+       i2 < 0 || i2 >= n || i3 < 0 || i3 >= n){
+      Msg::Warning("Bad vertex in tetrahedralization");
+      continue;
+    }
+    PointData *p[4] = {&points[i0], &points[i1],
+                       &points[i2], &points[i3]};
     int numComp = 0;
     std::vector<double> *vec = 0;
-    if((int)p[0]->v.size() == 3 + 9 * numSteps &&
-       (int)p[1]->v.size() == 3 + 9 * numSteps &&
-       (int)p[2]->v.size() == 3 + 9 * numSteps &&
-       (int)p[3]->v.size() == 3 + 9 * numSteps){
+    if((int)p[0]->val.size() == 9 * numSteps &&
+       (int)p[1]->val.size() == 9 * numSteps &&
+       (int)p[2]->val.size() == 9 * numSteps &&
+       (int)p[3]->val.size() == 9 * numSteps){
       numComp = 9; data2->NbTS++; vec = &data2->TS;
     }
-    else if((int)p[0]->v.size() == 3 + 3 * numSteps &&
-            (int)p[1]->v.size() == 3 + 3 * numSteps &&
-            (int)p[2]->v.size() == 3 + 3 * numSteps &&
-            (int)p[2]->v.size() == 3 + 3 * numSteps){
+    else if((int)p[0]->val.size() == 3 * numSteps &&
+            (int)p[1]->val.size() == 3 * numSteps &&
+            (int)p[2]->val.size() == 3 * numSteps &&
+            (int)p[3]->val.size() == 3 * numSteps){
       numComp = 3; data2->NbVS++; vec = &data2->VS;
     }
-    else{
+    else if((int)p[0]->val.size() == numSteps &&
+            (int)p[1]->val.size() == numSteps &&
+            (int)p[2]->val.size() == numSteps &&
+            (int)p[3]->val.size() == numSteps){
       numComp = 1; data2->NbSS++; vec = &data2->SS;
     }
-    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v[0]);
-    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v[1]);
-    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v[2]);
+    else{
+      Msg::Warning("Bad data in tetrahedralization");
+      continue;
+    }
+    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v->x());
+    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v->y());
+    for(int nod = 0; nod < 4; nod++) vec->push_back(p[nod]->v->z());
     for(int step = 0; step < numSteps; step++)
       for(int nod = 0; nod < 4; nod++)
         for(int comp = 0; comp < numComp; comp++)
-          vec->push_back(p[nod]->v[3 + numComp * step + comp]);
+          vec->push_back(p[nod]->val[numComp * step + comp]);
   }
 
+  for(unsigned int i = 0; i < tets.size(); i++)
+    delete tets[i];
+  for(unsigned int i = 0; i < vertices.size(); i++)
+    delete vertices[i];
+
   for(int i = 0; i < data1->getNumTimeSteps(); i++)
     data2->Time.push_back(data1->getTime(i));
   data2->setName(data1->getName() + "_Tetrahedralize");
   data2->setFileName(data1->getName() + "_Tetrahedralize.pos");
   data2->finalize();
-
   return v2;
 #endif
 }