Loading...
Searching...
No Matches
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
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::distributor
see the distributor page for the full documentation
Definition distributor.h:69
rheolef::distributor::dis_size
size_type dis_size() const
global and local sizes
Definition distributor.h:214
rheolef::distributor::decide
static const size_type decide
Definition distributor.h:83
rheolef::geo_base_rep
base class for M=sequential or distributed meshes representations
Definition geo.h:528
rheolef::geo_base_rep::_locator
geo_locate< T, M > _locator
Definition geo.h:692
rheolef::geo_base_rep::_tracer_ray_boundary
geo_trace_ray_boundary< T, M > _tracer_ray_boundary
Definition geo.h:693
rheolef::geo_base_rep::size_type
base::size_type size_type
Definition geo.h:533
rheolef::geo_base_rep::seq_trace_move
size_type seq_trace_move(const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
Definition geo_trace_move.cc:45
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_trace_move
size_type dis_trace_move(const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
Definition geo_trace_move.cc:83
rheolef::geo_base_rep< T, distributed >::comm
const communicator & comm() const
Definition geo.h:637
rheolef::geo_rep< T, distributed >::geo_rep
geo_rep()
Definition geo.cc:112
rheolef::geo_rep< T, distributed >::trace_ray_boundary
void trace_ray_boundary(const disarray< point_basic< T >, distributed > &x, const disarray< point_basic< T >, distributed > &v, disarray< size_type, distributed > &dis_ie, disarray< point_basic< T >, distributed > &y, bool do_check=false) const
Definition geo_trace_ray_boundary.cc:602
rheolef::geo_rep< T, distributed >::size_type
base::size_type size_type
Definition geo.h:934
rheolef::geo_rep< T, distributed >::locate
void locate(const disarray< point_basic< T >, distributed > &x, disarray< size_type, distributed > &dis_ie, bool do_check=false) const
Definition geo_locate.cc:441
rheolef::geo_rep< T, sequential >::geo_rep
geo_rep()
Definition geo.cc:85
rheolef::geo_rep< T, sequential >::size_type
base::size_type size_type
Definition geo.h:791
rheolef::geo_rep
sequential mesh representation
Definition geo.h:778
rheolef::point_basic
Definition point.h:86
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::distributed
Definition communicator.h:128
rheolef::sequential
Definition communicator.h:62