Loading...
Searching...
No Matches
geo_nearest.cc
Go to the documentation of this file.
//
// given x, search x* the closest point in the boundary of omega
// gives also K* in mesh such that x* in K*
//
// author: Pierre.Saramito@imag.fr
//
// date: 9 october 2014
//
// implementation note:
// use CGAL::Point_set_2 for 2d nearest point search
// http://doc.cgal.org/latest/Point_set_2/group__PkgPointSet2.html
// TODO: 3D: with CGAL::Neighbor_search::Tree
// http://doc.cgal.org/latest/Spatial_searching/index.html#Chapter_dD_Spatial_Searching
// TODO: distributed: try to store vertex index in Point_2 data structure ? How to do that ?
// http://doc.cgal.org/latest/Kernel_23/index.html#title20 : extensible kernel
// http://stackoverflow.com/questions/2418332/customizing-cgal-kernel-with-my-own-point-class
#include "rheolef/geo_nearest.h"
#include "rheolef/geo.h"
#include "rheolef/point_util.h"
// internal includes:
#ifdef _RHEOLEF_HAVE_CGAL
#include "rheolef/cgal_traits.h"
#include <CGAL/Point_set_2.h>
#endif // _RHEOLEF_HAVE_CGAL
namespace rheolef {
#ifdef _RHEOLEF_HAVE_CGAL
// ---------------------------------------------------------------------
// 1) the geo_nearest interface
// ---------------------------------------------------------------------
template <class T, class M>
geo_nearest<T,M>::~geo_nearest()
{
if (_ptr != 0) {
delete_macro(_ptr);
}
}
template <class T, class M>
typename geo_nearest<T,M>::size_type
geo_nearest<T,M>::seq_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
if (_ptr == 0) { _ptr = make_ptr(omega); }
return _ptr->seq_nearest (omega, x, x_nearest);
}
template <class T, class M>
typename geo_nearest<T,M>::size_type
geo_nearest<T,M>::dis_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
if (_ptr == 0) { _ptr = make_ptr(omega); }
return _ptr->dis_nearest (omega, x, x_nearest);
}
// ---------------------------------------------------------------------
// 2) the abstract data structure
// ---------------------------------------------------------------------
template <class T, class M>
class geo_nearest_abstract_rep {
public:
typedef typename disarray<T,M>::size_type size_type;
virtual ~geo_nearest_abstract_rep() {}
virtual void initialize (const geo_base_rep<T,M>& omega) const = 0;
virtual size_type seq_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const = 0;
virtual size_type dis_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const = 0;
};
// ---------------------------------------------------------------------
// 3) the concrete, and dimension dependent, data structure
// = cgal "Point_set_2" delaunay mesh for 2D
// = cgal "Neighbor_search::Tree" octree for 3D
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
struct geo_nearest_rep : public geo_nearest_abstract_rep<T,M> { };
// ---------------------------------------------------------------------
// 3a) the 2D implementation with cgal::point_set
// ---------------------------------------------------------------------
template <class T, class M>
class geo_nearest_rep<T,M,2> : public geo_nearest_abstract_rep<T,M> {
public:
// typedef:
typedef typename geo_base_rep<T,M>::size_type size_type;
// allocators:
geo_nearest_rep() : _point_set(), _ball() {}
geo_nearest_rep(const geo_base_rep<T,M>& omega) : _point_set(), _ball() { initialize(omega); }
~geo_nearest_rep() {}
void initialize (const geo_base_rep<T,M>& omega) const;
// accessors:
size_type seq_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const;
size_type dis_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const;
protected:
typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
typedef typename Kernel::Point_2 Point_2;
// data:
mutable CGAL::Point_set_2<Kernel> _point_set;
mutable disarray<index_set,M> _ball;
};
template<class T, class M>
static
void
build_ball (
const geo_base_rep<T,M>& omega,
disarray<index_set,M>& ball)
{
typedef typename geo_basic<T,M>::size_type size_type;
// 1) build ball(xi) := { K; xi is a vertex of K }
size_type map_dim = omega.map_dimension();
if (map_dim <= 0) return; // a set of points have no neighbours
index_set emptyset;
distributor vertex_ownership = omega.sizes().ownership_by_dimension[0];
ball.resize (vertex_ownership, emptyset);
for (geo::const_iterator iter = omega.begin(map_dim), last = omega.end(map_dim); iter != last; ++iter) {
const geo_element& K = *iter;
index_set dis_ie_set;
dis_ie_set += K.dis_ie();
for (size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
size_type dis_inod = K [loc_inod];
size_type dis_iv = omega.dis_inod2dis_iv (dis_inod);
ball.dis_entry (dis_iv) += dis_ie_set; // union with {dis_ie}
}
}
ball.dis_entry_assembly();
}
template <class T, class M>
void
geo_nearest_rep<T,M,2>::initialize(const geo_base_rep<T,M>& omega) const
{
size_type map_dim = omega.map_dimension();
distributor vertex_ownership = omega.sizes().ownership_by_dimension[0];
distributor side_ownership = omega.sizes().ownership_by_dimension[map_dim-1];
size_type first_dis_iv = vertex_ownership.first_index();
size_type first_dis_isid = side_ownership.first_index();
std::vector<bool> marked (omega.n_vertex(), false);
#ifdef TODO
// TODO: passer boundary_guard de geo_basic a geo_base_rep
boundary_guard (*this);
#endif // TODO
check_macro (omega.have_domain_indirect ("boundary"), "nearest: geo may have boundary domain defined");
const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
std::list<Point_2> Lr;
for (typename domain_indirect_basic<M>::const_iterator_ioige iter = boundary.ioige_begin(), last = boundary.ioige_end(); iter != last; ++iter) {
size_type dis_isid = (*iter).index();
const geo_element& S = omega.dis_get_geo_element(map_dim-1,dis_isid);
for (size_type iloc = 0, nloc = S.size(); iloc < nloc; ++iloc) {
size_type dis_iv = S[iloc];
if (!vertex_ownership.is_owned(dis_iv)) continue;
size_type iv = dis_iv - first_dis_iv;
if (marked[iv]) continue;
const point& xi = omega.node(iv);
Point_2 pi (xi[0],xi[1]);
Lr.push_back (pi);
marked[iv] = true;
}
}
_point_set.insert (Lr.begin(),Lr.end());
build_ball (omega, _ball);
omega.neighbour_guard();
}
template <class T, class M>
typename geo_nearest_rep<T,M,2>::size_type
geo_nearest_rep<T,M,2>::seq_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
// -----------------------------------------
// 1) nearest neighbor : get its coordinates
// -----------------------------------------
Point_2 x_cgal (x[0],x[1]);
typename CGAL::Point_set_2<Kernel>::Vertex_handle v_nearest = _point_set.nearest_neighbor(x_cgal);
point xv_nearest (v_nearest->point().x(), v_nearest->point().y());
x_nearest = xv_nearest;
// ------------------------------------------
// 2) nearest neighbor : get its vertex index
// ------------------------------------------
// TODO: try to store vertex index in Point_2 data structure ? How to do that ?
// => avoid a geo_locate(x) and facilitate the distributed version impementation
size_type map_dim = omega.map_dimension();
size_type dis_ie0 = 0;
#ifdef _RHEOLEF_HAVE_MPI
size_type nproc = omega.comm().size();
if (is_distributed<M>::value && nproc > 1) {
fatal_macro ("dis_nearest/dis_locate: not yet");
} else {
dis_ie0 = omega.seq_locate (xv_nearest);
}
#else // _RHEOLEF_HAVE_MPI
dis_ie0 = omega.seq_locate (xv_nearest);
#endif // _RHEOLEF_HAVE_MPI
check_macro (dis_ie0 != std::numeric_limits<size_type>::max(), "invalid element containing nearest point");
const geo_element& K0 = omega.dis_get_geo_element (map_dim, dis_ie0);
size_type dis_iv_nearest = std::numeric_limits<size_type>::max();
for (size_type iloc = 0, nloc = K0.size(); iloc < nloc; ++iloc) {
if (dist(omega.node(K0[iloc]), xv_nearest) == 0) {
dis_iv_nearest = K0[iloc];
break;
}
}
// --------------------------------------------
// 3) search in ball(K0) for nearest bdry point
// --------------------------------------------
index_set ball_nearest = _ball.dis_at (dis_iv_nearest);
check_macro (map_dim > 1, "invalid map_dim="<<map_dim);
size_type dis_ie_nearest = dis_ie0;
Float d_nearest = dist(x,xv_nearest);
distributor element_ownership = omega.sizes().ownership_by_dimension[map_dim];
distributor side_ownership = omega.sizes().ownership_by_dimension[map_dim-1];
for (index_set::const_iterator iter = ball_nearest.begin(), last = ball_nearest.end(); iter != last; ++iter) {
size_type dis_ie = *iter;
const geo_element& K = omega.dis_get_geo_element (map_dim, dis_ie);
for (size_type iloc_sid = 0, nloc_sid = K.n_subgeo(map_dim-1); iloc_sid < nloc_sid; ++iloc_sid) {
size_type dis_isid = (map_dim == 2) ? K.edge(iloc_sid) : K.face(iloc_sid);
if (!side_ownership.is_owned(dis_isid)) continue;
const geo_element& S = omega.dis_get_geo_element (map_dim-1, dis_isid);
if (S.master(1) != std::numeric_limits<size_type>::max()) continue; // not a bdry side
if (map_dim == 2) {
// 1) S=[a,b]: compute y = the intersection D(a,b) with D(x,n)
// that satisfies: ay.n = 0 and xy.t = 0
const point& a = omega.node(S[0]);
const point& b = omega.node(S[1]);
point t = b-a;
point n (t[1], -t[0]);
Float ab2 = dist2(a,b);
Float det = -ab2;
Float f0 = t[0]*x[0] + t[1]*x[1];
Float f1 = n[0]*a[0] + n[1]*a[1];
point f (f0,f1);
Float y0 = (n[1]*f[0] - t[1]*f[1])/det;
Float y1 = - (n[0]*f[0] - t[0]*f[1])/det;
point y (y0,y1);
check_macro (fabs(dot(y-a,n)) + fabs(dot(y-x,t)) < std::numeric_limits<Float>::epsilon(), "intersection problem");
if (dot(y-a,t) < 0) y = a; // y is outside of segment [a,b], at left
else if (dot(y-b,t) > 0) y = b; // y is outside of segment [a,b], at right
Float d = dist(x, y);
if (d >= d_nearest) continue;
x_nearest = y;
d_nearest = d;
dis_ie_nearest = dis_ie;
} else {
error_macro ("3D: not yet");
}
}
}
trace_macro ("point " << x << " nearest_point " << x_nearest);
return dis_ie_nearest;
}
template <class T, class M>
typename geo_nearest_rep<T,M,2>::size_type
geo_nearest_rep<T,M,2>::dis_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
size_type dis_ie = seq_nearest (omega, x, x_nearest);
#ifdef _RHEOLEF_HAVE_MPI
size_type nproc = omega.comm().size();
if (is_distributed<M>::value && nproc > 1) {
fatal_macro ("dis_nearest: not yet");
}
#endif // _RHEOLEF_HAVE_MPI
return dis_ie;
}
// ---------------------------------------------------------------------
// 3b) the 3D implementation with cgal::Tree
// ---------------------------------------------------------------------
// TODO
// ---------------------------------------------------------------------
// 3c) allocator
// ---------------------------------------------------------------------
template <class T, class M>
geo_nearest_abstract_rep<T,M>*
geo_nearest<T,M>::make_ptr (const geo_base_rep<T,M>& omega)
{
check_macro (omega.dimension() == omega.map_dimension(), "geo_nearest: map_dim < dim: not supported");
switch (omega.map_dimension()) {
case 2: return new_macro((geo_nearest_rep<T,M,2>)(omega));
//case 3: return new_macro((geo_nearest_rep<T,M,3>)(omega));
default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
}
}
// ---------------------------------------------------------------------
// 6) geo::nearest()
// ---------------------------------------------------------------------
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::seq_nearest (
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
return _nearestor.seq_nearest (*this, x, x_nearest);
}
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::dis_nearest (
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
return _nearestor.dis_nearest (*this, x, x_nearest);
}
template <class T, class M>
static
void
sequential_nearest (
const geo_base_rep<T,M>& omega,
const disarray<point_basic<T>,M>& x,
disarray<point_basic<T>,M>& x_nearest,
disarray<std::size_t, M>& dis_ie)
{
typedef typename geo_base_rep<T,M>::size_type size_type;
dis_ie.resize (x.ownership());
x_nearest.resize (x.ownership());
for (size_type i = 0, n = x.size(); i < n; ++i) {
dis_ie[i] = omega.seq_locate (x[i], dis_ie[i]);
if (dis_ie[i] != std::numeric_limits<size_type>::max()) {
x_nearest[i] = x[i];
continue;
}
dis_ie[i] = omega.seq_nearest (x[i], x_nearest[i]);
check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
"nearest: failed at x="<<ptos(x[i],omega.dimension()));
}
}
template <class T>
void
geo_rep<T,sequential>::nearest (
const disarray<point_basic<T>,sequential>& x,
disarray<point_basic<T>,sequential>& x_nearest,
disarray<size_type, sequential>& dis_ie) const
{
sequential_nearest (*this, x, x_nearest, dis_ie);
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::nearest (
const disarray<point_basic<T>,distributed>& x,
disarray<point_basic<T>,distributed>& x_nearest,
disarray<size_type, distributed>& dis_ie) const
{
size_type nproc = base::comm().size();
if (nproc == 1) {
sequential_nearest (*this, x, x_nearest, dis_ie);
return;
}
// TODO: also search for nearest when we go outside of the domain
// TODO: search seq the nearest in the current partition and then
// TODO: merge all nearests across partitions ?
locate (x, dis_ie);
x_nearest = x;
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// no CGAL: no nearest
// ----------------------------------------------------------------------------
#else // ! _RHEOLEF_HAVE_CGAL
template <class T, class M>
geo_nearest<T,M>::~geo_nearest()
{
}
template <class T, class M>
typename geo_nearest<T,M>::size_type
geo_nearest<T,M>::seq_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
return 0;
}
template <class T, class M>
typename geo_nearest<T,M>::size_type
geo_nearest<T,M>::dis_nearest (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
return 0;
}
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::seq_nearest (
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
return 0;
}
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::dis_nearest (
const point_basic<T>& x,
point_basic<T>& x_nearest) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
return 0;
}
template <class T>
void
geo_rep<T,sequential>::nearest (
const disarray<point_basic<T>, sequential>& x,
disarray<point_basic<T>, sequential>& x_nearest,
disarray<typename geo_base_rep<T,sequential>::size_type, sequential>& idx) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::nearest (
const disarray<point_basic<T>, distributed>& x,
disarray<point_basic<T>, distributed>& x_nearest,
disarray<typename geo_base_rep<T,distributed>::size_type, distributed>& idx) const
{
fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
}
#endif // _RHEOLEF_HAVE_MPI
#endif // _RHEOLEF_HAVE_CGAL
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template class geo_nearest<Float,M>; \
template class geo_base_rep<Float,M>; \
template class geo_rep<Float,M>;
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
_RHEOLEF_instanciation
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
size_type
field::size_type size_type
Definition branch.cc:430
Float
see the Float page for the full documentation
point
see the point page for the full documentation
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::disarray::size_type
rep::base::size_type size_type
Definition disarray.h:501
rheolef::geo_base_rep
base class for M=sequential or distributed meshes representations
Definition geo.h:528
rheolef::geo_base_rep::map_dimension
size_type map_dimension() const
Definition geo.h:564
rheolef::geo_base_rep::seq_locate
size_type seq_locate(const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
Definition geo_locate.cc:405
rheolef::geo_base_rep::size_type
base::size_type size_type
Definition geo.h:533
rheolef::geo_base_rep::dimension
size_type dimension() const
Definition geo.h:563
rheolef::geo_base_rep::seq_nearest
size_type seq_nearest(const point_basic< T > &x, point_basic< T > &x_nearest) const
Definition geo_nearest.cc:331
rheolef::geo_base_rep::dis_nearest
size_type dis_nearest(const point_basic< T > &x, point_basic< T > &x_nearest) const
Definition geo_nearest.cc:339
rheolef::geo_nearest::seq_nearest
size_type seq_nearest(const geo_base_rep< T, M > &omega, const point_basic< T > &x, point_basic< T > &x_nearest) const
Definition geo_nearest.cc:62
rheolef::geo_nearest::dis_nearest
size_type dis_nearest(const geo_base_rep< T, M > &omega, const point_basic< T > &x, point_basic< T > &x_nearest) const
Definition geo_nearest.cc:72
rheolef::geo_nearest::~geo_nearest
~geo_nearest()
Definition geo_nearest.cc:54
rheolef::geo_nearest::make_ptr
static geo_nearest_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
Definition geo_nearest.cc:317
rheolef::geo_nearest::size_type
disarray< T, M >::size_type size_type
Definition geo_nearest.h:43
rheolef::geo_rep< T, distributed >::size_type
base::size_type size_type
Definition geo.h:934
rheolef::geo_rep
sequential mesh representation
Definition geo.h:778
rheolef::point_basic
Definition point.h:86
rheolef::t
rheolef::space_base_rep< T, M > t
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.n
int n
Definition mkgeo_ball.sh:150
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
mkgeo_ball.map_dim
int map_dim
Definition mkgeo_ball.sh:337
mkgeo_ugrid.boundary
int boundary
Definition mkgeo_ugrid.sh:181
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition space_numbering.cc:177
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::boundary_guard
void boundary_guard(const geo_basic< T, M > &omega)
Definition geo.cc:600
rheolef::initialize
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
Definition field_expr_recursive.h:627
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition vec_expr_v2.h:415
rheolef::dist2
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
rheolef::ptos
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
rheolef::dist
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
d
size_t d
Definition neumann-laplace-lambda.cc:28
f
Definition cavity_dg.h:29
rheolef::distributed
Definition communicator.h:128
rheolef::is_distributed::value
static const bool value
Definition communicator.h:71
rheolef::sequential
Definition communicator.h:62
M
Expr1::memory_type M
Definition vec_expr_v2.h:416