[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Markus Bürg |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Thu, 23 May 2019 10:24:00 -0400 (EDT) |
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);
- [Getfem-commits] [getfem-commits] branch mb-Use_rtree_in_poly_composite created (now 22cd69d), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject),
Markus Bürg <=
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23
- [Getfem-commits] (no subject), Markus Bürg, 2019/05/23