From 9470acd6e85389f4d357b24bb17a7c5854b1fcb4 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Fri, 18 Mar 2016 11:01:18 +0000
Subject: [PATCH] option for old point insertion algo

---
 Common/Context.h                     |  2 +-
 Common/DefaultOptions.h              |  2 ++
 Common/Options.cpp                   |  7 +++++++
 Common/Options.h                     |  1 +
 Geo/MElement.cpp                     |  4 ++--
 Mesh/meshGRegion.cpp                 | 29 ++++++++++++++--------------
 Mesh/meshGRegionBoundaryRecovery.cpp |  3 ++-
 7 files changed, 30 insertions(+), 18 deletions(-)

diff --git a/Common/Context.h b/Common/Context.h
index a135bf6d42..5fc86f059e 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -32,7 +32,7 @@ struct contextMeshOptions {
   double anisoMax, smoothRatio;
   int lcFromPoints, lcFromCurvature, lcExtendFromBoundary;
   int dual, voronoi, drawSkinOnly, colorCarousel, labelSampling;
-  int fileFormat, nbSmoothing, algo2d, algo3d, algoSubdivide;
+  int fileFormat, nbSmoothing, algo2d, algo3d, algoSubdivide, oldRefinement;
   int algoRecombine, recombineAll, recombine3DAll, recombine3DLevel, recombine3DConformity;
   int flexibleTransfinite;
   //-- for recombination test (amaury) --
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index b811fe1d6b..409eb7dff2 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1141,6 +1141,8 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "NumSubEdges" , opt_mesh_num_sub_edges , 2. ,
     "Number of edge subdivisions when displaying high order elements" },
 
+  { F|O, "OldRefinement" , opt_mesh_old_refinement , 1 ,
+    "Use old 3D point insertion algorithm" },
   { F|O, "Optimize" , opt_mesh_optimize , 0. ,
     "Optimize the mesh to improve the quality of tetrahedral elements" },
   { F|O, "OptimizeNetgen" , opt_mesh_optimize_netgen , 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 52804dbdc1..1681ea9a63 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -4894,6 +4894,13 @@ double opt_mesh_optimize_netgen(OPT_ARGS_NUM)
   return CTX::instance()->mesh.optimizeNetgen;
 }
 
+double opt_mesh_old_refinement(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX::instance()->mesh.oldRefinement = (int)val;
+  return CTX::instance()->mesh.oldRefinement;
+}
+
 double opt_mesh_refine_steps(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET){
diff --git a/Common/Options.h b/Common/Options.h
index a2952a8479..a905474248 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -403,6 +403,7 @@ double opt_geometry_match_geom_and_mesh(OPT_ARGS_NUM);
 double opt_mesh_label_sampling(OPT_ARGS_NUM);
 double opt_mesh_optimize(OPT_ARGS_NUM);
 double opt_mesh_optimize_netgen(OPT_ARGS_NUM);
+double opt_mesh_old_refinement(OPT_ARGS_NUM);
 double opt_mesh_refine_steps(OPT_ARGS_NUM);
 double opt_mesh_normals(OPT_ARGS_NUM);
 double opt_mesh_num_sub_edges(OPT_ARGS_NUM);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 55c70c029c..352467f3b7 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1164,8 +1164,8 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
     double gamma = gammaShapeMeasure();
     for(int i = 0; i < n; i++){
       if(first) first = false; else fprintf(fp, ",");
-      //      fprintf(fp, "%g", gamma);  FIXME
-      fprintf(fp, "%d", getVertex(i)->getNum());
+      fprintf(fp, "%g", gamma);
+      //fprintf(fp, "%d", getVertex(i)->getNum());
     }
   }
   if(printRho){
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 7c53d7fb9c..4ee499d7f8 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -589,7 +589,7 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
      fprintf(fp, "};\n");
      fclose(fp);
    }
-   
+
 
   }
   catch(int err){
@@ -613,20 +613,21 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
   gr->set(faces);
 
   // now do insertion of points
-#if 1
-  insertVerticesInRegion(gr, 2000000000, true);
-#else
-  void edgeBasedRefinement(const int numThreads,
-                           const int nptsatonce,
-                           GRegion *gr);
-  // just to remove tets that are not to be meshed
-  insertVerticesInRegion(gr, 0);
-  for(unsigned int i = 0; i < regions.size(); i++){
-    Msg::Info("Refining volume %d",regions[i]->tag());
-    edgeBasedRefinement(1, 1, regions[i]);
+  if(CTX::instance()->mesh.oldRefinement){
+    insertVerticesInRegion(gr, 2000000000, true);
+  }
+  else{
+    void edgeBasedRefinement(const int numThreads,
+                             const int nptsatonce,
+                             GRegion *gr);
+    // just to remove tets that are not to be meshed
+    insertVerticesInRegion(gr, 0);
+    for(unsigned int i = 0; i < regions.size(); i++){
+      Msg::Info("Refining volume %d",regions[i]->tag());
+      edgeBasedRefinement(1, 1, regions[i]);
+    }
+    // RelocateVertices (regions,-1);
   }
-  // RelocateVertices (regions,-1);
-#endif
 }
 
 void MeshDelaunayVolume(std::vector<GRegion*> &regions)
diff --git a/Mesh/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp
index 01802c7527..9e249a695a 100644
--- a/Mesh/meshGRegionBoundaryRecovery.cpp
+++ b/Mesh/meshGRegionBoundaryRecovery.cpp
@@ -512,7 +512,8 @@ bool tetgenmesh::reconstructmesh(void *p)
 
   if ((dupverts > 0l) || (unuverts > 0l)) {
     // Remove hanging nodes.
-    //    jettisonnodes();
+    // cannot call this here due to 8 additional exterior vertices we inserted
+    // jettisonnodes();
   }
 
   long tetnumber, facenumber;
-- 
GitLab