diff --git a/Geo/MFace.h b/Geo/MFace.h
index f6cec8c283424d8478c29635a11b33cc34ee56d3..013e969553a7faa422a7a08ec0cb07ded97c77c8 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -20,11 +20,13 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+#include <functional>
 #include <vector>
 #include "MVertex.h"
 #include "SVector3.h"
 #include "Numeric.h"
 #include "Context.h"
+#include "Hash.h"
 
 extern Context_T CTX;
 
@@ -32,6 +34,8 @@ extern Context_T CTX;
 class MFace {
  private:
   MVertex *_v[4];
+  char _si[4];                          // sorted indices
+
  public:
   MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3=0) 
   {
@@ -49,16 +53,60 @@ class MFace {
     else{
       _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
     }
+    // This is simply an unrolled insertion sort (hopefully fast).  Note that if
+    // _v[3] == 0, _v[3] is not sorted.
+    if(_v[1] < _v[0]) {
+      _si[0] = 1;
+      _si[1] = 0;
+    }
+    else {
+      _si[0] = 0;
+      _si[1] = 1;
+    }
+    if(_v[2] < _v[int(_si[1])]) {
+      _si[2] = _si[1];
+      if(_v[2] < _v[int(_si[0])]) {
+        _si[1] = _si[0];
+        _si[0] = 2;
+      }
+      else
+        _si[1] = 2;
+    }
+    else
+      _si[2] = 2;
+    if( _v[3] && _v[3] < _v[int(_si[2])]) {
+      _si[3] = _si[2];
+      if(_v[3] < _v[int(_si[1])]) {
+        _si[2] = _si[1];
+        if(_v[3] < _v[int(_si[0])]) {
+          _si[1] = _si[0];
+          _si[0] = 3;
+        }
+        else
+          _si[1] = 3;
+      }
+      else
+        _si[2] = 3;
+    }
+    else
+      _si[3] = 3;
   }
   inline int getNumVertices() const { return _v[3] ? 4 : 3; }
-  inline MVertex *getVertex(int i) const { return _v[i]; }
-  void getOrderedVertices(std::vector<MVertex*> &verts)
+  inline MVertex *getVertex(const int i) const { return _v[i]; }
+  inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
+  void getOrderedVertices(std::vector<MVertex*> &verts) const
   {
     for(int i = 0; i < getNumVertices(); i++)
-      verts.push_back(getVertex(i));
-    std::sort(verts.begin(), verts.end());
+      verts.push_back(getSortedVertex(i));
+  }
+  void getOrderedVertices(const MVertex **const verts) const
+  {
+    verts[0] = getSortedVertex(0);
+    verts[1] = getSortedVertex(1);
+    verts[2] = getSortedVertex(2);
+    verts[3] = getSortedVertex(3);
   }
-  SVector3 normal()
+  SVector3 normal() const
   {
     double n[3];
     normal3points(_v[0]->x(), _v[0]->y(), _v[0]->z(),
@@ -66,7 +114,7 @@ class MFace {
 		  _v[2]->x(), _v[2]->y(), _v[2]->z(), n);
     return SVector3(n[0], n[1], n[2]);
   }
-  SPoint3 barycenter()
+  SPoint3 barycenter() const
   {
     SPoint3 p(0., 0., 0.);
     int n = getNumVertices();
@@ -83,4 +131,73 @@ class MFace {
   }
 };
 
+//--The following function objects compare the addresses of the mesh vertices.
+//--Equal, Less, and a Hash are defined.
+
+struct Equal_Face : public std::binary_function<MFace, MFace, bool> {
+  bool operator()(const MFace &f1, const MFace &f2) const
+  {
+    return (f1.getSortedVertex(0) == f2.getSortedVertex(0) &&
+            f1.getSortedVertex(1) == f2.getSortedVertex(1) &&
+            f1.getSortedVertex(2) == f2.getSortedVertex(2) &&
+            f1.getSortedVertex(3) == f2.getSortedVertex(3));
+  }
+};
+
+struct Equal_FacePtr : public std::binary_function<MFace*, MFace*, bool> {
+  bool operator()(const MFace *const f1, const MFace *const f2) const
+  {
+    return (f1->getSortedVertex(0) == f2->getSortedVertex(0) &&
+            f1->getSortedVertex(1) == f2->getSortedVertex(1) &&
+            f1->getSortedVertex(2) == f2->getSortedVertex(2) &&
+            f1->getSortedVertex(3) == f2->getSortedVertex(3));
+  }
+};
+
+struct Less_Face : public std::binary_function<MFace, MFace, bool> {
+  bool operator()(const MFace &f1, const MFace &f2) const
+  {
+    if(f1.getSortedVertex(0) < f2.getSortedVertex(0)) return true;
+    if(f1.getSortedVertex(0) > f2.getSortedVertex(0)) return false;
+    if(f1.getSortedVertex(1) < f2.getSortedVertex(1)) return true;
+    if(f1.getSortedVertex(1) > f2.getSortedVertex(1)) return false;
+    if(f1.getSortedVertex(2) < f2.getSortedVertex(2)) return true;
+    if(f1.getSortedVertex(2) > f2.getSortedVertex(2)) return false;
+    if(f1.getSortedVertex(3) < f2.getSortedVertex(3)) return true;
+    return false;
+  }
+};
+
+struct Less_FacePtr : public std::binary_function<MFace*, MFace*, bool> {
+  bool operator()(const MFace *const f1, const MFace *const f2) const
+  {
+    if(f1->getSortedVertex(0) < f2->getSortedVertex(0)) return true;
+    if(f1->getSortedVertex(0) > f2->getSortedVertex(0)) return false;
+    if(f1->getSortedVertex(1) < f2->getSortedVertex(1)) return true;
+    if(f1->getSortedVertex(1) > f2->getSortedVertex(1)) return false;
+    if(f1->getSortedVertex(2) < f2->getSortedVertex(2)) return true;
+    if(f1->getSortedVertex(2) > f2->getSortedVertex(2)) return false;
+    if(f1->getSortedVertex(3) < f2->getSortedVertex(3)) return true;
+    return false;
+  }
+};
+
+struct Hash_Face : public std::unary_function<MFace, size_t> {
+  size_t operator()(const MFace &f) const
+  {
+    const MVertex *v[4];
+    f.getOrderedVertices(v);
+    return HashFNV1a<sizeof(MVertex*[4])>::eval(v);
+  }
+};
+
+struct Hash_FacePtr : public std::unary_function<MFace*, size_t> {
+  size_t operator()(const MFace *const f) const
+  {
+    const MVertex *v[4];
+    f->getOrderedVertices(v);
+    return HashFNV1a<sizeof(MVertex*[4])>::eval(v);
+  }
+};
+
 #endif