diff --git a/Geo/MElement.h b/Geo/MElement.h
index f44d75d75a5a744bfae3e28f3c9d3f0689ce6e65..f63bf2769159bca665c5dba0f234452169d2f62c 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -141,21 +141,24 @@ class MElement
 
   // get an edge representation for drawing
   virtual int getNumEdgesRep(bool curved) = 0;
-  virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) = 0;
+  virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z,
+                          SVector3 *n) = 0;
 
   // get the faces
   virtual int getNumFaces() = 0;
   virtual MFace getFace(int num) = 0;
 
   // give an MFace as input and get its local number, sign and rotation
-  virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign, int &rot) const
+  virtual void getFaceInfo(const MFace & face, int &ithFace, int &sign,
+                           int &rot) const
   {
     Msg::Error("Face information not available for this element");
   }
 
   // get a face representation for drawing
   virtual int getNumFacesRep(bool curved) = 0;
-  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n) = 0;
+  virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z,
+                          SVector3 *n) = 0;
 
   // get all the vertices on a edge or a face
   virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
@@ -173,6 +176,7 @@ class MElement
   virtual int getNumChildren() const { return 0; }
   virtual MElement *getChild(int i) const { return NULL; }
   virtual bool ownsParent() const { return false; }
+
   // get base element in case of MSubElement
   virtual const MElement *getBaseElement() const { return this; }
   virtual MElement *getBaseElement() { return this; }
@@ -188,7 +192,7 @@ class MElement
   virtual double maxEdge();
   virtual double minEdge();
 
-  // Max. distance between curved and straight element among all high-order nodes
+  // max. distance between curved and straight element among all high-order nodes
   double maxDistToStraight() const;
 
   // get the quality measures
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 52ba59c1d58a65abdabda46d0a7038ab12a07c19..2e949a0be9b0c6acd37c187fdd06a1d7e39c253d 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -14,24 +14,6 @@
 #include "GmshMessage.h"
 #include "StringUtils.h"
 
-double MVertexLessThanLexicographic::tolerance = 1.e-6;
-
-bool MVertexLessThanLexicographic::operator()(const MVertex *v1, const MVertex *v2) const
-{
-  if(v1->x() - v2->x() >  tolerance) return true;
-  if(v1->x() - v2->x() < -tolerance) return false;
-  if(v1->y() - v2->y() >  tolerance) return true;
-  if(v1->y() - v2->y() < -tolerance) return false;
-  if(v1->z() - v2->z() >  tolerance) return true;
-  return false;
-}
-
-bool MVertexLessThanNum::operator()(const MVertex *v1, const MVertex *v2) const
-{
-  if(v1->getNum() < v2->getNum()) return true;
-  return false;
-}
-
 double angle3Vertices(const MVertex *p1, const MVertex *p2, const MVertex *p3)
 {
   SVector3 a(p1->x() - p2->x(), p1->y() - p2->y(), p1->z() - p2->z());
@@ -349,13 +331,24 @@ void MVertex::writeSU2(FILE *fp, int dim, double scalingFactor)
             z() * scalingFactor, _index - 1);
 }
 
-std::set<MVertex*, MVertexLessThanLexicographic>::iterator
-MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+bool MVertexLessThanNum::operator()(const MVertex *v1, const MVertex *v2) const
+{
+  if(v1->getNum() < v2->getNum()) return true;
+  return false;
+}
+
+double MVertexLessThanLexicographic::tolerance = 1.e-6;
+
+bool MVertexLessThanLexicographic::operator()(const MVertex *v1, const MVertex *v2) const
 {
-  for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
-      it != pos.end(); ++it)
-    if(distance(*it) < MVertexLessThanLexicographic::tolerance) return it;
-  return pos.end();
+  // you should not use this unless you know what you are doing; to create
+  // unique vertices, use MVertexRTree
+  if(v1->x() - v2->x() >  tolerance) return true;
+  if(v1->x() - v2->x() < -tolerance) return false;
+  if(v1->y() - v2->y() >  tolerance) return true;
+  if(v1->y() - v2->y() < -tolerance) return false;
+  if(v1->z() - v2->z() >  tolerance) return true;
+  return false;
 }
 
 static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params)
@@ -375,7 +368,6 @@ static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> &params
     for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++){
       if((*it)->isSeam(gf)) {
         Range<double> range = (*it)->parBounds(0);
-	//	printf("%d %d %d\n",gv->tag(),(*it)->getBeginVertex()->tag(),(*it)->getEndVertex()->tag());
         if (gv == (*it)->getBeginVertex()){
           params.push_back((*it)->reparamOnFace(gf, range.low(),-1));
           params.push_back((*it)->reparamOnFace(gf, range.low(), 1));
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index c4049fe74bd1a97b0a7f08787cfbc9ed47ea3dc0..037db31d13b932d42fbd38fff2efa911237c45ae 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -18,17 +18,6 @@ class GEdge;
 class GFace;
 class MVertex;
 
-class MVertexLessThanLexicographic{
- public:
-  static double tolerance;
-  bool operator()(const MVertex *v1, const MVertex *v2) const;
-};
-
-class MVertexLessThanNum{
- public:
-  bool operator()(const MVertex *v1, const MVertex *v2) const;
-};
-
 // A mesh vertex.
 class MVertex{
  protected:
@@ -99,10 +88,6 @@ class MVertex{
     return sqrt(dx * dx + dy * dy + dz * dz);
   }
 
-  // linear coordinate search for the vertex in a set
-  std::set<MVertex*, MVertexLessThanLexicographic>::iterator
-  linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos);
-
   // IO routines
   void writeMSH(FILE *fp, bool binary=false, bool saveParametric=false,
                 double scalingFactor=1.0);
@@ -159,6 +144,17 @@ class MFaceVertex : public MVertex{
   }
 };
 
+class MVertexLessThanLexicographic{
+ public:
+  static double tolerance;
+  bool operator()(const MVertex *v1, const MVertex *v2) const;
+};
+
+class MVertexLessThanNum{
+ public:
+  bool operator()(const MVertex *v1, const MVertex *v2) const;
+};
+
 bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
                            SPoint2 &param1, SPoint2 &param2);
 bool reparamMeshVertexOnFace(MVertex *v, const GFace *gf, SPoint2 &param,
@@ -167,7 +163,8 @@ bool reparamMeshVertexOnEdge(MVertex *v, const GEdge *ge, double &param);
 
 double angle3Vertices(const MVertex *p1, const MVertex *p2, const MVertex *p3);
 
-inline double distance (MVertex *v1, MVertex *v2){
+inline double distance (MVertex *v1, MVertex *v2)
+{
   const double dx = v1->x() - v2->x();
   const double dy = v1->y() - v2->y();
   const double dz = v1->z() - v2->z();
diff --git a/Geo/MVertexRTree.h b/Geo/MVertexRTree.h
index bfc9aae9090a01e8c99b031e1d0eff8f06adc480..1ec83a687d92eadb51a727c8fe4a4916b60fdefe 100644
--- a/Geo/MVertexRTree.h
+++ b/Geo/MVertexRTree.h
@@ -21,7 +21,7 @@ class MVertexRTree{
   {
     MVertex **out = static_cast<MVertex**>(ctx);
     *out = v;
-    return false;
+    return false; // we're done searching
   }
  public:
   MVertexRTree(double tolerance = 1.e-8)
@@ -48,14 +48,14 @@ class MVertexRTree{
                    "mesh with tolerance %g", v->getNum(),
                    v->x(), v->y(), v->z(), _tol);
     }
-    return 1;
+    return 1; // one vertex not inserted
   }
   int insert(std::vector<MVertex*> &v, bool warnIfExists=false)
   {
     int num = 0;
     for(unsigned int i = 0; i < v.size(); i++)
       num += insert(v[i], warnIfExists);
-    return num;
+    return num; // number of vertices not inserted
   }
   MVertex *find(double x, double y, double z)
   {