diff --git a/Common/Context.h b/Common/Context.h
index a135bf6d42a610be95cde585dec6504f09f1d49b..5fc86f059e5d9aa8869405fdb646d1c6ab373eec 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 b811fe1d6b2d7f00be9b1955080d2a1cbd691b5c..409eb7dff254b2dc99a5727b25669569a6a9acba 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 52804dbdc1523fd38620ac73c7abd5b2300f786f..1681ea9a63f1ece7fcd480554b2a11e96aee6292 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 a2952a8479d87e52058c309b1ca9a995db9f270d..a905474248d8b60bd8f0621b0e9a940cecf97bfe 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 55c70c029c72cf7e114dd5e1023f7ee93048a4f0..352467f3b776ba290211b805b7b37ab09f4ae7b4 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 7c53d7fb9ccbe68163a562377ee18041a305dd80..4ee499d7f81714d6f2b14564bf0fac06886adddc 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 01802c7527cd7ad5fea2cdcd946efc49628e25a7..9e249a695acd4608d81892577fc0a2297389cabb 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;