diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 1031141c8ae4f5d5b5969857485e1dd9651e1fc5..3092d3448a4d73053f312471e28ce3fd1143da49 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1077,7 +1077,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   }
 
   // build a set with all points of the boundaries
-  std::set<MVertex*> all_vertices, boundary;
+  std::set<MVertex*, MVertexLessThanNum> all_vertices, boundary;
   std::list<GEdge*>::iterator ite = edges.begin();
   while(ite != edges.end()){
     if((*ite)->isSeam(gf)) return false;
@@ -1152,7 +1152,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   if(all_vertices.size() == 3){
     MVertex *vv[3];
     int i = 0;
-    for(std::set<MVertex*>::iterator it = all_vertices.begin();
+    for(std::set<MVertex*, MVertexLessThanNum>::iterator it = all_vertices.begin();
 	it != all_vertices.end(); it++){
       vv[i++] = *it;
     }
@@ -1170,7 +1170,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   std::vector<BDS_Point*> points(all_vertices.size());
   SBoundingBox3d bbox;
   int count = 0;
-  for(std::set<MVertex*>::iterator it = all_vertices.begin();
+  for(std::set<MVertex*, MVertexLessThanNum>::iterator it = all_vertices.begin();
       it != all_vertices.end(); it++){
     MVertex *here = *it;
     GEntity *ge = here->onWhat();
@@ -1266,7 +1266,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
   {
     std::vector<MVertex*> v;
     std::vector<MTriangle*> result;
-    v.insert(v.end(),all_vertices.begin(),all_vertices.end());
+    v.insert(v.end(), all_vertices.begin(), all_vertices.end());
 
     std::map<MVertex*,SPoint3> pos;
     for(unsigned int i = 0; i < v.size(); i++) {
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 8b1754088bf8335d999b84938f4e30e03165f315..e20182df3a695de2a04707cdc4265c45c492dde1 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -479,6 +479,14 @@ static void midpointsphere(GFace *gf, double u1, double v1, double u2, double v2
 }
 */
 
+bool edges_sort(std::pair<double, BDS_Edge*> a, std::pair<double, BDS_Edge*> b)
+{
+  if (a.first == b.first)
+    return ((*a.second) < (*b.second));
+  else
+    return (a.first < b.first);
+}
+
 void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
 {
   std::list<BDS_Edge*>::iterator it = m.edges.begin();
@@ -494,7 +502,7 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
     ++it;
   }
 
-  std::sort(edges.begin(), edges.end());
+  std::sort(edges.begin(), edges.end(), edges_sort);
 
   for (unsigned int i = 0; i < edges.size(); ++i){
     BDS_Edge *e = edges[i].second;
@@ -547,7 +555,7 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
     ++it;
   }
 
-  std::sort(edges.begin(), edges.end());
+  std::sort(edges.begin(), edges.end(), edges_sort);
 
   for (unsigned int i = 0; i < edges.size(); i++){
     BDS_Edge *e = edges[i].second;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 647f464037bc75e7a8f1fc0f704eb948b5748cbc..600e343ffa980a50b6d3b54b96e43a996a500a81 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1395,7 +1395,8 @@ namespace nglib {
 }
 using namespace nglib;
 
-static void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
+static void getAllBoundingVertices(GRegion *gr,
+                                   std::set<MVertex*, MVertexLessThanNum> &allBoundingVertices)
 {
   std::list<GFace*> faces = gr->faces();
   std::list<GFace*>::iterator it = faces.begin();
@@ -1418,10 +1419,10 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
   Ng_Init();
   Ng_Mesh *ngmesh = Ng_NewMesh();
 
-  std::set<MVertex*> allBoundingVertices;
+  std::set<MVertex*, MVertexLessThanNum> allBoundingVertices;
   getAllBoundingVertices(gr, allBoundingVertices);
 
-  std::set<MVertex*>::iterator itv = allBoundingVertices.begin();
+  std::set<MVertex*, MVertexLessThanNum>::iterator itv = allBoundingVertices.begin();
   int I = 1;
   while(itv != allBoundingVertices.end()){
     double tmp[3];
@@ -1683,7 +1684,7 @@ void meshGRegion::operator() (GRegion *gr)
 
   // replace discreteFaces by their compounds
   {
-    std::set<GFace*> mySet;
+    std::set<GFace*, GEntityLessThan> mySet;
     std::list<GFace*>::iterator it = faces.begin();
     while(it != faces.end()){
       if((*it)->getCompound())