Dear Yves and Kostas,
Thank you for your hints.
> I would like to ask if someone has a ready piece of code
> to share which searches for a convex enclosing given
> geometric point?
After looking at the source code of the classes:
getfem::interpolator_on_mesh_fem
and
getfem::interpolated_fem
I have found in them a similarly looking method:
bool find_a_point(base_node pt, base_node &ptr, size_type &cv)
which implements the algorithm I needed.
With some copy and paste from other classes :-) (mostly mesh_fem)
I ended up with the class convext_locator for which the source code
I send in the appendix, for those interested with it.
(I am not sure if I maintain context dependencies properly with this class
but it seems to work).
Here is an example of the use of convex_locator to refine a mesh
around a fixed geometric location couple of times:
getfem::mesh myMesh;
/* ... initialize myMesh */
dal::bit_vector e;
getfem::convex_locator locator(myMesh);
bgeot::base_node pt;bgeot::sc(pt) = 0.5, 0.3;
for (int i=0; i<3; i++) {
e.clear();
int i = locator.find_all_convexes(pt, e);
std::cout << "Found " << i << "convexes\n";
myMesh.Bank_refine(e);
}
I think that a comment is needed here. The convex search in
convex_locator is probably not 100% bullet proof -- consider for
instance point on the convex facet. I haven't looked at
the implementation of bgeot::geotrans_inv_convex used in the
convex search algorithm but I suspect that as with most of the
geometric predicates based on standard floating point arithmetic
there is a chance of some erratic behaviour of the predicate implementation.
Most of the time we can live with that but for some applications (mesh
generation or contact zone search) this might be an important issue.
Regards,
Roman
PS. Yves, if you think it is worth to add convex_locator to GetFEM
I can add proper documentation for it and test code.
--
Roman Putanowicz, PhD < [email protected] >
Institute for Computational Civil Engng (L-5)
Dept. of Civil Engng, Cracow Univ. of Technology
www.l5.pk.edu.pl, tel. +48 12 628 2569, fax 2034
namespace getfem {
class convex_locator : public context_dependencies {
public:
explicit convex_locator(const mesh &me);
convex_locator();
virtual ~convex_locator();
void init_with_mesh(const mesh &me);
void update_from_context(void) const;
bool find_a_convex(base_node pt, size_type &cv, base_node *ptr=NULL);
size_type find_all_convexes(base_node pt, dal::bit_vector &bv);
protected:
void build_rtree();
const mesh *linked_mesh_;
bgeot::rtree boxtree;
size_type cv_stored;
bgeot::rtree::pbox_set boxlst;
bgeot::geotrans_inv_convex gic;
};
void convex_locator::init_with_mesh(const mesh &me) {
GMM_ASSERT1(linked_mesh_ == 0, "Mesh level set already initialized");
linked_mesh_ = &me;
this->add_dependency(me);
this->build_rtree();
}
convex_locator::convex_locator ( const getfem::mesh& me ) {
linked_mesh_ = 0;
init_with_mesh(me);
}
convex_locator::convex_locator() {
linked_mesh_ = 0;
}
convex_locator::~convex_locator() {
}
bool convex_locator::find_a_convex ( base_node pt, size_type& cv, base_node* ptr ) {
GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized convex_locator");
context_check();
base_node refpt(pt.size());
if (NULL == ptr) {
ptr = &refpt;
}
bool gt_invertible;
if (cv_stored != size_type(-1) && gic.invert(pt, *ptr, gt_invertible)) {
cv = cv_stored;
if (gt_invertible) {
return true;
}
}
boxtree.find_boxes_at_point(pt, boxlst);
bgeot::rtree::pbox_set::const_iterator it = boxlst.begin(),ite = boxlst.end();
for (; it != ite; ++it) {
gic = bgeot::geotrans_inv_convex(linked_mesh_->convex((*it)->id),
linked_mesh_->trans_of_convex((*it)->id));
cv_stored = (*it)->id;
if (gic.invert(pt, *ptr, gt_invertible)) {
cv = (*it)->id;
return true;
}
}
return false;
}
size_type convex_locator::find_all_convexes( base_node pt, dal::bit_vector &bv) {
GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized convex_locator");
context_check();
size_type counter=0;
base_node refpt(pt.size());
bool gt_invertible;
boxtree.find_boxes_at_point(pt, boxlst);
bgeot::rtree::pbox_set::const_iterator it = boxlst.begin(),ite = boxlst.end();
for (; it != ite; ++it) {
gic = bgeot::geotrans_inv_convex(linked_mesh_->convex((*it)->id),
linked_mesh_->trans_of_convex((*it)->id));
cv_stored = (*it)->id;
if (gic.invert(pt, refpt, gt_invertible)) {
bv.add((*it)->id);
counter++;
}
}
return counter;
}
void convex_locator::update_from_context ( void ) const {
convex_locator *self = const_cast<convex_locator*>(this);
self->cv_stored = size_type(-1);
self->build_rtree();
}
void convex_locator::build_rtree () {
GMM_ASSERT1(linked_mesh_ != 0, "Uninitialized convex_locator");
base_node min, max;
scalar_type EPS=1E-13;
cv_stored = size_type(-1);
boxtree.clear();
for (dal::bv_visitor cv(linked_mesh_->convex_index()); !cv.finished(); ++cv) {
bounding_box(min, max, linked_mesh_->points_of_convex(cv),
linked_mesh_->trans_of_convex(cv));
for (unsigned k=0; k < min.size(); ++k) { min[k]-=EPS; max[k]+=EPS; }
boxtree.add_box(min, max, cv);
}
}
} /* namespace getfem */
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
