diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 5aafe141e37b02933acc4b60ae5350eec84eddd6..6cb9df2a45844487715b5be901f90f846d87c637 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -631,7 +631,8 @@ void GEdge::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<doubl
   }
 }
 
-void GEdge::mesh(bool verbose){
+void GEdge::mesh(bool verbose)
+{
 #if defined(HAVE_MESH)
   meshGEdge mesher;
   mesher(this);
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 60df47ebe44e8ce25406ed2986206ab1e52a0edf..c7937afcdf0c840088b787eb79e9d318e1644f4c 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -254,27 +254,6 @@ bool GEdgeCompound::getLocalParameter(const double &t,
   return false;
 }
 
-void GEdgeCompound::getCompoundParameter(GEdge *ge,
-                                         const double &tLoc,
-                                         double &t) const
-{
-  if(_pars.empty()){
-    Msg::Error("Edge compound has no parametrization");
-    return;
-  }
-  for (int iEdge = 0; iEdge < (int)_compound.size(); iEdge++){
-    if (ge == _compound[iEdge]){
-      double tmin = _pars[iEdge];
-      double tmax = _pars[iEdge+1];
-      Range<double> b = _compound[iEdge]->parBounds(0);
-      t = _orientation[iEdge] ?
-        tmin + (tLoc - b.low())/(b.high()-b.low()) * (tmax-tmin):
-        tmax - (tLoc - b.low())/(b.high()-b.low()) * (tmax-tmin);
-      return;
-    }
-  }
-}
-
 void GEdgeCompound::parametrize()
 {
   _pars.clear();
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index 3efdc13a302dc420b3acd54e7249865bf1bcf3a3..390df7a35179f12bbb3320d8e27ae0f12d3b67b1 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -20,7 +20,6 @@ class GEdgeCompound : public GEdge {
  public:
   void parametrize();
   bool getLocalParameter(const double &t, int &iEdge, double & tLoc) const;
-  void getCompoundParameter(GEdge *ge, const double &tLoc, double &t) const;
   GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound);
   GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound,
                 std::vector<int> &orientation); // confidence in the input
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 739d988edc3a14edd895e19a4ae7f12d2dd57ba2..fe88a576fa15c34fc9ef0631ed42bd215b87f6d8 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1388,14 +1388,14 @@ bool GFace::fillPointCloud(double maxDist,
   return true;
 }
 
-void GFace::mesh(bool verbose) {
+void GFace::mesh(bool verbose)
+{
 #if defined(HAVE_MESH)
   meshGFace mesher;
-  mesher (this, verbose);
+  mesher(this, verbose);
 #endif
 }
 
-
 void GFace::lloyd(int nbiter, int infn)
 {
 #if defined(HAVE_MESH) && defined(HAVE_BFGS)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index e9e5bdfff12f4f21ec267565cdc98073be216d7d..0b09589ee0efacc92b2d8461764acfdfafaf07ab 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -956,7 +956,7 @@ bool GFaceCompound::parametrize() const
     _rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
 
     linearSystemPETSc<double> sys;
-	  printf("system 0 = %p\n", &sys);
+
     #if 1
     _rbf->RbfLapSurface_local_CPM_sparse(_ordered, false, _rbf->getXYZ(), _rbf->getN(), sys);
     _rbf->solveHarmonicMap_sparse(sys, _rbf->getXYZ().size1()* 3,_ordered, _coords, coordinates);
@@ -1432,63 +1432,82 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
 	gec->getLocalParameter(tGlob,iEdge,tLoc);
 	std::vector<GEdge*> gev = gec->getCompounds();
 	GEdge *ge = gev[iEdge];
-
-	// left and right vertex of the Edge
-	MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
-	MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
-	std::map<MVertex*,SPoint3>::iterator itL = coordinates.find(v0);
-	std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
-
-	// for the Edge, find the left and right vertices of the initial
-	// 1D mesh and interpolate to find (u,v)
-	//MVertex *vL = v0;
-	MVertex *vR = v1;
-	double tB = ge->parBounds(0).low();
-	double tE = ge->parBounds(0).high();
-	int j = 0;
-	tL=tB;
-	bool found = false;
-	while(j < (int)ge->mesh_vertices.size()){
-	  vR = ge->mesh_vertices[j];
-          if(vR->getPolynomialOrder() > 1){ j++; continue; }
-	  vR->getParameter(0,tR);
-	  if(!vR->getParameter(0,tR)) {
-	    Msg::Error("Vertex vr %p not an MEdgeVertex", vR);
-	    return SPoint2();
-	  }
-	  if(tLoc > tL && tLoc < tR){
-	    found = true;
-	    itR = coordinates.find(vR);
-	    if(itR == coordinates.end()){
-	      Msg::Error("Vertex %p (%g %g %g) not found", vR, vR->x(), vR->y(), vR->z());
-              return SPoint2(0,0);
-	    }
-	    break;
-	  }
-	  else{
-	    itL = coordinates.find(vR);
-	    //vL = vR;
-	    tL = tR;
-	  }
-	  j++;
-	}
-	if(!found) {
-	  vR = v1;
-	  tR = tE;
-	}
-
-	// linear interpolation between tL and tR
-	double uloc, vloc;
-	uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
-	vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
-	return SPoint2(uloc,vloc);
+        std::map<MVertex*,SPoint3>::iterator itL;
+        std::map<MVertex*,SPoint3>::iterator itR;
+
+        if(ge->geomType() == GEntity::DiscreteCurve){
+          discreteEdge *de = dynamic_cast<discreteEdge*>(ge);
+          int i = (int)tLoc;
+          MLine *l = de->lines[i];
+          MVertex *vL = l->getVertex(0);
+          MVertex *vR = l->getVertex(1);
+          itL = coordinates.find(vL);
+          if(itL == coordinates.end()){
+            Msg::Error("Vertex %p (%g %g %g) not found", vL, vL->x(), vL->y(), vL->z());
+            return SPoint2(0, 0);
+          }
+          itR = coordinates.find(vR);
+          if(itR == coordinates.end()){
+            Msg::Error("Vertex %p (%g %g %g) not found", vR, vR->x(), vR->y(), vR->z());
+            return SPoint2(0, 0);
+          }
+          double xi, uloc, vloc;
+          xi = tLoc - i;
+          uloc = (1 - xi) * itL->second.x() + xi * itR->second.x();
+          vloc = (1 - xi) * itL->second.y() + xi * itR->second.y();
+          return SPoint2(uloc,vloc);
+        }
+        else{
+          // left and right vertex of the Edge
+          MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
+          MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
+          itL = coordinates.find(v0);
+          itR = coordinates.find(v1);
+          // for the Edge, find the left and right vertices of the initial
+          // 1D mesh and interpolate to find (u,v)
+          //MVertex *vL = v0;
+          MVertex *vR = v1;
+          double tB = ge->parBounds(0).low();
+          double tE = ge->parBounds(0).high();
+          int j = 0;
+          tL=tB;
+          bool found = false;
+          while(j < (int)ge->mesh_vertices.size()){
+            vR = ge->mesh_vertices[j];
+            if(vR->getPolynomialOrder() > 1){ j++; continue; }
+            vR->getParameter(0, tR);
+            if(!vR->getParameter(0, tR)) {
+              Msg::Error("Vertex vr %p not an MEdgeVertex", vR);
+              return SPoint2();
+            }
+            if(tLoc > tL && tLoc < tR){
+              found = true;
+              itR = coordinates.find(vR);
+              if(itR == coordinates.end()){
+                Msg::Error("Vertex %p (%g %g %g) not found", vR, vR->x(), vR->y(), vR->z());
+                return SPoint2(0, 0);
+              }
+              break;
+            }
+            else{
+              itL = coordinates.find(vR);
+              //vL = vR;
+              tL = tR;
+            }
+            j++;
+          }
+          if(!found) {
+            vR = v1;
+            tR = tE;
+          }
+          // linear interpolation between tL and tR
+          double uloc, vloc;
+          uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
+          vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
+          return SPoint2(uloc,vloc);
+        }
       }
 
-    // }
-    // else if(v->onWhat()->dim() == 2){
-    //   Msg::Debug("Get coordinates of Compound Face cannot find vertex");
-    // }
-
   }
 
   // never here
@@ -1517,7 +1536,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
     lsys = new linearSystemPETSc<double>;
     lsys->setParameter("petscOptions",
     "-pc_type ilu -ksp_rtol 1.e-12 -pc_factor_levels 10");
-    
+
     //    lsys->setParameter("petscOptions",
     //    		       "-ksp_type preonly -pc_type lu -pc_factor_mat_solver_package umfpack");
 #elif defined(HAVE_GMM)
@@ -1585,7 +1604,6 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
     }
   }
   else if (_type == ALREADYFIXED){
-    printf("already fixed boundaries \n");
     for(unsigned int i = 0; i < _ordered.size(); i++){
       MVertex *v = _ordered[i];
       SPoint3 uv = coordinates[v];
@@ -2006,7 +2024,7 @@ SPoint2 GFaceCompound::parFromVertex(MVertex *v) const
     v->getParameter(1, V);
   }
   if (v->onWhat()->dim()==1 ||
-     (v->onWhat()->dim()==2 && 
+     (v->onWhat()->dim()==2 &&
       U == -1 && V==-1)){ //if MFaceVertex created on edge in bunin
     SPoint2 sp = getCoordinates(v);
     U = sp.x();
@@ -2035,7 +2053,6 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
 
   //create new octree with new mesh elements
   if(!octNew){
-    //printf("create new octree \n");
     std::vector<MElement *> myElems;
     for (unsigned int i = 0; i < triangles.size(); i++) myElems.push_back(triangles[i]);
     for (unsigned int i = 0; i < quadrangles.size(); i++) myElems.push_back(quadrangles[i]);
@@ -2067,11 +2084,9 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
       }
     }
 
-    //printf("size octrre new = %d \n", myParamElems.size());
     octNew = new MElementOctree(myParamElems);
 
     //build kdtree boundary nodes in parametric space
-    //printf("build bc kdtree \n");
     int nbBCNodes  = myBCNodes.size();
     ANNpointArray uv_nodes = annAllocPts(nbBCNodes, 3);
     std::set<SPoint2>::iterator itp = myBCNodes.begin();
@@ -2107,7 +2122,6 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
   }
   //if element not found in new octree find closest point
   else{
-    //printf("not found in new octree \n");
     GPoint gp(50,50,50);
     double pt[3] = {par1,par2,0.0};
     ANNidx index[2];
@@ -2163,15 +2177,12 @@ GPoint GFaceCompound::point(double par1, double par2) const
   GFaceCompoundTriangle *lt;
   getTriangle(par1, par2, &lt, U,V);
   if(!lt){
-    //printf("POINT (%g %g) NOT FOUND --> find closest \n", par1,par2);
     if (meshStatistics.status == GFace::DONE && triangles.size()+quadrangles.size() != 0) {
-      //printf("look in remeshed octree \n");
       GPoint gp = pointInRemeshedOctree(par1,par2);
       gp.setNoSuccess();
       return gp;
     }
     else{
-      //printf("look in kdtree \n");
       double pt[3] = {par1,par2, 0.0};
       ANNidx index[2];
       ANNdist dist[2];
@@ -2199,7 +2210,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
       sys2x2(M, R, X);
       U = X[0];
       V = X[1];
-      //printf("found closest point (%g %g) U V =%g %g \n", u,v, U,V);
     }
   }
 
@@ -2277,8 +2287,6 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 &param) const
   MTriangle *tri=NULL;
   if (lt) tri = lt->tri;
   else {
-    //printf("FIRSTDER POINT NOT FOUND --> kdtree \n");
-    //printf("uv=%g %g \n", param.x(), param.y());
     double pt[3] = {param.x(), param.y(), 0.0};
     ANNidx index[2];
     ANNdist dist[2];
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 7e8347dc917ff708c4251b46aca77375a89e346f..df2976c2ab405d30e3d1979a8e1a8e4b37d25fe7 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1276,9 +1276,10 @@ static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &el
   }
 }
 
-void GModel:: _createGeometryOfDiscreteEntities() {
-  if (CTX::instance()->meshDiscrete){
-    Msg::Info("creating the geometry of discrete curves");
+void GModel::_createGeometryOfDiscreteEntities(bool force)
+{
+  if (force || CTX::instance()->meshDiscrete){
+    Msg::Info("Creating the geometry of discrete curves");
     for(eiter it = firstEdge(); it != lastEdge(); ++it){
       if((*it)->geomType() == GEntity::DiscreteCurve) {
 	discreteEdge *de = dynamic_cast<discreteEdge*> (*it);
@@ -1286,7 +1287,10 @@ void GModel:: _createGeometryOfDiscreteEntities() {
 	de->createGeometry();
       }
     }
-    Msg::Info("creating the geometry of discrete surfaces");
+  }
+
+  if (CTX::instance()->meshDiscrete){
+    Msg::Info("Creating the geometry of discrete surfaces");
     for(fiter it = firstFace(); it != lastFace(); ++it){
       if((*it)->geomType() == GEntity::DiscreteSurface) {
 	discreteFace *df = dynamic_cast<discreteFace*> (*it);
@@ -2078,7 +2082,13 @@ void GModel::createTopologyFromMesh(int ignoreHoles)
   //create old format (necessary e.g. for old-style extruded boundary layers)
   exportDiscreteGEOInternals();
 
+  // FIXME: this whole thing will disappear, but for now we need this to make
+  // old compounds work:
+  if(!CTX::instance()->meshDiscrete)
+    _createGeometryOfDiscreteEntities(true);
+
   double t2 = Cpu();
+
   Msg::StatusBar(true, "Done creating topology from mesh (%g s)", t2 - t1);
 }
 
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 4329431d24423e8698685b28f32c84d26f19269e..9beb563afedc1cbc43fa1637d4408dddde11383b 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -120,7 +120,7 @@ class GModel
   void _storeElementsInEntities(std::map<int, std::vector<MElement*> > &map);
 
   // Discrete Entities have to have their mesh moved to a geometry container
-  void _createGeometryOfDiscreteEntities();
+  void _createGeometryOfDiscreteEntities(bool force=false);
 
   // loop over all vertices connected to elements and associate
   // geometrical entity
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index c6bf343a7b16698777c593942c89adf20d20cd0c..1d873ca5ced91f4a99806e648c2a7dcbc2af0d7c 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -467,8 +467,7 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
 bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 &param,
                              bool onSurface)
 {
-
-  if (gf->geomType() == GEntity::CompoundSurface ){
+  if(gf->geomType() == GEntity::CompoundSurface){
     GFaceCompound *gfc = (GFaceCompound*) gf;
     param = gfc->parFromVertex(v);
     return true;
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 3b0b037410f47b15ea234241e5bbd6f842158bd2..73ab50c8e4e93112568ddc6c27f6556f105290ed 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -36,7 +36,7 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
 
 void discreteEdge::createTopo()
 {
-  if(!createdTopo){ 
+  if(!createdTopo){
     orderMLines();
     setBoundVertices();
     createdTopo = true;
@@ -58,7 +58,7 @@ void discreteEdge::orderMLines()
   }
 
   // find a lonly MLine
-  for (std::list<MLine*>::iterator it = segments.begin(); 
+  for (std::list<MLine*>::iterator it = segments.begin();
        it != segments.end(); ++it){
     MVertex *vL = (*it)->getVertex(0);
     MVertex *vR = (*it)->getVertex(1);
@@ -74,7 +74,7 @@ void discreteEdge::orderMLines()
   MLine *firstLine;
   if (boundv.size() == 2){   // non periodic
     firstLine = (boundv.begin())->second;
-    for (std::list<MLine*>::iterator it = segments.begin(); 
+    for (std::list<MLine*>::iterator it = segments.begin();
          it != segments.end(); ++it){
       if (*it == firstLine){
         segments.erase(it);
@@ -87,7 +87,7 @@ void discreteEdge::orderMLines()
     segments.erase(segments.begin());
   }
   else{
-    Msg::Error("EdgeCompound %d is wrong (it has %d end points)", 
+    Msg::Error("EdgeCompound %d is wrong (it has %d end points)",
                tag(), boundv.size());
   }
 
@@ -99,7 +99,7 @@ void discreteEdge::orderMLines()
   while (first != last){
     if (segments.empty())break;
     bool found = false;
-    for (std::list<MLine*>::iterator it = segments.begin(); 
+    for (std::list<MLine*>::iterator it = segments.begin();
          it != segments.end(); ++it){
       MLine *e = *it;
       if (e->getVertex(0) == last){
@@ -135,15 +135,15 @@ void discreteEdge::orderMLines()
 
   //lines is now a list of ordered MLines
   lines = _m;
-  
+
   //mesh_vertices
   mesh_vertices.clear();
   for (unsigned int i = 0; i < lines.size(); ++i){
     MVertex *v1 = lines[i]->getVertex(0);
     MVertex *v2 = lines[i]->getVertex(1);
-    if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v1) ==  
+    if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v1) ==
         mesh_vertices.end()) mesh_vertices.push_back(v1);
-    if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v2) ==  
+    if (std::find(mesh_vertices.begin(), mesh_vertices.end(), v2) ==
         mesh_vertices.end()) mesh_vertices.push_back(v2);
   }
 
@@ -151,7 +151,7 @@ void discreteEdge::orderMLines()
   if (lines.size() < 2) return;
   if (_orientation[0] && lines[0]->getVertex(1) != lines[1]->getVertex(1)
       && lines[0]->getVertex(1) != lines[1]->getVertex(0)){
-    for (unsigned int i = 0; i < lines.size(); i++) 
+    for (unsigned int i = 0; i < lines.size(); i++)
       _orientation[i] = !_orientation[i];
   }
 }
@@ -160,11 +160,11 @@ void discreteEdge::setBoundVertices()
 {
   if (boundv.size() == 2){
     std::vector<GVertex*> bound_vertices;
-    for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); 
+    for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin();
          iter != boundv.end(); iter++){
       MVertex* vE = (iter)->first;
       bool existVertex  = false;
-      for(GModel::viter point = model()->firstVertex(); 
+      for(GModel::viter point = model()->firstVertex();
           point != model()->lastVertex(); point++){
         if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
           bound_vertices.push_back((*point));
@@ -174,7 +174,7 @@ void discreteEdge::setBoundVertices()
       }
       if(!existVertex){
         GVertex *gvB = new discreteVertex
-          (model(), model()->getMaxElementaryNumber(0) + 1);      
+          (model(), model()->getMaxElementaryNumber(0) + 1);
         //printf("*** Created discreteVertex %d\n", gvB->tag());
         bound_vertices.push_back(gvB);
         vE->setEntity(gvB);
@@ -182,22 +182,22 @@ void discreteEdge::setBoundVertices()
         gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
         model()->add(gvB);
       }
-      std::vector<MVertex*>::iterator itve = 
+      std::vector<MVertex*>::iterator itve =
         std::find(mesh_vertices.begin(), mesh_vertices.end(), vE);
       if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
 
       for(GModel::eiter edge = model()->firstEdge(); edge != model()->lastEdge(); edge++){
-        std::vector<MVertex*>::iterator itve = 
+        std::vector<MVertex*>::iterator itve =
           std::find((*edge)->mesh_vertices.begin(), (*edge)->mesh_vertices.end(), vE);
         if (itve != (*edge)->mesh_vertices.end()) (*edge)->mesh_vertices.erase(itve);
       }
       for(GModel::fiter face = model()->firstFace(); face != model()->lastFace(); face++){
-        std::vector<MVertex*>::iterator itve = 
+        std::vector<MVertex*>::iterator itve =
           std::find((*face)->mesh_vertices.begin(), (*face)->mesh_vertices.end(), vE);
         if (itve != (*face)->mesh_vertices.end()) (*face)->mesh_vertices.erase(itve);
       }
       for(GModel::riter reg = model()->firstRegion(); reg != model()->lastRegion(); reg++){
-        std::vector<MVertex*>::iterator itve = 
+        std::vector<MVertex*>::iterator itve =
           std::find((*reg)->mesh_vertices.begin(), (*reg)->mesh_vertices.end(), vE);
         if (itve != (*reg)->mesh_vertices.end()) (*reg)->mesh_vertices.erase(itve);
       }
@@ -212,7 +212,7 @@ void discreteEdge::setBoundVertices()
     std::vector<MLine*>::const_iterator it = lines.begin();
     MVertex* vE = (*it)->getVertex(0);
     bool existVertex  = false;
-    for(GModel::viter point = model()->firstVertex(); 
+    for(GModel::viter point = model()->firstVertex();
         point != model()->lastVertex(); point++){
       if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
         bound_vertex = (*point);
@@ -231,25 +231,25 @@ void discreteEdge::setBoundVertices()
       gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
       model()->add(gvB);
     }
-    std::vector<MVertex*>::iterator itve = 
+    std::vector<MVertex*>::iterator itve =
       std::find(mesh_vertices.begin(),mesh_vertices.end(), vE);
     if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
-   
-    for(GModel::eiter edge = model()->firstEdge(); 
+
+    for(GModel::eiter edge = model()->firstEdge();
         edge != model()->lastEdge(); edge++){
-      std::vector<MVertex*>::iterator itve = 
+      std::vector<MVertex*>::iterator itve =
         std::find((*edge)->mesh_vertices.begin(), (*edge)->mesh_vertices.end(), vE);
       if (itve != (*edge)->mesh_vertices.end()) (*edge)->mesh_vertices.erase(itve);
     }
-    for(GModel::fiter face = model()->firstFace(); 
+    for(GModel::fiter face = model()->firstFace();
         face != model()->lastFace(); face++){
-      std::vector<MVertex*>::iterator itve = 
+      std::vector<MVertex*>::iterator itve =
         std::find((*face)->mesh_vertices.begin(), (*face)->mesh_vertices.end(), vE);
       if (itve != (*face)->mesh_vertices.end()) (*face)->mesh_vertices.erase(itve);
     }
-    for(GModel::riter reg = model()->firstRegion(); 
+    for(GModel::riter reg = model()->firstRegion();
         reg != model()->lastRegion(); reg++){
-      std::vector<MVertex*>::iterator itve = 
+      std::vector<MVertex*>::iterator itve =
         std::find((*reg)->mesh_vertices.begin(), (*reg)->mesh_vertices.end(), vE);
       if (itve != (*reg)->mesh_vertices.end()) (*reg)->mesh_vertices.erase(itve);
     }
@@ -267,10 +267,10 @@ void discreteEdge::setBoundVertices()
     _pars[0]=0   _pars[1]=1    _pars[2]=2             _pars[N+1]=N+1
 */
 void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
-                               std::less<MVertex*> > > &face2Vert, 
+                               std::less<MVertex*> > > &face2Vert,
                                std::map<GRegion*, std::map<MVertex*, MVertex*,
                                std::less<MVertex*> > > &region2Vert)
-{ 
+{
   return;
   if (_pars.empty()){
     for (unsigned int i = 0; i < lines.size() + 1; i++){
@@ -278,11 +278,11 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
     }
   }
   //Replace MVertex by MedgeVertex
-  std::map<MVertex*, MVertex*, std::less<MVertex*> > old2new;   
+  std::map<MVertex*, MVertex*, std::less<MVertex*> > old2new;
   old2new.clear();
   std::vector<MVertex* > newVertices;
   std::vector<MLine*> newLines;
-  
+
   MVertex *vL = getBeginVertex()->mesh_vertices[0];
   int i = 0;
   for(i = 0; i < (int)lines.size() - 1; i++){
@@ -290,7 +290,7 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
     if (_orientation[i] == 1 ) vR = lines[i]->getVertex(1);
     else vR = lines[i]->getVertex(0);
     int param = i+1;
-    MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this, 
+    MVertex *vNEW = new MEdgeVertex(vR->x(),vR->y(),vR->z(), this,
                                     param, -1., vR->getNum());
     old2new.insert(std::make_pair(vR,vNEW));
     newVertices.push_back(vNEW);
@@ -308,15 +308,15 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
   mesh_vertices = newVertices;
   lines.clear();
   lines = newLines;
-  
+
   //  we need to recreate lines, triangles and tets
   //  that contain those new MEdgeVertices
-  
-   for(std::list<GFace*>::iterator iFace = l_faces.begin(); 
+
+   for(std::list<GFace*>::iterator iFace = l_faces.begin();
        iFace != l_faces.end(); ++iFace){
 
      // for each face find correspondane face2Vertex
-     std::map<GFace*, std::map<MVertex*, MVertex*, 
+     std::map<GFace*, std::map<MVertex*, MVertex*,
        std::less<MVertex*> > >::iterator itmap = face2Vert.find(*iFace);
      if (itmap == face2Vert.end()) {
        face2Vert.insert(std::make_pair(*iFace, old2new));
@@ -332,38 +332,38 @@ void discreteEdge::parametrize(std::map<GFace*, std::map<MVertex*, MVertex*,
      // do the same for regions
      for ( int j = 0; j < (*iFace)->numRegions(); j++){
        GRegion *r = (*iFace)->getRegion(j);
-       std::map<GRegion*,  std::map<MVertex*, MVertex*, 
+       std::map<GRegion*,  std::map<MVertex*, MVertex*,
          std::less<MVertex*> > >::iterator itmap = region2Vert.find(r);
        if (itmap == region2Vert.end()) {
 	 region2Vert.insert(std::make_pair(r, old2new));
        }
        else{
 	 std::map<MVertex*, MVertex*, std::less<MVertex*> > mapVert = itmap->second;
-	 for (std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator it = 
+	 for (std::map<MVertex*, MVertex*, std::less<MVertex*> >::iterator it =
                 old2new.begin(); it != old2new.end(); it++)
 	   mapVert.insert(*it);
 	 itmap->second = mapVert;
        }
      }
- 
+
    }
 }
 
 void discreteEdge::computeNormals () const
 {
   _normals.clear();
-  for (unsigned int i= 0; i < mesh_vertices.size(); i++) 
+  for (unsigned int i= 0; i < mesh_vertices.size(); i++)
     _normals[mesh_vertices[i]] = SVector3(0.0,0.0,0.0);
   _normals[getBeginVertex()->mesh_vertices[0]] = SVector3(0.0,0.0,0.0);
   _normals[getBeginVertex()->mesh_vertices[0]] = SVector3(0.0,0.0,0.0);
   double J[3][3];
 
-  for(std::list<GFace*>::const_iterator iFace = l_faces.begin(); 
+  for(std::list<GFace*>::const_iterator iFace = l_faces.begin();
       iFace != l_faces.end(); ++iFace){
     for (unsigned int i = 0; i < (*iFace)->triangles.size(); ++i){
       MTriangle *t = (*iFace)->triangles[i];
       for (int j = 0; j < 3; j++){
-        std::map<MVertex*, SVector3, std::less<MVertex*> >::iterator 
+        std::map<MVertex*, SVector3, std::less<MVertex*> >::iterator
           itn = _normals.find(t->getVertex(j));
         if (itn != _normals.end()){
           t->getJacobian(0, 0, 0, J);
@@ -371,7 +371,7 @@ void discreteEdge::computeNormals () const
           SVector3 d2(J[1][0], J[1][1], J[1][2]);
           SVector3 n = crossprod(d1, d2);
           n.normalize();
-          _normals[t->getVertex(j)] += n; 
+          _normals[t->getVertex(j)] += n;
         }
       }
     }
@@ -446,7 +446,7 @@ double discreteEdge::curvature(double par) const
   if( !Curvature::valueAlreadyComputed() ) {
     Msg::Warning("Need to compute discrete curvature (in discreteEdge)");
     Curvature::typeOfCurvature type = Curvature::RUSIN; //RUSIN; //RBF
-    curvature.computeCurvature(model(), type); 
+    curvature.computeCurvature(model(), type);
   }
 
   curvature.edgeNodalValues(lines[iEdge],c0, c1, 1);
@@ -475,6 +475,7 @@ Range<double> discreteEdge::parBounds(int i) const
 
 void discreteEdge::createGeometry(){
   if (discrete_lines.empty()){
+
     createTopo();
     // copy the mesh
     for (unsigned int i = 0; i < mesh_vertices.size(); i++){
@@ -493,6 +494,9 @@ void discreteEdge::createGeometry(){
 	discrete_lines.push_back(new MLine(v00,v01));
       }
       else{
+        lines[i]->setVertex(0, v1);
+        lines[i]->setVertex(1, v0);
+        _orientation[i] = 1;
 	discrete_lines.push_back(new MLine(v01,v00));
       }
     }
@@ -501,7 +505,7 @@ void discreteEdge::createGeometry(){
     for (unsigned int i = 1; i < discrete_lines.size(); i++){
       _pars.push_back((double)i);
       MVertex *newv = discrete_lines[i]->getVertex(0);
-      discrete_vertices.push_back(newv);      
+      discrete_vertices.push_back(newv);
     }
     _pars.push_back((double)discrete_lines.size());
   }
@@ -542,6 +546,6 @@ void discreteEdge::mesh(bool verbose){
 void discreteEdge::writeGEO(FILE *fp)
 {
   if (getBeginVertex() && getEndVertex())
-    fprintf(fp, "Discrete Line(%d) = {%d,%d};\n", tag(), 
-            getBeginVertex()->tag(), getEndVertex()->tag());  
+    fprintf(fp, "Discrete Line(%d) = {%d,%d};\n", tag(),
+            getBeginVertex()->tag(), getEndVertex()->tag());
 }
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 2e70b3caf90a122b7d5395790f2a0efb2abb2cb4..1211b07f71f84f23cbb9b55762323197fd6e25a3 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -129,7 +129,6 @@ Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 &param) const
 void discreteFace::secondDer(const SPoint2 &param,
                              SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
 {
-
   if (getCompound()){
     return getCompound()->secondDer(param, dudu, dvdv, dudv);
   }
@@ -140,7 +139,8 @@ void discreteFace::secondDer(const SPoint2 &param,
 }
 
 // FIXME PAB ----> really create an atlas !!!!!!!!!!!!!!!
-void discreteFace::createGeometry() {
+void discreteFace::createGeometry()
+{
 #if defined(HAVE_ANN) && defined(HAVE_SOLVER)
   if (!_atlas.empty())return;
   // parametrization is done here !!!
@@ -150,20 +150,24 @@ void discreteFace::createGeometry() {
 #endif
 }
 
-void discreteFace::gatherMeshes() {
+void discreteFace::gatherMeshes()
+{
 #if defined(HAVE_ANN) && defined(HAVE_SOLVER)
   for (unsigned int i=0;i<triangles.size();i++)delete triangles[i];
   for (unsigned int i=0;i<mesh_vertices.size();i++)delete mesh_vertices[i];
   triangles.clear();
   mesh_vertices.clear();
   for (unsigned int i=0;i<_atlas.size();i++){
-    triangles.insert(triangles.begin(), _atlas[i]->triangles.begin(), _atlas[i]->triangles.end());
-    mesh_vertices.insert(mesh_vertices.begin(), _atlas[i]->mesh_vertices.begin(), _atlas[i]->mesh_vertices.end());
+    triangles.insert(triangles.begin(), _atlas[i]->triangles.begin(),
+                     _atlas[i]->triangles.end());
+    mesh_vertices.insert(mesh_vertices.begin(), _atlas[i]->mesh_vertices.begin(),
+                         _atlas[i]->mesh_vertices.end());
   }
 #endif
 }
 
-void discreteFace::mesh(bool verbose) {
+void discreteFace::mesh(bool verbose)
+{
 #if defined(HAVE_ANN) && defined(HAVE_SOLVER) && defined(HAVE_MESH)
   if (!CTX::instance()->meshDiscrete) return;
   for (unsigned int i=0;i<_atlas.size();i++)
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index e1986c36a0793c80001bbc89a6b0c02218c5fdf5..42a7388a5544098eddbefcd289b873c0b48f5aeb 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -399,15 +399,19 @@ int Mesh2DWithBoundaryLayers(GModel *m)
 
   if(sourceEdges.empty() && sourceFaces.empty()) return 0;
 
-  // from Trevor Strickler -- Just in case ReplaceDuplicates() erases the ExtrudeParams::mesh.scaleLast
-  // flag, should check all bounding regions of this curve to see if scaleLast is set.
-  // if so, reset it in the extrudeParams (maybe this could be done in the TreeUtils....
-  // but I do not want to change the code too much and create a bug.
-  // The developers should decide that.
-  if( ExtrudeParams::calcLayerScaleFactor[0]  || ExtrudeParams::calcLayerScaleFactor[1] ){
-    unsigned int num_changed = FixErasedExtrScaleFlags(m, faceSkipScaleCalc, edgeSkipScaleCalc);
-    if( num_changed )
-      Msg::Warning("%d entities were changed from ScaleLast = false to ScaleLast = true", num_changed);
+  // from Trevor Strickler -- Just in case ReplaceDuplicates() erases the
+  // ExtrudeParams::mesh.scaleLast flag, should check all bounding regions of
+  // this curve to see if scaleLast is set.  if so, reset it in the
+  // extrudeParams (maybe this could be done in the TreeUtils....  but I do not
+  // want to change the code too much and create a bug.  The developers should
+  // decide that.
+  if(ExtrudeParams::calcLayerScaleFactor[0] ||
+     ExtrudeParams::calcLayerScaleFactor[1]){
+    unsigned int num_changed = FixErasedExtrScaleFlags(m, faceSkipScaleCalc,
+                                                       edgeSkipScaleCalc);
+    if(num_changed)
+      Msg::Warning("%d entities were changed from ScaleLast = false to ScaleLast = true",
+                   num_changed);
   }
   // compute mesh dependencies in source faces (so we can e.g. create
   // a boundary layer on an extruded mesh)
@@ -430,10 +434,13 @@ int Mesh2DWithBoundaryLayers(GModel *m)
       otherFaces.insert(*it);
 
   // mesh source surfaces (and dependencies)
-  for(int i = 0; i < 10; i++) // FIXME: should check PENDING status instead!
-    std::for_each(sourceFacesDependencies.begin(), sourceFacesDependencies.end(),
-                  meshGFace());
-  std::for_each(sourceFaces.begin(), sourceFaces.end(), meshGFace());
+  for(int i = 0; i < 10; i++){ // FIXME: should check PENDING status instead!
+    for(std::set<GFace*>::iterator it = sourceFacesDependencies.begin();
+        it != sourceFacesDependencies.end(); it++)
+      (*it)->mesh(false);
+  }
+  for(std::set<GFace*>::iterator it = sourceFaces.begin(); it != sourceFaces.end(); it++)
+    (*it)->mesh(false);
 
   // make sure the source surfaces for the boundary layers are
   // oriented correctly (normally we do this only after the 3D mesh is
@@ -479,7 +486,8 @@ int Mesh2DWithBoundaryLayers(GModel *m)
   // remesh non-source edges (since they might have been modified by
   // the change in boundary layer points)
   std::for_each(otherFaces.begin(), otherFaces.end(), deMeshGFace());
-  std::for_each(otherEdges.begin(), otherEdges.end(), meshGEdge());
+  for(std::set<GEdge*>::iterator it = otherEdges.begin(); it != otherEdges.end(); it++)
+    (*it)->mesh(false);
 
   // mesh the curves bounding the boundary layers by extrusion using
   // the smooth normal field
@@ -501,7 +509,8 @@ int Mesh2DWithBoundaryLayers(GModel *m)
   }
 
   // mesh non-source surfaces
-  std::for_each(otherFaces.begin(), otherFaces.end(), meshGFace());
+  for(std::set<GFace*>::iterator it = otherFaces.begin(); it != otherFaces.end(); it++)
+    (*it)->mesh(false);
 
   // mesh the surfaces bounding the boundary layers by extrusion using
   // the smooth normal field
@@ -517,4 +526,3 @@ int Mesh2DWithBoundaryLayers(GModel *m)
 
   return 1;
 }
-
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 54f8e06a976d093fa93e3f17e572686846a8dc7a..1f309a4cecb422ee36da0cd12274a16ce8ee2554 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -222,8 +222,9 @@ static void Mesh0D(GModel *m)
     GVertex *gv = *it;
     if (gv->meshMaster() != gv){
       if (gv->correspondingVertices.empty()){
-        GVertex *master = dynamic_cast<GVertex*> (gv->meshMaster());
-        if(master)gv->correspondingVertices[gv->mesh_vertices[0]] = (master->mesh_vertices[0]);
+        GVertex *master = dynamic_cast<GVertex*>(gv->meshMaster());
+        if(master)
+          gv->correspondingVertices[gv->mesh_vertices[0]] = master->mesh_vertices[0];
       }
     }
   }
@@ -295,7 +296,8 @@ static void PrintMesh2dStatistics(GModel *m)
     if((*it)->geomType() != GEntity::DiscreteSurface){
       worst = std::min((*it)->meshStatistics.worst_element_shape, worst);
       best = std::max((*it)->meshStatistics.best_element_shape, best);
-      avg += (*it)->meshStatistics.average_element_shape * (*it)->meshStatistics.nbTriangle;
+      avg += (*it)->meshStatistics.average_element_shape *
+        (*it)->meshStatistics.nbTriangle;
       e_avg += (*it)->meshStatistics.efficiency_index;
       e_long = std::max((*it)->meshStatistics.longest_edge_length, e_long);
       e_short = std::min((*it)->meshStatistics.smallest_edge_length, e_short);
@@ -309,9 +311,11 @@ static void PrintMesh2dStatistics(GModel *m)
     }
   }
 
-  Msg::Info("*** Efficiency index for surface mesh tau=%g ", 100*exp(e_avg/(double)nTotE));
+  Msg::Info("*** Efficiency index for surface mesh tau=%g ",
+            100*exp(e_avg/(double)nTotE));
 
-  fprintf(statreport,"\t%16s\t%d\t\t%d\t\t", m->getName().c_str(), numFaces, nUnmeshed);
+  fprintf(statreport,"\t%16s\t%d\t\t%d\t\t", m->getName().c_str(), numFaces,
+          nUnmeshed);
   fprintf(statreport,"%d\t\t%8.7f\t%8.7f\t%8.7f\t%d\t\t%8.7f\t",
           nTotT, avg / (double)nTotT, best, worst, nTotGoodQuality,
           (double)nTotGoodQuality / nTotT);
@@ -330,9 +334,9 @@ static void Mesh2D(GModel *m)
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
     (*it)->meshStatistics.status = GFace::PENDING;
 
-  // boundary layers are special: their generation (including vertices
-  // and curve meshes) is global as it depends on a smooth normal
-  // field generated from the surface mesh of the source surfaces
+  // boundary layers are special: their generation (including vertices and curve
+  // meshes) is global as it depends on a smooth normal field generated from the
+  // surface mesh of the source surfaces
   if(!Mesh2DWithBoundaryLayers(m)){
     std::set<GFace*, GEntityLessThan> cf, f;
     for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
@@ -354,7 +358,7 @@ static void Mesh2D(GModel *m)
       for(size_t K = 0 ; K < temp.size() ; K++){
         if (temp[K]->meshStatistics.status == GFace::PENDING){
           backgroundMesh::current()->unset();
-	  //          meshGFace mesher(true);
+	  // meshGFace mesher(true);
           temp[K]->mesh(true);
 #if defined(HAVE_BFGS)
           if(CTX::instance()->mesh.optimizeLloyd){
@@ -480,9 +484,8 @@ void fillv_(std::multimap<MVertex*, MElement*> &vertexToElement,
   }
 }
 
-
-int LaplaceSmoothing (GRegion *gr) {
-
+int LaplaceSmoothing (GRegion *gr)
+{
   std::multimap<MVertex*, MElement*> vertexToElement;
   fillv_(vertexToElement, (gr)->tetrahedra.begin(), (gr)->tetrahedra.end());
   fillv_(vertexToElement, (gr)->hexahedra.begin(),  (gr)->hexahedra.end());
@@ -510,13 +513,13 @@ int LaplaceSmoothing (GRegion *gr) {
     double minQual2 = 1.e22;
     for (it = it_low; it != it_up; ++it) {
       minQual2 = std::min(minQual2,it->second->minSICNShapeMeasure());
-      if (minQual2 < minQual){	
+      if (minQual2 < minQual){
 	v->setXYZ (xold,yold,zold);
 	break;
       }
     }
     if (minQual < minQual2) N++;
-  }  
+  }
   return N;
 }
 
@@ -548,8 +551,8 @@ void buildUniqueFaces (GRegion *gr, std::set<MFace,Less_Face> &bnd)
   }
 }
 
-bool MakeMeshConformal   (GModel *gm, int howto) {
-
+bool MakeMeshConformal(GModel *gm, int howto)
+{
   fs_cont search;
   buildFaceSearchStructure(gm, search);
   std::set<MFace,Less_Face> bnd;
@@ -558,7 +561,7 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
     buildUniqueFaces (gr,bnd);
   }
   // bnd2 contains non conforming faces
-  
+
    std::set<MFace,Less_Face> bnd2;
   for (std::set<MFace,Less_Face>::iterator itf = bnd.begin(); itf != bnd.end(); ++itf){
     GFace *gfound = findInFaceSearchStructure (*itf,search);
@@ -573,20 +576,24 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
   std::set<MFace,Less_Face> ncf;
   for (std::set<MFace,Less_Face>::iterator itf = bnd2.begin(); itf != bnd2.end(); ++itf){
     const MFace &f = *itf;
-    if (f.getNumVertices () == 4){ // quad face
-      std::set<MFace,Less_Face>::iterator it1 = bnd2.find (MFace(f.getVertex(0),f.getVertex(1),f.getVertex(2)));
-      std::set<MFace,Less_Face>::iterator it2 = bnd2.find (MFace(f.getVertex(2),f.getVertex(3),f.getVertex(0)));
+    if (f.getNumVertices() == 4){ // quad face
+      std::set<MFace,Less_Face>::iterator it1 =
+        bnd2.find (MFace(f.getVertex(0),f.getVertex(1),f.getVertex(2)));
+      std::set<MFace,Less_Face>::iterator it2 =
+        bnd2.find (MFace(f.getVertex(2),f.getVertex(3),f.getVertex(0)));
       if (it1 != bnd2.end() && it2 != bnd2.end()){
-	ncf.insert(MFace (f.getVertex(1),f.getVertex(2), f.getVertex(3),f.getVertex(0) )); 
+	ncf.insert(MFace(f.getVertex(1), f.getVertex(2),
+                         f.getVertex(3),f.getVertex(0)));
       }
       else {
 	it1 = bnd2.find (MFace(f.getVertex(0),f.getVertex(1),f.getVertex(3)));
 	it2 = bnd2.find (MFace(f.getVertex(3),f.getVertex(1),f.getVertex(2)));
 	if (it1 != bnd2.end() && it2 != bnd2.end()){
-	  ncf.insert(MFace (f.getVertex(0),f.getVertex(1), f.getVertex(2),f.getVertex(3) )); 
+	  ncf.insert(MFace(f.getVertex(0), f.getVertex(1),
+                           f.getVertex(2), f.getVertex(3)));
 	}
 	else {
-	  Msg::Error ("MakeMeshConformal: wrong mesh topology");
+	  Msg::Error("MakeMeshConformal: wrong mesh topology");
 	  return false;
 	}
       }
@@ -608,7 +615,7 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
 	}
 	else {
 	  faces.push_back(MFace(it->getVertex(0),it->getVertex(1),it->getVertex(3)));
-	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));	
+	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));
 	}
       }
       // HEX IS ONLY SURROUNED BY COMPATIBLE ELEMENTS
@@ -622,13 +629,15 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
 	for (unsigned int j=0;j<faces.size();j++){
 	  MFace &f = faces[j];
 	  if (f.getNumVertices() == 4){
-	    gr->pyramids.push_back(new MPyramid (f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), newv));  
+	    gr->pyramids.push_back(new MPyramid(f.getVertex(0), f.getVertex(1),
+                                                f.getVertex(2), f.getVertex(3), newv));
 	  }
 	  else {
-	    gr->tetrahedra.push_back(new MTetrahedron (f.getVertex(0), f.getVertex(1), f.getVertex(2), newv));  
+	    gr->tetrahedra.push_back(new MTetrahedron(f.getVertex(0), f.getVertex(1),
+                                                      f.getVertex(2), newv));
 	  }
 	}
-      }    
+      }
     }
     gr->hexahedra = remainingHexes;
     remainingHexes.clear();
@@ -644,7 +653,7 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
 	}
 	else {
 	  faces.push_back(MFace(it->getVertex(0),it->getVertex(1),it->getVertex(3)));
-	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));	
+	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));
 	}
       }
       // HEX IS ONLY SURROUNED BY COMPATIBLE ELEMENTS
@@ -658,13 +667,15 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
 	for (unsigned int j=0;j<faces.size();j++){
 	  MFace &f = faces[j];
 	  if (f.getNumVertices() == 4){
-	    gr->pyramids.push_back(new MPyramid (f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), newv));  
+	    gr->pyramids.push_back(new MPyramid(f.getVertex(0), f.getVertex(1),
+                                                f.getVertex(2), f.getVertex(3), newv));
 	  }
 	  else {
-	    gr->tetrahedra.push_back(new MTetrahedron (f.getVertex(0), f.getVertex(1), f.getVertex(2), newv));  
+	    gr->tetrahedra.push_back(new MTetrahedron(f.getVertex(0), f.getVertex(1),
+                                                      f.getVertex(2), newv));
 	  }
 	}
-      }    
+      }
     }
     gr->prisms = remainingPrisms;
   }
@@ -672,7 +683,8 @@ bool MakeMeshConformal   (GModel *gm, int howto) {
   return true;
 }
 
-void TestConformity   (GModel *gm) {
+void TestConformity(GModel *gm)
+{
   fs_cont search;
   buildFaceSearchStructure(gm, search);
   int count = 0;
@@ -691,7 +703,7 @@ void TestConformity   (GModel *gm) {
       }
     }
     printf("vol(%d) = %12.5E\n",gr->tag(),vol);
-    
+
     for (std::set<MFace,Less_Face>::iterator itf = bnd.begin(); itf != bnd.end(); ++itf){
       GFace *gfound = findInFaceSearchStructure (*itf,search);
       if (!gfound){
@@ -781,10 +793,12 @@ static void Mesh3D(GModel *m)
           sup.execute(gr);
         }
         PostOp post;
-        post.execute(gr,CTX::instance()->mesh.recombine3DLevel, CTX::instance()->mesh.recombine3DConformity); //0: no pyramid, 1: single-step, 2: two-steps (conforming), true: fill non-conformities with trihedra
-  
-	//	while(LaplaceSmoothing (gr)){
-	//	}
+        post.execute(gr,CTX::instance()->mesh.recombine3DLevel,
+                     CTX::instance()->mesh.recombine3DConformity);
+        // 0: no pyramid, 1: single-step, 2: two-steps (conforming),
+        // true: fill non-conformities with trihedra
+	// while(LaplaceSmoothing (gr)){
+	// }
         nb_elements_recombination += post.get_nb_elements();
         nb_hexa_recombination += post.get_nb_hexahedra();
         vol_element_recombination += post.get_vol_elements();
@@ -799,7 +813,7 @@ static void Mesh3D(GModel *m)
       }
     }
   }
-  
+
   if(CTX::instance()->mesh.recombine3DAll){
     Msg::Info("RECOMBINATION timing:");
     Msg::Info(" --- CUMULATIVE TIME RECOMBINATION : %g s.", time_recombination);
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 49c1d120f810197f53fb0a39ff030c0c64f6654c..59f2a8894ba02fedfd4476fba52914b6ce32a801 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -305,7 +305,8 @@ void copyMesh(GEdge *from, GEdge *to, int direction)
 
 void deMeshGEdge::operator() (GEdge *ge)
 {
-  if(ge->geomType() == GEntity::DiscreteCurve && !CTX::instance()->meshDiscrete)return;
+  if(ge->geomType() == GEntity::DiscreteCurve && !CTX::instance()->meshDiscrete)
+    return;
   ge->deleteMesh();
   ge->meshStatistics.status = GEdge::PENDING;
   ge->correspondingVertices.clear();
@@ -346,7 +347,7 @@ static void filterPoints (GEdge*ge) {
       CTX::instance()->mesh.flexibleTransfinite) &&
      CTX::instance()->mesh.algoRecombine != 0){
     if(CTX::instance()->mesh.recombineAll){
-      forceOdd = true;      
+      forceOdd = true;
     }
   }
 
@@ -368,7 +369,7 @@ static void filterPoints (GEdge*ge) {
     if (d < lc * .3) {
       lengths.push_back(std::make_pair(lc/d,v));
     }
-    else 
+    else
       v0=v;
   }
   std::sort(lengths.begin(),lengths.end());
@@ -376,7 +377,7 @@ static void filterPoints (GEdge*ge) {
   if (forceOdd) {
     while (last %2 != 0)last--;
   }
-  /*    
+  /*
 	if (CTX::instance()->mesh.algoRecombine == 2){
 	if (last < 4)last = 0;
 	while (last %4 != 0)last--;