diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 617ce1ffc6071d4262857896de7cf58861430c50..163daeb22fb88b3ad2fbde632050fdd9f7fbd5e2 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -423,17 +423,19 @@ void meshGEdge::operator() (GEdge *ge)
   // force odd number of points if blossom is used for recombination
   if((ge->meshAttributes.method != MESH_TRANSFINITE ||
       CTX::instance()->mesh.flexibleTransfinite) &&
-     CTX::instance()->mesh.algoRecombine == 1){
+     CTX::instance()->mesh.algoRecombine != 0){
     if(CTX::instance()->mesh.recombineAll){
       if (N % 2 == 0) N++;
-      N = increaseN(N);
+      if (CTX::instance()->mesh.algoRecombine == 2)
+	N = increaseN(N);
     }
     else{
       std::list<GFace*> faces = ge->faces();
       for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
         if((*it)->meshAttributes.recombine){
 	  if (N % 2 == 0) N ++;
-	  N = increaseN(N);
+	  if (CTX::instance()->mesh.algoRecombine == 2)
+	    N = increaseN(N);
           break;
         }
       }
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index daa31656aa64553c995febc28fa97852c7eb22cf..6d3a9c129964d31712e98c10ada5b61435c89cbc 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -47,10 +47,7 @@
 #include "filterElements.h"
 
 // define this to use the old initial delaunay
-#define OLD_TRI_CODE
-
-// define this to use th old quad code
-#define OLD_QUAD_CODE
+//#define OLD_CODE_DELAUNAY 1
 
 static void computeElementShapes(GFace *gf, double &worst, double &avg,
                                  double &best, int &nT, int &greaterThan)
@@ -88,7 +85,7 @@ public:
   // remove one point every two and remember middle points
   quadMeshRemoveHalfOfOneDMesh (GFace* gf) : _gf(gf){
     // only do it if a recombination has to be done
-    if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine){
+    if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && CTX::instance()->mesh.algoRecombine == 2){
       //      printf("GFace %d removing half of the points in the 1D mesh\n",gf->tag());
       std::list<GEdge*> edges = gf->edges();
       std::list<GEdge*>::iterator ite = edges.begin();
@@ -97,6 +94,7 @@ public:
 	  std::vector<MLine*> temp;
 	  (*ite)->mesh_vertices.clear();
 	  for(unsigned int i = 0; i< (*ite)->lines.size(); i+=2){
+	    if (i+1 >= (*ite)->lines.size())Msg::Fatal("1D mesh cannot be divided by 2");
 	    MVertex *v1 = (*ite)->lines[i]->getVertex(0);
 	    MVertex *v2 = (*ite)->lines[i]->getVertex(1);
 	    MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
@@ -189,7 +187,7 @@ public:
   }
   void finish (){
     backgroundMesh::setSizeFactor(1.0);
-    if(CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine){
+    if((CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine) && CTX::instance()->mesh.algoRecombine == 2){
       // recombine the elements on the half mesh
       recombineIntoQuads(_gf,true,true,.1);
       Msg::Info("subdividing");
@@ -1199,7 +1197,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   // use a divide & conquer type algorithm to create a triangulation.
   // We add to the triangulation a box with 4 points that encloses the
   // domain.
-#if defined(OLD_TRI_CODE)
+#ifdef OLD_CODE_DELAUNAY
   {
     // compute the bounding box in parametric space
     SVector3 dd(bbox.max(), bbox.min());
@@ -1675,11 +1673,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   // delete the mesh
   delete m;
 
-  if (1){
-    if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
-       !CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh)
-      recombineIntoQuads(gf);
-  }
+  if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
+     !CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh && CTX::instance()->mesh.algoRecombine != 2)
+    recombineIntoQuads(gf);
 
   computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
                        gf->meshStatistics.average_element_shape,
@@ -2051,7 +2047,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
   // Use a divide & conquer type algorithm to create a triangulation.
   // We add to the triangulation a box with 4 points that encloses the
   // domain.
-#if 1 //OLD_TRI_CODE
+#if 1 //OLD_CODE_DELAUNAY
   {
     DocRecord doc(nbPointsTotal + 4);
     int count = 0;
@@ -2464,11 +2460,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
   // delete the mesh
   delete m;
 
-  if (1){
-    if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
-       !CTX::instance()->mesh.optimizeLloyd)
-      recombineIntoQuads(gf);
-  }
+  if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
+     !CTX::instance()->mesh.optimizeLloyd && CTX::instance()->mesh.algoRecombine != 2)
+    recombineIntoQuads(gf);
 
   computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
                        gf->meshStatistics.average_element_shape,
@@ -2562,9 +2556,7 @@ void meshGFace::operator() (GFace *gf, bool print)
     return;
   }
 
-#if !defined(OLD_QUAD_CODE)
   quadMeshRemoveHalfOfOneDMesh halfmesh (gf);
-#endif
 
   if ((gf->getNativeType() != GEntity::AcisModel ||
        (!gf->periodic(0) && !gf->periodic(1))) &&
@@ -2582,9 +2574,7 @@ void meshGFace::operator() (GFace *gf, bool print)
   Msg::Debug("Type %d %d triangles generated, %d internal vertices",
              gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
 
-#if !defined(OLD_QUAD_CODE)
   halfmesh.finish();
-#endif
 }
 
 bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 560ca4c299d81cc32f0bd03fd8c05f91586815af..99a1d29a03bad96c2989fdf43e88dd0c28591c64 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1816,7 +1816,6 @@ static void initialSquare(std::vector<MVertex*> &v,
   box[1] = new MVertex (bbox.max().x(),bbox.min().y(),0);
   box[2] = new MVertex (bbox.max().x(),bbox.max().y(),0);
   box[3] = new MVertex (bbox.min().x(),bbox.max().y(),0);
-  std::vector<MTriangle*> t_box;
   MTriangle *t0 = new MTriangle (box[0],box[1],box[2]);
   MTriangle *t1 = new MTriangle (box[2],box[3],box[0]);
   t.push_back(new MTri3(t0,0.0));
@@ -1844,9 +1843,12 @@ MTri3 * getTriToBreak (MVertex *v, std::vector<MTri3*> &t, int &NB_GLOBAL_SEARCH
 }
 
 bool triOnBox (MTriangle *t, MVertex *box[4]){
-  for (size_t i = 0;i<3;i++)
-    for (size_t j = 0;j<4;j++)
-      if (t->getVertex(i) == box[j])return true;
+  for (size_t i = 0;i<3;i++){
+    for (size_t j = 0;j<4;j++){
+      if (t->getVertex(i) == box[j])
+	return true;
+    }
+  }
   return false;
 }
 
@@ -1869,13 +1871,16 @@ void delaunayMeshIn2D(std::vector<MVertex*> &v,
   initialSquare (v,box,t);
 
   int NB_GLOBAL_SEARCH = 0;
-  int AVG_ITER = 0;
+  double AVG_ITER = 0;
+  double AVG_CAVSIZE = 0;
 
   double t1 = Cpu();
 
+  Msg::Info("Delaunay 2D SORTING");
   if(hilbertSort) SortHilbert(v);
 
   double ta=0,tb=0,tc=0,td=0,T;
+  Msg::Info("Delaunay 2D INSERTING");
   for (size_t i=0;i<v.size();i++){
     MVertex *pv = v[i];
 
@@ -1883,7 +1888,7 @@ void delaunayMeshIn2D(std::vector<MVertex*> &v,
     T = Cpu();
     MTri3 * found = getTriToBreak (pv,t,NB_GLOBAL_SEARCH,NITER);
     ta += Cpu()-T;
-    AVG_ITER += NITER;
+    AVG_ITER += (double)NITER;
     if(!found) {
       Msg::Error("Cannot insert a point in 2D Delaunay");
       continue;
@@ -1893,6 +1898,7 @@ void delaunayMeshIn2D(std::vector<MVertex*> &v,
 
     T = Cpu();
     recurFindCavity(shell, cavity, pv, found);
+    AVG_CAVSIZE += (double)cavity.size();
     tb += Cpu()-T;
     //double V = 0.0;
     //for (unsigned int k=0;k<cavity.size();k++)V+=fabs(cavity[k]->tri()->getVolume());
@@ -1940,7 +1946,8 @@ void delaunayMeshIn2D(std::vector<MVertex*> &v,
   }
 
   double t2 = Cpu();
-  Msg::Debug("Delaunay 2D done for %d points : CPU = %g, %d global searches, AVG walk size %g",v.size(), t2-t1,NB_GLOBAL_SEARCH,1.+(double)AVG_ITER/v.size());
+  Msg::Info("Delaunay 2D done for %d points : CPU = %g, %d global searches, AVG walk size %g , AVG cavity size %g",
+	    v.size(), t2-t1,NB_GLOBAL_SEARCH,1.+AVG_ITER/v.size(),AVG_CAVSIZE/v.size());
   //  printf("%g %g %g %g --> %g(%g)\n",ta,tb,tc,td,t2-t1,ta+tb+tc+td);
 
   if (edgesToRecover)recoverEdges (t,*edgesToRecover);
@@ -1955,9 +1962,9 @@ void delaunayMeshIn2D(std::vector<MVertex*> &v,
     }
     delete t[i];
   }
-
-  if (removeBox)for (int i=0;i<4;i++)delete box[i];
-  else for (int i=0;i<4;i++)v.push_back(box[i]);
+  
+  if (removeBox){for (int i=0;i<4;i++)delete box[i];}
+  else {for (int i=0;i<4;i++)v.push_back(box[i]);}
 
   //  fprintf(f,"};\n");
   //  fclose(f);
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a91f8b73d4f59ea996a6b5ba53610f1b5c614d98..b3bdb9290c3d7fd35ce069b303b175fc8f72da3a 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -3165,7 +3165,7 @@ int recombineWithBlossom(GFace *gf, double dx, double dy,
   std::set<MElement*> touched;
   std::vector<std::pair<MElement*,MElement*> > toProcess;
 
-  if(CTX::instance()->mesh.algoRecombine == 1){
+  if(CTX::instance()->mesh.algoRecombine != 0){
 #if defined(HAVE_BLOSSOM)
     int ncount = gf->triangles.size();
     if (ncount % 2 == 0) {
@@ -3373,7 +3373,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
   std::set<MElement*> touched;
   std::vector<std::pair<MElement*,MElement*> > toProcess;
 
-  if(CTX::instance()->mesh.algoRecombine == 1){
+  if(CTX::instance()->mesh.algoRecombine != 0){
 #if defined(HAVE_BLOSSOM)
     int ncount = gf->triangles.size();
     if (ncount % 2 != 0) {
@@ -3619,7 +3619,7 @@ void recombineIntoQuads(GFace *gf,
     laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing);
   }
   // blossom-quad algo
-  if(success && CTX::instance()->mesh.algoRecombine == 1){
+  if(success && CTX::instance()->mesh.algoRecombine != 0){
     if(topologicalOpti){
       if(haveParam){
         if (saveAll) gf->model()->writeMSH("smoothed.msh");
diff --git a/Plugin/Triangulate.cpp b/Plugin/Triangulate.cpp
index 9ae174950b9c5642246312db70a2bdf956ba6cb3..aee7eb3934453425f62e85add0decd2cedfff81e 100644
--- a/Plugin/Triangulate.cpp
+++ b/Plugin/Triangulate.cpp
@@ -117,17 +117,32 @@ PView *GMSH_TriangulatePlugin::execute(PView *v)
   discreteFace *s = new discreteFace
     (GModel::current(), GModel::current()->getNumFaces() + 1);
   s->computeMeanPlane(points);
-  double plan[3][3];
-  s->getMeanPlaneData(plan);
-  for(unsigned int i = 0; i < points.size(); i++) project(points[i], plan);
+  double x, y, z, VX[3], VY[3];
+  s->getMeanPlaneData(VX, VY, x, y, z);
+  //  printf("mean plane %g %g %g (%g %g %g), (%g %g %g)\n",x,y,z,VX[0],VX[1],VX[2],VY[0],VY[1],VY[2]);
+  SBoundingBox3d bbox;
+  for(unsigned int i = 0; i < points.size(); i++) bbox += points[i]->point();
+  double lc = 10 * norm(SVector3(bbox.max(), bbox.min()));
+
+  std::map<MVertex*,SPoint3> temp;
+  for(unsigned int i = 0; i < points.size(); i++) {    
+    SPoint3 pp (points[i]->x(),points[i]->y(),points[i]->z());
+    temp[points[i]] = pp;
+    double u, v, vec[3] = {points[i]->x() - x, 
+			   points[i]->y() - y, 
+			   points[i]->z() - z};
+    prosca(vec, VX, &u);
+    prosca(vec, VY, &v);
+    points[i]->x() = u;
+    points[i]->y() = v;
+    points[i]->z() = 0;
+    //    printf("points[%d] = %g %g %g\n",i,points[i]->x() ,points[i]->y() ,points[i]->z()); 
+  }
   delete s;
 
 #if 0 // old code
 
   // get lc
-  SBoundingBox3d bbox;
-  for(unsigned int i = 0; i < points.size(); i++) bbox += points[i]->point();
-  double lc = 10 * norm(SVector3(bbox.max(), bbox.min()));
 
   // build a point record structure for the divide and conquer algorithm
   DocRecord doc(points.size());
@@ -191,7 +206,22 @@ PView *GMSH_TriangulatePlugin::execute(PView *v)
 #else // new code
   Msg::Info("Triangulating data points (new code)...");
   std::vector<MTriangle*> tris;
-  delaunayMeshIn2D(points, tris, true, 0, false);
+  for(unsigned int i = 0; i < points.size(); i++) {    
+    double XX = 1.e-12 * lc * (double)rand() / (double)RAND_MAX;
+    double YY = 1.e-12 * lc * (double)rand() / (double)RAND_MAX;
+    double ZZ = 1.e-17 * lc * (double)rand() / (double)RAND_MAX;
+    points[i]->x() += XX;
+    points[i]->y() += YY;
+    points[i]->z() += ZZ;
+  }
+  delaunayMeshIn2D(points, tris, true, 0, true);
+  for(unsigned int i = 0; i < points.size(); i++){
+    SPoint3 pp = temp[points[i]];
+    points[i]->x() = pp.x();
+    points[i]->y() = pp.y();
+    points[i]->z() = pp.z();
+  }
+
   PView *v2 = new PView();
   PViewDataList *data2 = getDataList(v2);
   for(unsigned int i = 0; i < tris.size(); i++){