[Getfem-commits] (no subject)

[Markus Bürg]

branch: mb-Use_rtree_in_poly_composite
commit 5b45591e62bc979a5eeb6357bbcead5b748d2955
Author: mb <address@hidden>
Date:   Thu May 23 16:22:51 2019 +0200

    Make sure that boxes in rtree are unique to obtain optimal performance.
---
 src/getfem/bgeot_rtree.h                        | 33 +++++++++++---
 src/getfem_generic_assembly_interpolation.cc    | 58 ++++++++++++++-----------
 src/getfem_interpolation_on_deformed_domains.cc | 32 +++++++-------
 3 files changed, 76 insertions(+), 47 deletions(-)

diff --git a/src/getfem/bgeot_rtree.h b/src/getfem/bgeot_rtree.h
index eac890a..c481b13 100644
--- a/src/getfem/bgeot_rtree.h
+++ b/src/getfem/bgeot_rtree.h
@@ -40,6 +40,7 @@
 
 #include <set>
 #include "bgeot_small_vector.h"
+#include "bgeot_node_tab.h"
 
 namespace bgeot {
 
@@ -48,12 +49,29 @@ namespace bgeot {
     const base_node *min, *max;
   };
 
-  struct box_index_compare {
+  struct box_index_id_compare {
     bool operator()(const box_index *plhs, const box_index *prhs) const {
       return plhs->id < prhs->id;
     }
   };
 
+  struct box_index_topology_compare {
+    bool is_less(const base_node &lhs, const base_node &rhs) const {
+      GMM_ASSERT2(lhs.size() == rhs.size(), "size mismatch");
+      const scalar_type EPS = 1e-13;
+      for (size_type i = 0; i < lhs.size(); ++i)
+        if (abs(lhs[i] - rhs[i]) > EPS) {
+          return lhs[i] < rhs[i];
+        }
+      return false;
+    }
+
+    bool operator()(const box_index &lhs, const box_index &rhs) const {
+      return is_less(*lhs.min, *rhs.min) ||
+             (!is_less(*rhs.min, *lhs.min) && is_less(*lhs.max, *rhs.max));
+    }
+  };
+
   struct rtree_elt_base {
     enum { RECTS_PER_LEAF=8 };
     bool isleaf_;
@@ -71,18 +89,20 @@ namespace bgeot {
    */
   class rtree {
   public:
-    using box_cont = std::deque<box_index> ;
+    using box_cont = std::set<box_index,box_index_topology_compare> ;
     using pbox_cont = std::vector<const box_index*>;
-    using pbox_set = std::set<const box_index *, box_index_compare>;
+    using pbox_set = std::set<const box_index *, box_index_id_compare>;
 
     rtree(scalar_type EPS = 1e-13);
     rtree(const rtree&) = delete;
     rtree& operator = (const rtree&) = delete;
 
-    void add_box(base_node min, base_node max, size_type id=size_type(-1)) {
-      box_index bi; bi.min = min; bi.max = max;
+    size_type add_box(const base_node &min, const base_node &max, size_type 
id=size_type(-1)) {
+      box_index bi;
+      bi.min = &nodes[nodes.add_node(min, EPS)];
+      bi.max = &nodes[nodes.add_node(max, EPS)];
       bi.id = (id + 1) ? id : boxes.size();
-      boxes.push_back(bi);
+      return boxes.emplace(std::move(bi)).first->id;
     }
     size_type nb_boxes() const { return boxes.size(); }
     void clear() { root = std::unique_ptr<rtree_elt_base>(); boxes.clear(); }
@@ -151,6 +171,7 @@ namespace bgeot {
     }
 
     const scalar_type EPS;
+    node_tab nodes;
     box_cont boxes;
     std::unique_ptr<rtree_elt_base> root;
     getfem::lock_factory locks_;
diff --git a/src/getfem_generic_assembly_interpolation.cc 
b/src/getfem_generic_assembly_interpolation.cc
index 561b788..3ea6c00 100644
--- a/src/getfem_generic_assembly_interpolation.cc
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -514,6 +514,9 @@ namespace getfem {
     mutable bool extract_variable_done;
     mutable bool extract_data_done;
 
+  private:
+    mutable std::map<size_type, std::vector<size_type>> box_to_convexes;
+
   public:
     void update_from_context() const {
       recompute_elt_boxes = true;
@@ -587,6 +590,7 @@ namespace getfem {
       // Element_boxes update (if necessary)
       if (recompute_elt_boxes) {
 
+        box_to_convexes.clear();
         element_boxes.clear();
         base_node bmin(N), bmax(N);
         const dal::bit_vector&
@@ -621,7 +625,7 @@ namespace getfem {
           for (auto &&val : bmin) val -= h*0.2;
           for (auto &&val : bmax) val += h*0.2;
 
-          element_boxes.add_box(bmin, bmax, cv);
+          box_to_convexes[element_boxes.add_box(bmin, bmax)].push_back(cv);
         }
         element_boxes.build_tree();
         recompute_elt_boxes = false;
@@ -703,32 +707,34 @@ namespace getfem {
       size_type best_cv(-1);
       base_node best_P_ref;
       for (size_type i = boxes.size(); i > 0; --i) {
-
-        cv = boxes[i-1]->id;
-        gic.init(target_mesh.points_of_convex(cv),
-                 target_mesh.trans_of_convex(cv));
-
-        bool converged;
-        bool is_in = gic.invert(P, P_ref, converged, 1E-4);
-        // cout << "cv = " << cv << " P = " << P << " P_ref = " << P_ref << 
endl;
-        // cout << " is_in = " << int(is_in) << endl;
-        // for (size_type iii = 0;
-        //     iii < target_mesh.points_of_convex(cv).size(); ++iii)
-        //  cout << target_mesh.points_of_convex(cv)[iii] << endl;
-
-        if (converged) {
-          if (is_in) {
-            face_num = short_type(-1); // Should detect potential faces ?
-            ret_type = 1;
-            break;
-          } else {
-            scalar_type dist
-              = target_mesh.trans_of_convex(cv)->convex_ref()->is_in(P_ref);
-            if (dist < best_dist) {
-              best_dist = dist;
-              best_cv = cv;
-              best_P_ref = P_ref;
+        for (auto convex : box_to_convexes.at(boxes[i-1]->id)) {
+          gic.init(target_mesh.points_of_convex(convex),
+                   target_mesh.trans_of_convex(convex));
+
+          bool converged;
+          bool is_in = gic.invert(P, P_ref, converged, 1E-4);
+          // cout << "cv = " << cv << " P = " << P << " P_ref = " << P_ref << 
endl;
+          // cout << " is_in = " << int(is_in) << endl;
+          // for (size_type iii = 0;
+          //     iii < target_mesh.points_of_convex(cv).size(); ++iii)
+          //  cout << target_mesh.points_of_convex(cv)[iii] << endl;
+
+          if (converged) {
+            cv = convex;
+            if (is_in) {
+              face_num = short_type(-1); // Should detect potential faces ?
+              ret_type = 1;
+              break;
+            } else {
+              scalar_type dist
+                = target_mesh.trans_of_convex(cv)->convex_ref()->is_in(P_ref);
+              if (dist < best_dist) {
+                best_dist = dist;
+                best_cv = cv;
+                best_P_ref = P_ref;
+              }
             }
+            break;
           }
         }
       }
diff --git a/src/getfem_interpolation_on_deformed_domains.cc 
b/src/getfem_interpolation_on_deformed_domains.cc
index c6f8610..fecd9ca 100644
--- a/src/getfem_interpolation_on_deformed_domains.cc
+++ b/src/getfem_interpolation_on_deformed_domains.cc
@@ -85,8 +85,8 @@ class  interpolate_transformation_on_deformed_domains
   contact_boundary slave; //also marked with a source or X prefix/suffix
 
   mutable bgeot::rtree element_boxes;
-  mutable std::vector<size_type> box_to_convex; //index to obtain
-                                                //a convex number from a box 
number
+  mutable std::map<size_type, std::vector<size_type>> box_to_convex; //index 
to obtain a convex
+                                                                     //number 
from a box number
   mutable bgeot::geotrans_inv_convex gic;
   mutable fem_precomp_pool fppool;
 
@@ -115,7 +115,6 @@ class  interpolate_transformation_on_deformed_domains
     dal::bit_vector points_already_interpolated;
     std::vector<base_node> transformed_points(m.nb_max_points());
     box_to_convex.clear();
-    box_to_convex.reserve(region.size());
 
     for (getfem::mr_visitor v(region, m); !v.finished(); ++v) {
       auto cv   = v.cv();
@@ -154,8 +153,7 @@ class  interpolate_transformation_on_deformed_domains
       }
 
       // Store the bounding box and additional information.
-      element_boxes.add_box(bmin, bmax, box_to_convex.size());
-      box_to_convex.push_back(cv);
+      box_to_convex[element_boxes.add_box(bmin, bmax)].push_back(cv);
     }
     element_boxes.build_tree();
   }
@@ -294,17 +292,21 @@ public:
     while (!bset.empty())
     {
       auto itmax = most_central_box(bset, pt_x);
-      cv = box_to_convex[(*itmax)->id];
-      auto deformed_nodes_y = deformed_master_nodes(cv);
-      gic.init(deformed_nodes_y, target_mesh.trans_of_convex(cv));
-      auto converged = true;
-      auto is_in = gic.invert(pt_x, P_ref, converged);
-      if (is_in && converged) {
-        face_num = static_cast<short_type>(-1);
-        transformation_success = true;
-        break;
+
+      for (auto i : box_to_convex.at((*itmax)->id))
+      {
+        auto deformed_nodes_y = deformed_master_nodes(i);
+        gic.init(deformed_nodes_y, target_mesh.trans_of_convex(i));
+        auto converged = true;
+        auto is_in = gic.invert(pt_x, P_ref, converged);
+        if (is_in && converged) {
+          cv = i;
+          face_num = static_cast<short_type>(-1);
+          transformation_success = true;
+          break;
+        }
       }
-      if (bset.size() == 1) break;
+      if (transformation_success || (bset.size() == 1)) break;
       bset.erase(itmax);
     }