diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 31033616d6d817613f273c3298a37ac1a896c752..2f8e145200f2cd74173db8addc510c25e83e2f82 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -743,13 +743,13 @@ MElement *GModel::getMeshElementByCoord(SPoint3 &p, int dim, bool strict)
   return _octree->find(p.x(), p.y(), p.z(), dim, strict);
 }
 
-std::vector<MElement*> GModel::getMeshElementsByCoord(SPoint3 &p, int dim)
+std::vector<MElement*> GModel::getMeshElementsByCoord(SPoint3 &p, int dim, bool strict)
 {
   if(!_octree){
     Msg::Debug("Rebuilding mesh element octree");
     _octree = new MElementOctree(this);
   }
-  return _octree->findAll(p.x(), p.y(), p.z(), dim);
+  return _octree->findAll(p.x(), p.y(), p.z(), dim, strict);
 }
 
 MVertex *GModel::getMeshVertexByTag(int n)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index cc7b3c3374430fa2e7631e7bef2bc8b8e2b802f4..8aa87fd874cc24cad759acb6f01b0bb8010beacf 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -313,7 +313,7 @@ class GModel
 
   // access a mesh element by coordinates (using an octree search)
   MElement *getMeshElementByCoord(SPoint3 &p, int dim=-1, bool strict=true);
-  std::vector<MElement*> getMeshElementsByCoord(SPoint3 &p, int dim=-1);
+  std::vector<MElement*> getMeshElementsByCoord(SPoint3 &p, int dim=-1, bool strict=true);
 
   // access a mesh element by tag, using the element cache
   MElement *getMeshElementByTag(int n);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 96085d86d9b17fc58b11b3576e1daf12e772c5c9..0127d3ea354d5b17f398937c7380ea75dad04b8b 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -615,7 +615,7 @@ int GModel::readMSH(const std::string &name)
             MElement *e = createElementMSH(this, num, type, physical, elementary,
                                            partition, vertices, elements, physicals,
                                            own, p);
-	    
+
 #if (FAST_ELEMENTS==1)
 	  elems[num] = e;
 	  elemreg[num] = elementary;
@@ -631,7 +631,7 @@ int GModel::readMSH(const std::string &name)
           delete [] data;
           numElementsPartial += numElms;
         }
-      } 
+      }
 #if (FAST_ELEMENTS==1)
 	for(int i = 0; i < 10; i++)
 	  elements[i].clear();
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index 7a614916c9f496f1440057e699f79134f964c58c..ff352299a01a2f3c8383984866cd971ee4d3e52f 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -123,8 +123,11 @@ MElementOctree::~MElementOctree()
 }
 
 
-std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, int dim)
+std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, int dim, bool strict)
 {
+  double maxTol = 1.;
+  double tolIncr = 10.;
+
   double P[3] = {x, y, z};
   std::vector<void*> v;
   std::vector<MElement*> e;
@@ -134,6 +137,53 @@ std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, in
     MElement *el = (MElement*) *it;
     if (dim == -1 || el->getDim() == dim)e.push_back(el);
   }
+  if (e.empty() && !strict && _gm) {
+    double initialTol = MElement::getTolerance();
+    double tol = initialTol;
+    while (tol < maxTol){
+      tol *= tolIncr;
+      MElement::setTolerance(tol);
+      std::vector<GEntity*> entities;
+      _gm->getEntities(entities);
+      for(unsigned int i = 0; i < entities.size(); i++){
+	for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+	  MElement* el = entities[i]->getMeshElement(j);
+	  if (dim == -1 ||  el->getDim() == dim){
+	    if (MElementInEle(el, P)){
+              e.push_back(el);
+	    }
+	  }
+	}
+      }
+      if(!e.empty()) {
+        MElement::setTolerance(initialTol);
+        return e;
+      }
+    }
+    MElement::setTolerance(initialTol);
+  }
+  else if (e.empty() && !strict && !_gm){
+    double initialTol = MElement::getTolerance();
+    double tol = initialTol;
+    while (tol < maxTol){
+      tol *= tolIncr;
+      MElement::setTolerance(tol);
+      for(unsigned int i = 0; i < _elems.size(); i++){
+        MElement* el = _elems[i];
+        if (dim == -1 ||  el->getDim() == dim){
+          if (MElementInEle(el, P)){
+            e.push_back(el);
+          }
+        }
+      }
+      if(!e.empty()) {
+        MElement::setTolerance(initialTol);
+        return e;
+      }
+    }
+    MElement::setTolerance(initialTol);
+    //Msg::Warning("Point %g %g %g not found",x,y,z);
+  }
   return e;
 }
 
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 86d226fa98efb7a88604af7615ef83c7e7dbdd6e..f93ce0299ab5ef62c4b5caeb89513c9beaed0ae8 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -23,7 +23,7 @@ class MElementOctree{
   ~MElementOctree();
   MElement *find(double x, double y, double z, int dim = -1, bool strict = false);
   Octree *getInternalOctree(){ return _octree; }
-  std::vector<MElement *> findAll(double x, double y, double z, int dim);
+  std::vector<MElement *> findAll(double x, double y, double z, int dim, bool strict = false);
 };
 
 #endif
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 8e60f50a6af679f62974b488c7a1ace7892755b5..f5231ff9f9dd24e06ef0171ae827ea325def9579 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -11,7 +11,7 @@
 /*
  * MLine
  *
- *   0----------1 --> u
+ *   0-----+-----1 --> u
  *
  */
 class MLine : public MElement {