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