[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sat, 25 May 2019 10:42:55 -0400 (EDT) |
branch: mb-Use_rtree_in_poly_composite
commit 8d4fe766d3a3c9b8512ffba9f7b8531784f787dc
Author: Yves Renard <address@hidden>
Date: Sat May 25 16:42:39 2019 +0200
some fixes for the existing utilizations of rtree
---
src/bgeot_rtree.cc | 109 ++++++++++++++---------
src/getfem/bgeot_poly_composite.h | 2 +-
src/getfem/bgeot_rtree.h | 14 ++-
src/getfem_contact_and_friction_common.cc | 3 +-
src/getfem_contact_and_friction_large_sliding.cc | 1 +
src/getfem_fem_global_function.cc | 1 +
src/getfem_global_function.cc | 1 +
src/getfem_interpolated_fem.cc | 1 +
src/getfem_mesh_im_level_set.cc | 1 +
src/getfem_mesh_slicers.cc | 1 +
10 files changed, 81 insertions(+), 53 deletions(-)
diff --git a/src/bgeot_rtree.cc b/src/bgeot_rtree.cc
index 2ba9aa9..aea8279 100644
--- a/src/bgeot_rtree.cc
+++ b/src/bgeot_rtree.cc
@@ -53,8 +53,7 @@ namespace bgeot {
const base_node& min2, const base_node& max2,
scalar_type EPS) {
for (size_type i=0; i < min1.size(); ++i)
- if ((abs(min1[i] - min2[i]) > EPS && min1[i] > min2[i]) ||
- (abs(max1[i] - max2[i]) > EPS && max1[i] < max2[i])) return false;
+ if ((min1[i] > min2[i]+EPS) || (max1[i] < max2[i]-EPS)) return false;
return true;
}
@@ -62,8 +61,7 @@ namespace bgeot {
const base_node& min2, const base_node& max2,
scalar_type EPS) {
for (size_type i=0; i < min1.size(); ++i)
- if ((abs(max1[i] - min2[i]) > EPS && max1[i] < min2[i]) ||
- (abs(min1[i] - max2[i]) > EPS && min1[i] > max2[i])) return false;
+ if ((max1[i] < min2[i]-EPS) || (min1[i] > max2[i]+EPS)) return false;
return true;
}
@@ -71,8 +69,8 @@ namespace bgeot {
struct intersection_p {
const base_node &min, &max;
const scalar_type EPS;
- intersection_p(const base_node& min_, const base_node& max_, scalar_type
EPS)
- : min(min_), max(max_), EPS{EPS} {}
+ intersection_p(const base_node& min_, const base_node& max_, scalar_type
EPS_)
+ : min(min_), max(max_), EPS(EPS_) {}
bool operator()(const base_node& min2, const base_node& max2) const
{ return r1_inter_r2(min,max,min2,max2,EPS); }
bool accept(const base_node& min2, const base_node& max2) const
@@ -83,8 +81,8 @@ namespace bgeot {
struct contains_p {
const base_node &min, &max;
const scalar_type EPS;
- contains_p(const base_node& min_, const base_node& max_, scalar_type EPS)
- : min(min_), max(max_), EPS{EPS} {}
+ contains_p(const base_node& min_, const base_node& max_, scalar_type EPS_)
+ : min(min_), max(max_), EPS(EPS_) {}
bool operator()(const base_node& min2, const base_node& max2) const
{ return r1_ge_r2(min2,max2,min,max,EPS); }
bool accept(const base_node& min2, const base_node& max2) const
@@ -95,8 +93,8 @@ namespace bgeot {
struct contained_p {
const base_node &min, &max;
const scalar_type EPS;
- contained_p(const base_node& min_, const base_node& max_, scalar_type EPS)
- : min(min_), max(max_), EPS{EPS} {}
+ contained_p(const base_node& min_, const base_node& max_, scalar_type EPS_)
+ : min(min_), max(max_), EPS(EPS_) {}
bool accept(const base_node& min2, const base_node& max2) const
{ return r1_inter_r2(min,max,min2,max2,EPS); }
bool operator()(const base_node& min2, const base_node& max2) const
@@ -107,11 +105,11 @@ namespace bgeot {
struct has_point_p {
const base_node &P;
const scalar_type EPS;
- has_point_p(const base_node& P_, scalar_type EPS) : P(P_), EPS{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;
+ if (P[i] < min2[i]-EPS) return false;
+ if (P[i] > max2[i]+EPS) return false;
}
return true;
}
@@ -157,8 +155,8 @@ namespace bgeot {
intersect_line_and_box(const base_node& org_,
const base_small_vector &dirv_,
const base_node& min_, const base_node& max_,
- scalar_type EPS)
- : org(org_), dirv(dirv_), min(min_), max(max_), EPS{EPS} {}
+ scalar_type EPS_)
+ : org(org_), dirv(dirv_), min(min_), max(max_), EPS(EPS_) {}
bool operator()(const base_node& min2, const base_node& max2) const {
size_type N = org.size();
GMM_ASSERT1(N == min2.size(), "Dimensions mismatch");
@@ -181,14 +179,29 @@ namespace bgeot {
{ return operator()(min2,max2); }
};
-
- rtree::rtree(scalar_type EPS) : EPS{EPS}
+ size_type rtree::add_box(const base_node &min, const base_node &max,
+ size_type id) {
+ box_index bi;
+ if (tree_built) {
+ GMM_WARNING3("Add a box when the tree is already built cancel the tree. "
+ "Unefficient operation.");
+ tree_built = false; root = std::unique_ptr<rtree_elt_base>();
+ }
+ bi.min = &nodes[nodes.add_node(min, EPS)];
+ bi.max = &nodes[nodes.add_node(max, EPS)];
+ bi.id = (id + 1) ? id : boxes.size();
+ return boxes.emplace(std::move(bi)).first->id;
+ }
+
+
+ rtree::rtree(scalar_type EPS_) : EPS(EPS_), tree_built(false)
{}
void rtree::clear() {
root = std::unique_ptr<rtree_elt_base>();
boxes.clear();
nodes.clear();
+ tree_built = false;
}
template <typename Predicate>
@@ -212,37 +225,45 @@ namespace bgeot {
void rtree::find_intersecting_boxes(const base_node& bmin,
const base_node& bmax,
pbox_set& boxlst) const {
+
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(),boxlst,intersection_p(bmin,bmax, EPS));
+ GMM_ASSERT1(tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(),boxlst,intersection_p(bmin,bmax, EPS));
}
void rtree::find_containing_boxes(const base_node& bmin,
- const base_node& bmax, pbox_set& boxlst)
const {
+ const base_node& bmax,
+ pbox_set& boxlst) const {
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(), boxlst, contains_p(bmin,bmax, EPS));
+ GMM_ASSERT1( tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(), boxlst, contains_p(bmin,bmax, EPS));
}
void rtree::find_contained_boxes(const base_node& bmin,
- const base_node& bmax, pbox_set& boxlst)
const {
+ const base_node& bmax,
+ pbox_set& boxlst) const {
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(), boxlst, contained_p(bmin,bmax, EPS));
+ GMM_ASSERT1(tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(), boxlst, contained_p(bmin,bmax, EPS));
}
void rtree::find_boxes_at_point(const base_node& P, pbox_set& boxlst) const {
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(), boxlst, has_point_p(P, EPS));
+ GMM_ASSERT1(tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(), boxlst, has_point_p(P, EPS));
}
void rtree::find_line_intersecting_boxes(const base_node& org,
const base_small_vector& dirv,
pbox_set& boxlst) const {
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(),boxlst,intersect_line(org, dirv));
+ GMM_ASSERT1(tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(),boxlst,intersect_line(org, dirv));
}
void rtree::find_line_intersecting_boxes(const base_node& org,
@@ -251,9 +272,10 @@ namespace bgeot {
const base_node& bmax,
pbox_set& boxlst) const {
boxlst.clear();
- GMM_ASSERT2(root, "Boxtree not initialised.");
- find_matching_boxes_(root.get(), boxlst,
- intersect_line_and_box(org, dirv, bmin, bmax, EPS));
+ GMM_ASSERT1(tree_built, "Boxtree not initialised.");
+ if (root)
+ find_matching_boxes_(root.get(), boxlst,
+ intersect_line_and_box(org, dirv, bmin, bmax, EPS));
}
/*
@@ -339,17 +361,20 @@ namespace bgeot {
}
void rtree::build_tree() {
- if (boxes.size() == 0) return;
- getfem::local_guard lock = locks_.get_lock();
- assert(root == 0);
- pbox_cont b(boxes.size());
- pbox_cont::iterator b_it = b.begin();
- base_node bmin(*boxes.begin()->min), bmax(*boxes.begin()->max);
- for (box_cont::const_iterator it=boxes.begin(); it != boxes.end(); ++it) {
- update_box(bmin,bmax,*(*it).min,*(*it).max);
- *b_it++ = &(*it);
+ if (!tree_built) {
+ if (boxes.size() == 0) { tree_built = true; return; }
+ getfem::local_guard lock = locks_.get_lock();
+ assert(root == 0);
+ pbox_cont b(boxes.size());
+ pbox_cont::iterator b_it = b.begin();
+ base_node bmin(*boxes.begin()->min), bmax(*boxes.begin()->max);
+ for (box_cont::const_iterator it=boxes.begin(); it != boxes.end(); ++it)
{
+ update_box(bmin,bmax,*(*it).min,*(*it).max);
+ *b_it++ = &(*it);
+ }
+ root = build_tree_(b, bmin, bmax, 0);
+ tree_built = true;
}
- root = build_tree_(b, bmin, bmax, 0);
}
static void dump_tree_(rtree_elt_base *p, int level, size_type& count) {
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index 38dce95..1d696d9 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -131,7 +131,7 @@ namespace bgeot {
const base_matrix &M, const ITER &it, const base_node &p1, base_node &p2) {
for (dim_type d = 0; d < p2.size(); ++d) {
p2[d] = 0;
- auto &col = mat_col(M, d);
+ auto col = mat_col(M, d);
for (dim_type i = 0; i < p1.size(); ++i) p2[d] += col[i] * (*(it + i) -
p1[i]);
}
}
diff --git a/src/getfem/bgeot_rtree.h b/src/getfem/bgeot_rtree.h
index c996c8e..5c3c993 100644
--- a/src/getfem/bgeot_rtree.h
+++ b/src/getfem/bgeot_rtree.h
@@ -60,7 +60,7 @@ namespace bgeot {
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) {
+ if (gmm::abs(lhs[i] - rhs[i]) > EPS) {
return lhs[i] < rhs[i];
}
return false;
@@ -93,17 +93,12 @@ namespace bgeot {
using pbox_cont = std::vector<const box_index*>;
using pbox_set = std::set<const box_index *, box_index_id_compare>;
- rtree(scalar_type EPS = 1e-13);
+ rtree(scalar_type EPS = 0);
rtree(const rtree&) = delete;
rtree& operator = (const rtree&) = delete;
- 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();
- return boxes.emplace(std::move(bi)).first->id;
- }
+ size_type add_box(const base_node &min, const base_node &max,
+ size_type id=size_type(-1));
size_type nb_boxes() const { return boxes.size(); }
void clear();
@@ -174,6 +169,7 @@ namespace bgeot {
node_tab nodes;
box_cont boxes;
std::unique_ptr<rtree_elt_base> root;
+ bool tree_built;
getfem::lock_factory locks_;
};
diff --git a/src/getfem_contact_and_friction_common.cc
b/src/getfem_contact_and_friction_common.cc
index e24eae2..d0a5326 100644
--- a/src/getfem_contact_and_friction_common.cc
+++ b/src/getfem_contact_and_friction_common.cc
@@ -649,6 +649,7 @@ namespace getfem {
element_boxes_info.push_back(influence_box(i, cv, v.f(), n_mean));
}
}
+ element_boxes.build_tree();
}
void multi_contact_frame::compute_potential_contact_pairs_influence_boxes() {
@@ -1458,7 +1459,7 @@ namespace getfem {
}
}
}
-
+ face_boxes.build_tree();
}
public:
diff --git a/src/getfem_contact_and_friction_large_sliding.cc
b/src/getfem_contact_and_friction_large_sliding.cc
index 16cdd7e..7fd947e 100644
--- a/src/getfem_contact_and_friction_large_sliding.cc
+++ b/src/getfem_contact_and_friction_large_sliding.cc
@@ -1479,6 +1479,7 @@ namespace getfem {
face_of_elements.push_back(v.f());
}
}
+ element_boxes.build_tree();
}
diff --git a/src/getfem_fem_global_function.cc
b/src/getfem_fem_global_function.cc
index 74ef34b..e97e96e 100644
--- a/src/getfem_fem_global_function.cc
+++ b/src/getfem_fem_global_function.cc
@@ -87,6 +87,7 @@ namespace getfem {
functions[i]->bounding_box(bmin, bmax);
boxtree.add_box(bmin, bmax, i);
}
+ boxtree.build_tree();
scalar_type EPS=1E-13;
size_type max_dof(0);
diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc
index 1e140d0..c571072 100644
--- a/src/getfem_global_function.cc
+++ b/src/getfem_global_function.cc
@@ -852,6 +852,7 @@ namespace getfem {
for (auto&& val : bmax) val += EPS;
boxtree.add_box(bmin, bmax, cv);
}
+ boxtree.build_tree();
}
bool interpolator_on_mesh_fem::find_a_point(const base_node &pt,
diff --git a/src/getfem_interpolated_fem.cc b/src/getfem_interpolated_fem.cc
index e58e3e7..42614d1 100644
--- a/src/getfem_interpolated_fem.cc
+++ b/src/getfem_interpolated_fem.cc
@@ -33,6 +33,7 @@ namespace getfem {
for (unsigned k=0; k < min.size(); ++k) { min[k]-=EPS; max[k]+=EPS; }
boxtree.add_box(min, max, cv);
}
+ boxtree.build_tree();
}
bool interpolated_fem::find_a_point(base_node pt, base_node &ptr,
diff --git a/src/getfem_mesh_im_level_set.cc b/src/getfem_mesh_im_level_set.cc
index ef47b6c..ca36d85 100644
--- a/src/getfem_mesh_im_level_set.cc
+++ b/src/getfem_mesh_im_level_set.cc
@@ -667,6 +667,7 @@ namespace getfem {
size_type is = global_intersection.add_segment(i1, i2);
rtree_seg.add_box(min, max, is);
+ rtree_seg.build_tree(); // Not efficient !
const base_node &PE1
diff --git a/src/getfem_mesh_slicers.cc b/src/getfem_mesh_slicers.cc
index 446008b..2d09625 100644
--- a/src/getfem_mesh_slicers.cc
+++ b/src/getfem_mesh_slicers.cc
@@ -306,6 +306,7 @@ namespace getfem {
bgeot::bounding_box(min,max,slm.points_of_convex(cv),slm.trans_of_convex(cv));
tree.add_box(min, max, cv);
}
+ tree.build_tree();
}
void slicer_mesh_with_mesh::exec(mesh_slicer &ms) {