Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
geo_nearest.cc
Go to the documentation of this file.
1
21//
22// given x, search x* the closest point in the boundary of omega
23// gives also K* in mesh such that x* in K*
24//
25// author: Pierre.Saramito@imag.fr
26//
27// date: 9 october 2014
28//
29// implementation note:
30// use CGAL::Point_set_2 for 2d nearest point search
31// http://doc.cgal.org/latest/Point_set_2/group__PkgPointSet2.html
32// TODO: 3D: with CGAL::Neighbor_search::Tree
33// http://doc.cgal.org/latest/Spatial_searching/index.html#Chapter_dD_Spatial_Searching
34// TODO: distributed: try to store vertex index in Point_2 data structure ? How to do that ?
35// http://doc.cgal.org/latest/Kernel_23/index.html#title20 : extensible kernel
36// http://stackoverflow.com/questions/2418332/customizing-cgal-kernel-with-my-own-point-class
37#include "rheolef/geo_nearest.h"
38#include "rheolef/geo.h"
39#include "rheolef/point_util.h"
40
41// internal includes:
42#ifdef _RHEOLEF_HAVE_CGAL
43#include "rheolef/cgal_traits.h"
44#include <CGAL/Point_set_2.h>
45#endif // _RHEOLEF_HAVE_CGAL
46
47namespace rheolef {
49#ifdef _RHEOLEF_HAVE_CGAL
50// ---------------------------------------------------------------------
51// 1) the geo_nearest interface
52// ---------------------------------------------------------------------
53template <class T, class M>
55{
56 if (_ptr != 0) {
57 delete_macro(_ptr);
58 }
59}
60template <class T, class M>
63 const geo_base_rep<T,M>& omega,
64 const point_basic<T>& x,
65 point_basic<T>& x_nearest) const
66{
67 if (_ptr == 0) { _ptr = make_ptr(omega); }
68 return _ptr->seq_nearest (omega, x, x_nearest);
69}
70template <class T, class M>
73 const geo_base_rep<T,M>& omega,
74 const point_basic<T>& x,
75 point_basic<T>& x_nearest) const
76{
77 if (_ptr == 0) { _ptr = make_ptr(omega); }
78 return _ptr->dis_nearest (omega, x, x_nearest);
79}
80// ---------------------------------------------------------------------
81// 2) the abstract data structure
82// ---------------------------------------------------------------------
83template <class T, class M>
84class geo_nearest_abstract_rep {
85public:
86 typedef typename disarray<T,M>::size_type size_type;
87 virtual ~geo_nearest_abstract_rep() {}
88 virtual void initialize (const geo_base_rep<T,M>& omega) const = 0;
89 virtual size_type seq_nearest (
90 const geo_base_rep<T,M>& omega,
91 const point_basic<T>& x,
92 point_basic<T>& x_nearest) const = 0;
93 virtual size_type dis_nearest (
94 const geo_base_rep<T,M>& omega,
95 const point_basic<T>& x,
96 point_basic<T>& x_nearest) const = 0;
97};
98// ---------------------------------------------------------------------
99// 3) the concrete, and dimension dependent, data structure
100// = cgal "Point_set_2" delaunay mesh for 2D
101// = cgal "Neighbor_search::Tree" octree for 3D
102// ---------------------------------------------------------------------
103template <class T, class M, size_t D>
104struct geo_nearest_rep : public geo_nearest_abstract_rep<T,M> { };
105
106// ---------------------------------------------------------------------
107// 3a) the 2D implementation with cgal::point_set
108// ---------------------------------------------------------------------
109template <class T, class M>
110class geo_nearest_rep<T,M,2> : public geo_nearest_abstract_rep<T,M> {
111public:
112
113// typedef:
114
116
117// allocators:
118
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;
123
124// accessors:
125
126 size_type seq_nearest (
127 const geo_base_rep<T,M>& omega,
128 const point_basic<T>& x,
129 point_basic<T>& x_nearest) const;
130 size_type dis_nearest (
131 const geo_base_rep<T,M>& omega,
132 const point_basic<T>& x,
133 point_basic<T>& x_nearest) const;
134
135protected:
136 typedef typename geo_cgal_traits<T,2>::Kernel Kernel;
137 typedef typename Kernel::Point_2 Point_2;
138
139// data:
140 mutable CGAL::Point_set_2<Kernel> _point_set;
141 mutable disarray<index_set,M> _ball;
142};
143template<class T, class M>
144static
145void
146build_ball (
147 const geo_base_rep<T,M>& omega,
148 disarray<index_set,M>& ball)
149{
150 typedef typename geo_basic<T,M>::size_type size_type;
151 // 1) build ball(xi) := { K; xi is a vertex of K }
152 size_type map_dim = omega.map_dimension();
153 if (map_dim <= 0) return; // a set of points have no neighbours
154 index_set emptyset;
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++) {
162 size_type dis_inod = K [loc_inod];
163 size_type dis_iv = omega.dis_inod2dis_iv (dis_inod);
164 ball.dis_entry (dis_iv) += dis_ie_set; // union with {dis_ie}
165 }
166 }
167 ball.dis_entry_assembly();
168}
169template <class T, class M>
170void
171geo_nearest_rep<T,M,2>::initialize(const geo_base_rep<T,M>& omega) const
172{
173 size_type map_dim = omega.map_dimension();
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);
179#ifdef TODO
180 // TODO: passer boundary_guard de geo_basic a geo_base_rep
181 boundary_guard (*this);
182#endif // TODO
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) {
187 size_type dis_isid = (*iter).index();
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) {
190 size_type dis_iv = S[iloc];
191 if (!vertex_ownership.is_owned(dis_iv)) continue;
192 size_type iv = dis_iv - first_dis_iv;
193 if (marked[iv]) continue;
194 const point& xi = omega.node(iv);
195 Point_2 pi (xi[0],xi[1]);
196 Lr.push_back (pi);
197 marked[iv] = true;
198 }
199 }
200 _point_set.insert (Lr.begin(),Lr.end());
201 build_ball (omega, _ball);
202 omega.neighbour_guard();
203}
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
210{
211 // -----------------------------------------
212 // 1) nearest neighbor : get its coordinates
213 // -----------------------------------------
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;
218 // ------------------------------------------
219 // 2) nearest neighbor : get its vertex index
220 // ------------------------------------------
221 // TODO: try to store vertex index in Point_2 data structure ? How to do that ?
222 // => avoid a geo_locate(x) and facilitate the distributed version impementation
223 size_type map_dim = omega.map_dimension();
224 size_type dis_ie0 = 0;
225#ifdef _RHEOLEF_HAVE_MPI
226 size_type nproc = omega.comm().size();
227 if (is_distributed<M>::value && nproc > 1) {
228 fatal_macro ("dis_nearest/dis_locate: not yet");
229 } else {
230 dis_ie0 = omega.seq_locate (xv_nearest);
231 }
232#else // _RHEOLEF_HAVE_MPI
233 dis_ie0 = omega.seq_locate (xv_nearest);
234#endif // _RHEOLEF_HAVE_MPI
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];
241 break;
242 }
243 }
244 // --------------------------------------------
245 // 3) search in ball(K0) for nearest bdry point
246 // --------------------------------------------
247 index_set ball_nearest = _ball.dis_at (dis_iv_nearest);
248 check_macro (map_dim > 1, "invalid map_dim="<<map_dim);
249 size_type dis_ie_nearest = dis_ie0;
250 Float d_nearest = dist(x,xv_nearest);
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) {
254 size_type dis_ie = *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) {
257 size_type dis_isid = (map_dim == 2) ? K.edge(iloc_sid) : K.face(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; // not a bdry side
261 if (map_dim == 2) {
262 // 1) S=[a,b]: compute y = the intersection D(a,b) with D(x,n)
263 // that satisfies: ay.n = 0 and xy.t = 0
264 const point& a = omega.node(S[0]);
265 const point& b = omega.node(S[1]);
266 point t = b-a;
267 point n (t[1], -t[0]);
268 Float ab2 = dist2(a,b);
269 Float det = -ab2;
270 Float f0 = t[0]*x[0] + t[1]*x[1];
271 Float f1 = n[0]*a[0] + n[1]*a[1];
272 point f (f0,f1);
273 Float y0 = (n[1]*f[0] - t[1]*f[1])/det;
274 Float y1 = - (n[0]*f[0] - t[0]*f[1])/det;
275 point y (y0,y1);
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; // y is outside of segment [a,b], at left
278 else if (dot(y-b,t) > 0) y = b; // y is outside of segment [a,b], at right
279 Float d = dist(x, y);
280 if (d >= d_nearest) continue;
281 x_nearest = y;
282 d_nearest = d;
283 dis_ie_nearest = dis_ie;
284 } else {
285 error_macro ("3D: not yet");
286 }
287 }
288 }
289 trace_macro ("point " << x << " nearest_point " << x_nearest);
290 return dis_ie_nearest;
291}
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
298{
299 size_type dis_ie = seq_nearest (omega, x, x_nearest);
300#ifdef _RHEOLEF_HAVE_MPI
301 size_type nproc = omega.comm().size();
302 if (is_distributed<M>::value && nproc > 1) {
303 fatal_macro ("dis_nearest: not yet");
304 }
305#endif // _RHEOLEF_HAVE_MPI
306 return dis_ie;
307}
308// ---------------------------------------------------------------------
309// 3b) the 3D implementation with cgal::Tree
310// ---------------------------------------------------------------------
311// TODO
312// ---------------------------------------------------------------------
313// 3c) allocator
314// ---------------------------------------------------------------------
315template <class T, class M>
316geo_nearest_abstract_rep<T,M>*
318{
319 check_macro (omega.dimension() == omega.map_dimension(), "geo_nearest: map_dim < dim: not supported");
320 switch (omega.map_dimension()) {
321 case 2: return new_macro((geo_nearest_rep<T,M,2>)(omega));
322 //case 3: return new_macro((geo_nearest_rep<T,M,3>)(omega));
323 default: error_macro ("unsupported dimension d=" << omega.dimension()); return 0;
324 }
325}
326// ---------------------------------------------------------------------
327// 6) geo::nearest()
328// ---------------------------------------------------------------------
329template <class T, class M>
332 const point_basic<T>& x,
333 point_basic<T>& x_nearest) const
334{
335 return _nearestor.seq_nearest (*this, x, x_nearest);
336}
337template <class T, class M>
340 const point_basic<T>& x,
341 point_basic<T>& x_nearest) const
342{
343 return _nearestor.dis_nearest (*this, x, x_nearest);
344}
345template <class T, class M>
346static
347void
348sequential_nearest (
349 const geo_base_rep<T,M>& omega,
350 const disarray<point_basic<T>,M>& x,
351 disarray<point_basic<T>,M>& x_nearest,
353{
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()) {
360 x_nearest[i] = x[i];
361 continue;
362 }
363 dis_ie[i] = omega.seq_nearest (x[i], x_nearest[i]);
364 check_macro (dis_ie[i] != std::numeric_limits<size_type>::max(),
365 "nearest: failed at x="<<ptos(x[i],omega.dimension()));
366 }
367}
368template <class T>
369void
374{
375 sequential_nearest (*this, x, x_nearest, dis_ie);
376}
377#ifdef _RHEOLEF_HAVE_MPI
378template <class T>
379void
384{
385 size_type nproc = base::comm().size();
386 if (nproc == 1) {
387 sequential_nearest (*this, x, x_nearest, dis_ie);
388 return;
389 }
390 // TODO: also search for nearest when we go outside of the domain
391 // TODO: search seq the nearest in the current partition and then
392 // TODO: merge all nearests across partitions ?
393 locate (x, dis_ie);
394 x_nearest = x;
395}
396#endif // _RHEOLEF_HAVE_MPI
397// ----------------------------------------------------------------------------
398// no CGAL: no nearest
399// ----------------------------------------------------------------------------
400#else // ! _RHEOLEF_HAVE_CGAL
401template <class T, class M>
403{
404}
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
411{
412 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
413 return 0;
414}
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
421{
422 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
423 return 0;
424}
425template <class T, class M>
428 const point_basic<T>& x,
429 point_basic<T>& x_nearest) const
430{
431 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
432 return 0;
433}
434template <class T, class M>
437 const point_basic<T>& x,
438 point_basic<T>& x_nearest) const
439{
440 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
441 return 0;
442}
443template <class T>
444void
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
449{
450 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
451}
452#ifdef _RHEOLEF_HAVE_MPI
453template <class T>
454void
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
459{
460 fatal_macro ("geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
461}
462#endif // _RHEOLEF_HAVE_MPI
463#endif // _RHEOLEF_HAVE_CGAL
464// ----------------------------------------------------------------------------
465// instanciation in library
466// ----------------------------------------------------------------------------
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>;
471
473#ifdef _RHEOLEF_HAVE_MPI
474_RHEOLEF_instanciation(Float,distributed)
475#endif // _RHEOLEF_HAVE_MPI
476
477} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the point page for the full documentation
see the disarray page for the full documentation
Definition disarray.h:497
rep::base::size_type size_type
Definition disarray.h:501
base class for M=sequential or distributed meshes representations
Definition geo.h:528
size_type map_dimension() const
Definition geo.h:564
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
Definition geo.h:533
size_type dimension() const
Definition geo.h:563
size_type seq_nearest(const point_basic< T > &x, point_basic< T > &x_nearest) const
size_type dis_nearest(const point_basic< T > &x, point_basic< T > &x_nearest) const
size_type seq_nearest(const geo_base_rep< T, M > &omega, const point_basic< T > &x, point_basic< T > &x_nearest) const
size_type dis_nearest(const geo_base_rep< T, M > &omega, const point_basic< T > &x, point_basic< T > &x_nearest) const
static geo_nearest_abstract_rep< T, M > * make_ptr(const geo_base_rep< T, M > &omega)
disarray< T, M >::size_type size_type
Definition geo_nearest.h:43
base::size_type size_type
Definition geo.h:934
sequential mesh representation
Definition geo.h:778
rheolef::space_base_rep< T, M > t
#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 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 boundary_guard(const geo_basic< T, M > &omega)
Definition geo.cc:600
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
T dist2(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:292
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
T dist(const point_basic< T > &x, const point_basic< T > &y)
Definition point.h:298
Definition cavity_dg.h:29
static const bool value
Expr1::memory_type M