diff --git a/Common/OctreeInternals.cpp b/Common/OctreeInternals.cpp
index d76f27b86ac002f90b8130c7b64ab336d9f0ebd8..4e45e79ff5e203d56f43e20bb00ff9fc7efccae1 100644
--- a/Common/OctreeInternals.cpp
+++ b/Common/OctreeInternals.cpp
@@ -295,7 +295,8 @@ void *searchElement(octantBucket *_buckets_head, double *_pt, globalInfo *_globa
      
   ptrBucket = findElementBucket(_buckets_head, _pt);
   if (ptrBucket == NULL) {
-    Msg::Error("The point is not in the domain");
+    // this is not an error
+    Msg::Debug("The point is not in the domain");
     return NULL;
   }     
 
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 1fa18ff896806ab428623435f798f3330c43ec4c..57243bdafb345b825d600618da4126a964d1fa2d 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -28,6 +28,7 @@
 extern Context_T CTX;
 
 int MElement::_globalNum = 0;
+double MElement::_isInsideTolerance = 1.e-6;
 double MElementLessThanLexicographic::tolerance = 1.e-6;
 
 void MElement::_getEdgeRep(MVertex *v0, MVertex *v1, 
diff --git a/Geo/MElement.h b/Geo/MElement.h
index ab3475377bcc3edc528ddde04611b475969b22d3..3bfd5002f596ba58afbee2b25542572e2a5dbce8 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -33,6 +33,9 @@ class MElement
   // a visibility flag
   char _visible;
  protected:
+  // the tolerance used to determine if a point is inside an element,
+  // in parametric coordinates
+  static double _isInsideTolerance;
   void _getEdgeRep(MVertex *v0, MVertex *v1, 
                    double *x, double *y, double *z, SVector3 *n,
                    int faceIndex=-1);
@@ -55,6 +58,10 @@ class MElement
   // reset the global node number
   static void resetGlobalNumber(){ _globalNum = 0; }
 
+  // set/get the tolerance for isInside() test
+  static void setTolerance (const double tol){ _isInsideTolerance = tol; }
+  static double getTolerance () { return _isInsideTolerance; }
+
   // return the tag of the element
   virtual int getNum(){ return _num; }
 
@@ -187,7 +194,7 @@ class MElement
 
   // test if a point, given in parametric coordinates, belongs to the
   // element
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8) = 0;
+  virtual bool isInside(double u, double v, double w) = 0;
 
   // interpolate the given nodal data (resp. its gradient, curl and
   // divergence) at point (u,v,w) in parametric coordinates
@@ -296,7 +303,7 @@ class MPoint : public MElement {
   {
     s[0][0] = s[0][1] = s[0][2] = 0.;
   }
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
     return true;
   }
@@ -357,8 +364,9 @@ class MLine : public MElement {
     MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
   }
   virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < -(1. + tol) || u > (1. + tol))
       return false;
     return true;
@@ -569,8 +577,9 @@ class MTriangle : public MElement {
     MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
   }
   virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
       return false; 
     return true;
@@ -862,8 +871,9 @@ class MQuadrangle : public MElement {
     s[2][0] =  0.25 * (1. + v); s[2][1] =  0.25 * (1. + u); s[2][2] = 0.;
     s[3][0] = -0.25 * (1. + v); s[3][1] =  0.25 * (1. - u); s[3][2] = 0.;
   }
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < -(1. + tol) || v < -(1. + tol) || u > (1. + tol) || v > (1. + tol))
       return false;
     return true;
@@ -1190,8 +1200,9 @@ class MTetrahedron : public MElement {
   virtual double etaShapeMeasure();
   void xyz2uvw(double xyz[3], double uvw[3]);
   virtual const gmshFunctionSpace* getFunctionSpace(int o=-1) const;
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < (-tol) || v < (-tol) || w < (-tol) || u > ((1. + tol) - v - w))
       return false;
     return true;
@@ -1657,8 +1668,9 @@ class MHexahedron : public MElement {
     s[7][1] =  0.125 * (1. - u) * (1. + w);
     s[7][2] =  0.125 * (1. - u) * (1. + v);
   }
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < -(1. + tol) || v < -(1. + tol) || w < -(1. + tol) || 
        u > (1. + tol) || v > (1. + tol) || w > (1. + tol))
       return false;
@@ -2131,8 +2143,9 @@ class MPrism : public MElement {
     s[5][1] =  0.5 * (1. + w)    ;
     s[5][2] =  0.5 * v           ;
   }
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(w > (1. + tol) || w < -(1. + tol) || u < (1. + tol)
        || v < (1. + tol) || u > ((1. + tol) - v))
       return false;
@@ -2577,8 +2590,9 @@ class MPyramid : public MElement {
     s[4][1] = 0.;
     s[4][2] = 1.;
   }
-  virtual bool isInside(double u, double v, double w, double tol=1.e-8)
+  virtual bool isInside(double u, double v, double w)
   {
+    double tol = _isInsideTolerance;
     if(u < (w - (1. + tol)) || u > ((1. + tol) - w) || v < (w - (1. + tol)) ||
        v > ((1. + tol) - w) || w < (-tol) || w > (1. + tol))
       return false;
diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp
index 087ca754269d4d236ff587dc031c4fb92c80dcb2..a7731aa4afd831362d654e3ade6d22c6fc7f53d1 100644
--- a/Post/OctreePost.cpp
+++ b/Post/OctreePost.cpp
@@ -417,11 +417,14 @@ bool OctreePost::searchScalar(double x, double y, double z, double *values,
                               int step, double *size)
 {
   bool a = _searchScalar(x, y, z, values, step, size);
-  if (!a){
-    double oldeps = element::getTolerance();
+  if(!a){
+    double oldeps1 = element::getTolerance();
+    double oldeps2 = MElement::getTolerance();
     element::setTolerance(10.);
+    MElement::setTolerance(10.);
     a = _searchScalar(x, y, z, values, step, size);
-    element::setTolerance(oldeps);
+    element::setTolerance(oldeps1);
+    MElement::setTolerance(oldeps2);
   }    
   if (!a) Msg::Debug("No element found containing point (%g,%g,%g)", x, y, z);
   return a;