37 #include "rheolef/geo_nearest.h"
38 #include "rheolef/geo.h"
39 #include "rheolef/point_util.h"
42 #ifdef _RHEOLEF_HAVE_CGAL
43 #include "rheolef/cgal_traits.h"
44 #include <CGAL/Point_set_2.h>
49 #ifdef _RHEOLEF_HAVE_CGAL
53 template <
class T,
class M>
60 template <
class T,
class M>
67 if (_ptr == 0) { _ptr = make_ptr(omega); }
68 return _ptr->seq_nearest (omega, x, x_nearest);
70 template <
class T,
class M>
77 if (_ptr == 0) { _ptr = make_ptr(omega); }
78 return _ptr->dis_nearest (omega, x, x_nearest);
83 template <
class T,
class M>
84 class 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,
94 const geo_base_rep<T,M>& omega,
103 template <
class T,
class M,
size_t D>
104 struct geo_nearest_rep :
public geo_nearest_abstract_rep<T,M> { };
109 template <
class T,
class M>
110 class 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,
131 const geo_base_rep<T,M>& omega,
137 typedef typename Kernel::Point_2 Point_2;
140 mutable CGAL::Point_set_2<Kernel> _point_set;
141 mutable disarray<index_set,M> _ball;
143 template<
class T,
class M>
147 const geo_base_rep<T,M>& omega,
148 disarray<index_set,M>& ball)
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++) {
164 ball.dis_entry (dis_iv) += dis_ie_set;
167 ball.dis_entry_assembly();
169 template <
class T,
class M>
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();
204 template <
class T,
class M>
206 geo_nearest_rep<T,M,2>::seq_nearest (
207 const geo_base_rep<T,M>& omega,
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);
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];
274 Float y1 = - (
n[0]*
f[0] -
t[0]*
f[1])/det;
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;
292 template <
class T,
class M>
294 geo_nearest_rep<T,M,2>::dis_nearest (
295 const geo_base_rep<T,M>& omega,
299 size_type dis_ie = seq_nearest (omega, x, x_nearest);
300 #ifdef _RHEOLEF_HAVE_MPI
315 template <
class T,
class M>
316 geo_nearest_abstract_rep<T,M>*
321 case 2:
return new_macro((geo_nearest_rep<T,M,2>)(omega));
329 template <
class T,
class M>
335 return _nearestor.seq_nearest (*
this, x, x_nearest);
337 template <
class T,
class M>
343 return _nearestor.dis_nearest (*
this, x, x_nearest);
345 template <
class T,
class M>
355 dis_ie.resize (x.ownership());
356 x_nearest.resize (x.ownership());
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);
401 template <
class T,
class M>
405 template <
class T,
class M>
408 const geo_base_rep<T,M>& omega,
412 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
415 template <
class T,
class M>
418 const geo_base_rep<T,M>& omega,
422 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
425 template <
class T,
class M>
431 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
434 template <
class T,
class M>
440 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
450 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
452 #ifdef _RHEOLEF_HAVE_MPI
460 fatal_macro (
"geo: nearest not available (HINT: recompile Rheolef with the CGAL library)");
467 #define _RHEOLEF_instanciation(T,M) \
468 template class geo_nearest<Float,M>; \
469 template class geo_base_rep<Float,M>; \
470 template class geo_rep<Float,M>;
473 #ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the point page for the full documentation
see the disarray page for the full documentation
base class for M=sequential or distributed meshes representations
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 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
void nearest(const disarray< point_basic< T >, distributed > &x, disarray< point_basic< T >, distributed > &x_nearest, disarray< size_type, distributed > &dis_ie) const
base::size_type size_type
void nearest(const disarray< point_basic< T >, sequential > &x, disarray< point_basic< T >, sequential > &x_nearest, disarray< size_type, sequential > &dis_ie) const
sequential mesh representation
rheolef::space_base_rep< T, M > t
#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)")
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.
T dist(const point_basic< T > &x, const point_basic< T > &y)
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)
dot(x,y): see the expression page for the full documentation
void boundary_guard(const geo_basic< T, M > &omega)
T dist2(const point_basic< T > &x, const point_basic< T > &y)
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