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>
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;
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;
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,
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>());
403 template <
class T,
class M>
409 return _locator.seq_locate (*
this, x, dis_ie_guest);
411 template <
class T,
class M>
417 return _locator.dis_locate (*
this, x, dis_ie_guest);
429 dis_ie.resize (x.ownership());
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;
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));
527 template <
class T,
class M>
531 template <
class T,
class M>
534 const geo_base_rep<T,M>& omega,
538 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
541 template <
class T,
class M>
544 const geo_base_rep<T,M>& omega,
548 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
551 template <
class T,
class M>
557 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
560 template <
class T,
class M>
566 fatal_macro (
"geo: locator not available (HINT: recompile Rheolef with the CGAL library)");
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
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) \
595 template class geo_locate<Float,M>; \
596 template class geo_base_rep<Float,M>; \
597 template class geo_rep<Float,M>;
600 #ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the disarray page for the full documentation
see the distributor page for the full documentation
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
const communicator_type & comm() const
geo_element_hack::size_type size_type
base class for M=sequential or distributed meshes representations
size_type dis_locate(const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
const node_type & dis_node(size_type dis_inod) const
size_type map_dimension() const
size_type seq_locate(const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
base::size_type size_type
size_type dimension() const
size_type size(size_type dim) const
see the geo_element page for the full documentation
size_type seq_locate(const geo_base_rep< T, M > &omega, const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
size_type dis_locate(const geo_base_rep< T, M > &omega, const point_basic< T > &x, size_type dis_ie_guest=std::numeric_limits< size_type >::max()) const
static geo_locate_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
disarray< T, M >::size_type size_type
base::size_type size_type
void locate(const disarray< point_basic< T >, distributed > &x, disarray< size_type, distributed > &dis_ie, bool do_check=false) const
void locate(const disarray< point_basic< T >, sequential > &x, disarray< size_type, sequential > &dis_ie, bool do_check=false) const
sequential mesh representation
#define trace_macro(message)
#define error_macro(message)
#define fatal_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
bool contains(const geo_element &K, const disarray< point_basic< T >, M > &node, const point_basic< T > &x)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
This file is part of Rheolef.
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
void compute_bbox(const geo_base_rep< T, M > &omega, const geo_element &K, point_basic< T > &xmin, point_basic< T > &xmax)
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
t operator()(const t &a, const t &b)
std::string ptos(const point_basic< T > &x, int d=3)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
CGAL::Filtered_kernel_adaptor< custom_cgal::kernel_2d< T > > Kernel
CGAL::Filtered_kernel_adaptor< custom_cgal::kernel_2d< T > > Kernel
CGAL::Filtered_kernel_adaptor< custom_cgal::kernel_3d< T > > Kernel