49#ifdef _RHEOLEF_HAVE_CGAL
53template <
class T,
class M>
60template <
class T,
class M>
67 if (_ptr == 0) { _ptr = make_ptr(omega); }
68 return _ptr->seq_nearest (omega, x, x_nearest);
70template <
class T,
class M>
77 if (_ptr == 0) { _ptr = make_ptr(omega); }
78 return _ptr->dis_nearest (omega, x, x_nearest);
83template <
class T,
class M>
84class geo_nearest_abstract_rep {
87 virtual ~geo_nearest_abstract_rep() {}
88 virtual void initialize (
const geo_base_rep<T,M>& omega)
const = 0;
90 const geo_base_rep<T,M>& omega,
91 const point_basic<T>& x,
92 point_basic<T>& x_nearest)
const = 0;
94 const geo_base_rep<T,M>& omega,
95 const point_basic<T>& x,
96 point_basic<T>& x_nearest)
const = 0;
103template <
class T,
class M,
size_t D>
104struct geo_nearest_rep :
public geo_nearest_abstract_rep<T,M> { };
109template <
class T,
class M>
110class geo_nearest_rep<
T,
M,2> :
public geo_nearest_abstract_rep<T,M> {
119 geo_nearest_rep() : _point_set(), _ball() {}
120 geo_nearest_rep(
const geo_base_rep<T,M>& omega) : _point_set(), _ball() {
initialize(omega); }
121 ~geo_nearest_rep() {}
122 void initialize (
const geo_base_rep<T,M>& omega)
const;
127 const geo_base_rep<T,M>& omega,
128 const point_basic<T>& x,
129 point_basic<T>& x_nearest)
const;
131 const geo_base_rep<T,M>& omega,
132 const point_basic<T>& x,
133 point_basic<T>& x_nearest)
const;
136 typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
137 typedef typename Kernel::Point_2 Point_2;
140 mutable CGAL::Point_set_2<Kernel> _point_set;
141 mutable disarray<index_set,M> _ball;
143template<
class T,
class M>
147 const geo_base_rep<T,M>& omega,
148 disarray<index_set,M>& ball)
150 typedef typename geo_basic<T,M>::size_type
size_type;
153 if (map_dim <= 0)
return;
155 distributor vertex_ownership = omega.sizes().ownership_by_dimension[0];
156 ball.resize (vertex_ownership, emptyset);
157 for (geo::const_iterator iter = omega.begin(map_dim), last = omega.end(map_dim); iter != last; ++iter) {
158 const geo_element& K = *iter;
159 index_set dis_ie_set;
160 dis_ie_set += K.dis_ie();
161 for (
size_type loc_inod = 0, loc_nnod = K.size(); loc_inod < loc_nnod; loc_inod++) {
163 size_type dis_iv = omega.dis_inod2dis_iv (dis_inod);
164 ball.dis_entry (dis_iv) += dis_ie_set;
167 ball.dis_entry_assembly();
169template <
class T,
class M>
171geo_nearest_rep<T,M,2>::initialize(
const geo_base_rep<T,M>& omega)
const
174 distributor vertex_ownership = omega.sizes().ownership_by_dimension[0];
175 distributor side_ownership = omega.sizes().ownership_by_dimension[
map_dim-1];
176 size_type first_dis_iv = vertex_ownership.first_index();
177 size_type first_dis_isid = side_ownership.first_index();
178 std::vector<bool> marked (omega.n_vertex(),
false);
183 check_macro (omega.have_domain_indirect (
"boundary"),
"nearest: geo may have boundary domain defined");
184 const domain_indirect_basic<M>&
boundary = omega.get_domain_indirect (
"boundary");
185 std::list<Point_2> Lr;
186 for (
typename domain_indirect_basic<M>::const_iterator_ioige iter =
boundary.ioige_begin(), last =
boundary.ioige_end(); iter != last; ++iter) {
188 const geo_element& S = omega.dis_get_geo_element(map_dim-1,dis_isid);
189 for (
size_type iloc = 0, nloc = S.size(); iloc < nloc; ++iloc) {
191 if (!vertex_ownership.is_owned(dis_iv))
continue;
193 if (marked[iv])
continue;
194 const point& xi = omega.node(iv);
195 Point_2 pi (xi[0],xi[1]);
200 _point_set.insert (Lr.begin(),Lr.end());
201 build_ball (omega, _ball);
202 omega.neighbour_guard();
204template <
class T,
class M>
205typename geo_nearest_rep<T,M,2>::size_type
206geo_nearest_rep<T,M,2>::seq_nearest (
207 const geo_base_rep<T,M>& omega,
208 const point_basic<T>& x,
209 point_basic<T>& x_nearest)
const
214 Point_2 x_cgal (x[0],x[1]);
215 typename CGAL::Point_set_2<Kernel>::Vertex_handle v_nearest = _point_set.nearest_neighbor(x_cgal);
216 point xv_nearest (v_nearest->point().x(), v_nearest->point().y());
217 x_nearest = xv_nearest;
225#ifdef _RHEOLEF_HAVE_MPI
230 dis_ie0 = omega.seq_locate (xv_nearest);
233 dis_ie0 = omega.seq_locate (xv_nearest);
235 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
"invalid element containing nearest point");
236 const geo_element& K0 = omega.dis_get_geo_element (map_dim, dis_ie0);
237 size_type dis_iv_nearest = std::numeric_limits<size_type>::max();
238 for (
size_type iloc = 0, nloc = K0.size(); iloc < nloc; ++iloc) {
239 if (
dist(omega.node(K0[iloc]), xv_nearest) == 0) {
240 dis_iv_nearest = K0[iloc];
247 index_set ball_nearest = _ball.dis_at (dis_iv_nearest);
248 check_macro (map_dim > 1,
"invalid map_dim="<<map_dim);
251 distributor element_ownership = omega.sizes().ownership_by_dimension[
map_dim];
252 distributor side_ownership = omega.sizes().ownership_by_dimension[
map_dim-1];
253 for (index_set::const_iterator iter = ball_nearest.begin(), last = ball_nearest.end(); iter != last; ++iter) {
255 const geo_element& K = omega.dis_get_geo_element (map_dim, dis_ie);
256 for (
size_type iloc_sid = 0, nloc_sid = K.n_subgeo(map_dim-1); iloc_sid < nloc_sid; ++iloc_sid) {
258 if (!side_ownership.is_owned(dis_isid))
continue;
259 const geo_element& S = omega.dis_get_geo_element (map_dim-1, dis_isid);
260 if (S.master(1) != std::numeric_limits<size_type>::max())
continue;
264 const point&
a = omega.node(S[0]);
265 const point&
b = omega.node(S[1]);
270 Float f0 = t[0]*x[0] + t[1]*x[1];
273 Float y0 = (
n[1]*
f[0] - t[1]*
f[1])/det;
274 Float y1 = - (
n[0]*
f[0] - t[0]*
f[1])/det;
276 check_macro (fabs(
dot(y-a,n)) + fabs(
dot(y-x,t)) < std::numeric_limits<Float>::epsilon(),
"intersection problem");
277 if (
dot(y-a,t) < 0) y =
a;
278 else if (
dot(y-b,t) > 0) y =
b;
280 if (
d >= d_nearest)
continue;
283 dis_ie_nearest = dis_ie;
289 trace_macro (
"point " << x <<
" nearest_point " << x_nearest);
290 return dis_ie_nearest;
292template <
class T,
class M>
293typename geo_nearest_rep<T,M,2>::size_type
294geo_nearest_rep<T,M,2>::dis_nearest (
295 const geo_base_rep<T,M>& omega,
296 const point_basic<T>& x,
297 point_basic<T>& x_nearest)
const
299 size_type dis_ie = seq_nearest (omega, x, x_nearest);
300#ifdef _RHEOLEF_HAVE_MPI
315template <
class T,
class M>
316geo_nearest_abstract_rep<T,M>*
321 case 2:
return new_macro((geo_nearest_rep<T,M,2>)(omega));
329template <
class T,
class M>
335 return _nearestor.seq_nearest (*
this, x, x_nearest);
337template <
class T,
class M>
343 return _nearestor.dis_nearest (*
this, x, x_nearest);
345template <
class T,
class M>
355 dis_ie.resize (x.ownership());
356 x_nearest.resize (x.ownership());
357 for (
size_type i = 0, n = x.size(); i < n; ++i) {
358 dis_ie[i] = omega.
seq_locate (x[i], dis_ie[i]);
359 if (dis_ie[i] != std::numeric_limits<size_type>::max()) {
363 dis_ie[i] = omega.
seq_nearest (x[i], x_nearest[i]);
364 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
375 sequential_nearest (*
this, x, x_nearest, dis_ie);
377#ifdef _RHEOLEF_HAVE_MPI
387 sequential_nearest (*
this, x, x_nearest, dis_ie);
401template <
class T,
class M>
405template <
class T,
class M>
408 const geo_base_rep<T,M>& omega,
409 const point_basic<T>& x,
410 point_basic<T>& x_nearest)
const
412 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
415template <
class T,
class M>
418 const geo_base_rep<T,M>& omega,
419 const point_basic<T>& x,
420 point_basic<T>& x_nearest)
const
422 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
425template <
class T,
class M>
428 const point_basic<T>& x,
429 point_basic<T>& x_nearest)
const
431 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
434template <
class T,
class M>
437 const point_basic<T>& x,
438 point_basic<T>& x_nearest)
const
440 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
445geo_rep<T,sequential>::nearest (
446 const disarray<point_basic<T>, sequential>& x,
447 disarray<point_basic<T>, sequential>& x_nearest,
448 disarray<
typename geo_base_rep<T,sequential>::size_type, sequential>& idx)
const
450 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
452#ifdef _RHEOLEF_HAVE_MPI
455geo_rep<T,distributed>::nearest (
456 const disarray<point_basic<T>, distributed>& x,
457 disarray<point_basic<T>, distributed>& x_nearest,
458 disarray<
typename geo_base_rep<T,distributed>::size_type, distributed>& idx)
const
460 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
467#define _RHEOLEF_instanciation(T,M) \
468template class geo_nearest<Float,M>; \
469template class geo_base_rep<Float,M>; \
470template class geo_rep<Float,M>;
473#ifdef _RHEOLEF_HAVE_MPI