use absolute geometricl tolerance when creating bboxes in octree; better

isInside tolerance in ref coordinates

Geo/MElement.h

@@ -264,7 +264,8 @@ class MElement
   // integration routines
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   {
-    Msg::Error("No integration points defined for this type of element: %d",this->getType());
+    Msg::Error("No integration points defined for this type of element: %d",
+               this->getType());
   }
 
   // IO routines

Geo/MElementOctree.cpp

@@ -7,6 +7,7 @@
 #include "MElement.h"
 #include "MElementOctree.h"
 #include "Octree.h"
+#include "Context.h"
 
 static void MElementBB(void *a, double *min, double *max)
 {
@@ -21,6 +22,13 @@ static void MElementBB(void *a, double *min, double *max)
     min[1] = std::min(min[1], v->y()); max[1] = std::max(max[1], v->y());
     min[2] = std::min(min[2], v->z()); max[2] = std::max(max[2], v->z());
   }
+
+  // make bounding boxes larger up to (absolute) geometrical tolerance
+  double eps = CTX::instance()->geom.tolerance;
+  for(int i = 0; i < 3; i++){
+    min[i] -= eps;
+    max[i] += eps;
+  }
 }
 
 static void MElementCentroid(void *a, double *x)
@@ -50,10 +58,15 @@ static int MElementInEle(void *a, double *x)
 MElementOctree::MElementOctree(GModel *m)
 {
   SBoundingBox3d bb = m->bounds();
-  double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
-  double size[3] = {bb.max().x() - bb.min().x(),
-                    bb.max().y() - bb.min().y(),
-                    bb.max().z() - bb.min().z()};
+  // make bounding box larger up to (absolute) geometrical tolerance
+  double eps = CTX::instance()->geom.tolerance;
+  SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
+  bbmin -= bbeps;
+  bbmax += bbeps;
+  double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
+  double size[3] = {bbmax.x() - bbmin.x(),
+                    bbmax.y() - bbmin.y(),
+                    bbmax.z() - bbmin.z()};
   const int maxElePerBucket = 100; // memory vs. speed trade-off
   _octree = Octree_Create(maxElePerBucket, min, size,
                           MElementBB, MElementCentroid, MElementInEle);
@@ -75,10 +88,15 @@ MElementOctree::MElementOctree(std::vector<MElement*> &v)
                     v[i]->getVertex(j)->z());
     }
   }
-  double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
-  double size[3] = {bb.max().x() - bb.min().x(),
-                    bb.max().y() - bb.min().y(),
-                    bb.max().z() - bb.min().z()};
+  // make bounding box larger up to (absolute) geometrical tolerance
+  double eps = CTX::instance()->geom.tolerance;
+  SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
+  bbmin -= bbeps;
+  bbmax += bbeps;
+  double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
+  double size[3] = {bbmax.x() - bbmin.x(),
+                    bbmax.y() - bbmin.y(),
+                    bbmax.z() - bbmin.z()};
   const int maxElePerBucket = 100; // memory vs. speed trade-off
   _octree = Octree_Create(maxElePerBucket, min, size,
                           MElementBB, MElementCentroid, MElementInEle);

Mesh/Field.cpp

@@ -1146,7 +1146,9 @@ class PostViewField : public Field
       update_needed = false;
     }
     double l = 0.;
-    if(!octree->searchScalarWithTol(x, y, z, &l, 0, 0, 10.))
+    // use large tolerance (in element reference coordinates) to
+    // maximize chance of finding an element
+    if(!octree->searchScalarWithTol(x, y, z, &l, 0, 0, 1.))
       Msg::Info("No element found containing point (%g,%g,%g)", x, y, z);
     if(l <= 0 && crop_negative_values) return MAX_LC;
     return l;

Post/OctreePost.cpp

@@ -13,6 +13,7 @@
 #include "shapeFunctions.h"
 #include "GModel.h"
 #include "MElement.h"
+#include "Context.h"
 
 // helper routines for list-based views
 
@@ -33,6 +34,13 @@ static void minmax(int n, double *X, double *Y, double *Z,
     max[1] = (Y[i] > max[1]) ? Y[i] : max[1];
     max[2] = (Z[i] > max[2]) ? Z[i] : max[2];
   }
+
+  // make bounding boxes larger up to (absolute) geometrical tolerance
+  double eps = CTX::instance()->geom.tolerance;
+  for(int i = 0; i < 3; i++){
+    min[i] -= eps;
+    max[i] += eps;
+  }
 }
 
 static void centroid(int n, double *X, double *Y, double *Z, double *c)
@@ -233,10 +241,15 @@ OctreePost::OctreePost(PView *v)
     }
 
     SBoundingBox3d bb = l->getBoundingBox();
-    double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
-    double size[3] = {bb.max().x() - bb.min().x(),
-                      bb.max().y() - bb.min().y(),
-                      bb.max().z() - bb.min().z()};                   
+    // make bounding box larger up to (absolute) geometrical tolerance
+    double eps = CTX::instance()->geom.tolerance;
+    SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
+    bbmin -= bbeps;
+    bbmax += bbeps;
+    double min[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
+    double size[3] = {bbmax.x() - bbmin.x(),
+                      bbmax.y() - bbmin.y(),
+                      bbmax.z() - bbmin.z()};
     const int maxElePerBucket = 100; // memory vs. speed trade-off
     
     _SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);

Post/PViewData.h

@@ -213,6 +213,9 @@ class PViewData {
   virtual bool isRemote(){ return false; }
   virtual int fillRemoteVertexArrays(std::string &options){ return 0; }
 
+  // get MElement (if view supports it)
+  virtual MElement *getElement(int step, int entity, int element);
+
   // I/O routines
   virtual bool writeSTL(std::string fileName);
   virtual bool writeTXT(std::string fileName);
@@ -220,8 +223,7 @@ class PViewData {
                         bool append=false);
   virtual bool writeMSH(std::string fileName, bool binary=false);
   virtual bool writeMED(std::string fileName);
-  //
-  virtual MElement *getElement (int step, int entity, int element);
+
   static void registerBindings(binding *b);
 };
 

Post/PViewDataGModel.h

@@ -227,6 +227,8 @@ class PViewDataGModel : public PViewData {
   bool getValueByIndex(int step, int dataIndex, int node, int comp, double &val);
   // get underlying model
   GModel* getModel(int step){ return _steps[step]->getModel(); }
+  // get MElement
+  MElement *getElement(int step, int entity, int element);
 
   // Add some data "on the fly" (data is stored in a map, indexed by
   // node or element number depending on the type of dataset)
@@ -240,7 +242,6 @@ class PViewDataGModel : public PViewData {
   bool writeMSH(std::string fileName, bool binary=false);
   bool readMED(std::string fileName, int fileIndex);
   bool writeMED(std::string fileName);
-  MElement *getElement (int step, int entity, int element);
 };
 
 #endif