[Getfem-commits] (no subject)

[Markus Bürg]

branch: mb-Use_rtree_in_poly_composite
commit 29ab6b9a79ab8775b4a50cf42c53473b39895dc2
Author: mb <address@hidden>
Date:   Thu May 23 16:16:05 2019 +0200

    Provide tolerance to correctly capture points with round-off errors.
---
 src/bgeot_rtree.cc       | 17 ++++++++++++-----
 src/getfem/bgeot_rtree.h |  3 ++-
 2 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/bgeot_rtree.cc b/src/bgeot_rtree.cc
index 1412dbd..ae03ce5 100644
--- a/src/bgeot_rtree.cc
+++ b/src/bgeot_rtree.cc
@@ -99,10 +99,13 @@ namespace bgeot {
   /* match boxes containing P */
   struct has_point_p {
     const base_node P;
-    has_point_p(const base_node& P_) : P(P_) {}
-    bool operator()(const base_node& min2, const base_node& max2) {
-      for (size_type i=0; i < P.size(); ++i)
-        if (P[i] < min2[i] || P[i] > max2[i]) return false;
+    const scalar_type EPS;
+    has_point_p(const base_node& P_, scalar_type EPS) : P(P_), EPS{EPS} {}
+    bool operator()(const base_node& min2, const base_node& max2) const {
+      for (size_type i = 0; i < P.size(); ++i) {
+        if ((abs(P[i] - min2[i]) > EPS) && (P[i] < min2[i])) return false;
+        if ((abs(max2[i] - P[i]) > EPS) && (P[i] > max2[i])) return false;
+      }
       return true;
     }
     bool accept(const base_node& min2, const base_node& max2) const
@@ -170,6 +173,9 @@ namespace bgeot {
   };
 
 
+  rtree::rtree(scalar_type EPS) : EPS{EPS}
+  {}
+
   template <typename Predicate>
   static void find_matching_boxes_(rtree_elt_base *n, rtree::pbox_set& boxlst,
                                    const Predicate &p) {
@@ -209,7 +215,8 @@ namespace bgeot {
 
   void rtree::find_boxes_at_point(const base_node& P, pbox_set& boxlst) {
     boxlst.clear(); if (!root) build_tree();
-    if (root) find_matching_boxes_(root.get(), boxlst, has_point_p(P));
+    GMM_ASSERT1(root, "Boxtree not initialised.");
+    find_matching_boxes_(root.get(), boxlst, has_point_p(P, EPS));
   }
 
   void rtree::find_line_intersecting_boxes(const base_node& org,
diff --git a/src/getfem/bgeot_rtree.h b/src/getfem/bgeot_rtree.h
index d26f5f8..dd7a496 100644
--- a/src/getfem/bgeot_rtree.h
+++ b/src/getfem/bgeot_rtree.h
@@ -75,7 +75,7 @@ namespace bgeot {
     using pbox_cont = std::vector<const box_index*>;
     using pbox_set = std::set<const box_index *, box_index_compare>;
 
-    rtree() = default;
+    rtree(scalar_type EPS = 1e-13);
     rtree(const rtree&) = delete;
     rtree& operator = (const rtree&) = delete;
 
@@ -150,6 +150,7 @@ namespace bgeot {
         idvec.push_back((*it)->id);
     }
 
+    const scalar_type EPS;
     box_cont boxes;
     std::unique_ptr<rtree_elt_base> root;
     getfem::lock_factory locks_;