[Getfem-commits] (no subject)

[Markus Bürg]

branch: mb-Use_rtree_in_poly_composite
commit 4f307a386dcd4276cd27d5dd0f6b39851e904220
Author: mb <address@hidden>
Date:   Thu May 23 16:11:11 2019 +0200

    Use rtree in poly_composite to speed up evaluation.
---
 src/bgeot_poly_composite.cc         | 14 +++++++-
 src/getfem/bgeot_poly_composite.h   | 69 +++++++++++++++++++------------------
 src/getfem_fem_composite.cc         |  6 ++--
 src/getfem_integration_composite.cc |  6 ++--
 4 files changed, 54 insertions(+), 41 deletions(-)

diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 3d504ad..97c49ad 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -63,6 +63,10 @@ namespace bgeot {
   }
 
   mesh_precomposite::mesh_precomposite(const basic_mesh &m) {
+    initialise(m);
+  }
+
+  void mesh_precomposite::initialise(const basic_mesh &m) {
     msh = &m;
     det.resize(m.nb_convex());
     orgs.resize(m.nb_convex());
@@ -70,6 +74,8 @@ namespace bgeot {
     for (size_type i = 0; i <= m.points().index().last_true(); ++i) {
       vertexes.add(m.points()[i]);
     }
+
+    base_node min, max;
     for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
 
       pgeometric_trans pgt = m.trans_of_convex(cv);
@@ -89,8 +95,14 @@ namespace bgeot {
       gmm::mult(gmm::transposed(pc), gmm::transposed(G), B0);
       det[cv] = gmm::lu_inverse(B0);
       gtrans[cv] = B0;
-      orgs[cv] = m.points_of_convex(cv)[0];
+
+      auto points_of_convex = m.points_of_convex(cv);
+      orgs[cv] = points_of_convex[0];
+      bounding_box(min, max, points_of_convex);
+      box_to_convexes_map[box_tree.add_box(min, max)].push_back(cv);
     }
+
+    box_tree.build_tree();
   }
 
   DAL_TRIPLE_KEY(base_poly_key, short_type, short_type, 
std::vector<opt_long_scalar_type>);
diff --git a/src/getfem/bgeot_poly_composite.h 
b/src/getfem/bgeot_poly_composite.h
index 5fb79af..2546bba 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -43,6 +43,7 @@
 
 #include "bgeot_poly.h"
 #include "bgeot_mesh.h"
+#include "bgeot_rtree.h"
 
 // TODO : Use of rtree instead of dal::dynamic_tree_sorted<base_node,
 //        imbricated_box_less>
@@ -74,6 +75,8 @@ namespace bgeot {
 
     const basic_mesh *msh;
     PTAB vertexes;
+    rtree box_tree;
+    std::map<size_type, std::vector<size_type>> box_to_convexes_map;
     std::vector<base_matrix> gtrans;
     std::vector<scalar_type> det;
     std::vector<base_node> orgs;
@@ -83,6 +86,7 @@ namespace bgeot {
     dim_type dim(void) const { return msh->dim(); }
     pgeometric_trans trans_of_convex(size_type ic) const
     { return msh->trans_of_convex(ic); }
+    void initialise(const basic_mesh &m);
     
     mesh_precomposite(const basic_mesh &m);
     mesh_precomposite(void) : msh(0) {}
@@ -142,52 +146,49 @@ namespace bgeot {
       return poly_of_subelt(l).eval(p.begin());
     }
 
-    base_node p0(mp->dim());
+    auto dim = mp->dim();
+    base_node p0(dim);
     std::copy(it, it + mp->dim(), p0.begin());
-    mesh_structure::ind_cv_ct::const_iterator itc, itce;
 
-    mesh_precomposite::PTAB::const_sorted_iterator
-      it1 = mp->vertexes.sorted_ge(p0), it2 = it1;
-    size_type i1 = it1.index(), i2;
-
-    --it2; i2 = it2.index();
-
-
-    while (i1 != size_type(-1) || i2 != size_type(-1))
+    auto &box_tree = mp->box_tree;
+    rtree::pbox_set box_list;
+    box_tree.find_boxes_at_point(p0, box_list);
+    while (box_list.size())
     {
-      if (i1 != size_type(-1))
+      auto pmax_box = *box_list.begin();
+
+      if (box_list.size() > 1)
       {
-        const mesh_structure::ind_cv_ct &tc
-          = mp->linked_mesh().convex_to_point(i1);
-        itc = tc.begin(); itce = tc.end();
-        for (; itc != itce; ++itc)
+        auto max_rate = -1.0;
+        for (auto &&pbox : box_list)
         {
-          size_type ii = *itc;
-          mult_diff_transposed(mp->gtrans[ii], it, mp->orgs[ii], p0);
-          if (mp->trans_of_convex(ii)->convex_ref()->is_in(p0) < 1E-10) {
-            return to_scalar(poly_of_subelt(ii).eval(local_coordinate ? 
p0.begin() : it));
+          auto box_rate = 1.0;
+          for (size_type i = 0; i < dim; ++i)
+          {
+            auto h = pbox->max->at(i) - pbox->min->at(i);
+            if (h > 0)
+            {
+              auto rate = std::min(pbox->max->at(i) - p0[i], p0[i] - 
pbox->min->at(i)) / h;
+              box_rate = std::min(rate, box_rate);
+            }
+          }
+          if (box_rate > max_rate)
+          {
+            pmax_box = pbox;
+            max_rate = box_rate;
           }
         }
-        ++it1; i1 = it1.index();
       }
 
-      if (i2 != size_type(-1))
-      {
-        const mesh_structure::ind_cv_ct &tc
-          = mp->linked_mesh().convex_to_point(i2);
-        itc = tc.begin(); itce = tc.end();
-        for (; itc != itce; ++itc)
-        {
-          size_type ii = *itc;
-          mult_diff_transposed(mp->gtrans[ii], it, mp->orgs[ii], p0);
-          if (mp->trans_of_convex(ii)->convex_ref()->is_in(p0) < 1E-10) {
-            return to_scalar(poly_of_subelt(ii).eval(local_coordinate ? 
p0.begin() : it));
-          }
+      for (auto cv : mp->box_to_convexes_map.at(pmax_box->id)) {
+        mult_diff_transposed(mp->gtrans[cv], it, mp->orgs[cv], p0);
+        if (mp->trans_of_convex(cv)->convex_ref()->is_in(p0) < 1E-10) {
+          return to_scalar(poly_of_subelt(cv).eval(local_coordinate ? 
p0.begin() : it));
         }
-        --it2; i2 = it2.index();
       }
+      if (box_list.size() == 1) break;
+      box_list.erase(pmax_box);
     }
-
     GMM_ASSERT1(
       false, "Element not found in composite polynomial: " << base_node(*it, 
*it + mp->dim()));
   }
diff --git a/src/getfem_fem_composite.cc b/src/getfem_fem_composite.cc
index 2241b98..8236347 100644
--- a/src/getfem_fem_composite.cc
+++ b/src/getfem_fem_composite.cc
@@ -170,7 +170,7 @@ namespace getfem {
     m.add_triangle(i0, i2, i3);
     m.add_triangle(i0, i3, i1);
     m.add_triangle(i0, i1, i2);
-    mp = bgeot::mesh_precomposite(m);
+    mp.initialise(m);
 
     std::stringstream s("1-x-y;1-x-y;1-x-y;x;x;x;y;y;y;3-3*x-3*y;3*x;3*y;");
 
@@ -299,7 +299,7 @@ namespace getfem {
     m.add_triangle(i0, i2, i3);
     m.add_triangle(i0, i3, i1);
     m.add_triangle(i0, i1, i2);
-    mp = bgeot::mesh_precomposite(m);
+    mp.initialise(m);
 
 
  std::stringstream s
@@ -538,7 +538,7 @@ namespace getfem {
     m.add_triangle(i2, i0, i4);
     m.add_triangle(i3, i2, i4);
     m.add_triangle(i0, i1, i4);
-    mp = bgeot::mesh_precomposite(m);
+    mp.initialise(m);
 
  std::stringstream s
    ("2 - 3*x - 3*y + 6*x*y + x^3 - 3*x^2*y;"
diff --git a/src/getfem_integration_composite.cc 
b/src/getfem_integration_composite.cc
index 9562343..3da1f7b 100644
--- a/src/getfem_integration_composite.cc
+++ b/src/getfem_integration_composite.cc
@@ -130,7 +130,7 @@ namespace getfem {
     jfs.m.add_triangle(i0, i2, i3);
     jfs.m.add_triangle(i0, i3, i1);
     jfs.m.add_triangle(i0, i1, i2);
-    jfs.mp = bgeot::mesh_precomposite(jfs.m);
+    jfs.mp.initialise(jfs.m);
 
     mesh_im mi(jfs.m);
     mi.set_integration_method(jfs.m.convex_index(), pim);
@@ -169,7 +169,7 @@ namespace getfem {
     jfs.m.add_triangle(i2, i0, i4);
     jfs.m.add_triangle(i3, i2, i4);
     jfs.m.add_triangle(i0, i1, i4);
-    jfs.mp = bgeot::mesh_precomposite(jfs.m);
+    jfs.mp.initialise(jfs.m);
 
     mesh_im mi(jfs.m);
     mi.set_integration_method(jfs.m.convex_index(), pim);
@@ -205,7 +205,7 @@ namespace getfem {
     size_type i4 = jfs.m.add_point(base_node( 0.0,  0.0, 1.0));
     jfs.m.add_tetrahedron(i0, i1, i2, i4);
     jfs.m.add_tetrahedron(i1, i3, i2, i4);
-    jfs.mp = bgeot::mesh_precomposite(jfs.m);
+    jfs.mp.initialise(jfs.m);
 
     mesh_im mi(jfs.m);
     mi.set_integration_method(jfs.m.convex_index(), pim);