diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 49d2af7e6b2c077bd87577d567b6ff0ce6f325e5..e4bbf7da66c9a01ec75ee2b035d4fda9d7418e9b 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -24,7 +24,6 @@ CenterlineField.cpp
     BoundaryLayers.cpp 
     BDS.cpp 
     HighOrder.cpp 
-      highOrderSmoother.cpp 
       highOrderTools.cpp 
     meshPartition.cpp
     meshRefine.cpp
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 1a9049d7bec2fec732d106907768f094a7f29a1f..619ef0b109628e54369a1f9df1b4080daf7182e4 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -8,7 +8,6 @@
 //
 
 #include "HighOrder.h"
-#include "highOrderSmoother.h"
 #include "highOrderTools.h"
 #include "MLine.h"
 #include "MTriangle.h"
@@ -55,7 +54,7 @@ static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N
 		p0.y() + i * du * (p1.y()-p0.y()),
 		p0.z() + i * du * (p1.z()-p0.z()));
     double t;
-    GPoint gp = ge->closestPoint(pi, t);		
+    GPoint gp = ge->closestPoint(pi, t);
     u[i] = gp.u();
     //    printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
   }  
@@ -257,8 +256,7 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
 
 static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
                             edgeContainer &edgeVertices, bool linear,
-                            int nPts = 1, highOrderSmoother *displ2D = 0,
-                            highOrderSmoother *displ3D = 0)
+                            int nPts = 1)
 {
   if(ge->geomType() == GEntity::DiscreteCurve ||
      ge->geomType() == GEntity::BoundaryLayerCurve)
@@ -280,7 +278,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
       bool reparamOK = true;
       if(!linear) {
         reparamOK &= reparamMeshVertexOnEdge(v0, ge, u0);
-        if(ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0])
+        if(ge->periodic(0) && ge->getEndVertex()->getNumMeshVertices() > 0 && v1 == ge->getEndVertex()->mesh_vertices[0])
           u1 = ge->parBounds(0).high();
         else
           reparamOK &= reparamMeshVertexOnEdge(v1, ge, u1);
@@ -321,11 +319,6 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
 	  GPoint pc = ge->point(US[count]);
           v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,pc.u()); 
 	  //	  printf("Edge %d(%g) %d(%g) new vertex %g %g at %g\n",v0->getNum(),u0,v1->getNum(),u1,v->x(),v->y(), US[count]);
-          if (displ2D || displ3D){
-            SPoint3 pc2 = edge.interpolate(t);          
-            if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-            if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-          }
         }
         temp.push_back(v);
         // this destroys the ordering of the mesh vertices on the edge
@@ -343,9 +336,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
 
 static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
                             edgeContainer &edgeVertices, bool linear,
-                            int nPts = 1, 
-                            highOrderSmoother *displ2D = 0,
-                            highOrderSmoother *displ3D = 0)
+                            int nPts = 1)
 {
   if(gf->geomType() == GEntity::DiscreteSurface ||
      gf->geomType() == GEntity::BoundaryLayerSurface)
@@ -383,11 +374,6 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
         else{
           GPoint pc = gf->point(US[j + 1], VS[j + 1]);
           v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]);
-          if (displ2D || displ3D){
-            SPoint3 pc2 = edge.interpolate(t);          
-            if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-            if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-          }
         }
         temp.push_back(v);
         gf->mesh_vertices.push_back(v);
@@ -403,8 +389,7 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
 
 static void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
                             edgeContainer &edgeVertices, bool linear,
-                            int nPts = 1, highOrderSmoother *displ2D = 0,
-                            highOrderSmoother *displ3D = 0)
+                            int nPts = 1)
 {
   for(int i = 0; i < ele->getNumEdges(); i++){
     MEdge edge = ele->getEdge(i);
@@ -435,8 +420,7 @@ static void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
 
 static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele, 
                             std::vector<MVertex*> &vf, faceContainer &faceVertices, 
-                            bool linear, int nPts = 1, highOrderSmoother *displ2D = 0,
-                            highOrderSmoother *displ3D = 0)
+                            bool linear, int nPts = 1)
 {
   if(gf->geomType() == GEntity::DiscreteSurface ||
      gf->geomType() == GEntity::BoundaryLayerSurface)
@@ -492,11 +476,6 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
           else{
             v = new MVertex(X, Y, Z, gf);
           }
-          if(displ3D || displ2D){
-            SPoint3 pc2 = face.interpolate(t1, t2);
-            if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-            if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
-          }       
         }
         // should be expensive -> induces a new search each time
         vtcs.push_back(v);
@@ -775,6 +754,25 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
     }
     start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ;
     break;
+  case TYPE_PYR:
+    printf("aaaaaand...\n");
+    switch (nPts){
+    case 0: 
+    case 1: return;
+    case 2: points = polynomialBases::find(MSH_PYR_30)->points; break;
+    case 3: points = polynomialBases::find(MSH_PYR_55)->points; break;
+    case 4: points = polynomialBases::find(MSH_PYR_91)->points; break;
+    case 5: points = polynomialBases::find(MSH_PYR_140)->points; break;
+    case 6: points = polynomialBases::find(MSH_PYR_204)->points; break;
+    case 7: points = polynomialBases::find(MSH_PYR_285)->points; break;
+    case 8: points = polynomialBases::find(MSH_PYR_385)->points; break;
+    default :  
+      Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
+      break;
+    }
+    start = ( nPts+2) * ( (nPts+2) + 1) * (2*(nPts+2) + 1) / 6  - (nPts-1) * ( (nPts-1) + 1) * (2*(nPts-1) + 1) / 6;
+    break;
+
   }
 
   for(int k = start; k < points.size1(); k++){
@@ -791,14 +789,13 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
 }
 
 static void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear, 
-                         int nbPts = 1, highOrderSmoother *displ2D = 0,
-                         highOrderSmoother *displ3D = 0)
+                         int nbPts = 1)
 {
   std::vector<MLine*> lines2;
   for(unsigned int i = 0; i < ge->lines.size(); i++){
     MLine *l = ge->lines[i];
     std::vector<MVertex*> ve;
-    getEdgeVertices(ge, l, ve, edgeVertices, linear, nbPts, displ2D, displ3D);
+    getEdgeVertices(ge, l, ve, edgeVertices, linear, nbPts);
     if(nbPts == 1)
       lines2.push_back(new MLine3(l->getVertex(0), l->getVertex(1), ve[0]));
     else
@@ -812,12 +809,10 @@ static void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
 MTriangle *setHighOrder(MTriangle *t, GFace *gf, 
                         edgeContainer &edgeVertices, 
                         faceContainer &faceVertices, 
-                        bool linear, bool incomplete, int nPts, 
-                        highOrderSmoother *displ2D,
-                        highOrderSmoother *displ3D)
+                        bool linear, bool incomplete, int nPts)
 {
   std::vector<MVertex*> ve, vf;
-  getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+  getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts);
   if(nPts == 1){
     return new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),
                           ve[0], ve[1], ve[2]);
@@ -828,17 +823,9 @@ MTriangle *setHighOrder(MTriangle *t, GFace *gf,
                             ve, nPts + 1);
     }
     else{
-      if (displ2D && gf->geomType() == GEntity::Plane){
-        MTriangle incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2));
-        getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts, displ2D, 
-                        displ3D);
-      }
-      else{
-        MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
-                         ve, nPts + 1);
-        getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts, displ2D, 
-                        displ3D);
-      }
+      MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+		       ve, nPts + 1);
+      getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
       ve.insert(ve.end(), vf.begin(), vf.end());
       return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
                             ve, nPts + 1);
@@ -849,12 +836,10 @@ MTriangle *setHighOrder(MTriangle *t, GFace *gf,
 static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
                                  edgeContainer &edgeVertices, 
                                  faceContainer &faceVertices, 
-                                 bool linear, bool incomplete, int nPts, 
-                                 highOrderSmoother *displ2D,
-                                 highOrderSmoother *displ3D)
+                                 bool linear, bool incomplete, int nPts)
 {
   std::vector<MVertex*> ve, vf;
-  getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+  getEdgeVertices(gf, q, ve, edgeVertices, linear, nPts);
   if(incomplete){
     if(nPts == 1){
       return new MQuadrangle8(q->getVertex(0), q->getVertex(1), q->getVertex(2), 
@@ -866,15 +851,9 @@ static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
     }
   } 
   else {
-    if (displ2D && gf->geomType() == GEntity::Plane){
-      MQuadrangle incpl(q->getVertex(0), q->getVertex(1), q->getVertex(2), 
-                        q->getVertex(3));
-      getFaceVertices(gf, &incpl, q, vf, faceVertices, linear, nPts, displ2D, displ3D);
-    }else{
-      MQuadrangleN incpl(q->getVertex(0), q->getVertex(1), q->getVertex(2), 
-                         q->getVertex(3), ve, nPts + 1);
-      getFaceVertices(gf, &incpl, q, vf, faceVertices, linear, nPts, displ2D, displ3D);
-    }
+    MQuadrangleN incpl(q->getVertex(0), q->getVertex(1), q->getVertex(2), 
+		       q->getVertex(3), ve, nPts + 1);
+    getFaceVertices(gf, &incpl, q, vf, faceVertices, linear, nPts);
     ve.insert(ve.end(), vf.begin(), vf.end());
     if(nPts == 1){
       return new MQuadrangle9(q->getVertex(0), q->getVertex(1), q->getVertex(2),
@@ -889,14 +868,13 @@ static MQuadrangle *setHighOrder(MQuadrangle *q, GFace *gf,
 
 static void setHighOrder(GFace *gf, edgeContainer &edgeVertices, 
                          faceContainer &faceVertices, bool linear, bool incomplete,
-                         int nPts = 1, highOrderSmoother *displ2D = 0,
-                         highOrderSmoother *displ3D = 0)
+                         int nPts = 1)
 {
   std::vector<MTriangle*> triangles2;
   for(unsigned int i = 0; i < gf->triangles.size(); i++){
     MTriangle *t = gf->triangles[i];
     MTriangle *tNew = setHighOrder(t, gf,edgeVertices, faceVertices, linear, incomplete,
-                                   nPts, displ2D, displ3D);
+                                   nPts);
     triangles2.push_back(tNew);
     delete t;
   }
@@ -905,7 +883,7 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
   for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
     MQuadrangle *q = gf->quadrangles[i];
     MQuadrangle *qNew = setHighOrder(q, gf, edgeVertices, faceVertices, linear,
-                                     incomplete, nPts, displ2D, displ3D);
+                                     incomplete, nPts);
     quadrangles2.push_back(qNew);
     delete q;
   }
@@ -915,14 +893,13 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
 
 static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices, 
                          faceContainer &faceVertices, bool linear, bool incomplete,
-                         int nPts = 1, highOrderSmoother *displ2D = 0,
-                         highOrderSmoother *displ3D = 0)
+                         int nPts = 1)
 {
   std::vector<MTetrahedron*> tetrahedra2;
   for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
     MTetrahedron *t = gr->tetrahedra[i];
     std::vector<MVertex*> ve, vf, vr;
-    getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+    getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
     if(nPts == 1){
       tetrahedra2.push_back
         (new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2), 
@@ -947,7 +924,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
   for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
     MHexahedron *h = gr->hexahedra[i];
     std::vector<MVertex*> ve, vf, vr;
-    getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+    getEdgeVertices(gr, h, ve, edgeVertices, linear, nPts);
     if(nPts == 1){
       if(incomplete){
 	hexahedra2.push_back
@@ -990,8 +967,8 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
   std::vector<MPrism*> prisms2;
   for(unsigned int i = 0; i < gr->prisms.size(); i++){
     MPrism *p = gr->prisms[i];
-    std::vector<MVertex*> ve, vf;
-    getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+    std::vector<MVertex*> ve, vf, vr;
+    getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
     if(incomplete){
       prisms2.push_back
         (new MPrism15(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
@@ -999,12 +976,31 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
                       ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8]));
     }
     else{
-      getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
-      prisms2.push_back
-        (new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
-                      p->getVertex(3), p->getVertex(4), p->getVertex(5), 
-                      ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
-                      vf[0], vf[1], vf[2]));
+      if (nPts == 1) {
+	getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+	prisms2.push_back
+	  (new MPrism18(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+			p->getVertex(3), p->getVertex(4), p->getVertex(5), 
+			ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
+			vf[0], vf[1], vf[2]));
+      } else {
+
+	//getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+	//ve.insert(ve.end(), vf.begin(), vf.end());
+	//MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2), p->getVertex(3),
+	//		   p->getVertex(4), p->getVertex(5), ve, nPts + 1);
+	//getRegionVertices(gr, &incpl, p, vr, linear, nPts);
+	//if (nPts == 0) {
+	//   printf("Screwed\n");
+	//  }
+	//ve.insert(ve.end(), vr.begin(), vr.end());
+	//MPrism* n = new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+	//				  p->getVertex(3), p->getVertex(4), p->getVertex(5),
+	//				  ve, nPts+1);
+	//if (!mappingIsInvertible(n))
+	//  Msg::Warning("Found invalid curved volume element (# %d in list)", i);
+	//prisms2.push_back(n);
+      }
     }
     delete p;
   }
@@ -1013,20 +1009,35 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
   std::vector<MPyramid*> pyramids2;
   for(unsigned int i = 0; i < gr->pyramids.size(); i++){
     MPyramid *p = gr->pyramids[i];
-    std::vector<MVertex*> ve, vf;
-    getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts, displ2D, displ3D);
-    if(incomplete){
-      pyramids2.push_back
-        (new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
-                        p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
-                        ve[3], ve[4], ve[5], ve[6], ve[7]));
+    std::vector<MVertex*> ve, vf, vr;
+    getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
+    if (nPts == 1) {
+      if(incomplete){
+	pyramids2.push_back
+	  (new MPyramid13(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+			  p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
+			  ve[3], ve[4], ve[5], ve[6], ve[7]));
+      }
+      else{
+	getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+	pyramids2.push_back
+	  (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+			  p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
+			  ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
+      }
     }
-    else{
+    else {
       getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
-      pyramids2.push_back
-        (new MPyramid14(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
-                        p->getVertex(3), p->getVertex(4), ve[0], ve[1], ve[2], 
-                        ve[3], ve[4], ve[5], ve[6], ve[7], vf[0]));
+      ve.insert(ve.end(), vf.begin(), vf.end());     
+      MPyramidN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+                         p->getVertex(3), p->getVertex(4), ve, nPts + 1);
+      getRegionVertices(gr, &incpl, p, vr, linear, nPts); 
+      ve.insert(ve.end(), vr.begin(), vr.end());
+      printf("Hmmm, let's see : %d\n", nPts);
+      MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1), p->getVertex(2), 
+				  p->getVertex(3), p->getVertex(4), ve, nPts + 1);
+      pyramids2.push_back(n);
+
     }
     delete p;
   }
@@ -1247,25 +1258,30 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
 
   // m->writeMSH("BEFORE.msh");
 
-  highOrderSmoother *displ2D = 0; 
-  highOrderSmoother *displ3D = 0; 
+  /*
   if(CTX::instance()->mesh.smoothInternalEdges){
     displ2D = new highOrderSmoother(2);
     displ3D = new highOrderSmoother(3);
   }
+  */
 
   // then create new second order vertices/elements
   edgeContainer edgeVertices;
   faceContainer faceVertices;
   
-  Msg::StatusBar(2, true, "Meshing curves order %d...", order);
-  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
-    setHighOrder(*it, edgeVertices, linear, nPts, displ2D, displ3D);
+  int counter = 1;
+  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
+    Msg::StatusBar(2, true, "Meshing curves order %d (%i/%i)...", order, counter, m->getNumEdges());
+    counter++;
+    setHighOrder(*it, edgeVertices, linear, nPts);
+  }
 
-  Msg::StatusBar(2, true, "Meshing surfaces order %d...", order);
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
-    setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,
-                 displ2D, displ3D);
+  counter = 1;
+  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
+    Msg::StatusBar(2, true, "Meshing surfaces order %d (%i/%i)...", order, counter, m->getNumFaces());
+    counter++;
+    setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
+  }
 
   highOrderTools hot(m);
 
@@ -1277,18 +1293,22 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
   // printJacobians(m, "smoothness_b.pos");
   // m->writeMSH("RAW.msh");
 
-  if (displ2D){
+  /*if (displ2D){
     checkHighOrderTriangles("Before optimization", m, bad, worst);
     for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
       displ2D->optimize(*it,edgeVertices,faceVertices);
     checkHighOrderTriangles("After optimization", m, bad, worst);
-  }
+    }
+  */
 
-  Msg::StatusBar(2, true, "Meshing volumes order %d...", order);
-  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
-    setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,
-                 displ2D, displ3D);
+  counter = 1;
+  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
+    Msg::StatusBar(2, true, "Meshing volumes order %d (%i/%i)...", order, counter, m->getNumRegions());
+    counter++;
+    setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
+  }
 
+/*
   // smooth the 3D regions
   if (displ3D){
     checkHighOrderTetrahedron("Before optimization", m, bad, worst);
@@ -1296,9 +1316,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
       displ3D->smooth(*it);
     checkHighOrderTetrahedron("After optimization", m, bad, worst);
   }
-
-  if(displ2D) delete displ2D;
-  if(displ3D) delete displ3D;
+*/
 
   double t2 = Cpu();
 
@@ -1306,7 +1324,13 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
   // m->writeMSH("SMOOTHED.msh");
   // FIXME !!
   if (!linear){
-  //hot.ensureMinimumDistorsion(0.1);
+    for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
+      std::vector<MElement*> v;
+      v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
+      v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
+      hot.applySmoothingTo(v, -10000000., true);
+    }
+    hot.ensureMinimumDistorsion(0.1);
     checkHighOrderTriangles("Final mesh", m, bad, worst);
     checkHighOrderTetrahedron("Final mesh", m, bad, worst);
   }
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index e886e91601ee39e216a08ba49d900f564dc56e76..a0c0e1a068b130fa06b20f02db9bb0b2e054da7c 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -19,16 +19,12 @@ typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MVertex*> > edgeCont
 // are the high order representation of the face
 typedef std::map<MFace, std::vector<MVertex*>, Less_Face> faceContainer;
 
-class highOrderSmoother;
-
 void SetOrder1(GModel *m);
 void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
 MTriangle* setHighOrder(MTriangle *t, GFace *gf, 
                         edgeContainer &edgeVertices, 
                         faceContainer &faceVertices, 
-                        bool linear, bool incomplete, int nPts = 1, 
-                        highOrderSmoother *displ2D = 0,
-                        highOrderSmoother *displ3D = 0);
+                        bool linear, bool incomplete, int nPts = 1);
 void checkHighOrderTriangles(const char* cc, GModel *m, 
                              std::vector<MElement*> &bad, double &minJGlob);