geo_trace_move.cc

Go to the documentation of this file.

//  y = trace_move(x,v);
//  given x and v, search K in mesh such that y=x+v in K
//  and when x+v goes outside, ray trace the y point on the boundary
//
// author: Pierre.Saramito@imag.fr
//
// date: 15 march 2012
//
// TODO: avoid the 2nd locate in trace_move ?
//  en utilisant une table de liens S -> K (face-element: meme procs)
//  => trace_ray renvoie un numero d'element tout de suite
// les liens S -> sont dans le maillage des la lecture : geo::get
// ils sont aussi utiles pour Galerkin-discontinu
//
#include "rheolef/geo.h"
 
namespace rheolef {
 
// --------------------------------------------------------------------------
// 1) one point x
// --------------------------------------------------------------------------
// TODO: add dis_ie_guest
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::seq_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
{
  y = x+v;
  size_type dis_ie = _locator.seq_locate (*this,y);
  if (dis_ie != std::numeric_limits<size_type>::max()) {
    // y = x+v falls inside Omega
    return dis_ie;
  } 
  // y = x+v falls outside Omega: compute y = ray(x,v) intersection with the Omega boundary 
  bool hit = _tracer_ray_boundary.seq_trace_ray_boundary (*this, x, v, y);
  if (hit) {
    // TODO: hit knows the face boundary: could get the element without locator
    dis_ie = _locator.seq_locate (*this,y);
    check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
    return dis_ie;
  }
  // x is on the boundary and also y (tangential velocity)
  dis_ie = _locator.seq_locate (*this,y);
  if (dis_ie != std::numeric_limits<size_type>::max()) {
    return dis_ie;
  } 
  // x was outside Omega ? most of the time, y is exactly on the boundary:
  // boundary projection in the multi-connected case with: y = omega.nearest(x+v)
  point_basic<T> y_nearest;
  dis_ie = seq_nearest (y, y_nearest);
  if (dis_ie != std::numeric_limits<size_type>::max()) {
    y = y_nearest;
    return dis_ie;
  } 
  // machine precision problem ?
  error_macro ("invalid seq_trace_move");
  return std::numeric_limits<size_type>::max();
}
template <class T, class M>
typename geo_base_rep<T,M>::size_type
geo_base_rep<T,M>::dis_trace_move (
                const point_basic<T>&     x,
                const point_basic<T>&     v,
                      point_basic<T>&     y) const
{
  y = x+v;
  size_type dis_ie = _locator.dis_locate (*this,y);
  if (dis_ie != std::numeric_limits<size_type>::max()) {
    // y = x+v falls inside Omega
    return dis_ie;
  } 
  // y = x+v falls outside Omega: compute y = ray(x,v) intersection with the Omega boundary 
  bool hit = _tracer_ray_boundary.dis_trace_ray_boundary (*this, x, v, y);
  if (hit) {
    // TODO: hit knows the face boundary: could get the element without locator
    dis_ie = _locator.dis_locate (*this,y);
    check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
    return dis_ie;
  }
  // x was outside Omega ?
  error_macro ("invalid dis_trace_move");
  return std::numeric_limits<size_type>::max();
}
// --------------------------------------------------------------------------
// 2) disarrays x & v : groups comms in the distributed case
// --------------------------------------------------------------------------
template <class T>
void
geo_rep<T,sequential>::trace_move (
                const disarray<point_basic<T>,sequential>&     x,
                const disarray<point_basic<T>,sequential>&     v,
                      disarray<size_type, sequential>&         dis_ie,
                      disarray<point_basic<T>,sequential>&     y) const
{
  check_macro (x.ownership() == v.ownership(), "invalid ranges for x and v argumenst");
  dis_ie.resize (x.ownership());
  y.resize      (x.ownership());
  for (size_type i = 0, n = x.size(); i < n; i++) {
    // TODO: use dis_ie[i] as a guest
    dis_ie[i] = base::seq_trace_move (x[i], v[i], y[i]);
  }
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::trace_move (
                const disarray<point_basic<T>,distributed>&     x,
                const disarray<point_basic<T>,distributed>&     v,
                      disarray<size_type, distributed>&         dis_ie, // in: guest; out: localized
                      disarray<point_basic<T>,distributed>&     y) const
{
  trace_macro("trace_move...");
  check_macro (x.ownership() == v.ownership(), "invalid ranges for x and v argumenst");
  if (dis_ie.ownership() != x.ownership()) dis_ie.resize (x.ownership());
  y.resize      (x.ownership());
  // ---------------------------------------------------------------
  // 1) locate: solve all x(i)+y(i) that remains inside the domain: 
  // ---------------------------------------------------------------
  for (size_type i = 0, n = x.size(); i < n; i++) {
    y[i] = x[i] + v[i];
  }
  bool do_stop_when_failed = false;
  trace_macro("trace_move/locate...");
  locate (y, dis_ie, do_stop_when_failed);
  // -----------------------------------------------------------------
  // 2) for all x+v that goes outside
  // compact all failed points and run trace_ray_boundary (fld_x,fld_v)
  // ---------------------------------------------------------------
  std::list<size_type>  failed;
  for (size_type i = 0, n = x.size(); i < n; i++) {
    if (dis_ie[i] != std::numeric_limits<size_type>::max()) continue;
    // here x[i)+v[i] cross the boundary:
    failed.push_back (i);
  }
  distributor fld_ownership (distributor::decide, base::comm(), failed.size());
  if (fld_ownership.dis_size() == 0) {
    // all are solved: nothing more to do
    trace_macro("trace_move done(1)");
    return;
  }
  disarray<point_basic<T> > fld_x (fld_ownership);
  disarray<point_basic<T> > fld_v (fld_ownership);
  disarray<point_basic<T> > fld_y (fld_ownership);
  disarray<size_type>  fld_dis_ie (fld_ownership, std::numeric_limits<size_type>::max());
  typename std::list<size_type>::const_iterator iter = failed.begin();
  for (size_type fld_i = 0, fld_n = fld_x.size(); fld_i < fld_n; ++fld_i, ++iter) {
    size_type i = *iter;
    fld_x [fld_i] = x[i];
    fld_v [fld_i] = v[i];
  }
  do_stop_when_failed = true;
  trace_macro("trace_move/ray_boundary..."); 
  trace_ray_boundary (fld_x, fld_v, fld_dis_ie, fld_y, do_stop_when_failed);
  trace_macro("trace_move/ray_boundary = " 
        << fld_x.dis_size() << "/" << x.dis_size() << " = "
        << 1.0*fld_x.dis_size()/x.dis_size());
  // -----------------------------------------------------------------
  // 3) unpack the result
  // ---------------------------------------------------------------
  iter = failed.begin();
  for (size_type fld_i = 0, fld_n = fld_x.size(); fld_i < fld_n; ++fld_i, ++iter) {
    size_type i = *iter;
    dis_ie[i] = fld_dis_ie [fld_i];
         y[i] = fld_y      [fld_i];
  }
  trace_macro("trace_move done(2)"); 
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_base_rep<Float,sequential>;
template class geo_rep<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class geo_base_rep<Float,distributed>;
template class geo_rep<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
 
} // namespace rheolef