diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 929ce4354804c7942626c8fa3f33220374e0defd..b811fe1d6b2d7f00be9b1955080d2a1cbd691b5c 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -948,7 +948,7 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "Algorithm" , opt_mesh_algo2d , ALGO_2D_AUTO ,
     "2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)" },
   { F|O, "Algorithm3D" , opt_mesh_algo3d , ALGO_3D_DELAUNAY ,
-    "3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)" },
+    "3D mesh algorithm (1=Delaunay, 2=New Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)" },
   { F|O, "AngleSmoothNormals" , opt_mesh_angle_smooth_normals , 30.0 ,
     "Threshold angle below which normals are not smoothed" },
   { F|O, "AnisoMax" , opt_mesh_aniso_max, 1.e33,
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index 85592a19fbc0c5a1280db8b56978f6b05ce6f593..2abd8dd9a4d37e31ec16944d81bbea4d29c6956f 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -241,6 +241,7 @@
 
 // 3D meshing algorithms (numbers should not be changed)
 #define ALGO_3D_DELAUNAY       1
+#define ALGO_3D_DELAUNAY_NEW   2
 #define ALGO_3D_FRONTAL        4
 #define ALGO_3D_FRONTAL_DEL    5
 #define ALGO_3D_FRONTAL_HEX    6
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 903e7bdc6e4b2b7de42626a284427e2e1ea1c4f9..c8e0758cc7d2a8c8d7ecd0e62763a5991b6ad619 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -38,7 +38,6 @@
 #include "discreteFace.h"
 #include "filterElements.h"
 
-
 #if defined(HAVE_ANN)
 #include "ANN/ANN.h"
 #endif
@@ -421,31 +420,6 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
   }
 }
 
-#endif
-
-static void addOrRemove(const MFace &f,
-			MElement *e,
-			std::map<MFace,MElement*,Less_Face> & bfaces,
-			splitQuadRecovery &sqr)
-{
-  {
-    std::map<MFace, MVertex*, Less_Face>::const_iterator it = sqr._invmap.find(f);
-    if (it != sqr._invmap.end()){
-      addOrRemove (MFace(it->second, f.getVertex(0),f.getVertex(1)),e,bfaces,sqr);
-      addOrRemove (MFace(it->second, f.getVertex(1),f.getVertex(2)),e,bfaces,sqr);
-      addOrRemove (MFace(it->second, f.getVertex(2),f.getVertex(3)),e,bfaces,sqr);
-      addOrRemove (MFace(it->second, f.getVertex(3),f.getVertex(0)),e,bfaces,sqr);
-      return;
-    }
-  }
-
-  std::map<MFace,MElement*,Less_Face>::iterator it = bfaces.find(f);
-  if (it == bfaces.end())bfaces.insert(std::make_pair(f,e));
-  else bfaces.erase(it);
-}
-
-
-#if defined(HAVE_TETGEN)
 bool CreateAnEmptyVolumeMesh(GRegion *gr)
 {
   printf("creating an empty volume mesh\n");
@@ -470,24 +444,10 @@ bool CreateAnEmptyVolumeMesh(GRegion *gr)
   return true;
 }
 
-#else
-
-bool CreateAnEmptyVolumeMesh(GRegion *gr)
-{
-  Msg::Error("You should compile with TETGEN in order to create an empty volume mesh");
-  return false;
-}
-
-#endif
-
 void MeshDelaunayVolumeTetgen(std::vector<GRegion*> &regions)
 {
   if(regions.empty()) return;
 
-#if !defined(HAVE_TETGEN)
-  Msg::Error("Tetgen is not compiled in this version of Gmsh");
-#else
-
   for(unsigned int i = 0; i < regions.size(); i++)
     Msg::Info("Meshing volume %d (Delaunay)", regions[i]->tag());
 
@@ -588,12 +548,16 @@ void MeshDelaunayVolumeTetgen(std::vector<GRegion*> &regions)
  if (sqr.buildPyramids (gr->model())){
    RelocateVertices (regions);
  }
-#endif
 }
 
-// uncomment this to test the new code
-#if !defined(HAVE_TETGEN)
-#define NEW_CODE
+#else
+
+bool CreateAnEmptyVolumeMesh(GRegion *gr)
+{
+  Msg::Error("You should compile with TETGEN in order to create an empty volume mesh");
+  return false;
+}
+
 #endif
 
 static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> &regions)
@@ -647,11 +611,16 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
 {
   if(regions.empty()) return;
 
-#if !defined(NEW_CODE) && defined(HAVE_TETGEN)
-  MeshDelaunayVolumeTetgen(regions);
+  if(CTX::instance()->mesh.algo3d == ALGO_3D_DELAUNAY_NEW){
+    MeshDelaunayVolumeNewCode(regions);
+  }
+  else{
+#if defined(HAVE_TETGEN)
+    MeshDelaunayVolumeTetgen(regions);
 #else
-  MeshDelaunayVolumeNewCode(regions);
+    Msg::Error("Old Delaunay algorithm requires Tetgen");
 #endif
+  }
   return;
 }
 
@@ -987,7 +956,7 @@ void meshGRegion::operator() (GRegion *gr)
   }
   else if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL){
 #if !defined(HAVE_NETGEN)
-    Msg::Error("Netgen is not compiled in this version of Gmsh");
+    Msg::Error("Frontal algorithm requires Netgen");
 #else
     Msg::Info("Meshing volume %d (Frontal)", gr->tag());
     // orient the triangles of with respect to this region
@@ -1015,7 +984,7 @@ void optimizeMeshGRegionNetgen::operator() (GRegion *gr)
   if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY) return;
 
 #if !defined(HAVE_NETGEN)
-  Msg::Error("Netgen is not compiled in this version of Gmsh");
+  Msg::Error("Netgen optimizer is not compiled in this version of Gmsh");
 #else
   Msg::Info("Optimizing volume %d", gr->tag());
   // import mesh into netgen, including volume tets
diff --git a/Mesh/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp
index 724358da0d858ddc64056cdeb60be0190e290554..7eab024953688066e348838cca51ec1c0580c266 100644
--- a/Mesh/meshGRegionBoundaryRecovery.cpp
+++ b/Mesh/meshGRegionBoundaryRecovery.cpp
@@ -17,15 +17,8 @@
 #include "MTetrahedron.h"
 #include "OS.h"
 
-#if defined(HAVE_TETGEN)
-
-bool meshGRegionBoundaryRecovery(GRegion *gr)
+namespace tetgenBR
 {
-  Msg::Error("Should not call meshGRegionBoundaryRecovery with Tetgen");
-  return false;
-}
-
-#else
 
 #define REAL double
 
@@ -94,25 +87,13 @@ public:
 };
 
 // redefinition of predicates using our own
-#define exactinit robustPredicates::exactinit
-#define incircle robustPredicates::incircle
-#define insphere robustPredicates::insphere
-#define orient2d robustPredicates::orient2d
 #define orient3d robustPredicates::orient3d
-double orient4d(double*, double *, double *, double *, double *,
-                double, double, double, double, double){ return 0.;}
+static double orient4d(double*, double *, double *, double *, double *,
+                double, double, double, double, double){ return 0.; }
 
 #include "tetgenBR.h"
 #include "tetgenBR.cxx"
 
-bool meshGRegionBoundaryRecovery(GRegion *gr)
-{
-  tetgenmesh *m = new tetgenmesh();
-  bool ret = m->reconstructmesh((void*)gr);
-  delete m;
-  return ret;
-}
-
 bool tetgenmesh::reconstructmesh(void *p)
 {
   GRegion *_gr = (GRegion*)p;
@@ -1050,4 +1031,12 @@ void tetgenmesh::outmesh2medit(const char* mfilename)
   fclose(outfile);
 }
 
-#endif
+} // end namespace
+
+bool meshGRegionBoundaryRecovery(GRegion *gr)
+{
+  tetgenBR::tetgenmesh *m = new tetgenBR::tetgenmesh();
+  bool ret = m->reconstructmesh((void*)gr);
+  delete m;
+  return ret;
+}