geo_trace_move.cc

Go to the documentation of this file.

//  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
#include "rheolef/geo.h"

namespace rheolef {

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()) {
    return dis_ie;
  } 
  bool hit = _tracer_ray_boundary.seq_trace_ray_boundary (*this, x, v, y);
  if (hit) {
    dis_ie = _locator.seq_locate (*this,y);
    check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
    return dis_ie;
  }
  for (size_t i = 0; i < _dimension; i++) {
    y[i] = (v[i] < 0) ? _xmin[i] : _xmax[i];
  }
  dis_ie = _locator.seq_locate (*this,y);
  if (dis_ie != std::numeric_limits<size_type>::max()) {
    return dis_ie;
  } 
  error_macro ("invalid dis_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()) {
    return dis_ie;
  } 
  bool hit = _tracer_ray_boundary.dis_trace_ray_boundary (*this, x, v, y);
  if (hit) {
    dis_ie = _locator.dis_locate (*this,y);
    check_macro (dis_ie != std::numeric_limits<size_type>::max(), "invalid ray computation");
    return dis_ie;
  }
  error_macro ("invalid dis_trace_move");
  return std::numeric_limits<size_type>::max();
}
template <class T>
void
geo_rep<T,sequential>::trace_move (
                const array<point_basic<T>,sequential>&     x,
                const array<point_basic<T>,sequential>&     v,
                      array<size_type, sequential>&         dis_ie,
                      array<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++) {
    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 array<point_basic<T>,distributed>&     x,
                const array<point_basic<T>,distributed>&     v,
                      array<size_type, distributed>&         dis_ie, // in: guest; out: localized
                      array<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());
  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);
  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;
    failed.push_back (i);
  }
  distributor fld_ownership (distributor::decide, base::comm(), failed.size());
  if (fld_ownership.dis_size() == 0) {
    trace_macro("trace_move done(1)");
    return;
  }
  array<point_basic<T> > fld_x (fld_ownership);
  array<point_basic<T> > fld_v (fld_ownership);
  array<point_basic<T> > fld_y (fld_ownership);
  array<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());
  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];
  }
#ifdef TO_CLEAN
  for (size_type i = 0, n = x.size(); i < n; i++) {
    check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),    
        "trace_move: localization failed");
  }
#endif // TO_CLEAN
  trace_macro("trace_move done(2)"); 
}
#endif // _RHEOLEF_HAVE_MPI
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