32 #include "rheolef/geo_trace_ray_boundary.h"
33 #include "rheolef/geo.h"
36 #include "rheolef/cgal_traits.h"
37 #include "rheolef/point_util.h"
44 #ifdef _RHEOLEF_HAVE_CGAL
45 #ifdef _RHEOLEF_HAVE_INTEL_CXX
46 #include <CGAL/config.h>
48 #ifndef CGAL_CFG_NO_CPP0X_TUPLE
49 #define CGAL_CFG_NO_CPP0X_TUPLE
52 #ifndef CGAL_CFG_NO_TR1_TUPLE
53 #define CGAL_CFG_NO_TR1_TUPLE
56 #ifndef CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
57 #define CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
61 #include <CGAL/AABB_tree.h>
62 #include <CGAL/AABB_traits.h>
63 #include <CGAL/AABB_triangle_primitive.h>
68 #ifdef _RHEOLEF_HAVE_CGAL
77 struct hit_t : std::pair<bool,point_basic<T> > {
78 typedef std::pair<bool,point_basic<T> > base;
81 template<
class Archive>
89 #ifdef _RHEOLEF_HAVE_MPI
97 struct is_mpi_datatype<
rheolef::hit_t<double> > : mpl::true_ { };
105 struct nearest_hit :
public std::binary_function<hit_t<T>, hit_t<T>, hit_t<T> > {
108 const hit_t<T>&
operator() (
const hit_t<T>&
a,
const hit_t<T>&
b) {
109 if (!
b.first)
return a;
110 if (!
a.first)
return b;
114 if (da < db)
return a;
115 if (db < da)
return b;
126 dis_compute_nearest (
132 hit_t<T> hit_y (hit, y);
133 #ifdef _RHEOLEF_HAVE_MPI
134 hit_y = mpi::all_reduce (comm, hit_y, nearest_hit<T>(x));
142 template <
class T,
class M>
149 template <
class T,
class M>
157 if (_ptr == 0) { _ptr = make_ptr(omega); }
158 return _ptr->seq_trace_ray_boundary (omega, x, v, y);
160 template <
class T,
class M>
168 if (_ptr == 0) { _ptr = make_ptr(omega); }
169 return _ptr->dis_trace_ray_boundary (omega, x, v, y);
174 template <
class T,
class M>
175 class geo_trace_ray_boundary_abstract_rep {
178 virtual ~geo_trace_ray_boundary_abstract_rep() {}
179 virtual void initialize (
const geo_base_rep<T,M>& omega)
const = 0;
180 virtual bool seq_trace_ray_boundary (
181 const geo_base_rep<T,M>& omega,
185 virtual bool dis_trace_ray_boundary (
186 const geo_base_rep<T,M>& omega,
195 template <
class T,
class M,
size_t D>
196 struct geo_trace_ray_boundary_rep :
public geo_trace_ray_boundary_abstract_rep<T,M> { };
201 template <
class T,
class M>
202 class geo_trace_ray_boundary_rep<
T,
M,1> :
public geo_trace_ray_boundary_abstract_rep<T,M> {
211 geo_trace_ray_boundary_rep() {}
212 geo_trace_ray_boundary_rep(
const geo_base_rep<T,M>& omega) {
initialize(omega); }
213 ~geo_trace_ray_boundary_rep() {}
214 void initialize (
const geo_base_rep<T,M>& omega)
const;
218 bool seq_trace_ray_boundary (
219 const geo_base_rep<T,M>& omega,
223 bool dis_trace_ray_boundary (
224 const geo_base_rep<T,M>& omega,
229 template <
class T,
class M>
234 template <
class T,
class M>
236 geo_trace_ray_boundary_rep<T,M,1>::seq_trace_ray_boundary (
237 const geo_base_rep<T,M>& omega,
246 T d_min = std::numeric_limits<T>::max();
247 const domain_indirect_basic<M>&
boundary = omega.get_domain_indirect (
"boundary");
249 for (
size_t ioige = 0, noige =
boundary.size(); ioige < noige; ioige++) {
250 const geo_element_indirect& oige =
boundary.oige(ioige);
252 const geo_element& S = omega.get_geo_element(0,ie);
256 if ((v0 > 0 && a0 > x0) || (v0 < 0 && a0 < x0) || (v0 == 0 && a0 == x0)) {
268 template <
class T,
class M>
270 geo_trace_ray_boundary_rep<T,M,1>::dis_trace_ray_boundary (
271 const geo_base_rep<T,M>& omega,
276 bool hit = seq_trace_ray_boundary (omega, x, v, y);
278 dis_compute_nearest (omega.comm(), x, hit, y);
284 template <
class T,
class M>
285 class geo_trace_ray_boundary_rep<
T,
M,2> :
public geo_trace_ray_boundary_abstract_rep<T,M> {
294 geo_trace_ray_boundary_rep() {}
295 geo_trace_ray_boundary_rep(
const geo_base_rep<T,M>& omega) {
initialize(omega); }
296 ~geo_trace_ray_boundary_rep() {}
297 void initialize (
const geo_base_rep<T,M>& omega)
const;
301 bool seq_trace_ray_boundary (
302 const geo_base_rep<T,M>& omega,
306 bool dis_trace_ray_boundary (
307 const geo_base_rep<T,M>& omega,
312 template <
class T,
class M>
317 template <
class T,
class M>
319 geo_trace_ray_boundary_rep<T,M,2>::seq_trace_ray_boundary (
320 const geo_base_rep<T,M>& omega,
325 typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<T> > Kernel;
326 typedef typename Kernel::Point_2 Point;
327 typedef typename Kernel::Ray_2 Ray;
328 typedef typename Kernel::Vector_2 Vector2d;
329 typedef typename Kernel::Segment_2 Segment;
331 Point x1 (x[0], x[1]);
332 Vector2d v1 (v[0], v[1]);
333 Ray ray_query(x1,v1);
336 T d_min = std::numeric_limits<T>::max();
337 size_t n_intersect = 0;
338 const domain_indirect_basic<M>&
boundary = omega.get_domain_indirect (
"boundary");
340 for (
size_t ioige = 0, noige =
boundary.size(); ioige < noige; ioige++) {
341 const geo_element_indirect& oige =
boundary.oige(ioige);
343 const geo_element& S = omega.get_geo_element(1,ie);
347 Point a1 (
a[0],
a[1]);
348 Point b1 (
b[0],
b[1]);
350 CGAL::Object obj = intersection (s, ray_query);
351 const Point* ptr_xo = 0;
352 const Segment* ptr_so = 0;
353 if ((ptr_xo = CGAL::object_cast<Point> (&obj))) {
355 const Point& xo = *ptr_xo;
363 }
else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
365 const Segment& so = *ptr_so;
386 template <
class T,
class M>
388 geo_trace_ray_boundary_rep<T,M,2>::dis_trace_ray_boundary (
389 const geo_base_rep<T,M>& omega,
394 bool hit = seq_trace_ray_boundary (omega, x, v, y);
396 dis_compute_nearest (omega.comm(), x, hit, y);
402 template <
class T,
class M>
403 class geo_trace_ray_boundary_rep<
T,
M,3> :
public geo_trace_ray_boundary_abstract_rep<T,M> {
410 typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<T> > Kernel;
413 typedef typename Kernel::Point_3 Point;
414 typedef typename Kernel::Ray_3 Ray;
415 typedef typename Kernel::Vector_3 Vector3d;
416 typedef typename Kernel::Segment_3 Segment;
417 typedef typename Kernel::Triangle_3 Triangle;
418 typedef typename std::list<Triangle>::iterator Iterator;
419 typedef CGAL::AABB_triangle_primitive<Kernel,Iterator> Primitive;
420 typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
421 typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
422 typedef typename Tree::Object_and_primitive_id Object_and_primitive_id;
423 typedef typename Tree::Primitive_id Primitive_id;
427 geo_trace_ray_boundary_rep() : _tree() {}
428 geo_trace_ray_boundary_rep(
const geo_base_rep<T,M>& omega) : _tree() {
initialize(omega); }
429 ~geo_trace_ray_boundary_rep() {}
430 void initialize (
const geo_base_rep<T,M>& omega)
const;
434 bool seq_trace_ray_boundary (
435 const geo_base_rep<T,M>& omega,
439 bool dis_trace_ray_boundary (
440 const geo_base_rep<T,M>& omega,
449 template <
class T,
class M>
455 std::list<Triangle> triangles;
456 const domain_indirect_basic<M>&
boundary = omega.get_domain_indirect (
"boundary");
458 for (
size_t ioige = 0, noige =
boundary.size(); ioige < noige; ioige++) {
459 const geo_element_indirect& oige =
boundary.oige(ioige);
461 const geo_element& S = omega.get_geo_element(2,ie);
462 switch (S.variant()) {
467 Point a1 (
a[0],
a[1],
a[2]);
468 Point b1 (
b[0],
b[1],
b[2]);
469 Point c1 (
c[0],
c[1],
c[2]);
470 triangles.push_back(Triangle(a1,b1,c1));
478 Point a1 (
a[0],
a[1],
a[2]);
479 Point b1 (
b[0],
b[1],
b[2]);
480 Point c1 (
c[0],
c[1],
c[2]);
481 Point d1 (
d[0],
d[1],
d[2]);
482 triangles.push_back(Triangle(a1,b1,c1));
483 triangles.push_back(Triangle(a1,c1,d1));
486 default:
error_macro (
"unexpected element type: "<<S.name());
490 _tree.rebuild (triangles.begin(), triangles.end());
491 trace_macro (
"create the 3d CGAL::AABB_tree done");
493 template <
class T,
class M>
495 geo_trace_ray_boundary_rep<T,M,3>::seq_trace_ray_boundary (
496 const geo_base_rep<T,M>& omega,
501 Point x1 (x[0], x[1], x[2]);
502 Vector3d v1 (v[0], v[1], v[2]);
503 Ray ray_query (x1,v1);
506 std::list<Object_and_primitive_id> intersections;
507 _tree.all_intersections (ray_query, std::back_inserter(intersections));
508 T d_min = std::numeric_limits<T>::max();
510 for (
typename std::list<Object_and_primitive_id>::iterator i = intersections.begin(), last = intersections.end(); i != last; i++) {
511 CGAL::Object obj = (*i).first;
512 Iterator
id = (*i).second;
513 const Point* ptr_yo = 0;
514 const Segment* ptr_so = 0;
515 if ((ptr_yo = CGAL::object_cast<Point> (&obj))) {
516 const Point& yo = *ptr_yo;
524 }
else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
525 const Segment& so = *ptr_so;
546 template <
class T,
class M>
548 geo_trace_ray_boundary_rep<T,M,3>::dis_trace_ray_boundary (
549 const geo_base_rep<T,M>& omega,
554 bool hit = seq_trace_ray_boundary (omega, x, v, y);
556 dis_compute_nearest (omega.comm(), x, hit, y);
562 template <
class T,
class M>
563 geo_trace_ray_boundary_abstract_rep<T,M>*
568 case 1:
return new_macro((geo_trace_ray_boundary_rep<T,M,1>)(omega));
569 case 2:
return new_macro((geo_trace_ray_boundary_rep<T,M,2>)(omega));
570 case 3:
return new_macro((geo_trace_ray_boundary_rep<T,M,3>)(omega));
587 bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*
this, x[i], v[i], y[i]);
589 dis_ie[i] = base::_locator.seq_locate (*
this, y[i]);
591 dis_ie[i] = std::numeric_limits<size_type>::max();
594 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
595 "trace_ray_boundary: failed at x="<<
ptos(x[i],base::_dimension)<<
",v="<<
ptos(v[i],base::_dimension));
599 #ifdef _RHEOLEF_HAVE_MPI
610 const size_type large = std::numeric_limits<size_type>::max();
611 const T infty = std::numeric_limits<T>::max();
615 bool is_seq = (x.comm().size() == 1);
616 std::list<size_type> failed;
618 bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*
this, x[i], v[i], y[i]);
621 dis_ie[i] = base::_locator.seq_locate (*
this, y[i]);
622 if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
630 for (
size_t k = 0; k < base::_dimension; k++) {
631 y[i][k] = std::max (std::min (y[i][k]+eps, base::_xmax[k]-eps), base::_xmin[k]);
633 dis_ie[i] = base::_locator.seq_locate (*
this, y[i]);
635 std::cerr << std::setprecision (17);
636 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
637 "trace_ray_boundary: locator failed at y="<<
ptos(y[i],base::_dimension));
640 failed.push_back (i);
649 if (fld_dis_size == 0) {
655 size_type last_fld_dis_i = fld_ownership. last_index();
657 std::vector<pt2_t<T> > massive_failed (fld_dis_size, unset2);
658 typename std::list<size_type>::iterator iter = failed.begin();
659 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
661 massive_failed [fld_dis_i] =
pt2_t<T>(x[i],v[i]);
663 std::vector<pt2_t<T> > massive_query (fld_dis_size, unset2);
665 fld_ownership.
comm(),
666 massive_failed.begin().operator->(),
667 massive_failed.size(),
668 massive_query.begin().operator->(),
673 std::vector<id_pt_t<T> > massive_result (fld_dis_size, unset);
675 for (
size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
676 bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*
this,
677 massive_query [fld_dis_i].first,
678 massive_query [fld_dis_i].second,
679 massive_result[fld_dis_i].second);
681 massive_result [fld_dis_i].first = base::_locator.seq_locate (*
this, massive_result[fld_dis_i].second);
685 for (
size_type fld_dis_i = last_fld_dis_i; fld_dis_i < fld_dis_size; ++fld_dis_i) {
686 bool hit = base::_tracer_ray_boundary.seq_trace_ray_boundary (*
this,
687 massive_query [fld_dis_i].first,
688 massive_query [fld_dis_i].second,
689 massive_result[fld_dis_i].second);
691 massive_result [fld_dis_i].first = base::_locator.seq_locate (*
this, massive_result[fld_dis_i].second);
695 std::vector<id_pt_t<T> > massive_merged (fld_dis_size, unset);
697 fld_ownership.
comm(),
698 massive_result.begin().operator->(),
699 massive_result.size(),
700 massive_merged.begin().operator->(),
704 iter = failed.begin();
705 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
707 dis_ie[i] = massive_merged [fld_dis_i].first;
708 y[i] = massive_merged [fld_dis_i].second;
710 if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
715 dis_ie[i] = base::_locator.seq_locate (*
this, y[i]);
718 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
719 "trace_ray_boundary: failed at x="
720 << std::setprecision(17)
721 <<
ptos(x[i],base::_dimension)<<
", v="
722 <<
ptos(v[i],base::_dimension)<<
", y="
723 <<
ptos(y[i],base::_dimension));
733 template <
class T,
class M>
737 template <
class T,
class M>
740 const geo_base_rep<T,M>& omega,
745 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
748 template <
class T,
class M>
751 const geo_base_rep<T,M>& omega,
756 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
764 disarray<size_type, sequential>& dis_ie,
768 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
770 #ifdef _RHEOLEF_HAVE_MPI
776 disarray<size_type, distributed>& dis_ie,
780 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
787 template class geo_rep<Float,sequential>;
788 template class geo_trace_ray_boundary<Float,sequential>;
790 #ifdef _RHEOLEF_HAVE_MPI
791 template class geo_rep<Float,distributed>;
792 template class geo_trace_ray_boundary<Float,distributed>;
field::size_type size_type
see the communicator page for the full documentation
see the disarray page for the full documentation
see the distributor page for the full documentation
size_type dis_size() const
global and local sizes
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
static const size_type decide
const communicator_type & comm() const
geo_element_hack::size_type size_type
base class for M=sequential or distributed meshes representations
size_type map_dimension() const
base::size_type size_type
size_type dimension() const
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
base::size_type size_type
void 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=false) const
sequential mesh representation
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
~geo_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
static geo_trace_ray_boundary_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
static const variant_type q
static const variant_type e
static const variant_type p
static const variant_type t
#define trace_macro(message)
#define error_macro(message)
#define fatal_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void serialize(Archive &ar, class rheolef::point_basic< T > &x, const unsigned int version)
This file is part of Rheolef.
T dist(const point_basic< T > &x, const point_basic< T > &y)
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
t operator()(const t &a, const t &b)
T dist2(const point_basic< T > &x, const point_basic< T > &y)
std::string typeid_name(const char *name, bool do_indent)
std::string ptos(const point_basic< T > &x, int d=3)