diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 96bf32dde147274f38a27a0160daab485abaf581..fb88fe0c134cdb9d2495ed29abce82c24d8253a3 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -485,13 +485,14 @@ bool GFaceCompound::parametrize() const
   //-----------------
   if (_mapping == HARMONIC){
     Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag()); 
-    fillNeumannBCS();
+    //fillNeumannBCS();
     parametrize(ITERU,HARMONIC); 
     parametrize(ITERV,HARMONIC);
   }
   //Multiscale Laplace parametrization
   //-----------------
   else if (_mapping == MULTISCALE){
+    printf("multiscale exit \n");
     std::vector<MElement*> _elements;
     for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
       for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 346edff76f3bb98bea5af64f71b0eb1f7d516276..6c81e9ef1ab2a2acee0a605480489dee4b6db14c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1081,6 +1081,7 @@ void GModel::createTopologyFromMesh()
 {
 
   Msg::Info("Creating topology from mesh");
+  double t1 = Cpu();
 
   // for each discreteRegion, create topology
   std::vector<discreteRegion*> discRegions;
@@ -1102,8 +1103,12 @@ void GModel::createTopologyFromMesh()
   
   //EMI TODO
   //check for closed faces
-  
+ 
   createTopologyFromFaces(discFaces);
+
+  double t2 = Cpu();
+  Msg::Info("Creating topology from mesh done (%g s)", t2-t1);
+ 
   
   //create old format (necessaray for boundary layers)
   exportDiscreteGEOInternals();
@@ -1131,7 +1136,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
   // create reverse map, for each face find set of MEdges that are
   // candidate for new discrete Edges
   std::map<int, std::vector<int> > face2Edges;
-  
+
   while (!map_edges.empty()){
 
     std::vector<MEdge> myEdges;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 644564efd0844fcf455fe953c1b3f95c1a7fcdb0..08b7c47ca3392a0fbd047f15929818687272b96e 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -11,6 +11,8 @@
 #include "MEdge.h"
 #include "Geo.h"
 #include "GFaceCompound.h"
+#include "Context.h"
+#include "Os.h"
 
 discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
 {
@@ -22,22 +24,19 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
 void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_edges)
 {
 
-  // find the boundary edges
-  std::list<MEdge> bound_edges;
+  std::set<MEdge, Less_Edge> bound_edges;
   for (unsigned int iFace = 0; iFace  < getNumMeshElements() ; iFace++) {
     MElement *e = getMeshElement(iFace);
     for (int iEdge = 0; iEdge < e->getNumEdges(); iEdge++) {
       MEdge tmp_edge =  e->getEdge(iEdge);
-      if (std::find(bound_edges.begin(), bound_edges.end(), tmp_edge) == 
-          bound_edges.end())
-        bound_edges.push_back(tmp_edge);
-      else
-        bound_edges.erase(std::find(bound_edges.begin(), bound_edges.end(), tmp_edge));
+      std::set<MEdge, Less_Edge >::iterator itset = bound_edges.find(tmp_edge);
+      if (itset == bound_edges.end()) 	bound_edges.insert(tmp_edge);
+      else bound_edges.erase(itset);
     }
   }
 
   // for the boundary edges, associate the tag of the current discrete face
-  for (std::list<MEdge>::iterator itv = bound_edges.begin(); 
+  for (std::set<MEdge, Less_Edge>::iterator itv = bound_edges.begin(); 
        itv != bound_edges.end(); ++itv){
     std::map<MEdge, std::vector<int>, Less_Edge >::iterator itmap = map_edges.find(*itv);
     if (itmap == map_edges.end()){
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 6d173753154de23a995cd6854e90235f67e4d7bd..dad3b0777edb8121e7ded7c0d92ae2be9909e00f 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -359,7 +359,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
 		nextLevel->recur,nextLevel->region, genus, AR, nbParts);  
       partition(*nextLevel, nbParts, MULTILEVEL);
     }
-    else if (genus == 0  &&  AR > 3 || genus == 0  &&  NB > 1){
+    else if (genus == 0  &&  AR > 3  || genus == 0  &&  NB > 1){
       int nbParts = 2;
       Msg::Info("Mesh partition: level (%d-%d)  is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
  		nextLevel->recur,nextLevel->region, AR, NB, nbParts);  
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 370bb2508bca56f2a2b2bebc99b47a9f704e5ace..2ff9fd997803ee835381fec8973e5140a8be6e2d 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -1,3 +1,4 @@
+
 #include "dgConservationLaw.h"
 #include "dgConservationLawWaveEquation.h"
 #include "function.h"
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index aef8830e523f9b44cf43212e1ef4393014f90763..a3e871fc9a99ed173b5b3808572fe5b2994a48a4 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1078,8 +1078,10 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
   recur_cut_ (1.0, M_PI, 0.0, root,left,right);
   connected_left_right(left, right);
 
-  if ( elements.size() != left.size()+right.size()) 
+  if ( elements.size() != left.size()+right.size()) {
     Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)",  elements.size(), left.size()+right.size());
+    exit(1);
+  }
 
   elements.clear();
   elements.insert(elements.end(),left.begin(),left.end());