geo_locate.cc

Go to the documentation of this file.

//
// given x, search K in mesh such that x in K
//
// author: Pierre.Saramito@imag.fr
//
// date: 12 march 2012
//
// implementation note:
//  use CGAL::Segment_tree for element location
//
#include "rheolef/geo_locate.h"
#include "rheolef/geo.h"
#include "rheolef/geo_element_contains.h"
#include "rheolef/point_util.h"
 
// internal includes:
#ifdef _RHEOLEF_HAVE_CGAL
#include "rheolef/cgal_traits.h"
#include <CGAL/Segment_tree_k.h>
#include <CGAL/Range_segment_tree_traits.h>
#endif // _RHEOLEF_HAVE_CGAL
 
#ifdef _RHEOLEF_HAVE_CGAL
namespace rheolef { namespace details {
// CGAL::Segment_tree_1::make_tree is broken in recent CGAL (6.x):
// it passes a Tree_anchor pointer instead of a reference to the
// Segment_tree_d constructor. Rebuild the 1D tree by hand.
template <class Tree, class Iterator>
inline void cgal_make_tree (Tree& tree, Iterator first, Iterator last)
{
  tree.make_tree (first, last);
}
template <class Traits, class Iterator>
inline void cgal_make_tree (CGAL::Segment_tree_1<Traits>& tree, Iterator first, Iterator last)
{
  typedef typename CGAL::Segment_tree_1<Traits>::Segment_tree_1_type tree_type;
  delete tree.segment_tree_1;
  tree.segment_tree_1 = new tree_type (*tree.anchor);
  tree.segment_tree_1->make_tree (first, last);
}
}} // namespace rheolef::details
#endif // _RHEOLEF_HAVE_CGAL
 
namespace rheolef {
 
// ---------------------------------------------------------------------
// 1) utility: bbox of a geo_element
// ---------------------------------------------------------------------
template <class T, class M>
void
compute_bbox (const geo_base_rep<T,M>& omega, const geo_element& K, point_basic<T>& xmin, point_basic<T>& xmax) 
{
  // TODO: omega.order == 1 only: K.size == K.n_node
  typedef typename geo_base_rep<T,M>::size_type size_type;
  size_type d = omega.dimension();
  for (size_type j = 0; j < d; j++) {
    xmin[j] =  std::numeric_limits<T>::max();
    xmax[j] = -std::numeric_limits<T>::max();
  }
  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
    size_type dis_inod = K[iloc];
    const point_basic<T>& x = omega.dis_node(dis_inod);
    for (size_type j = 0; j < d; j++) {
      xmin[j] = std::min(x[j], xmin[j]);
      xmax[j] = std::max(x[j], xmax[j]);
    }
  }
}
#ifdef _RHEOLEF_HAVE_CGAL
// -----------------------------------------------------------------------------------
// 2) missing in cgal: inspirated from "examples/RangeSegmentTrees/include/Tree_Traits.h"
// -----------------------------------------------------------------------------------
template <class Kernel, class Val>
class Segment_tree_map_traits_1 {
 public:
  typedef  Val                        Value;
  typedef typename Kernel::FT         Point_1 ;
  typedef Point_1                     Key;
  typedef Point_1                     Key_1;
  typedef std::pair<Key,Key>          Pure_interval;
  typedef std::pair<Pure_interval, Val> Interval;
 
  class C_Low_1{
  public:
    Key_1 operator()(const Interval& i)
    { return i.first.first;}
  };
 
  class C_High_1{
  public:
    Key_1 operator()(const Interval& i)
    { return i.first.second;}
  };
 
  class C_Key_1{
  public:
    Key_1 operator()(const Key& k)
    { return k;}
  };
 
  class C_Compare_1{
  public:
    bool operator()(Key_1 k1, Key_1 k2)
    {
      return std::less<Float>()(k1,k2);
    }
  };
 
  typedef C_Compare_1 compare_1;
  typedef C_Low_1 low_1;
  typedef C_High_1 high_1;
  typedef C_Key_1 key_1;
};
// -----------------------------------------------------------------------------------
// 3) helper a generic dimension=D implementation with cgal
// -----------------------------------------------------------------------------------
 
template <class T, size_t D> struct cgal_locate_traits {};
 
template <class T>
struct cgal_locate_traits<T,1> {
 
// typedef:
 
  typedef typename geo_base_rep<T,sequential>::size_type                 size_type;
 
  typedef typename geo_cgal_traits<T,1>::Kernel                       Kernel;
  typedef Segment_tree_map_traits_1<Kernel,size_type>                 Traits;
  typedef ::CGAL::Segment_tree_1<Traits>                              Segment_tree_type;
  typedef typename Traits::Interval                                   Interval;
  typedef typename Traits::Pure_interval                              Pure_interval;
  typedef typename Traits::Key                                        Key;
 
  static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
  {
    return Pure_interval(Key(xmin[0]),
                         Key(xmax[0]));
  }
  static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
  {
    return Pure_interval (Key(x[0]-eps),
                          Key(x[0]+eps));
  }
  static void put (std::ostream& out, const Pure_interval& b) {
    out << "[" << b.first << "," << b.second << "["; 
  }
};
template <class T>
struct cgal_locate_traits<T,2> {
 
// typedef:
 
  typedef typename geo_base_rep<T,sequential>::size_type                 size_type;
 
  typedef typename geo_cgal_traits<T,2>::Kernel                       Kernel;
  typedef ::CGAL::Segment_tree_map_traits_2<Kernel,size_type>         Traits;
  typedef ::CGAL::Segment_tree_2<Traits>                              Segment_tree_type;
  typedef typename Traits::Interval                                   Interval;
  typedef typename Traits::Pure_interval                              Pure_interval;
  typedef typename Traits::Key                                        Key;
 
  static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
  {
    return Pure_interval(Key(xmin[0], xmin[1]),
                         Key(xmax[0], xmax[1]));
  }
  static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
  {
    return Pure_interval (Key(x[0]-eps, x[1]-eps),
                          Key(x[0]+eps, x[1]+eps));
  }
  static void put (std::ostream& out, const Pure_interval& b) {
    out << "[" << b.first.x() << "," << b.second.x() << "[x[" 
               << b.first.y() << "," << b.second.y() << "["; 
  }
};
template <class T>
struct cgal_locate_traits<T,3> {
 
// typedef:
 
  typedef typename geo_base_rep<T,sequential>::size_type                 size_type;
 
  typedef typename geo_cgal_traits<T,3>::Kernel                       Kernel;
  typedef ::CGAL::Segment_tree_map_traits_3<Kernel,size_type>         Traits;
  typedef ::CGAL::Segment_tree_3<Traits>                              Segment_tree_type;
  typedef typename Traits::Interval                                   Interval;
  typedef typename Traits::Pure_interval                              Pure_interval;
  typedef typename Traits::Key                                        Key;
 
  static Pure_interval make_cgal_bbox (const point_basic<T>& xmin, const point_basic<T>& xmax)
  {
    return Pure_interval(Key(xmin[0], xmin[1], xmin[2]),
                         Key(xmax[0], xmax[1], xmax[2]));
  }
  static Pure_interval make_cgal_point_window (const point_basic<T>& x, const T& eps)
  {
    return Pure_interval (Key(x[0]-eps, x[1]-eps, x[2]-eps),
                          Key(x[0]+eps, x[1]+eps, x[2]+eps));
  }
  static void put (std::ostream& out, const Pure_interval& b) {
    out << "[" << b.first.x() << "," << b.second.x() << "[x[" 
               << b.first.y() << "," << b.second.y() << "[x["
               << b.first.z() << "," << b.second.z() << "["; 
  }
};
// ---------------------------------------------------------------------
// 4) the geo_locate interface
// ---------------------------------------------------------------------
template <class T, class M>
geo_locate<T,M>::~geo_locate()
{
  if (_ptr != 0) {
    delete_macro(_ptr);
  }
}
template <class T, class M>
typename geo_locate<T,M>::size_type
geo_locate<T,M>::seq_locate (
    const geo_base_rep<T,M>& omega,
    const point_basic<T>&    x,
    size_type                dis_ie_guest) const
{
  if (_ptr == 0) { _ptr = make_ptr(omega); }
  return _ptr->seq_locate (omega, x, dis_ie_guest);
}
template <class T, class M>
typename geo_locate<T,M>::size_type
geo_locate<T,M>::dis_locate (
    const geo_base_rep<T,M>& omega,
    const point_basic<T>&    x,
    size_type                dis_ie_guest) const
{
  if (_ptr == 0) { _ptr = make_ptr(omega); }
  return _ptr->dis_locate (omega, x, dis_ie_guest);
}
// ---------------------------------------------------------------------
// 5) the tree box abstract data structure
// ---------------------------------------------------------------------
template <class T, class M>
class geo_locate_abstract_rep {
public:
  typedef typename disarray<T,M>::size_type size_type;
  virtual ~geo_locate_abstract_rep() {}
  virtual void initialize (const geo_base_rep<T,M>& omega) const = 0;
  virtual size_type seq_locate (
        const geo_base_rep<T,M>& omega,
        const point_basic<T>&    x,
        size_type                dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
  virtual size_type dis_locate (
        const geo_base_rep<T,M>& omega,
        const point_basic<T>&    x,
        size_type                dis_ie_guest = std::numeric_limits<size_type>::max()) const = 0;
};
// ---------------------------------------------------------------------
// 5) the tree box concrete, and dimension dependent, data structure
//    = cgal "segment" tree
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
class geo_locate_rep : public geo_locate_abstract_rep<T,M> {
public:
 
// typedef:
 
  typedef typename geo_base_rep<T,M>::size_type                          size_type;
 
  typedef typename cgal_locate_traits<T,D>::Segment_tree_type         Segment_tree_type;
  typedef typename cgal_locate_traits<T,D>::Interval                  Interval;
 
// allocators:
 
  geo_locate_rep() : _tree() {}
  geo_locate_rep(const geo_base_rep<T,M>& omega) : _tree() { initialize(omega); } 
  ~geo_locate_rep() {}
  void initialize (const geo_base_rep<T,M>& omega) const; 
 
// accessors:
 
  size_type seq_locate (
        const geo_base_rep<T,M>& omega,
        const point_basic<T>&    x,
        size_type                dis_ie_guest = std::numeric_limits<size_type>::max()) const;
  size_type dis_locate (
        const geo_base_rep<T,M>& omega,
        const point_basic<T>&    x,
        size_type                dis_ie_guest = std::numeric_limits<size_type>::max()) const;
 
// data:
protected:
  mutable Segment_tree_type _tree; // cgal window query requires a mutable _tree (why?)
};
// ---------------------------------------------------------------------
// 5a) allocator
// ---------------------------------------------------------------------
template <class T, class M>
geo_locate_abstract_rep<T,M>*
geo_locate<T,M>::make_ptr (const geo_base_rep<T,M>& omega)
{
  check_macro (omega.dimension() == omega.map_dimension(), "geo_locate: map_dim < dim: not supported");
  switch (omega.dimension()) {
   case 1: return new_macro((geo_locate_rep<T,M,1>)(omega));
   case 2: return new_macro((geo_locate_rep<T,M,2>)(omega));
   case 3: return new_macro((geo_locate_rep<T,M,3>)(omega));
   default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
  }
}
// ---------------------------------------------------------------------
// 5b) initialize
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
void
geo_locate_rep<T,M,D>::initialize (const geo_base_rep<T,M>& omega) const
{
  // create the corresponding vector of bounding boxes
  trace_macro ("geo::locate initialize...");
  std::list<Interval> boxes;
  point_basic<T> xmin, xmax;
  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
    const geo_element& K = omega[ie];
    compute_bbox (omega, K, xmin, xmax);
    boxes.push_back (Interval(cgal_locate_traits<T,D>::make_cgal_bbox (xmin,xmax),ie));
  }
  // create tree:
  // _tree.clear(); // TODO : when e.g. nodes update, rebuild boxes & tree
  details::cgal_make_tree (_tree, boxes.begin(), boxes.end());
  trace_macro ("geo::locate initialize done");
}
// ---------------------------------------------------------------------
// 5c) locate
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
typename geo_locate_rep<T,M,D>::size_type
geo_locate_rep<T,M,D>::seq_locate (
  const geo_base_rep<T,M>&  omega,
  const point_basic<T>&     x,
  size_type                 dis_ie_guest) const
{
  if (dis_ie_guest != std::numeric_limits<size_type>::max()) {
    // have a guest:
    const distributor& ownership = omega.ownership();
    if (ownership.is_owned (dis_ie_guest)) {
      size_type first_dis_ie = ownership.first_index();
      size_type ie_guest = dis_ie_guest - first_dis_ie;
      const geo_element& K_guest = omega[ie_guest];
      bool intersect = details::contains (K_guest, omega.get_nodes(), x);
      if (intersect) {
        return K_guest.dis_ie();
      }
    }
  }
  typedef typename geo_base_rep<T,M>::size_type size_type;
  size_type dis_ie = std::numeric_limits<size_type>::max();
  // epsilon is only for bbox: then, predicate on element K is exact
  // TODO: compute epsilon with omega.hmin scale ?
  // static const T eps = 1e5*std::numeric_limits<T>::epsilon();
  // static const T eps = 1000*std::numeric_limits<T>::epsilon();
  static const T eps = sqrt(std::numeric_limits<T>::epsilon());
  Interval xe = Interval (cgal_locate_traits<T,D>::make_cgal_point_window (x, eps), 0);
  std::list<Interval> intersected_boxes;
  // point query = inverse range query ; from ::CGAL documentation:
  // "In order to perform an inverse range query, a range query of epsilon width has to be performed.
  // We preferred not to offer an extra function for this sort of query, since the inverse range
  // query is a special case of the range query (window_query)"
  _tree.window_query (xe, std::back_inserter(intersected_boxes));
  for (typename std::list<Interval>::iterator j = intersected_boxes.begin(); j != intersected_boxes.end(); j++) {
    size_type ie = (*j).second;
    const geo_element& K = omega[ie];
    bool intersect = details::contains (K, omega.get_nodes(), x); 
    if (intersect) {
      dis_ie = K.dis_ie();
      break;
    }
  }
  return dis_ie;
}
// ---------------------------------------------------------------------
// 5c) dis_locate = one distributed computation
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
typename geo_locate_rep<T,M,D>::size_type
geo_locate_rep<T,M,D>::dis_locate (
  const geo_base_rep<T,M>&  omega,
  const point_basic<T>&     x,
  size_type                 dis_ie_guest) const
{
  // mpi::reduce with std::minimum handles unsigned=size_t as signed=long
  // so use here long to reduce, otherwise leads to a bug
  using signed_size_type = long; 
  const signed_size_type signed_large = -1;
  size_type dis_ie = seq_locate (omega, x, dis_ie_guest);
#ifdef _RHEOLEF_HAVE_MPI
  if (omega.comm().size() > 1 && is_distributed<M>::value) {
    // fix for numeric_limits<size_t> == -1 with long int
    dis_ie = mpi::all_reduce (omega.comm(), signed_size_type(dis_ie), mpi::maximum<signed_size_type>());
  }
#endif // _RHEOLEF_HAVE_MPI
  return dis_ie;
}
// ---------------------------------------------------------------------
// 6) geo::locate()
// ---------------------------------------------------------------------
// 6a) scalar case:
// ----------------
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::seq_locate (
  const point_basic<T>& x,
  size_type             dis_ie_guest) const
{
  return _locator.seq_locate (*this, x, dis_ie_guest);
}
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::dis_locate (
  const point_basic<T>& x,
  size_type             dis_ie_guest) const
{
  return _locator.dis_locate (*this, x, dis_ie_guest);
}
// ----------------
// 6b) disarray case:
// ----------------
template <class T>
void
geo_rep<T,sequential>::locate (
  const disarray<point_basic<T>, sequential>& x,
        disarray<size_type, sequential>&      dis_ie,
  bool do_check) const
{
  dis_ie.resize (x.ownership());
  for (size_type i = 0, n = x.size(); i < n; i++) {
    dis_ie[i] = base::_locator.seq_locate (*this, x[i], dis_ie[i]);
    if (do_check) {
      check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
        "locate: failed at x="<<ptos(x[i],base::_dimension));
    }
  }
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::locate (
  const disarray<point_basic<T>, distributed>& x,
        disarray<size_type, distributed>&      dis_ie,
        // dis_ie(ix) can contains a guest for x in K=omega[dis_ie(ix)]
  bool do_check) const
{
  // mpi::reduce with std::minimum handles unsigned=size_t as signed=long
  // so use here long to reduce, otherwise leads to a bug
  using signed_size_type = long;
  const signed_size_type signed_large = -1;
  const        size_type        large = std::numeric_limits<size_type>::max();
  const T                       infty = std::numeric_limits<T>::max();
  // 1) scan x and locate into the local mesh partition ; list when failed
  distributor ownership = x.ownership();
  if (dis_ie.ownership() != ownership) dis_ie.resize (ownership);
  size_type first_dis_i = ownership.first_index();
  std::list<id_pt_t<T> >  failed;
  for (size_type i = 0, n = x.size(); i < n; i++) {
    dis_ie[i] = base::_locator.seq_locate (*this, x[i], dis_ie[i]);
    if (dis_ie[i] == large) {
      size_type dis_i = first_dis_i + i;
      failed.push_back (id_pt_t<T>(dis_i, x[i]));
    }
  }
  // 2) merge the failed list into a massive disarray, then distributed to all
  communicator comm = ownership.comm();
  distributor fld_ownership (distributor::decide, comm, failed.size());
  size_type   fld_dis_size = fld_ownership.dis_size();
  if (comm.size() == 1 || fld_dis_size == 0) { 
    // no unsolved, on any procs: nothing to do !
trace_macro ("locate done(1)");
    return;
  }
  size_type first_fld_dis_i    = fld_ownership.first_index();
  size_type  last_fld_dis_i    = fld_ownership. last_index();
  id_pt_t<T> unset (large, point_basic<T>(infty,infty,infty));
  std::vector<id_pt_t<T> > massive_failed (fld_dis_size, unset);
  typename std::list<id_pt_t<T> >::iterator iter = failed.begin();
  for (size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
    massive_failed [fld_dis_i] = *iter;
  }
  std::vector<id_pt_t<T> > massive_query (fld_dis_size, unset);
  mpi::all_reduce (
        comm,
        massive_failed.begin().operator->(),
        massive_failed.size(),
        massive_query.begin().operator->(),
        id_pt_minimum<T>());
 
  // 3) run the locator on all failed points ON ALL PROCS, skipping local queries (already failed)
  std::vector<signed_size_type> massive_result (fld_dis_size, signed_large);
  // 3a) range [0:first[
  for (size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
    massive_result [fld_dis_i] = base::_locator.seq_locate (*this, massive_query[fld_dis_i].second);
  }
  // 3b) range [last,dis_size[
  for (size_type fld_dis_i = last_fld_dis_i; fld_dis_i < fld_dis_size; ++fld_dis_i) {
    massive_result [fld_dis_i] = base::_locator.seq_locate (*this, massive_query[fld_dis_i].second);
  }
  // 4) send & merge the results to all
  std::vector<signed_size_type> massive_merged (fld_dis_size, signed_large);
  mpi::all_reduce (
        comm,
        massive_result.begin().operator->(),
        massive_result.size(),
        massive_merged.begin().operator->(),
        mpi::maximum<signed_size_type>());
 
  // 5) store the local range into the distributed disarray dis_ie:
  for (size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
    size_type dis_i  = massive_query [fld_dis_i].first;
    check_macro (dis_i >= first_dis_i, "invalid index");
    size_type i = dis_i - first_dis_i;
    dis_ie[i] = massive_merged [fld_dis_i];
    if (do_check) {
      check_macro (dis_ie[i] != large,
        "dis_locate: failed at x="<<ptos(x[i],base::_dimension));
    }
  }
trace_macro ("locate done(2)");
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// no CGAL: no locators
// ----------------------------------------------------------------------------
#else // _RHEOLEF_HAVE_CGAL
template <class T, class M>
geo_locate<T,M>::~geo_locate()
{
}
template <class T, class M>
typename geo_locate<T,M>::size_type
geo_locate<T,M>::seq_locate (
    const geo_base_rep<T,M>& omega,
    const point_basic<T>&    x,
    size_type                dis_ie_guest) const
{
  fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
  return 0;
}
template <class T, class M>
typename geo_locate<T,M>::size_type
geo_locate<T,M>::dis_locate (
    const geo_base_rep<T,M>& omega,
    const point_basic<T>&    x,
    size_type                dis_ie_guest) const
{
  fatal_macro ("geo: locator 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_locate (
  const point_basic<T>& x,
  size_type             dis_ie_guest) const
{
  fatal_macro ("geo: locator 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_locate (
  const point_basic<T>& x,
  size_type             dis_ie_guest) const
{
  fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
  return 0;
}
template <class T>
void
geo_rep<T,sequential>::locate (
  const disarray<point_basic<T>, sequential>& x,
        disarray<size_type, sequential>&      dis_ie,
  bool do_check) const
{
  fatal_macro ("geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::locate (
  const disarray<point_basic<T>, distributed>& x,
        disarray<size_type, distributed>&      dis_ie,
        // dis_ie(ix) can contains a guest for x in K=omega[dis_ie(ix)]
  bool do_check) const
{
  fatal_macro ("geo: locator 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_locate<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