Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_trace_ray_boundary.cc
Go to the documentation of this file.
1
21//
22// given x and v, search S on boudary mesh such that ray(x,v) hits S
23// and returns the closest hit to x
24//
25// author: Pierre.Saramito@imag.fr
26//
27// date: 12 march 2012
28//
29// implementation note:
30// use CGAL::AABB_Tree for 3D surfacic element intersections
31//
32#include "rheolef/geo_trace_ray_boundary.h"
33#include "rheolef/geo.h"
34
35// internal includes:
36#include "rheolef/cgal_traits.h"
37#include "rheolef/point_util.h"
38
39// ------------------------------------------
40// assume that CGAL is configured for GNU C++
41// intel C++ has less 2011 standard features
42// thus re-config it!
43// ------------------------------------------
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
50#endif
51
52#ifndef CGAL_CFG_NO_TR1_TUPLE
53#define CGAL_CFG_NO_TR1_TUPLE
54#endif
55
56#ifndef CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
57#define CGAL_CFG_NO_CPP0X_VARIADIC_TEMPLATES
58#endif
59#endif // _RHEOLEF_HAVE_INTEL_CXX
60
61#include <CGAL/AABB_tree.h>
62#include <CGAL/AABB_traits.h>
63#include <CGAL/AABB_triangle_primitive.h>
64#endif // _RHEOLEF_HAVE_CGAL
65
66namespace rheolef {
67
68#ifdef _RHEOLEF_HAVE_CGAL
69// ---------------------------------------------------------------------
70// 0) utility: nearest hit over processors
71// ---------------------------------------------------------------------
72// formaly, we want to do in distributed ernvironment:
73// [hit,y] = mpi::all_reduce (omega.comm(), [hit,y], mpi::minimum<pair<bool,point> >());
74// where mpi::minimum<pair<bool,point> >(a,b) acts as the nearest point to x
75//
76template <class T>
77struct hit_t : std::pair<bool,point_basic<T> > {
78 typedef std::pair<bool,point_basic<T> > base;
79 hit_t () : base() {}
80 hit_t (bool hit, const point_basic<T>& y) : base(hit,y) {}
81 template<class Archive>
82 void serialize (Archive& ar, const unsigned int version) {
83 ar & base::first;
84 ar & base::second;
85 }
86};
87} // namespace rheolef
88
89#ifdef _RHEOLEF_HAVE_MPI
90// Some serializable types have a fixed amount of data stored at fixed field positions.
91// When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous copy operations.
92// To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
93namespace boost {
94 namespace mpi {
95 // TODO: when hit_t<T> is not a simple type, such as T=bigfloat or T=gmp, etc
96 template <>
97 struct is_mpi_datatype<rheolef::hit_t<double> > : mpl::true_ { };
98 } // namespace mpi
99} // namespace boost
100#endif // _RHEOLEF_HAVE_MPI
101
102namespace rheolef {
103
104template <class T>
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) {}
107 point_basic<T> _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;
111 // here: a & b have a valid ray intersecion
112 T da = dist2(a.second, _x);
113 T db = dist2(b.second, _x);
114 if (da < db) return a;
115 if (db < da) return b;
116 // here: a & b have an equi-distant ray intersecion:
117 // => it is the same intersection, at an inter-element boundary
118 return a;
119 }
120};
121// compute over processors the nearest hit y to x
122// (hit,y) are local to proc as input and global over procs as output
123template <class T>
124static
125void
126dis_compute_nearest (
127 const communicator& comm,
128 const point_basic<T>& x,
129 bool& hit,
130 point_basic<T>& y)
131{
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));
135#endif // _RHEOLEF_HAVE_MPI
136 hit = hit_y.first;
137 y = hit_y.second;
138}
139// ---------------------------------------------------------------------
140// 1) the ray tracer interface
141// ---------------------------------------------------------------------
142template <class T, class M>
144{
145 if (_ptr != 0) {
146 delete_macro(_ptr);
147 }
148}
149template <class T, class M>
150bool
152 const geo_base_rep<T,M>& omega,
153 const point_basic<T>& x,
154 const point_basic<T>& v,
155 point_basic<T>& y) const
156{
157 if (_ptr == 0) { _ptr = make_ptr(omega); }
158 return _ptr->seq_trace_ray_boundary (omega, x, v, y);
159}
160template <class T, class M>
161bool
163 const geo_base_rep<T,M>& omega,
164 const point_basic<T>& x,
165 const point_basic<T>& v,
166 point_basic<T>& y) const
167{
168 if (_ptr == 0) { _ptr = make_ptr(omega); }
169 return _ptr->dis_trace_ray_boundary (omega, x, v, y);
170}
171// ---------------------------------------------------------------------
172// 2) the ray tracer abtract data structure
173// ---------------------------------------------------------------------
174template <class T, class M>
175class geo_trace_ray_boundary_abstract_rep {
176public:
177 typedef typename disarray<T,M>::size_type size_type;
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;
190};
191// ---------------------------------------------------------------------
192// 2) the ray tracer concrete, and dimension dependent, data structure
193// = CGAL::AABB_Tree in 3D and empty otherwise
194// ---------------------------------------------------------------------
195template <class T, class M, size_t D>
196struct geo_trace_ray_boundary_rep : public geo_trace_ray_boundary_abstract_rep<T,M> { };
197
198// ---------------------------------------------------------------------
199// 2a) the 1D implementation
200// ---------------------------------------------------------------------
201template <class T, class M>
202class geo_trace_ray_boundary_rep<T,M,1> : public geo_trace_ray_boundary_abstract_rep<T,M> {
203public:
204
205// typedef:
206
208
209// allocators:
210
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;
215
216// accessors:
217
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;
228};
229template <class T, class M>
230void
231geo_trace_ray_boundary_rep<T,M,1>::initialize (const geo_base_rep<T,M>& omega) const
232{
233}
234template <class T, class M>
235bool
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
241{
242 T x0 = x[0];
243 T v0 = v[0];
244 T y0 = 0;
245 bool hit = false;
246 T d_min = std::numeric_limits<T>::max();
247 const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
248 check_macro (boundary.map_dimension() == 0, "unexpected boundary domain");
249 for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
250 const geo_element_indirect& oige = boundary.oige(ioige);
251 size_type ie = oige.index();
252 const geo_element& S = omega.get_geo_element(0,ie);
253 check_macro (S.variant() == reference_element::p, "unexpected element type: "<<S.name());
254 const point_basic<T>& a = omega.dis_node(S[0]);
255 T a0 = a[0];
256 if ((v0 > 0 && a0 > x0) || (v0 < 0 && a0 < x0) || (v0 == 0 && a0 == x0)) {
257 T di = fabs(x0-a0);
258 if (di < d_min) {
259 hit = true;
260 y0 = a0;
261 d_min = di;
262 }
263 }
264 }
265 y = point_basic<T>(y0);
266 return hit;
267}
268template <class T, class M>
269bool
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
275{
276 bool hit = seq_trace_ray_boundary (omega, x, v, y);
277 if (is_sequential<M>::value) return hit;
278 dis_compute_nearest (omega.comm(), x, hit, y);
279 return hit;
280}
281// ---------------------------------------------------------------------
282// 2b) the 2D implementation
283// ---------------------------------------------------------------------
284template <class T, class M>
285class geo_trace_ray_boundary_rep<T,M,2> : public geo_trace_ray_boundary_abstract_rep<T,M> {
286public:
287
288// typedef:
289
291
292// allocators:
293
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;
298
299// accessors:
300
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;
311};
312template <class T, class M>
313void
314geo_trace_ray_boundary_rep<T,M,2>::initialize (const geo_base_rep<T,M>& omega) const
315{
316}
317template <class T, class M>
318bool
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
324{
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;
330
331 Point x1 (x[0], x[1]);
332 Vector2d v1 (v[0], v[1]);
333 Ray ray_query(x1,v1);
334
335 bool hit = false;
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");
339 check_macro (boundary.map_dimension() == 1, "unexpected boundary domain");
340 for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
341 const geo_element_indirect& oige = boundary.oige(ioige);
342 size_type ie = oige.index();
343 const geo_element& S = omega.get_geo_element(1,ie);
344 check_macro (S.variant() == reference_element::e, "unexpected element type: "<<S.name());
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]);
349 Segment s (a1, b1);
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))) {
354 n_intersect++;
355 const Point& xo = *ptr_xo;
356 point_basic<T> y0 (xo.x(), xo.y());
357 T d0 = dist(y0,x);
358 if (d0 < d_min) {
359 hit = true;
360 y = y0;
361 d_min = d0;
362 }
363 } else if ((ptr_so = CGAL::object_cast<Segment> (&obj))) {
364 n_intersect += 2;
365 const Segment& so = *ptr_so;
366 point_basic<T> y0 (so[0].x(), so[0].y());
367 T d0 = dist(y0,x);
368 if (d0 < d_min) {
369 hit = true;
370 y = y0;
371 d_min = d0;
372 }
373 point_basic<T> y1 (so[1].x(), so[1].y());
374 T d1 = dist(y1,x);
375 if (d1 < d_min) {
376 hit = true;
377 y = y1;
378 d_min = d1;
379 }
380 } else {
381 // no intersection
382 }
383 }
384 return hit;
385}
386template <class T, class M>
387bool
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
393{
394 bool hit = seq_trace_ray_boundary (omega, x, v, y);
395 if (is_sequential<M>::value) return hit;
396 dis_compute_nearest (omega.comm(), x, hit, y);
397 return hit;
398}
399// ---------------------------------------------------------------------
400// 2c) the 3D implementation
401// ---------------------------------------------------------------------
402template <class T, class M>
403class geo_trace_ray_boundary_rep<T,M,3> : public geo_trace_ray_boundary_abstract_rep<T,M> {
404public:
405
406// typedef:
407
409
410 typedef CGAL::Filtered_kernel<CGAL::Simple_cartesian<T> > Kernel;
411 // typedef CGAL::Filtered_kernel_adaptor<mycgal::my_kernel_3d<T> > Kernel; // missing some 3d features
412
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;
424
425// allocators:
426
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;
431
432// accessors:
433
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;
444
445// data:
446protected:
447 mutable Tree _tree; // for cgal ray intersection query
448};
449template <class T, class M>
450void
451geo_trace_ray_boundary_rep<T,M,3>::initialize (const geo_base_rep<T,M>& omega) const
452{
453 // create the corresponding tree of bounding box
454 trace_macro ("create the 3d CGAL::AABB_tree...");
455 std::list<Triangle> triangles;
456 const domain_indirect_basic<M>& boundary = omega.get_domain_indirect ("boundary");
457 check_macro (boundary.map_dimension() == 2, "unexpected boundary domain");
458 for (size_t ioige = 0, noige = boundary.size(); ioige < noige; ioige++) {
459 const geo_element_indirect& oige = boundary.oige(ioige);
460 size_type ie = oige.index();
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));
471 break;
472 }
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));
484 break;
485 }
486 default: error_macro ("unexpected element type: "<<S.name());
487 }
488 }
489 // constructs the AABB tree
490 _tree.rebuild (triangles.begin(), triangles.end());
491 trace_macro ("create the 3d CGAL::AABB_tree done");
492}
493template <class T, class M>
494bool
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
500{
501 Point x1 (x[0], x[1], x[2]);
502 Vector3d v1 (v[0], v[1], v[2]);
503 Ray ray_query (x1,v1);
504
505 // computes all intersections with segment query (as pairs object - primitive_id)
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();
509 bool hit = false;
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());
518 T di = dist(yi,x);
519 if (di < d_min) {
520 hit = true;
521 y = yi;
522 d_min = di;
523 }
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());
527 T d0 = dist(y0,x);
528 if (d0 < d_min) {
529 hit = true;
530 y = y0;
531 d_min = d0;
532 }
533 point_basic<T> y1 (so[1].x(), so[1].y());
534 T d1 = dist(y1,x);
535 if (d1 < d_min) {
536 hit = true;
537 y = y1;
538 d_min = d1;
539 }
540 } else {
541 error_macro ("unexpected intersection type: " << typeid_name(obj.type().name(), false));
542 }
543 }
544 return hit;
545}
546template <class T, class M>
547bool
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
553{
554 bool hit = seq_trace_ray_boundary (omega, x, v, y);
555 if (is_sequential<M>::value) return hit;
556 dis_compute_nearest (omega.comm(), x, hit, y);
557 return hit;
558}
559// ---------------------------------------------------------------------
560// 3) the pointer allocator
561// ---------------------------------------------------------------------
562template <class T, class M>
563geo_trace_ray_boundary_abstract_rep<T,M>*
565{
566 check_macro (omega.dimension() == omega.map_dimension(), "geo_trace_ray_boundary: map_dim < dim: not supported");
567 switch (omega.dimension()) {
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));
571 default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
572 }
573}
574// ---------------------------------------------------------------------
575// 4) disarray interface (more efficient in the distributed case)
576// ---------------------------------------------------------------------
577template <class T>
578void
584 bool do_check) const
585{
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]);
588 if (hit) {
589 dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
590 } else {
591 dis_ie[i] = std::numeric_limits<size_type>::max();
592 }
593 if (do_check) {
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));
596 }
597 }
598}
599#ifdef _RHEOLEF_HAVE_MPI
600template <class T>
601void
607 bool do_check) const
608{
609 trace_macro ("ray_boundary...");
610 const size_type large = std::numeric_limits<size_type>::max();
611 const T infty = std::numeric_limits<T>::max();
612 // -----------------------------------------------------------------
613 // 1) first try a ray(x[i],v[i]) on the local boundary partition
614 // -----------------------------------------------------------------
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]);
619 if (hit) {
620 // here x[i)+v[i] cross the local boundary partition
621 dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
622 if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
623 // y was outside Omega ? most of the time, y is exactly on the boundary:
624 // dirty hack: project y on the bbox of the boundary:
625 // - this works in 1d when Omega is simply connected
626 // - this works in 2d & 3d when Omega is a rectangle or a paralellotope...
627 // TODO: boundary projection in the multi-connected case with: y = omega.nearest(x+v)
628 // T eps = 1e3*std::numeric_limits<T>::epsilon();
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]);
632 }
633 dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
634 }
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));
638 } else {
639 // here x[i)+v[i] cross the boundary, on a boundary part that is managed by another proc
640 failed.push_back (i);
641 }
642 }
643 // -----------------------------------------------------------------
644 // 3) for all x+v that goes outside and in another boundary partition
645 // then solve a distributed ray trace
646 // ---------------------------------------------------------------
647 distributor fld_ownership (distributor::decide, base::comm(), failed.size());
648 size_type fld_dis_size = fld_ownership.dis_size();
649 if (fld_dis_size == 0) {
650 // all are solved: nothing more to do
651 trace_macro ("ray_boundary done(1)");
652 return;
653 }
654 size_type first_fld_dis_i = fld_ownership.first_index();
655 size_type last_fld_dis_i = fld_ownership. last_index();
656 pt2_t<T> unset2 (point_basic<T>(infty,infty,infty), point_basic<T>(infty,infty,infty));
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) {
660 size_type i = *iter;
661 massive_failed [fld_dis_i] = pt2_t<T>(x[i],v[i]);
662 }
663 std::vector<pt2_t<T> > massive_query (fld_dis_size, unset2);
664 mpi::all_reduce (
665 fld_ownership.comm(),
666 massive_failed.begin().operator->(),
667 massive_failed.size(),
668 massive_query.begin().operator->(),
670
671 // 3) run the locator on all failed points ON ALL PROCS, skipping local queries (already failed)
672 id_pt_t<T> unset (large, point_basic<T>(infty,infty,infty));
673 std::vector<id_pt_t<T> > massive_result (fld_dis_size, unset);
674 // 3a) range [0:first[
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);
680 if (hit) {
681 massive_result [fld_dis_i].first = base::_locator.seq_locate (*this, massive_result[fld_dis_i].second);
682 }
683 }
684 // 3b) range [last,dis_size[
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);
690 if (hit) {
691 massive_result [fld_dis_i].first = base::_locator.seq_locate (*this, massive_result[fld_dis_i].second);
692 }
693 }
694 // 4) send & merge the results to all
695 std::vector<id_pt_t<T> > massive_merged (fld_dis_size, unset);
696 mpi::all_reduce (
697 fld_ownership.comm(),
698 massive_result.begin().operator->(),
699 massive_result.size(),
700 massive_merged.begin().operator->(),
702
703 // 5) store the local range into the distributed disarray dis_ie:
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) {
706 size_type i = *iter;
707 dis_ie[i] = massive_merged [fld_dis_i].first;
708 y[i] = massive_merged [fld_dis_i].second;
709
710 if (dis_ie[i] == std::numeric_limits<size_type>::max()) {
711 // then x belongs to the boundary and v.n > 0 on the boundary
712 // => the ray(x,v) goes outside
713 // => takes y=x
714 y[i] = x[i];
715 dis_ie[i] = base::_locator.seq_locate (*this, y[i]);
716 }
717 if (do_check) {
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));
724 }
725 }
726 trace_macro ("ray_boundary done(2)");
727}
728#endif // _RHEOLEF_HAVE_MPI
729// ----------------------------------------------------------------------------
730// CGAL lib is missing
731// ----------------------------------------------------------------------------
732#else // _RHEOLEF_HAVE_CGAL
733template <class T, class M>
735{
736}
737template <class T, class M>
738bool
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
744{
745 fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
746 return false;
747}
748template <class T, class M>
749bool
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
755{
756 fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
757 return false;
758}
759template <class T>
760void
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,
766 bool do_check) const
767{
768 fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
769}
770#ifdef _RHEOLEF_HAVE_MPI
771template <class T>
772void
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,
778 bool do_check) const
779{
780 fatal_macro ("geo: trace_ray_boundary not available (HINT: recompile Rheolef with the CGAL library)");
781}
782#endif // _RHEOLEF_HAVE_MPI
783#endif // _RHEOLEF_HAVE_CGAL
784// ----------------------------------------------------------------------------
785// instanciation in library
786// ----------------------------------------------------------------------------
787template class geo_rep<Float,sequential>;
788template class geo_trace_ray_boundary<Float,sequential>;
789
790#ifdef _RHEOLEF_HAVE_MPI
791template class geo_rep<Float,distributed>;
792template class geo_trace_ray_boundary<Float,distributed>;
793#endif // _RHEOLEF_HAVE_MPI
794
795} // namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the communicator page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
rep::base::size_type size_type
Definition disarray.h:501
see the distributor page for the full documentation
Definition distributor.h:69
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
Definition distributor.h:83
const communicator_type & comm() const
base class for M=sequential or distributed meshes representations
Definition geo.h:528
size_type map_dimension() const
Definition geo.h:564
base::size_type size_type
Definition geo.h:533
size_type dimension() const
Definition geo.h:563
base::size_type size_type
Definition geo.h:934
base::size_type size_type
Definition geo.h:791
sequential mesh representation
Definition geo.h:778
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
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)
Definition dis_macros.h:111
#define error_macro(message)
Definition dis_macros.h:49
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
std::string typeid_name(const char *name, bool do_indent)
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
static const bool value
Expr1::memory_type M