31 #include "rheolef/geo_locate.h"
32 #include "rheolef/geo.h"
33 #include "rheolef/geo_element_contains.h"
34 #include "rheolef/point_util.h"
37 #ifdef _RHEOLEF_HAVE_CGAL
38 #include "rheolef/cgal_traits.h"
39 #include <CGAL/Segment_tree_k.h>
40 #include <CGAL/Range_segment_tree_traits.h>
41 #endif // _RHEOLEF_HAVE_CGAL
48 template <
class T,
class M>
56 xmin[j] = std::numeric_limits<T>::max();
57 xmax[j] = -std::numeric_limits<T>::max();
59 for (
size_type iloc = 0, nloc = K.
size(); iloc < nloc; iloc++) {
63 xmin[j] = std::min(x[j], xmin[j]);
64 xmax[j] = std::max(x[j], xmax[j]);
68 #ifdef _RHEOLEF_HAVE_CGAL
72 template <
class Kernel,
class Val>
73 class Segment_tree_map_traits_1 {
76 typedef typename Kernel::FT Point_1 ;
78 typedef Point_1 Key_1;
79 typedef std::pair<Key,Key> Pure_interval;
80 typedef std::pair<Pure_interval, Val> Interval;
84 Key_1 operator()(
const Interval& i)
85 {
return i.first.first;}
90 Key_1 operator()(
const Interval& i)
91 {
return i.first.second;}
96 Key_1 operator()(
const Key& k)
102 bool operator()(Key_1 k1, Key_1 k2)
104 return std::less<Float>()(k1,k2);
108 typedef C_Compare_1 compare_1;
109 typedef C_Low_1 low_1;
110 typedef C_High_1 high_1;
111 typedef C_Key_1 key_1;
117 template <
class T,
size_t D>
struct cgal_locate_traits {};
120 struct cgal_locate_traits<
T,1> {
127 typedef Segment_tree_map_traits_1<Kernel,size_type> Traits;
128 typedef ::CGAL::Segment_tree_1<Traits> Segment_tree_type;
129 typedef typename Traits::Interval Interval;
130 typedef typename Traits::Pure_interval Pure_interval;
131 typedef typename Traits::Key Key;
135 return Pure_interval(Key(xmin[0]),
138 static Pure_interval make_cgal_point_window (
const point_basic<T>& x,
const T& eps)
140 return Pure_interval (Key(x[0]-eps),
143 static void put (std::ostream&
out,
const Pure_interval&
b) {
144 out <<
"[" <<
b.first <<
"," <<
b.second <<
"[";
148 struct cgal_locate_traits<
T,2> {
155 typedef ::CGAL::Segment_tree_map_traits_2<Kernel,size_type> Traits;
156 typedef ::CGAL::Segment_tree_2<Traits> Segment_tree_type;
157 typedef typename Traits::Interval Interval;
158 typedef typename Traits::Pure_interval Pure_interval;
159 typedef typename Traits::Key Key;
163 return Pure_interval(Key(xmin[0], xmin[1]),
164 Key(xmax[0], xmax[1]));
166 static Pure_interval make_cgal_point_window (
const point_basic<T>& x,
const T& eps)
168 return Pure_interval (Key(x[0]-eps, x[1]-eps),
169 Key(x[0]+eps, x[1]+eps));
171 static void put (std::ostream&
out,
const Pure_interval&
b) {
172 out <<
"[" <<
b.first.x() <<
"," <<
b.second.x() <<
"[x["
173 <<
b.first.y() <<
"," <<
b.second.y() <<
"[";
177 struct cgal_locate_traits<
T,3> {
184 typedef ::CGAL::Segment_tree_map_traits_3<Kernel,size_type> Traits;
185 typedef ::CGAL::Segment_tree_3<Traits> Segment_tree_type;
186 typedef typename Traits::Interval Interval;
187 typedef typename Traits::Pure_interval Pure_interval;
188 typedef typename Traits::Key Key;
192 return Pure_interval(Key(xmin[0], xmin[1], xmin[2]),
193 Key(xmax[0], xmax[1], xmax[2]));
195 static Pure_interval make_cgal_point_window (
const point_basic<T>& x,
const T& eps)
197 return Pure_interval (Key(x[0]-eps, x[1]-eps, x[2]-eps),
198 Key(x[0]+eps, x[1]+eps, x[2]+eps));
200 static void put (std::ostream&
out,
const Pure_interval&
b) {
201 out <<
"[" <<
b.first.x() <<
"," <<
b.second.x() <<
"[x["
202 <<
b.first.y() <<
"," <<
b.second.y() <<
"[x["
203 <<
b.first.z() <<
"," <<
b.second.z() <<
"[";
209 template <
class T,
class M>
216 template <
class T,
class M>
223 if (_ptr == 0) { _ptr = make_ptr(omega); }
224 return _ptr->seq_locate (omega, x, dis_ie_guest);
226 template <
class T,
class M>
233 if (_ptr == 0) { _ptr = make_ptr(omega); }
234 return _ptr->dis_locate (omega, x, dis_ie_guest);
239 template <
class T,
class M>
240 class geo_locate_abstract_rep {
243 virtual ~geo_locate_abstract_rep() {}
244 virtual void initialize (
const geo_base_rep<T,M>& omega)
const = 0;
246 const geo_base_rep<T,M>& omega,
248 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const = 0;
250 const geo_base_rep<T,M>& omega,
252 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const = 0;
258 template <
class T,
class M,
size_t D>
259 class geo_locate_rep :
public geo_locate_abstract_rep<T,M> {
266 typedef typename cgal_locate_traits<T,D>::Segment_tree_type Segment_tree_type;
267 typedef typename cgal_locate_traits<T,D>::Interval Interval;
271 geo_locate_rep() : _tree() {}
272 geo_locate_rep(
const geo_base_rep<T,M>& omega) : _tree() {
initialize(omega); }
274 void initialize (
const geo_base_rep<T,M>& omega)
const;
279 const geo_base_rep<T,M>& omega,
281 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const;
283 const geo_base_rep<T,M>& omega,
285 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const;
289 mutable Segment_tree_type _tree;
294 template <
class T,
class M>
295 geo_locate_abstract_rep<T,M>*
300 case 1:
return new_macro((geo_locate_rep<T,M,1>)(omega));
301 case 2:
return new_macro((geo_locate_rep<T,M,2>)(omega));
302 case 3:
return new_macro((geo_locate_rep<T,M,3>)(omega));
309 template <
class T,
class M,
size_t D>
315 std::list<Interval> boxes;
320 boxes.push_back (Interval(cgal_locate_traits<T,D>::make_cgal_bbox (xmin,xmax),ie));
324 _tree.make_tree (boxes.begin(), boxes.end());
330 template <
class T,
class M,
size_t D>
332 geo_locate_rep<T,M,D>::seq_locate (
333 const geo_base_rep<T,M>& omega,
337 if (dis_ie_guest != std::numeric_limits<size_type>::max()) {
339 const distributor& ownership = omega.ownership();
340 if (ownership.is_owned (dis_ie_guest)) {
341 size_type first_dis_ie = ownership.first_index();
342 size_type ie_guest = dis_ie_guest - first_dis_ie;
343 const geo_element& K_guest = omega[ie_guest];
346 return K_guest.dis_ie();
351 size_type dis_ie = std::numeric_limits<size_type>::max();
357 Interval xe = Interval (cgal_locate_traits<T,D>::make_cgal_point_window (x, eps), 0);
358 std::list<Interval> intersected_boxes;
363 _tree.window_query (xe, std::back_inserter(intersected_boxes));
364 for (
typename std::list<Interval>::iterator j = intersected_boxes.begin(); j != intersected_boxes.end(); j++) {
366 const geo_element& K = omega[ie];
378 template <
class T,
class M,
size_t D>
380 geo_locate_rep<T,M,D>::dis_locate (
381 const geo_base_rep<T,M>& omega,
385 size_type dis_ie = seq_locate (omega, x, dis_ie_guest);
386 #ifdef _RHEOLEF_HAVE_MPI
388 dis_ie = mpi::all_reduce (omega.comm(), dis_ie, mpi::minimum<size_type>());
390 #endif // _RHEOLEF_HAVE_MPI
398 template <
class T,
class M>
404 return _locator.seq_locate (*
this, x, dis_ie_guest);
406 template <
class T,
class M>
412 return _locator.dis_locate (*
this, x, dis_ie_guest);
424 dis_ie.resize (x.ownership());
426 dis_ie[i] = base::_locator.seq_locate (*
this, x[i], dis_ie[i]);
428 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
429 "locate: failed at x="<<
ptos(x[i],base::_dimension));
433 #ifdef _RHEOLEF_HAVE_MPI
443 const size_type large = std::numeric_limits<size_type>::max();
444 const T infty = std::numeric_limits<T>::max();
447 if (dis_ie.ownership() != ownership) dis_ie.resize (ownership);
449 std::list<id_pt_t<T> > failed;
451 dis_ie[i] = base::_locator.seq_locate (*
this, x[i], dis_ie[i]);
452 if (dis_ie[i] == large) {
458 communicator comm = ownership.
comm();
461 if (comm.size() == 1 || fld_dis_size == 0) {
466 size_type first_fld_dis_i = fld_ownership.first_index();
467 size_type last_fld_dis_i = fld_ownership. last_index();
469 std::vector<id_pt_t<T> > massive_failed (fld_dis_size, unset);
470 typename std::list<id_pt_t<T> >
::iterator iter = failed.begin();
471 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
472 massive_failed [fld_dis_i] = *iter;
474 std::vector<id_pt_t<T> > massive_query (fld_dis_size, unset);
477 massive_failed.begin().operator->(),
478 massive_failed.size(),
479 massive_query.begin().operator->(),
483 std::vector<size_type> massive_result (fld_dis_size, large);
485 for (
size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
486 massive_result [fld_dis_i] = base::_locator.seq_locate (*
this, massive_query[fld_dis_i].second);
489 for (
size_type fld_dis_i = last_fld_dis_i; fld_dis_i < fld_dis_size; ++fld_dis_i) {
490 massive_result [fld_dis_i] = base::_locator.seq_locate (*
this, massive_query[fld_dis_i].second);
493 std::vector<size_type> massive_merged (fld_dis_size, large);
496 massive_result.begin().operator->(),
497 massive_result.size(),
498 massive_merged.begin().operator->(),
499 mpi::minimum<size_type>());
502 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
503 size_type dis_i = massive_query [fld_dis_i].first;
504 check_macro (dis_i >= first_dis_i,
"invalid index");
506 dis_ie[i] = massive_merged [fld_dis_i];
509 "dis_locate: failed at x="<<
ptos(x[i],base::_dimension));
514 #endif // _RHEOLEF_HAVE_MPI
518 #else // _RHEOLEF_HAVE_CGAL
519 template <
class T,
class M>
523 template <
class T,
class M>
526 const geo_base_rep<T,M>& omega,
530 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
533 template <
class T,
class M>
536 const geo_base_rep<T,M>& omega,
540 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
543 template <
class T,
class M>
549 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
552 template <
class T,
class M>
558 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
565 disarray<size_type, sequential>& dis_ie,
568 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
570 #ifdef _RHEOLEF_HAVE_MPI
575 disarray<size_type, distributed>& dis_ie,
579 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
581 #endif // _RHEOLEF_HAVE_MPI
582 #endif // _RHEOLEF_HAVE_CGAL
586 #define _RHEOLEF_instanciation(T,M) \
587 template class geo_locate<Float,M>; \
588 template class geo_base_rep<Float,M>; \
589 template class geo_rep<Float,M>;
592 #ifdef _RHEOLEF_HAVE_MPI
594 #endif // _RHEOLEF_HAVE_MPI