From eb0f613adfdf8b55a317591738f9f08b384c7595 Mon Sep 17 00:00:00 2001
From: Nicolas Tardieu <nicolas.tardieu@edf.fr>
Date: Wed, 29 Jun 2005 09:57:32 +0000
Subject: [PATCH] New fonction to use Tetgen 's algorithm Attention : several
 problems   - the method deinitialized of Tetgen is not seen by the compiler  
 - it seems that the data structure tetgenio is destroyed after its creation,
 so that all meshing procedure is dine in the constructor

---
 Mesh/3D_Mesh_Tetgen.cpp | 276 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 276 insertions(+)
 create mode 100644 Mesh/3D_Mesh_Tetgen.cpp

diff --git a/Mesh/3D_Mesh_Tetgen.cpp b/Mesh/3D_Mesh_Tetgen.cpp
new file mode 100644
index 0000000000..b0326a2a10
--- /dev/null
+++ b/Mesh/3D_Mesh_Tetgen.cpp
@@ -0,0 +1,276 @@
+
+#include "Gmsh.h"
+#include "Geo.h"
+#include "Mesh.h"
+#include "Create.h"
+#include "Numeric.h"
+#include "Context.h"
+
+extern Context_T CTX;
+extern Mesh *THEM;
+
+#if !defined(HAVE_TETGEN)
+
+int Tetgen(Volume * v)
+{
+  if(CTX.mesh.algo3d == DELAUNAY_TETGEN)
+    Msg(GERROR, "Tetgen is not compiled in this version of Gmsh");
+  return 0;
+}
+
+#else
+
+#include "tetgen.h"
+
+
+class Tetgen{
+ private:
+  List_T *_surverts, *_volverts;
+  Volume *_vol;
+  tetgenio *in, *out;
+ public:
+  Tetgen(Volume *vol);
+  ~Tetgen();
+  void MeshVolume();
+  void TransferVolumeMesh();
+};
+
+Tetgen::Tetgen(Volume *vol)
+  : _volverts(0), _vol(vol)
+{
+  // ===================================
+  //     TRANSFER GMSH => TETGEN
+  // ===================================
+  
+  // creates Tetgen mesh structure
+  tetgenio *in = new tetgenio;
+  tetgenio *out = new tetgenio;
+  tetgenio::facet *f;
+  tetgenio::polygon *p;
+
+  // All indices start from 1.
+  in->firstnumber = 1;
+  
+  // Get all surface vertices (the same vertex can belong to several
+  // surfaces...)
+  Tree_T *tree = Tree_Create(sizeof(Vertex*), compareVertex);
+  for(int i = 0; i < List_Nbr(_vol->Surfaces); i++) {
+    Surface *s;
+    List_Read(_vol->Surfaces, i, &s);
+    Tree_Unit(tree, s->Vertices);
+  }
+  _surverts = Tree2List(tree);
+  List_Sort(_surverts, compareVertex);
+
+  Tree_Delete(tree);
+  
+  // Tetgen points input data
+  in->numberofpoints = List_Nbr(_surverts);
+  in->pointlist = new REAL[in->numberofpoints * 3];
+
+  // Transfer the vertices
+  for(int i = 0; i < List_Nbr(_surverts); i++){
+    Vertex *v;
+    List_Read(_surverts, i, &v);
+    in->pointlist[3*i + 0]  = v->Pos.X;
+    in->pointlist[3*i + 1]  = v->Pos.Y;
+    in->pointlist[3*i + 2]  = v->Pos.Z;
+  }
+
+  // Count number of facets
+  int NbOfFacets=0;
+  for(int i = 0; i < List_Nbr(_vol->Surfaces); i++) {
+    Surface *s;
+    List_Read(_vol->Surfaces, i, &s);
+    List_T *simplist = Tree2List(s->Simplexes);
+    NbOfFacets = NbOfFacets + List_Nbr(simplist);
+  }
+  // Tetgen elements input data
+  in->numberoffacets = NbOfFacets;
+  in->facetlist = new tetgenio::facet[in->numberoffacets];
+  in->facetmarkerlist = new int[in->numberoffacets];
+
+  // Transfert all surface simplices
+  int currentFacet=0;
+  for(int i = 0; i < List_Nbr(_vol->Surfaces); i++) {
+    Surface *s;
+    List_Read(_vol->Surfaces, i, &s);
+    int sign;
+    List_Read(_vol->SurfacesOrientations, i, &sign);
+    List_T *simplist = Tree2List(s->Simplexes);
+    for(int j = 0; j < List_Nbr(simplist); j++) {
+      Simplex *simp;
+      List_Read(simplist, j, &simp);
+      int tmp[3], index[3];
+      if(sign > 0){
+	index[0] = 0;
+	index[1] = 1;
+	index[2] = 2;
+      }
+      else{
+	index[0] = 0;
+	index[1] = 2;
+	index[2] = 1;
+      }
+      tmp[0] = 1 + List_ISearch(_surverts, &simp->V[index[0]], compareVertex);
+      tmp[1] = 1 + List_ISearch(_surverts, &simp->V[index[1]], compareVertex);
+      tmp[2] = 1 + List_ISearch(_surverts, &simp->V[index[2]], compareVertex);
+      
+      f = &in->facetlist[currentFacet];
+      f->numberofpolygons = 1;
+      f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+      f->numberofholes = 0;
+      f->holelist = NULL;
+      p = &f->polygonlist[0];
+      p->numberofvertices = 3;
+      p->vertexlist = new int[p->numberofvertices];
+      p->vertexlist[0] = tmp[0];
+      p->vertexlist[1] = tmp[1];
+      p->vertexlist[2] = tmp[2];
+      in->facetmarkerlist[currentFacet] = 0;
+      currentFacet=currentFacet+1;
+    }
+    List_Delete(simplist);
+  }
+  // In order to check the input...
+  in->save_nodes("barin");
+  in->save_poly("barin");
+
+
+  
+  
+  // ===================================
+  //     MESHING...
+  // ===================================
+  tetrahedralize("pqV", in, out);
+  
+  
+  
+  
+  // ===================================
+  //     TRANSFER TETGEN => GMSH
+  // ===================================
+  // Gets total number of vertices of Tetgen's mesh
+  int nbv = out->numberofpoints;
+  Msg(INFO, " NIKO : nbv=%i",nbv);
+
+  for(int i = 0; i < List_Nbr(_vol->Surfaces); i++) {
+    Surface *s;
+    List_Read(_vol->Surfaces, i, &s);
+    Tree_Delete(s->Simplexes);
+  }
+  Tree_Delete(_vol->Vertices);
+  Tree_Delete(THEM->Vertices);
+  _vol->Vertices = Tree_Create(sizeof(Simplex *), compareSimplex);
+  THEM->Vertices = Tree_Create(sizeof(Simplex *), compareSimplex);
+  THEM->MaxPointNum=0;
+
+  // Create new vertices
+  Vertex **vtable = (Vertex **)Malloc(nbv * sizeof(Vertex*));
+  for(int i = 0; i < nbv; i++) {
+    vtable[i] = Create_Vertex(++(THEM->MaxPointNum), out->pointlist[3*i + 0], out->pointlist[3*i + 1],out->pointlist[3*i + 2], 1., 0);
+    Tree_Add(_vol->Vertices, &vtable[i]);
+    Tree_Add(THEM->Vertices, &vtable[i]);
+  }
+
+  Msg(INFO, "<NIKO>: nbTriFace=%i",out->numberoftrifaces);
+  if(!nbv) {
+    Msg(STATUS3, "<NIKO>: AUCUN ELEMENT");
+    return;}
+
+  
+  // Get total number of simplices of Tetgen's mesh
+  int nbe = out->numberoftetrahedra;
+
+  // Create new volume simplices
+  for(int i = 0; i < nbe; i++) {
+    Simplex *simp = Create_Simplex(vtable[out->tetrahedronlist[4*i + 0]-1], vtable[out->tetrahedronlist[4*i + 1]-1],
+				   vtable[out->tetrahedronlist[4*i + 2]-1], vtable[out->tetrahedronlist[4*i + 3]-1]);
+    simp->iEnt = _vol->Num;
+    Tree_Add(_vol->Simplexes, &simp);
+    // also add a copy in the global simplex tree
+    Tree_Add(THEM->Simplexes, &simp);
+  }
+  
+  Free(vtable);
+}
+
+
+Tetgen::~Tetgen()
+{
+  List_Delete(_surverts);
+  List_Delete(_volverts);
+  // Strange : this method is not seen by the compiler 
+  // tetgenio::deinitialize();
+}
+
+void Tetgen::MeshVolume()
+{
+  tetrahedralize("pqV", in, out);
+}
+
+void Tetgen::TransferVolumeMesh()
+{
+  // Gets total number of vertices of Tetgen's mesh
+  int nbv = out->numberofpoints;
+  Msg(INFO, " NIKO : nbv=%i",nbv);
+
+  if(!nbv) {
+    Msg(STATUS3, " NIKO : AUCUN ELEMENT");
+    return;}
+
+  Vertex **vtable = (Vertex **)Malloc(nbv * sizeof(Vertex*));
+  
+  // Get existing unmodified surface vertices
+  for(int i = 0; i < List_Nbr(_surverts); i++){
+    List_Read(_surverts, i, &vtable[i]);
+    Tree_Insert(_vol->Vertices, &vtable[i]);
+    Tree_Insert(THEM->Vertices, &vtable[i]);
+  }
+
+  // Create new volume vertices
+  for(int i = List_Nbr(_surverts); i < nbv; i++) {
+    vtable[i] = Create_Vertex(++(THEM->MaxPointNum), out->pointlist[3*i + 0], out->pointlist[3*i + 1],out->pointlist[3*i + 2], 1., 0);
+    Tree_Add(_vol->Vertices, &vtable[i]);
+    Tree_Add(THEM->Vertices, &vtable[i]);
+  }
+
+  // Get total number of simplices of Tetgen's mesh
+  int nbe = out->numberoftetrahedra;
+
+  // Create new volume simplices
+  for(int i = 0; i < nbe; i++) {
+    Msg(INFO, " NIKO out->tetrahedronlist: =%i",out->tetrahedronlist[4*i + 0]);
+    Simplex *simp = Create_Simplex(vtable[out->tetrahedronlist[4*i + 0]], vtable[out->tetrahedronlist[4*i + 1]],
+				   vtable[out->tetrahedronlist[4*i + 2]], vtable[out->tetrahedronlist[4*i + 3]]);
+    simp->iEnt = _vol->Num;
+    Tree_Add(_vol->Simplexes, &simp);
+    // also add a copy in the global simplex tree
+    Tree_Add(THEM->Simplexes, &simp);
+  }
+  
+  Free(vtable);
+}
+
+
+int Mesh_Tetgen(Volume * v)
+{
+  if(CTX.mesh.algo3d != DELAUNAY_TETGEN)
+    return 0;
+
+  if(THEM->BGM.Typ == ONFILE){
+    Msg(GERROR, "Tetgen is not ready to be used with a background mesh");
+    return 0;
+  }
+
+  Msg(STATUS3, "Meshing volume %d", v->Num);
+  Tetgen tg(v);
+  //tg.MeshVolume();
+  //tg.TransferVolumeMesh();
+  
+  return 1;
+}
+
+
+
+#endif // !HAVE_NETGEN
-- 
GitLab