Loading...
Searching...
No Matches
geo_trace_ray_boundary.cc
Go to the documentation of this file.
//
// given x and v, search S on boudary mesh such that ray(x,v) hits S
// and returns the closest hit to x
//
// author: Pierre.Saramito@imag.fr
//
// date: 12 march 2012
//
// implementation note:
// use CGAL::AABB_Tree for 3D surfacic element intersections
//
#include "rheolef/geo_trace_ray_boundary.h"
#include "rheolef/geo.h"
// internal includes:
#include "rheolef/cgal_traits.h"
#include "rheolef/point_util.h"
// ------------------------------------------
// assume that CGAL is configured for GNU C++
// intel C++ has less 2011 standard features
// thus re-config it!
// ------------------------------------------
#ifdef _RHEOLEF_HAVE_CGAL
#ifdef _RHEOLEF_HAVE_INTEL_CXX
#include <CGAL/config.h>
#ifndef CGAL_CFG_NO_CPP0X_TUPLE
#define CGAL_CFG_NO_CPP0X_TUPLE
#endif
#ifndef CGAL_CFG_NO_TR1_TUPLE
#define CGAL_CFG_NO_TR1_TUPLE
#endif
#ifndef CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
#define CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
#endif
#endif // _RHEOLEF_HAVE_INTEL_CXX
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>
#endif // _RHEOLEF_HAVE_CGAL
namespace rheolef {
#ifdef _RHEOLEF_HAVE_CGAL
// ---------------------------------------------------------------------
// 0) utility: nearest hit over processors
// ---------------------------------------------------------------------
// formaly, we want to do in distributed ernvironment:
// [hit,y] = mpi::all_reduce (omega.comm(), [hit,y], mpi::minimum<pair<bool,point> >());
// where mpi::minimum<pair<bool,point> >(a,b) acts as the nearest point to x
//
template <class T>
struct hit_t : std::pair<bool,point_basic<T> > {
typedef std::pair<bool,point_basic<T> > base;
hit_t () : base() {}
hit_t (bool hit, const point_basic<T>& y) : base(hit,y) {}
template<class Archive>
void serialize (Archive& ar, const unsigned int version) {
ar & base::first;
ar & base::second;
}
};
} // namespace rheolef
#ifdef _RHEOLEF_HAVE_MPI
// Some serializable types have a fixed amount of data stored at fixed field positions.
// When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
// To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
namespace boost {
namespace mpi {
// TODO: when hit_t<T> is not a simple type, such as T=bigfloat or T=gmp, etc
template <>
struct is_mpi_datatype<rheolef::hit_t<double> > : mpl::true_ { };
} // namespace mpi
} // namespace boost
#endif // _RHEOLEF_HAVE_MPI
namespace rheolef {
template <class T>
struct nearest_hit : public std::binary_function<hit_t<T>, hit_t<T>, hit_t<T> > {
nearest_hit (const point_basic<T>& x = point_basic<T>()) : _x(x) {}
point_basic<T> _x;
const hit_t<T>& operator() (const hit_t<T>& a, const hit_t<T>& b) {
if (! b.first) return a;
if (! a.first) return b;
// here: a & b have a valid ray intersecion
T da = dist2(a.second, _x);
T db = dist2(b.second, _x);
if (da < db) return a;
if (db < da) return b;
// here: a & b have an equi-distant ray intersecion:
// => it is the same intersection, at an inter-element boundary
return a;
}
};
// compute over processors the nearest hit y to x
// (hit,y) are local to proc as input and global over procs as output
template <class T>
static
void
dis_compute_nearest (
const communicator& comm,
const point_basic<T>& x,
bool& hit,
point_basic<T>& y)
{
hit_t<T> hit_y (hit, y);
#ifdef _RHEOLEF_HAVE_MPI
hit_y = mpi::all_reduce (comm, hit_y, nearest_hit<T>(x));
#endif // _RHEOLEF_HAVE_MPI
hit = hit_y.first;
y = hit_y.second;
}
// ---------------------------------------------------------------------
// 1) the ray tracer interface
// ---------------------------------------------------------------------
template <class T, class M>
geo_trace_ray_boundary<T,M>::~geo_trace_ray_boundary()
{
if (_ptr != 0) {
delete_macro(_ptr);
}
}
template <class T, class M>
bool
geo_trace_ray_boundary<T,M>::seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
if (_ptr == 0) { _ptr = make_ptr(omega); }
return _ptr->seq_trace_ray_boundary (omega, x, v, y);
}
template <class T, class M>
bool
geo_trace_ray_boundary<T,M>::dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
if (_ptr == 0) { _ptr = make_ptr(omega); }
return _ptr->dis_trace_ray_boundary (omega, x, v, y);
}
// ---------------------------------------------------------------------
// 2) the ray tracer abtract data structure
// ---------------------------------------------------------------------
template <class T, class M>
class geo_trace_ray_boundary_abstract_rep {
public:
typedef typename disarray<T,M>::size_type size_type;
virtual ~geo_trace_ray_boundary_abstract_rep() {}
virtual void initialize (const geo_base_rep<T,M>& omega) const = 0;
virtual bool seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const = 0;
virtual bool dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const = 0;
};
// ---------------------------------------------------------------------
// 2) the ray tracer concrete, and dimension dependent, data structure
// = CGAL::AABB_Tree in 3D and empty otherwise
// ---------------------------------------------------------------------
template <class T, class M, size_t D>
struct geo_trace_ray_boundary_rep : public geo_trace_ray_boundary_abstract_rep<T,M> { };
// ---------------------------------------------------------------------
// 2a) the 1D implementation
// ---------------------------------------------------------------------
template <class T, class M>
class geo_trace_ray_boundary_rep<T,M,1> : public geo_trace_ray_boundary_abstract_rep<T,M> {
public:
// typedef:
typedef typename geo_base_rep<T,M>::size_type size_type;
// allocators:
geo_trace_ray_boundary_rep() {}
geo_trace_ray_boundary_rep(const geo_base_rep<T,M>& omega) { initialize(omega); }
~geo_trace_ray_boundary_rep() {}
void initialize (const geo_base_rep<T,M>& omega) const;
// accessors:
bool seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
bool dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
};
template <class T, class M>
void
geo_trace_ray_boundary_rep<T,M,1>::initialize (const geo_base_rep<T,M>& omega) const
{
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,1>::seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
T x0 = x[0];
T v0 = v[0];
T y0 = 0;
bool hit = false;
T d_min = std::numeric_limits<T>::max();
const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
check_macro (boundary.map_dimension() == 0, "unexpected boundary domain");
for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
const geo_element_indirect& oige = boundary.oige(ioige);
size_type ie = oige.index();
const geo_element& S = omega.get_geo_element(0,ie);
check_macro (S.variant() == reference_element::p, "unexpected element type: "<<S.name());
const point_basic<T>& a = omega.dis_node(S[0]);
T a0 = a[0];
if ((v0 > 0 && a0 > x0) || (v0 < 0 && a0 < x0) || (v0 == 0 && a0 == x0)) {
T di = fabs(x0-a0);
if (di < d_min) {
hit = true;
y0 = a0;
d_min = di;
}
}
}
y = point_basic<T>(y0);
return hit;
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,1>::dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
bool hit = seq_trace_ray_boundary (omega, x, v, y);
if (is_sequential<M>::value) return hit;
dis_compute_nearest (omega.comm(), x, hit, y);
return hit;
}
// ---------------------------------------------------------------------
// 2b) the 2D implementation
// ---------------------------------------------------------------------
template <class T, class M>
class geo_trace_ray_boundary_rep<T,M,2> : public geo_trace_ray_boundary_abstract_rep<T,M> {
public:
// typedef:
typedef typename geo_base_rep<T,M>::size_type size_type;
// allocators:
geo_trace_ray_boundary_rep() {}
geo_trace_ray_boundary_rep(const geo_base_rep<T,M>& omega) { initialize(omega); }
~geo_trace_ray_boundary_rep() {}
void initialize (const geo_base_rep<T,M>& omega) const;
// accessors:
bool seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
bool dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
};
template <class T, class M>
void
geo_trace_ray_boundary_rep<T,M,2>::initialize (const geo_base_rep<T,M>& omega) const
{
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,2>::seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<T> > Kernel;
typedef typename Kernel::Point_2 Point;
typedef typename Kernel::Ray_2 Ray;
typedef typename Kernel::Vector_2 Vector2d;
typedef typename Kernel::Segment_2 Segment;
Point x1 (x[0], x[1]);
Vector2d v1 (v[0], v[1]);
Ray ray_query(x1,v1);
bool hit = false;
T d_min = std::numeric_limits<T>::max();
size_t n_intersect = 0;
const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
check_macro (boundary.map_dimension() == 1, "unexpected boundary domain");
for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
const geo_element_indirect& oige = boundary.oige(ioige);
size_type ie = oige.index();
const geo_element& S = omega.get_geo_element(1,ie);
check_macro (S.variant() == reference_element::e, "unexpected element type: "<<S.name());
const point_basic<T>& a = omega.dis_node(S[0]);
const point_basic<T>& b = omega.dis_node(S[1]);
Point a1 (a[0], a[1]);
Point b1 (b[0], b[1]);
Segment s (a1, b1);
CGAL::Object obj = intersection (s, ray_query);
const Point* ptr_xo = 0;
const Segment* ptr_so = 0;
if ((ptr_xo = CGAL::object_cast<Point> (&obj))) {
n_intersect++;
const Point& xo = *ptr_xo;
point_basic<T> y0 (xo.x(), xo.y());
T d0 = dist(y0,x);
if (d0 < d_min) {
hit = true;
y = y0;
d_min = d0;
}
} else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
n_intersect += 2;
const Segment& so = *ptr_so;
point_basic<T> y0 (so[0].x(), so[0].y());
T d0 = dist(y0,x);
if (d0 < d_min) {
hit = true;
y = y0;
d_min = d0;
}
point_basic<T> y1 (so[1].x(), so[1].y());
T d1 = dist(y1,x);
if (d1 < d_min) {
hit = true;
y = y1;
d_min = d1;
}
} else {
// no intersection
}
}
return hit;
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,2>::dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
bool hit = seq_trace_ray_boundary (omega, x, v, y);
if (is_sequential<M>::value) return hit;
dis_compute_nearest (omega.comm(), x, hit, y);
return hit;
}
// ---------------------------------------------------------------------
// 2c) the 3D implementation
// ---------------------------------------------------------------------
template <class T, class M>
class geo_trace_ray_boundary_rep<T,M,3> : public geo_trace_ray_boundary_abstract_rep<T,M> {
public:
// typedef:
typedef typename geo_base_rep<T,M>::size_type size_type;
typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<T> > Kernel;
// typedef CGAL::Filtered_kernel_adaptor<mycgal::my_kernel_3d<T> > Kernel; // missing some 3d features
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Ray_3 Ray;
typedef typename Kernel::Vector_3 Vector3d;
typedef typename Kernel::Segment_3 Segment;
typedef typename Kernel::Triangle_3 Triangle;
typedef typename std::list<Triangle>::iterator Iterator;
typedef CGAL::AABB_triangle_primitive<Kernel,Iterator> Primitive;
typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
typedef typename Tree::Object_and_primitive_id Object_and_primitive_id;
typedef typename Tree::Primitive_id Primitive_id;
// allocators:
geo_trace_ray_boundary_rep() : _tree() {}
geo_trace_ray_boundary_rep(const geo_base_rep<T,M>& omega) : _tree() { initialize(omega); }
~geo_trace_ray_boundary_rep() {}
void initialize (const geo_base_rep<T,M>& omega) const;
// accessors:
bool seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
bool dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const;
// data:
protected:
mutable Tree _tree; // for cgal ray intersection query
};
template <class T, class M>
void
geo_trace_ray_boundary_rep<T,M,3>::initialize (const geo_base_rep<T,M>& omega) const
{
// create the corresponding tree of bounding box
trace_macro ("create the 3d CGAL::AABB_tree...");
std::list<Triangle> triangles;
const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
check_macro (boundary.map_dimension() == 2, "unexpected boundary domain");
for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
const geo_element_indirect& oige = boundary.oige(ioige);
size_type ie = oige.index();
const geo_element& S = omega.get_geo_element(2,ie);
switch (S.variant()) {
case reference_element::t: {
const point_basic<T>& a = omega.dis_node(S[0]);
const point_basic<T>& b = omega.dis_node(S[1]);
const point_basic<T>& c = omega.dis_node(S[2]);
Point a1 (a[0], a[1], a[2]);
Point b1 (b[0], b[1], b[2]);
Point c1 (c[0], c[1], c[2]);
triangles.push_back(Triangle(a1,b1,c1));
break;
}
case reference_element::q: {
const point_basic<T>& a = omega.dis_node(S[0]);
const point_basic<T>& b = omega.dis_node(S[1]);
const point_basic<T>& c = omega.dis_node(S[2]);
const point_basic<T>& d = omega.dis_node(S[3]);
Point a1 (a[0], a[1], a[2]);
Point b1 (b[0], b[1], b[2]);
Point c1 (c[0], c[1], c[2]);
Point d1 (d[0], d[1], d[2]);
triangles.push_back(Triangle(a1,b1,c1));
triangles.push_back(Triangle(a1,c1,d1));
break;
}
default: error_macro ("unexpected element type: "<<S.name());
}
}
// constructs the AABB tree
_tree.rebuild (triangles.begin(), triangles.end());
trace_macro ("create the 3d CGAL::AABB_tree done");
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,3>::seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
Point x1 (x[0], x[1], x[2]);
Vector3d v1 (v[0], v[1], v[2]);
Ray ray_query (x1,v1);
// computes all intersections with segment query (as pairs object - primitive_id)
std::list<Object_and_primitive_id> intersections;
_tree.all_intersections (ray_query, std::back_inserter(intersections));
T d_min = std::numeric_limits<T>::max();
bool hit = false;
for (typename std::list<Object_and_primitive_id>::iterator i = intersections.begin(), last = intersections.end(); i != last; i++) {
CGAL::Object obj = (*i).first;
Iterator id = (*i).second;
const Point* ptr_yo = 0;
const Segment* ptr_so = 0;
if ((ptr_yo = CGAL::object_cast<Point> (&obj))) {
const Point& yo = *ptr_yo;
point_basic<T> yi = point_basic<T>(yo.x(), yo.y(), yo.z());
T di = dist(yi,x);
if (di < d_min) {
hit = true;
y = yi;
d_min = di;
}
} else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
const Segment& so = *ptr_so;
point_basic<T> y0 (so[0].x(), so[0].y());
T d0 = dist(y0,x);
if (d0 < d_min) {
hit = true;
y = y0;
d_min = d0;
}
point_basic<T> y1 (so[1].x(), so[1].y());
T d1 = dist(y1,x);
if (d1 < d_min) {
hit = true;
y = y1;
d_min = d1;
}
} else {
error_macro ("unexpected intersection type: " << typeid_name(obj.type().name(), false));
}
}
return hit;
}
template <class T, class M>
bool
geo_trace_ray_boundary_rep<T,M,3>::dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
bool hit = seq_trace_ray_boundary (omega, x, v, y);
if (is_sequential<M>::value) return hit;
dis_compute_nearest (omega.comm(), x, hit, y);
return hit;
}
// ---------------------------------------------------------------------
// 3) the pointer allocator
// ---------------------------------------------------------------------
template <class T, class M>
geo_trace_ray_boundary_abstract_rep<T,M>*
geo_trace_ray_boundary<T,M>::make_ptr (const geo_base_rep<T,M>& omega)
{
check_macro (omega.dimension() == omega.map_dimension(), "geo_trace_ray_boundary: map_dim < dim: not supported");
switch (omega.dimension()) {
case 1: return new_macro((geo_trace_ray_boundary_rep<T,M,1>)(omega));
case 2: return new_macro((geo_trace_ray_boundary_rep<T,M,2>)(omega));
case 3: return new_macro((geo_trace_ray_boundary_rep<T,M,3>)(omega));
default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
}
}
// ---------------------------------------------------------------------
// 4) disarray interface (more efficient in the distributed case)
// ---------------------------------------------------------------------
template <class T>
void
geo_rep<T,sequential>::trace_ray_boundary (
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,
bool do_check) const
{
for (size_type i = 0, n = x.size(); i < n; i++) {
bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*this, x[i], v[i], y[i]);
if (hit) {
dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
} else {
dis_ie[i] = std::numeric_limits<size_type>::max();
}
if (do_check) {
check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
"trace_ray_boundary: failed at x="<<ptos(x[i],base::_dimension)<<",v="<<ptos(v[i],base::_dimension));
}
}
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::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) const
{
trace_macro ("ray_boundary...");
const size_type large = std::numeric_limits<size_type>::max();
const T infty = std::numeric_limits<T>::max();
// -----------------------------------------------------------------
// 1) first try a ray(x[i],v[i]) on the local boundary partition
// -----------------------------------------------------------------
bool is_seq = (x.comm().size() == 1);
std::list<size_type> failed;
for (size_type i = 0, n = x.size(); i < n; i++) {
bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*this, x[i], v[i], y[i]);
if (hit) {
// here x[i)+v[i] cross the local boundary partition
dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
// y was outside Omega ? most of the time, y is exactly on the boundary:
// dirty hack: project y on the bbox of the boundary:
// - this works in 1d when Omega is simply connected
// - this works in 2d & 3d when Omega is a rectangle or a paralellotope...
// TODO: boundary projection in the multi-connected case with: y = omega.nearest(x+v)
// T eps = 1e3*std::numeric_limits<T>::epsilon();
T eps = sqrt(std::numeric_limits<T>::epsilon());
for (size_t k = 0; k < base::_dimension; k++) {
y[i][k] = std::max (std::min (y[i][k]+eps, base::_xmax[k]-eps), base::_xmin[k]);
}
dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
}
std::cerr << std::setprecision (17);
check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
"trace_ray_boundary: locator failed at y="<<ptos(y[i],base::_dimension));
} else {
// here x[i)+v[i] cross the boundary, on a boundary part that is managed by another proc
failed.push_back (i);
}
}
// -----------------------------------------------------------------
// 3) for all x+v that goes outside and in another boundary partition
// then solve a distributed ray trace
// ---------------------------------------------------------------
distributor fld_ownership (distributor::decide, base::comm(), failed.size());
size_type fld_dis_size = fld_ownership.dis_size();
if (fld_dis_size == 0) {
// all are solved: nothing more to do
trace_macro ("ray_boundary done(1)");
return;
}
size_type first_fld_dis_i = fld_ownership.first_index();
size_type last_fld_dis_i = fld_ownership. last_index();
pt2_t<T> unset2 (point_basic<T>(infty,infty,infty), point_basic<T>(infty,infty,infty));
std::vector<pt2_t<T> > massive_failed (fld_dis_size, unset2);
typename std::list<size_type>::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) {
size_type i = *iter;
massive_failed [fld_dis_i] = pt2_t<T>(x[i],v[i]);
}
std::vector<pt2_t<T> > massive_query (fld_dis_size, unset2);
mpi::all_reduce (
fld_ownership.comm(),
massive_failed.begin().operator->(),
massive_failed.size(),
massive_query.begin().operator->(),
pt2_minimum<T>());
// 3) run the locator on all failed points ON ALL PROCS, skipping local queries (already failed)
id_pt_t<T> unset (large, point_basic<T>(infty,infty,infty));
std::vector<id_pt_t<T> > massive_result (fld_dis_size, unset);
// 3a) range [0:first[
for (size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*this,
massive_query [fld_dis_i].first,
massive_query [fld_dis_i].second,
massive_result[fld_dis_i].second);
if (hit) {
massive_result [fld_dis_i].first = base::_locator.seq_locate (*this, massive_result[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) {
bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*this,
massive_query [fld_dis_i].first,
massive_query [fld_dis_i].second,
massive_result[fld_dis_i].second);
if (hit) {
massive_result [fld_dis_i].first = base::_locator.seq_locate (*this, massive_result[fld_dis_i].second);
}
}
// 4) send & merge the results to all
std::vector<id_pt_t<T> > massive_merged (fld_dis_size, unset);
mpi::all_reduce (
fld_ownership.comm(),
massive_result.begin().operator->(),
massive_result.size(),
massive_merged.begin().operator->(),
id_pt_minimum<T>());
// 5) store the local range into the distributed disarray dis_ie:
iter = failed.begin();
for (size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
size_type i = *iter;
dis_ie[i] = massive_merged [fld_dis_i].first;
y[i] = massive_merged [fld_dis_i].second;
if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
// then x belongs to the boundary and v.n > 0 on the boundary
// => the ray(x,v) goes outside
// => takes y=x
y[i] = x[i];
dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
}
if (do_check) {
check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
"trace_ray_boundary: failed at x="
<< std::setprecision(17)
<<ptos(x[i],base::_dimension)<<", v="
<<ptos(v[i],base::_dimension)<<", y="
<<ptos(y[i],base::_dimension));
}
}
trace_macro ("ray_boundary done(2)");
}
#endif // _RHEOLEF_HAVE_MPI
// ----------------------------------------------------------------------------
// CGAL lib is missing
// ----------------------------------------------------------------------------
#else // _RHEOLEF_HAVE_CGAL
template <class T, class M>
geo_trace_ray_boundary<T,M>::~geo_trace_ray_boundary()
{
}
template <class T, class M>
bool
geo_trace_ray_boundary<T,M>::seq_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
return false;
}
template <class T, class M>
bool
geo_trace_ray_boundary<T,M>::dis_trace_ray_boundary (
const geo_base_rep<T,M>& omega,
const point_basic<T>& x,
const point_basic<T>& v,
point_basic<T>& y) const
{
fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
return false;
}
template <class T>
void
geo_rep<T,sequential>::trace_ray_boundary (
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,
bool do_check) const
{
fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
}
#ifdef _RHEOLEF_HAVE_MPI
template <class T>
void
geo_rep<T,distributed>::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) const
{
fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
}
#endif // _RHEOLEF_HAVE_MPI
#endif // _RHEOLEF_HAVE_CGAL
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template class geo_rep<Float,sequential>;
template class geo_trace_ray_boundary<Float,sequential>;
#ifdef _RHEOLEF_HAVE_MPI
template class geo_rep<Float,distributed>;
template class geo_trace_ray_boundary<Float,distributed>;
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
rheolef::disarray
see the disarray page for the full documentation
Definition disarray.h:497
rheolef::disarray::size_type
rep::base::size_type size_type
Definition disarray.h:501
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::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition distributor.h:158
rheolef::distributor::decide
static const size_type decide
Definition distributor.h:83
rheolef::distributor::comm
const communicator_type & comm() const
Definition distributor.h:152
rheolef::domain_indirect_basic
the finite element boundary domain
Definition domain_indirect.h:335
rheolef::geo_base_rep
base class for M=sequential or distributed meshes representations
Definition geo.h:528
rheolef::geo_base_rep< T, sequential >::_dimension
size_type _dimension
Definition geo.h:685
rheolef::geo_base_rep< T, sequential >::_locator
geo_locate< T, sequential > _locator
Definition geo.h:692
rheolef::geo_base_rep< T, sequential >::_tracer_ray_boundary
geo_trace_ray_boundary< T, sequential > _tracer_ray_boundary
Definition geo.h:693
rheolef::geo_base_rep::map_dimension
size_type map_dimension() const
Definition geo.h:564
rheolef::geo_base_rep::size_type
base::size_type size_type
Definition geo.h:533
rheolef::geo_base_rep< T, distributed >::_xmax
node_type _xmax
Definition geo.h:688
rheolef::geo_base_rep::dimension
size_type dimension() const
Definition geo.h:563
rheolef::geo_base_rep< T, distributed >::_xmin
node_type _xmin
Definition geo.h:687
rheolef::geo_base_rep< T, distributed >::comm
const communicator & comm() const
Definition geo.h:637
rheolef::geo_element_indirect
Definition geo_element_indirect.h:32
rheolef::geo_element
see the geo_element page for the full documentation
Definition geo_element.h:102
rheolef::geo_rep< T, distributed >::geo_rep
geo_rep()
Definition geo.cc:112
rheolef::geo_rep< T, distributed >::size
size_type size(size_type dim) const
Definition geo.h:998
rheolef::geo_rep< T, distributed >::size_type
base::size_type size_type
Definition geo.h:934
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::geo_trace_ray_boundary
Definition geo_trace_ray_boundary.h:41
rheolef::geo_trace_ray_boundary::_ptr
geo_trace_ray_boundary_abstract_rep< T, M > * _ptr
Definition geo_trace_ray_boundary.h:61
rheolef::geo_trace_ray_boundary::seq_trace_ray_boundary
bool seq_trace_ray_boundary(const geo_base_rep< T, M > &omega, const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
Definition geo_trace_ray_boundary.cc:151
rheolef::geo_trace_ray_boundary::~geo_trace_ray_boundary
~geo_trace_ray_boundary()
Definition geo_trace_ray_boundary.cc:143
rheolef::geo_trace_ray_boundary::dis_trace_ray_boundary
bool dis_trace_ray_boundary(const geo_base_rep< T, M > &omega, const point_basic< T > &x, const point_basic< T > &v, point_basic< T > &y) const
Definition geo_trace_ray_boundary.cc:162
rheolef::geo_trace_ray_boundary::make_ptr
static geo_trace_ray_boundary_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
Definition geo_trace_ray_boundary.cc:564
rheolef::point_basic
Definition point.h:86
rheolef::reference_element::q
static const variant_type q
Definition reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition reference_element.h:76
rheolef::reference_element::p
static const variant_type p
Definition reference_element.h:75
rheolef::reference_element::t
static const variant_type t
Definition reference_element.h:77
trace_macro
#define trace_macro(message)
Definition dis_macros.h:111
error_macro
#define error_macro(message)
Definition dis_macros.h:49
fatal_macro
#define fatal_macro(message)
Definition dis_macros.h:33
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
boost::mpi
Definition mpi_pair_datatype.h:52
boost
Definition mpi_pair_datatype.h:51
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition mkgeo_ball.sh:153
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
mkgeo_ugrid.boundary
int boundary
Definition mkgeo_ugrid.sh:181
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::typeid_name
std::string typeid_name(const char *name, bool do_indent)
Definition pretty_name.cc:78
rheolef::dist2
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
rheolef::ptos
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
rheolef::dist
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
d
size_t d
Definition neumann-laplace-lambda.cc:28
rheolef::distributed
Definition communicator.h:128
rheolef::Kernel
rheolef::id_pt_minimum
Definition point_util.h:67
rheolef::id_pt_t
Definition point_util.h:39
rheolef::is_sequential::value
static const bool value
Definition communicator.h:67
rheolef::pt2_minimum
Definition point_util.h:113
rheolef::pt2_t
Definition point_util.h:85
rheolef::sequential
Definition communicator.h:62
M
Expr1::memory_type M
Definition vec_expr_v2.h:416