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
72template <
class Kernel,
class Val>
73class 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;
85 {
return i.first.first;}
91 {
return i.first.second;}
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;
117template <
class T,
size_t D>
struct cgal_locate_traits {};
120struct cgal_locate_traits<
T,1> {
124 typedef typename geo_base_rep<T,sequential>::size_type
size_type;
126 typedef typename geo_cgal_traits<T,1>::Kernel Kernel;
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;
133 static Pure_interval make_cgal_bbox (
const point_basic<T>& xmin,
const point_basic<T>& xmax)
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 <<
"[";
148struct cgal_locate_traits<
T,2> {
152 typedef typename geo_base_rep<T,sequential>::size_type
size_type;
154 typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
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;
161 static Pure_interval make_cgal_bbox (
const point_basic<T>& xmin,
const point_basic<T>& xmax)
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() <<
"[";
177struct cgal_locate_traits<
T,3> {
181 typedef typename geo_base_rep<T,sequential>::size_type
size_type;
183 typedef typename geo_cgal_traits<T,3>::Kernel Kernel;
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;
190 static Pure_interval make_cgal_bbox (
const point_basic<T>& xmin,
const point_basic<T>& xmax)
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() <<
"[";
209template <
class T,
class M>
216template <
class T,
class M>
223 if (_ptr == 0) { _ptr = make_ptr(omega); }
224 return _ptr->seq_locate (omega, x, dis_ie_guest);
226template <
class T,
class M>
233 if (_ptr == 0) { _ptr = make_ptr(omega); }
234 return _ptr->dis_locate (omega, x, dis_ie_guest);
239template <
class T,
class M>
240class 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,
247 const point_basic<T>& x,
248 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const = 0;
250 const geo_base_rep<T,M>& omega,
251 const point_basic<T>& x,
252 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const = 0;
258template <
class T,
class M,
size_t D>
259class 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,
280 const point_basic<T>& x,
281 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const;
283 const geo_base_rep<T,M>& omega,
284 const point_basic<T>& x,
285 size_type dis_ie_guest = std::numeric_limits<size_type>::max())
const;
289 mutable Segment_tree_type _tree;
294template <
class T,
class M>
295geo_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));
309template <
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());
330template <
class T,
class M,
size_t D>
331typename geo_locate_rep<T,M,D>::size_type
332geo_locate_rep<T,M,D>::seq_locate (
333 const geo_base_rep<T,M>& omega,
334 const point_basic<T>& x,
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();
356 static const T eps = sqrt(std::numeric_limits<T>::epsilon());
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];
378template <
class T,
class M,
size_t D>
379typename geo_locate_rep<T,M,D>::size_type
380geo_locate_rep<T,M,D>::dis_locate (
381 const geo_base_rep<T,M>& omega,
382 const point_basic<T>& x,
387 using signed_size_type = long;
388 const signed_size_type signed_large = -1;
389 size_type dis_ie = seq_locate (omega, x, dis_ie_guest);
390#ifdef _RHEOLEF_HAVE_MPI
393 dis_ie = mpi::all_reduce (omega.comm(), signed_size_type(dis_ie), mpi::maximum<signed_size_type>());
403template <
class T,
class M>
409 return _locator.seq_locate (*
this, x, dis_ie_guest);
411template <
class T,
class M>
417 return _locator.dis_locate (*
this, x, dis_ie_guest);
429 dis_ie.resize (x.ownership());
430 for (
size_type i = 0, n = x.size(); i < n; i++) {
431 dis_ie[i] = base::_locator.seq_locate (*
this, x[i], dis_ie[i]);
433 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
434 "locate: failed at x="<<
ptos(x[i],base::_dimension));
438#ifdef _RHEOLEF_HAVE_MPI
449 using signed_size_type = long;
450 const signed_size_type signed_large = -1;
451 const size_type large = std::numeric_limits<size_type>::max();
452 const T infty = std::numeric_limits<T>::max();
455 if (dis_ie.ownership() != ownership) dis_ie.resize (ownership);
457 std::list<id_pt_t<T> > failed;
458 for (
size_type i = 0, n = x.size(); i < n; i++) {
459 dis_ie[i] = base::_locator.seq_locate (*
this, x[i], dis_ie[i]);
460 if (dis_ie[i] == large) {
466 communicator comm = ownership.
comm();
469 if (comm.size() == 1 || fld_dis_size == 0) {
475 size_type last_fld_dis_i = fld_ownership. last_index();
477 std::vector<id_pt_t<T> > massive_failed (fld_dis_size, unset);
478 typename std::list<id_pt_t<T> >::iterator iter = failed.begin();
479 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
480 massive_failed [fld_dis_i] = *iter;
482 std::vector<id_pt_t<T> > massive_query (fld_dis_size, unset);
485 massive_failed.begin().operator->(),
486 massive_failed.size(),
487 massive_query.begin().operator->(),
491 std::vector<signed_size_type> massive_result (fld_dis_size, signed_large);
493 for (
size_type fld_dis_i = 0; fld_dis_i < first_fld_dis_i; ++fld_dis_i) {
494 massive_result [fld_dis_i] = base::_locator.seq_locate (*
this, massive_query[fld_dis_i].second);
497 for (
size_type fld_dis_i = last_fld_dis_i; fld_dis_i < fld_dis_size; ++fld_dis_i) {
498 massive_result [fld_dis_i] = base::_locator.seq_locate (*
this, massive_query[fld_dis_i].second);
501 std::vector<signed_size_type> massive_merged (fld_dis_size, signed_large);
504 massive_result.begin().operator->(),
505 massive_result.size(),
506 massive_merged.begin().operator->(),
507 mpi::maximum<signed_size_type>());
510 for (
size_type fld_dis_i = first_fld_dis_i; fld_dis_i < last_fld_dis_i; ++fld_dis_i, ++iter) {
511 size_type dis_i = massive_query [fld_dis_i].first;
512 check_macro (dis_i >= first_dis_i,
"invalid index");
514 dis_ie[i] = massive_merged [fld_dis_i];
517 "dis_locate: failed at x="<<
ptos(x[i],base::_dimension));
527template <
class T,
class M>
531template <
class T,
class M>
534 const geo_base_rep<T,M>& omega,
535 const point_basic<T>& x,
538 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
541template <
class T,
class M>
544 const geo_base_rep<T,M>& omega,
545 const point_basic<T>& x,
548 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
551template <
class T,
class M>
554 const point_basic<T>& x,
557 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
560template <
class T,
class M>
563 const point_basic<T>& x,
566 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
571geo_rep<T,sequential>::locate (
572 const disarray<point_basic<T>, sequential>& x,
573 disarray<size_type, sequential>& dis_ie,
576 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
578#ifdef _RHEOLEF_HAVE_MPI
581geo_rep<T,distributed>::locate (
582 const disarray<point_basic<T>, distributed>& x,
583 disarray<size_type, distributed>& dis_ie,
587 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
594#define _RHEOLEF_instanciation(T,M) \
595template class geo_locate<Float,M>; \
596template class geo_base_rep<Float,M>; \
597template class geo_rep<Float,M>;
600#ifdef _RHEOLEF_HAVE_MPI