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
77struct hit_t : std::pair<bool,point_basic<T> > {
78 typedef std::pair<bool,point_basic<T> > base;
81 template<
class Archive>
82 void serialize (Archive& ar,
const unsigned int version) {
89#ifdef _RHEOLEF_HAVE_MPI
97 struct is_mpi_datatype<
rheolef::hit_t<double> > : mpl::true_ { };
105struct nearest_hit :
public std::binary_function<hit_t<T>, hit_t<T>, hit_t<T> > {
106 nearest_hit (
const point_basic<T>& x = point_basic<T>()) : _x(x) {}
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;
128 const point_basic<T>& x,
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));
142template <
class T,
class M>
149template <
class T,
class M>
157 if (_ptr == 0) { _ptr = make_ptr(omega); }
158 return _ptr->seq_trace_ray_boundary (omega, x, v, y);
160template <
class T,
class M>
168 if (_ptr == 0) { _ptr = make_ptr(omega); }
169 return _ptr->dis_trace_ray_boundary (omega, x, v, y);
174template <
class T,
class M>
175class 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,
182 const point_basic<T>& x,
183 const point_basic<T>& v,
184 point_basic<T>& y)
const = 0;
185 virtual bool dis_trace_ray_boundary (
186 const geo_base_rep<T,M>& omega,
187 const point_basic<T>& x,
188 const point_basic<T>& v,
189 point_basic<T>& y)
const = 0;
195template <
class T,
class M,
size_t D>
196struct geo_trace_ray_boundary_rep :
public geo_trace_ray_boundary_abstract_rep<T,M> { };
201template <
class T,
class M>
202class 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,
220 const point_basic<T>& x,
221 const point_basic<T>& v,
222 point_basic<T>& y)
const;
223 bool dis_trace_ray_boundary (
224 const geo_base_rep<T,M>& omega,
225 const point_basic<T>& x,
226 const point_basic<T>& v,
227 point_basic<T>& y)
const;
229template <
class T,
class M>
231geo_trace_ray_boundary_rep<T,M,1>::initialize (
const geo_base_rep<T,M>& omega)
const
234template <
class T,
class M>
236geo_trace_ray_boundary_rep<T,M,1>::seq_trace_ray_boundary (
237 const geo_base_rep<T,M>& omega,
238 const point_basic<T>& x,
239 const point_basic<T>& v,
240 point_basic<T>& y)
const
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);
254 const point_basic<T>&
a = omega.dis_node(S[0]);
256 if ((v0 > 0 && a0 > x0) || (v0 < 0 && a0 < x0) || (v0 == 0 && a0 == x0)) {
265 y = point_basic<T>(y0);
268template <
class T,
class M>
270geo_trace_ray_boundary_rep<T,M,1>::dis_trace_ray_boundary (
271 const geo_base_rep<T,M>& omega,
272 const point_basic<T>& x,
273 const point_basic<T>& v,
274 point_basic<T>& y)
const
276 bool hit = seq_trace_ray_boundary (omega, x, v, y);
278 dis_compute_nearest (omega.comm(), x, hit, y);
284template <
class T,
class M>
285class 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,
303 const point_basic<T>& x,
304 const point_basic<T>& v,
305 point_basic<T>& y)
const;
306 bool dis_trace_ray_boundary (
307 const geo_base_rep<T,M>& omega,
308 const point_basic<T>& x,
309 const point_basic<T>& v,
310 point_basic<T>& y)
const;
312template <
class T,
class M>
314geo_trace_ray_boundary_rep<T,M,2>::initialize (
const geo_base_rep<T,M>& omega)
const
317template <
class T,
class M>
319geo_trace_ray_boundary_rep<T,M,2>::seq_trace_ray_boundary (
320 const geo_base_rep<T,M>& omega,
321 const point_basic<T>& x,
322 const point_basic<T>& v,
323 point_basic<T>& y)
const
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);
345 const point_basic<T>&
a = omega.dis_node(S[0]);
346 const point_basic<T>&
b = omega.dis_node(S[1]);
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;
356 point_basic<T> y0 (xo.x(), xo.y());
363 }
else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
365 const Segment& so = *ptr_so;
366 point_basic<T> y0 (so[0].x(), so[0].y());
373 point_basic<T> y1 (so[1].x(), so[1].y());
386template <
class T,
class M>
388geo_trace_ray_boundary_rep<T,M,2>::dis_trace_ray_boundary (
389 const geo_base_rep<T,M>& omega,
390 const point_basic<T>& x,
391 const point_basic<T>& v,
392 point_basic<T>& y)
const
394 bool hit = seq_trace_ray_boundary (omega, x, v, y);
396 dis_compute_nearest (omega.comm(), x, hit, y);
402template <
class T,
class M>
403class 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,
436 const point_basic<T>& x,
437 const point_basic<T>& v,
438 point_basic<T>& y)
const;
439 bool dis_trace_ray_boundary (
440 const geo_base_rep<T,M>& omega,
441 const point_basic<T>& x,
442 const point_basic<T>& v,
443 point_basic<T>& y)
const;
449template <
class T,
class M>
451geo_trace_ray_boundary_rep<T,M,3>::initialize (
const geo_base_rep<T,M>& omega)
const
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()) {
464 const point_basic<T>&
a = omega.dis_node(S[0]);
465 const point_basic<T>&
b = omega.dis_node(S[1]);
466 const point_basic<T>&
c = omega.dis_node(S[2]);
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));
474 const point_basic<T>&
a = omega.dis_node(S[0]);
475 const point_basic<T>&
b = omega.dis_node(S[1]);
476 const point_basic<T>&
c = omega.dis_node(S[2]);
477 const point_basic<T>&
d = omega.dis_node(S[3]);
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");
493template <
class T,
class M>
495geo_trace_ray_boundary_rep<T,M,3>::seq_trace_ray_boundary (
496 const geo_base_rep<T,M>& omega,
497 const point_basic<T>& x,
498 const point_basic<T>& v,
499 point_basic<T>& y)
const
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;
517 point_basic<T> yi = point_basic<T>(yo.x(), yo.y(), yo.z());
524 }
else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
525 const Segment& so = *ptr_so;
526 point_basic<T> y0 (so[0].x(), so[0].y());
533 point_basic<T> y1 (so[1].x(), so[1].y());
546template <
class T,
class M>
548geo_trace_ray_boundary_rep<T,M,3>::dis_trace_ray_boundary (
549 const geo_base_rep<T,M>& omega,
550 const point_basic<T>& x,
551 const point_basic<T>& v,
552 point_basic<T>& y)
const
554 bool hit = seq_trace_ray_boundary (omega, x, v, y);
556 dis_compute_nearest (omega.comm(), x, hit, y);
562template <
class T,
class M>
563geo_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));
586 for (
size_type i = 0, n = x.size(); i < n; i++) {
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;
617 for (
size_type i = 0, n = x.size(); i < n; i++) {
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()) {
629 T eps = sqrt(std::numeric_limits<T>::epsilon());
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));
733template <
class T,
class M>
737template <
class T,
class M>
740 const geo_base_rep<T,M>& omega,
741 const point_basic<T>& x,
742 const point_basic<T>& v,
743 point_basic<T>& y)
const
745 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
748template <
class T,
class M>
751 const geo_base_rep<T,M>& omega,
752 const point_basic<T>& x,
753 const point_basic<T>& v,
754 point_basic<T>& y)
const
756 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
761geo_rep<T,sequential>::trace_ray_boundary (
762 const disarray<point_basic<T>,sequential>& x,
763 const disarray<point_basic<T>,sequential>& v,
764 disarray<size_type, sequential>& dis_ie,
765 disarray<point_basic<T>,sequential>& y,
768 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
770#ifdef _RHEOLEF_HAVE_MPI
773geo_rep<T,distributed>::trace_ray_boundary (
774 const disarray<point_basic<T>,distributed>& x,
775 const disarray<point_basic<T>,distributed>& v,
776 disarray<size_type, distributed>& dis_ie,
777 disarray<point_basic<T>,distributed>& y,
780 fatal_macro (
"geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
787template class geo_rep<Float,sequential>;
788template class geo_trace_ray_boundary<Float,sequential>;
790#ifdef _RHEOLEF_HAVE_MPI
791template class geo_rep<Float,distributed>;
792template class geo_trace_ray_boundary<Float,distributed>;