Rheolef  7.2
an efficient C++ finite element environment
geo_trace_ray_boundary.cc
Go to the documentation of this file.
1 //
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>
47 
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 
66 namespace 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 //
76 template <class T>
77 struct 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.:
93 namespace 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 
102 namespace rheolef {
103 
104 template <class T>
105 struct 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
123 template <class T>
124 static
125 void
126 dis_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 // ---------------------------------------------------------------------
142 template <class T, class M>
144 {
145  if (_ptr != 0) {
146  delete_macro(_ptr);
147  }
148 }
149 template <class T, class M>
150 bool
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 }
160 template <class T, class M>
161 bool
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 // ---------------------------------------------------------------------
174 template <class T, class M>
175 class geo_trace_ray_boundary_abstract_rep {
176 public:
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 // ---------------------------------------------------------------------
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> { };
197 
198 // ---------------------------------------------------------------------
199 // 2a) the 1D implementation
200 // ---------------------------------------------------------------------
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> {
203 public:
204 
205 // typedef:
206 
207  typedef typename geo_base_rep<T,M>::size_type size_type;
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 };
229 template <class T, class M>
230 void
231 geo_trace_ray_boundary_rep<T,M,1>::initialize (const geo_base_rep<T,M>& omega) const
232 {
233 }
234 template <class T, class M>
235 bool
236 geo_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 }
268 template <class T, class M>
269 bool
270 geo_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 // ---------------------------------------------------------------------
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> {
286 public:
287 
288 // typedef:
289 
290  typedef typename geo_base_rep<T,M>::size_type size_type;
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 };
312 template <class T, class M>
313 void
314 geo_trace_ray_boundary_rep<T,M,2>::initialize (const geo_base_rep<T,M>& omega) const
315 {
316 }
317 template <class T, class M>
318 bool
319 geo_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 }
386 template <class T, class M>
387 bool
388 geo_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 // ---------------------------------------------------------------------
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> {
404 public:
405 
406 // typedef:
407 
408  typedef typename geo_base_rep<T,M>::size_type size_type;
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:
446 protected:
447  mutable Tree _tree; // for cgal ray intersection query
448 };
449 template <class T, class M>
450 void
451 geo_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()) {
463  case reference_element::t: {
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  }
473  case reference_element::q: {
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 }
493 template <class T, class M>
494 bool
495 geo_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 }
546 template <class T, class M>
547 bool
548 geo_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 // ---------------------------------------------------------------------
562 template <class T, class M>
563 geo_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 // ---------------------------------------------------------------------
577 template <class T>
578 void
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
600 template <class T>
601 void
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->(),
669  pt2_minimum<T>());
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->(),
701  id_pt_minimum<T>());
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
733 template <class T, class M>
735 {
736 }
737 template <class T, class M>
738 bool
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 }
748 template <class T, class M>
749 bool
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 }
759 template <class T>
760 void
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
771 template <class T>
772 void
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 // ----------------------------------------------------------------------------
787 template class geo_rep<Float,sequential>;
788 template class geo_trace_ray_boundary<Float,sequential>;
789 
790 #ifdef _RHEOLEF_HAVE_MPI
791 template class geo_rep<Float,distributed>;
792 template 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
see the distributor page for the full documentation
Definition: distributor.h:69
size_type dis_size() const
global and local sizes
Definition: distributor.h:214
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:158
static const size_type decide
Definition: distributor.h:83
const communicator_type & comm() const
Definition: distributor.h:152
geo_element_hack::size_type size_type
Definition: geo.h:260
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
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
Definition: geo.h:934
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
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
size_t size_type
Definition: basis_get.cc:76
distributed
Definition: asr.cc:228
point_basic< T >
Definition: piola_fem.h:135
#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)")
void serialize(Archive &ar, class rheolef::point_basic< T > &x, const unsigned int version)
Definition: point.h:579
This file is part of Rheolef.
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:298
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
t operator()(const t &a, const t &b)
Definition: space.cc:386
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition: point.h:292
std::string typeid_name(const char *name, bool do_indent)
Definition: pretty_name.cc:78
std::string ptos(const point_basic< T > &x, int d=3)
Definition: point.h:413
static const bool value
Definition: communicator.h:67
Float epsilon
Expr1::memory_type M
Definition: vec_expr_v2.h:416