diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index c52b708871c0d1171229d2f33565bf6ced71c77d..87efec18eea866133f83201ea1ecba86ec92de55 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2840,6 +2840,11 @@ GEntity *GModel::addCone(std::vector<double> p1, std::vector<double> p2,
   return 0;
 }
 
+void GModel::healGeometry(double tolerance)
+{
+  if(_factory) _factory->healGeometry(this, tolerance);
+}
+
 GModel *GModel::computeBooleanUnion(GModel *tool, int createNewModel)
 {
   if(_factory)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index ccd36d4a601da672e281ccb89a0cf83b433cb27e..18d7d44760003016371336c7c8f644f4485b2cbc 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -524,6 +524,9 @@ class GModel
   GEntity *addCone(std::vector<double> p1, std::vector<double> p2, double radius1,
                    double radius2);
 
+  // heal geometry using the factory
+  void healGeometry(double tolerance = -1);
+
   // boolean operators acting on 2 models
   GModel *computeBooleanUnion(GModel *tool, int createNewModel=0);
   GModel *computeBooleanIntersection(GModel *tool, int createNewModel=0);
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 158208870803cde67b08df74117bd0986729305f..6bf5b365780642f6534ab15e39d03689fad9cc99 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -427,6 +427,16 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e,
 
 };
 
+void GeoFactory::healGeometry(GModel *gm, double tolerance)
+{
+  GModel *current = GModel::current();
+  GModel::setCurrent(gm);
+  ReplaceAllDuplicatesNew(tolerance);
+  gm->destroy();
+  gm->importGEOInternals();
+  GModel::setCurrent(current);
+}
+
 #if defined(HAVE_OCC)
 #include "OCCIncludes.h"
 #include "GModelIO_OCC.h"
@@ -1571,6 +1581,18 @@ GEntity *OCCFactory::addPipe(GModel *gm, GEntity *base, std::vector<GEdge *> wir
   return ret;
 }
 
+void OCCFactory::healGeometry(GModel *gm, double tolerance)
+{
+  if (tolerance < 0.)
+    tolerance = CTX::instance()->geom.tolerance;
+  if (!gm || !gm->_occ_internals)
+    return;
+  //gm->_occ_internals->healGeometry(tolerance, false, false, false, true, false, false);
+  gm->_occ_internals->healGeometry(tolerance, true, true, true, true, true, true);
+  gm->_occ_internals->buildLists();
+  gm->_occ_internals->buildGModel(gm);
+}
+
 //Prepare SGEOM integration
 #if defined(HAVE_SGEOM) && defined(HAVE_OCC)
 
@@ -1610,6 +1632,12 @@ GRegion* SGEOMFactory::addVolume(GModel *gm, std::vector<std::vector<GFace *> >
   return 0;
 }
 
+void SGEOMFactory::healGeometry(GModel *gm, double tolerance)
+{
+  Msg::Error("healGeometry not implemented yet for SGEOMFactory");
+  return 0;
+}
+
 #endif
 
 #endif
diff --git a/Geo/GModelFactory.h b/Geo/GModelFactory.h
index ba9aa2b6a4e06922087826e8c04d2788bd292738..ffab382cc8c9dd8765e9698246dc6f8e8f5114ee 100644
--- a/Geo/GModelFactory.h
+++ b/Geo/GModelFactory.h
@@ -217,6 +217,8 @@ class GModelFactory {
     Msg::Error("setPhysicalNumToEntitiesInBox not implemented yet");
   }
 
+  virtual void healGeometry(GModel *gm, double tolerance = -1.) = 0;
+
 };
 
 class GeoFactory : public GModelFactory {
@@ -230,6 +232,7 @@ class GeoFactory : public GModelFactory {
   std::vector<GFace *> addRuledFaces(GModel *gm, std::vector<std::vector<GEdge *> > edges);
   std::vector<GEntity*> extrudeBoundaryLayer(GModel *gm, GEntity *e, int nbLayers,
                                              double hLayers, int dir, int view);
+  void healGeometry(GModel *gm, double tolerance = -1.);
 };
 
 #if defined(HAVE_OCC)
@@ -288,6 +291,7 @@ class OCCFactory : public GModelFactory {
                                      std::vector<double> p1, std::vector<double> p2);
 
   void fillet(GModel *gm, std::vector<int> edges, double radius);
+  void healGeometry(GModel *gm, double tolerance = -1.);
 };
 
 #endif
@@ -301,6 +305,7 @@ class SGEOMFactory : public GModelFactory {
   GEdge *addLine(GModel *gm,GVertex *v1, GVertex *v2);
   GFace *addPlanarFace(GModel *gm, std::vector<std::vector<GEdge *> > edges);
   GRegion *addVolume(GModel *gm, std::vector<std::vector<GFace *> > faces);
+  void healGeometry(GModel *gm, double tolerance = -1.);
 };
 
 #endif
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 54d6c7a611e977c0a53765ce4dd2a3a7a1d55e16..ee81b72af36aa91bcc56e3b91dda991e8e3766fe 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2356,11 +2356,11 @@ static void MaxNumSurface(void *a, void *b)
     std::max(GModel::current()->getGEOInternals()->MaxSurfaceNum, s->Num);
 }
 
-static void ReplaceDuplicatePointsNew()
+static void ReplaceDuplicatePointsNew(double tol = -1.)
 {
   Msg::Info("New Coherence...");
-
-  double tol = CTX::instance()->geom.tolerance * CTX::instance()->lc;
+  if (tol < 0) 
+    tol = CTX::instance()->geom.tolerance * CTX::instance()->lc;
 
   // create kdtree
   std::map<MVertex*, Vertex*> v2V;
@@ -2927,6 +2927,15 @@ void ReplaceAllDuplicates()
   ReplaceAllDuplicates(report);
 }
 
+void ReplaceAllDuplicatesNew(double tol)
+{
+  if (tol < 0) 
+    tol = CTX::instance()->geom.tolerance * CTX::instance()->lc;
+  ReplaceDuplicatePointsNew(tol);
+  ReplaceDuplicateCurves(NULL);
+  ReplaceDuplicateSurfaces(NULL);
+}
+
 // Extrusion routines
 
 void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e)
diff --git a/Geo/Geo.h b/Geo/Geo.h
index 4097c296c0b95f3c4f29c452fc7f783659ecbf3f..f11a0edfd0d571461c372216d88fbe16553c6caf 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -392,6 +392,7 @@ int Extrude_ProtudeSurface(int type, int is,
 void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e);
 
 void ReplaceAllDuplicates();
+void ReplaceAllDuplicatesNew(double tol = -1.);
 
 bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2]);