OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
closest_point.hpp
Go to the documentation of this file.
1
17#ifndef __CLOSEST_POINT_HPP__
18#define __CLOSEST_POINT_HPP__
19
20#include "Grid/grid_dist_key.hpp"
21#include "algoim_hocp.hpp"
22
23// Width of extra padding around each grid patch needed to correctly construct kDTree in Algoim.
24constexpr int algoim_padding = 4;
25
34template<size_t wrapping_field, typename grid_type, typename wrapping_field_type = typename boost::mpl::at<typename grid_type::value_type::type,boost::mpl::int_<wrapping_field>>::type>
36{
37 const static unsigned int dim = grid_type::dims;
38 grid_type &gd;
39 int patch_id;
40 AlgoimWrapper(grid_type& ls_grid, const int pid) : gd(ls_grid), patch_id(pid) {}
41
43 double operator() (const blitz::TinyVector<int, dim> idx) const
44 {
45 long int local_key[dim];
46
47 auto ghost_offset = gd.getLocalGridsInfo().get(patch_id).Dbox.getKP1();
48
49 for (int d = 0; d < dim; ++d)
50 local_key[d] = idx(d) - algoim_padding;
51
52 // Generate OpenFPM grid_key object from local grid indices
53 grid_dist_key_dx<dim> grid_key(patch_id, grid_key_dx<dim> (local_key) + ghost_offset);
54
55 return gd.template get<wrapping_field>(grid_key);
56 }
57
58 template<size_t extend_field_temp, int poly_order, typename coord_type, typename dx_type, typename pos_type, typename key_type>
59 void extend(coord_type coord, dx_type dx, pos_type pos, key_type key) {
60
61 using Poly = typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
62
63 Poly field_poly = Poly(coord, *this, dx);
64 // Extension is first done to the temporary field. Otherwise interpolation will be affected.
65 gd.template get<extend_field_temp>(key) = field_poly(pos);
66 }
67
68};
69
70template<size_t wrapping_field, typename grid_type, typename wrapping_field_type, size_t N1>
71struct AlgoimWrapper<wrapping_field,grid_type,wrapping_field_type[N1]>
72{
73 const static unsigned int dim = grid_type::dims;
74 grid_type &gd;
75 int patch_id;
76 size_t comp_i;
77 AlgoimWrapper(grid_type& ls_grid, const int pid) : gd(ls_grid), patch_id(pid) {}
78
80 double operator() (const blitz::TinyVector<int, dim> idx) const
81 {
82 long int local_key[dim];
83
84 auto ghost_offset = gd.getLocalGridsInfo().get(patch_id).Dbox.getKP1();
85
86 for (int d = 0; d < dim; ++d)
87 local_key[d] = idx(d) - algoim_padding;
88
89 // Generate OpenFPM grid_key object from local grid indices
90 grid_dist_key_dx<dim> grid_key(patch_id, grid_key_dx<dim> (local_key) + ghost_offset);
91
92 return gd.template get<wrapping_field>(grid_key)[comp_i];
93 }
94
95 template<size_t extend_field_temp, int poly_order, typename coord_type, typename dx_type, typename pos_type, typename key_type>
96 void extend(coord_type coord, dx_type dx, pos_type pos, key_type key) {
97
98 using Poly = typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
99
100 for (int i = 0; i < N1; ++i) {
101 comp_i = i;
102 Poly field_poly = Poly(coord, *this, dx);
103 // Extension is first done to the temporary field. Otherwise interpolation will be affected.
104 gd.template get<extend_field_temp>(key)[i] = field_poly(pos);
105 }
106 }
107};
108
109template<size_t wrapping_field, typename grid_type, typename wrapping_field_type, size_t N1, size_t N2>
110struct AlgoimWrapper<wrapping_field,grid_type,wrapping_field_type[N1][N2]>
111{
112 const static unsigned int dim = grid_type::dims;
113 grid_type &gd;
114 int patch_id;
115 size_t comp_i, comp_j;
116 AlgoimWrapper(grid_type& ls_grid, const int pid) : gd(ls_grid), patch_id(pid) {}
117
119 double operator() (const blitz::TinyVector<int, dim> idx) const
120 {
121 long int local_key[dim];
122
123 auto ghost_offset = gd.getLocalGridsInfo().get(patch_id).Dbox.getKP1();
124
125 for (int d = 0; d < dim; ++d)
126 local_key[d] = idx(d) - algoim_padding;
127
128 // Generate OpenFPM grid_key object from local grid indices
129 grid_dist_key_dx<dim> grid_key(patch_id, grid_key_dx<dim> (local_key) + ghost_offset);
130
131 return gd.template get<wrapping_field>(grid_key)[comp_i][comp_j];
132 }
133
134 template<size_t extend_field_temp, int poly_order, typename coord_type, typename dx_type, typename pos_type, typename key_type>
135 void extend(coord_type coord, dx_type dx, pos_type pos, key_type key) {
136
137 using Poly = typename Algoim::StencilPoly<grid_type::dims, poly_order>::T_Poly;
138
139 for (int i = 0; i < N1; ++i) {
140 for (int j = 0; j < N2; ++j) {
141 comp_i = i;
142 comp_j = j;
143 Poly field_poly = Poly(coord, *this, dx);
144 // Extension is first done to the temporary field. Otherwise interpolation will be affected.
145 gd.template get<extend_field_temp>(key)[i][j] = field_poly(pos);
146 }
147 }
148 }
149};
150
151template<size_t wrapping_field, typename grid_type, typename wrapping_field_type, size_t N1, size_t N2, size_t N3>
152struct AlgoimWrapper<wrapping_field,grid_type,wrapping_field_type[N1][N2][N3]>
153{
154 const static unsigned int dim = grid_type::dims;
155 grid_type &gd;
156 int patch_id;
157 size_t comp_i, comp_j, comp_k;
158 AlgoimWrapper(grid_type& ls_grid, const int pid) : gd(ls_grid), patch_id(pid) {}
159
161 double operator() (const blitz::TinyVector<int, dim> idx) const
162 {
163 long int local_key[dim];
164
165 auto ghost_offset = gd.getLocalGridsInfo().get(patch_id).Dbox.getKP1();
166
167 for (int d = 0; d < dim; ++d)
168 local_key[d] = idx(d) - algoim_padding;
169
170 // Generate OpenFPM grid_key object from local grid indices
171 grid_dist_key_dx<dim> grid_key(patch_id, grid_key_dx<dim> (local_key) + ghost_offset);
172
173 return gd.template get<wrapping_field>(grid_key)[comp_i][comp_j][comp_k];
174 }
175
176 template<size_t extend_field_temp, int poly_order, typename coord_type, typename dx_type, typename pos_type, typename key_type>
177 void extend(coord_type coord, dx_type dx, pos_type pos, key_type key) {
178
179 using Poly = typename Algoim::StencilPoly<grid_type::dims, poly_order>::T_Poly;
180
181 for (int i = 0; i < N1; ++i) {
182 for (int j = 0; j < N2; ++j) {
183 for (int k = 0; k < N3; ++k) {
184 comp_i = i;
185 comp_j = j;
186 comp_k = k;
187 Poly field_poly = Poly(coord, *this, dx);
188 // Extension is first done to the temporary field. Otherwise interpolation will be affected.
189 gd.template get<extend_field_temp>(key)[i][j][k] = field_poly(pos);
190 }
191 }
192 }
193 }
194};
195
196
208template<size_t phi_field, size_t cp_field, int poly_order, typename grid_type>
209void estimateClosestPoint(grid_type &gd, const double nb_gamma)
210{
211 const unsigned int dim = grid_type::dims;
212 // Update the phi field in ghosts
213 gd.template ghost_get<phi_field>(KEEP_PROPERTIES);
214
215 // Stencil polynomial type
216 using Poly = typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
217
218 // Grid spacing along each dimension
219 blitz::TinyVector<double,dim> dx;
220 for(int d = 0; d < dim; ++d)
221 dx(d) = gd.spacing(d);
222
223 auto &patches = gd.getLocalGridsInfo();
224
225 grid_key_dx<dim> p_lo;
226 grid_key_dx<dim> p_hi;
227
228 for(int i = 0; i < patches.size();i++)
229 {
230 for(int d = 0; d < dim; ++d)
231 {
232 p_lo.set_d(d, patches.get(i).Dbox.getLow(d) + patches.get(i).origin[d]);
233 p_hi.set_d(d, patches.get(i).Dbox.getHigh(d) + patches.get(i).origin[d]);
234 }
235
237
238 // Find all cells containing the interface and construct the high-order polynomials
239 std::vector<Algoim::detail::CellPoly<dim,Poly>> cells;
240
241 blitz::TinyVector<int,dim> ext;
242
243 for(int d = 0; d < dim; ++d)
244 ext(d) = static_cast<int>(p_hi.get(d) - p_lo.get(d) + 1 + 2*algoim_padding);
245
246 Algoim::detail::createCellPolynomials(ext, phiwrap, dx, false, cells);
247
248 std::vector<blitz::TinyVector<double,dim>> points;
249 std::vector<int> pointcells;
250 Algoim::detail::samplePolynomials<dim,Poly>(cells, 2, dx, 0.0, points, pointcells);
251
252 Algoim::KDTree<double,dim> kdtree(points);
253
254 // In order to ensure that CP is estimated for all points in the narrowband, we add a buffer to the distance check.
255 double nb_gamma_plus_dx = nb_gamma + gd.spacing(0);
256 // Pass everything to the closest point computation engine
257 Algoim::ComputeHighOrderCP<dim,Poly> hocp(nb_gamma_plus_dx < std::numeric_limits<double>::max() ? nb_gamma_plus_dx*nb_gamma_plus_dx : std::numeric_limits<double>::max(), // squared bandradius
258 0.5*blitz::max(dx), // amount that each polynomial overlaps / size of the bounding ball in Newton's method
259 Algoim::sqr(std::max(1.0e-14, std::pow(blitz::max(dx), Poly::order))), // tolerance to determine convergence
260 cells, kdtree, points, pointcells, dx, 0.0);
261
262 auto it = gd.getSubDomainIterator(p_lo, p_hi);
263 while(it.isNext())
264 {
265 auto key = it.get();
266 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
267 {
268 auto key_g = gd.getGKey(key);
269 // NOTE: This is not the real grid coordinates, but internal coordinates for algoim
270 blitz::TinyVector<double,dim> patch_pos, cp;
271 for(int d = 0; d < dim; ++d)
272 patch_pos(d) = (key_g.get(d) - p_lo.get(d) + algoim_padding) * dx(d);
273
274 if (hocp.compute(patch_pos, cp))
275 {
276 for(int d = 0; d < dim; ++d)
277 gd.template get<cp_field>(key)[d] = cp(d);
278 }
279 else
280 {
281 std::cout<<"WARN: Closest point computation fails at : ";
282 for(int d = 0; d < dim; ++d)
283 {
284 std::cout<<key_g.get(d)<<" ";
285 gd.template get<cp_field>(key)[d] = -100.0;
286 }
287 std::cout<<"\n";
288 }
289 }
290 ++it;
291 }
292 }
293 return;
294}
295
308template<size_t phi_field, size_t cp_field, size_t extend_field, size_t extend_field_temp, int poly_order, typename grid_type>
309void extendLSField(grid_type &gd, const double nb_gamma)
310{
311 const unsigned int dim = grid_type::dims;
312 // Update the phi and cp fields in ghost
313 gd.template ghost_get<phi_field, cp_field, extend_field>(KEEP_PROPERTIES);
314
315 // Stencil polynomial object
316 using Poly = typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
317 auto &patches = gd.getLocalGridsInfo();
318 blitz::TinyVector<double,dim> dx;
319 for(int d = 0; d < dim; ++d)
320 dx(d) = gd.spacing(d);
321
322 grid_key_dx<dim> p_lo;
323 grid_key_dx<dim> p_hi;
324
325 for(int i = 0; i < patches.size();i++)
326 {
327 for(int d = 0; d < dim; ++d)
328 {
329 p_lo.set_d(d, patches.get(i).Dbox.getLow(d) + patches.get(i).origin[d]);
330 p_hi.set_d(d, patches.get(i).Dbox.getHigh(d) + patches.get(i).origin[d]);
331 }
332
333 auto it = gd.getSubDomainIterator(p_lo, p_hi);
334
335 while(it.isNext())
336 {
337 auto key = it.get();
338 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
339 {
340 blitz::TinyVector<int,dim> coord;
341 blitz::TinyVector<double,dim> pos;
342
343 for(int d = 0; d < dim; ++d)
344 {
345 double cp_d = gd.template get<cp_field>(key)[d];
346 coord(d) = static_cast<int>(floor(cp_d / gd.spacing(d)));
347 pos(d) = cp_d - coord(d)*gd.spacing(d);
348 }
349
351 fieldwrap.template extend<extend_field_temp,poly_order>(coord,dx,pos,key);
352 // Poly field_poly = Poly(coord, fieldwrap, dx);
353 // // Extension is first done to the temporary field. Otherwise interpolation will be affected.
354 // gd.template get<extend_field_temp>(key) = field_poly(pos);
355 }
356 ++it;
357 }
358 }
359
360 // Copy the results to the actual variable
361 typedef typename boost::mpl::at<typename grid_type::value_type::type,boost::mpl::int_<extend_field>>::type type_to_copy;
362 auto it = gd.getDomainIterator();
363 while(it.isNext())
364 {
365 auto key = it.get();
366 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
367 meta_copy<type_to_copy>::meta_copy_(gd.template get<extend_field_temp>(key),gd.template get<extend_field>(key));
368 ++it;
369 }
370}
371
381template<size_t phi_field, size_t cp_field, typename grid_type>
382void reinitializeLS(grid_type &gd, const double nb_gamma)
383{
384 const unsigned int dim = grid_type::dims;
385
386 // Update the cp_field in ghost
387 gd.template ghost_get<cp_field>(KEEP_PROPERTIES);
388
389 auto &patches = gd.getLocalGridsInfo();
390 blitz::TinyVector<double,dim> dx;
391 for(int d = 0; d < dim; ++d)
392 dx(d) = gd.spacing(d);
393
394 grid_key_dx<dim> p_lo;
395 grid_key_dx<dim> p_hi;
396
397 for(int i = 0; i < patches.size();i++)
398 {
399 for(int d = 0; d < dim; ++d)
400 {
401 p_lo.set_d(d, patches.get(i).Dbox.getLow(d) + patches.get(i).origin[d]);
402 p_hi.set_d(d, patches.get(i).Dbox.getHigh(d) + patches.get(i).origin[d]);
403 }
404
405 auto it = gd.getSubDomainIterator(p_lo, p_hi);
406
407 while(it.isNext())
408 {
409 auto key = it.get();
410 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
411 {
412 // Preserve the current sign of the SDF
413 double sign_fn = (gd.template get<phi_field>(key) >= 0.0)?1.0:-1.0;
414 auto key_g = gd.getGKey(key);
415
416 // Compute the Euclidean distance from gird coordinate to closest point coordinate
417 double distance = 0.0;
418 for(int d = 0; d < dim; ++d)
419 {
420 // NOTE: This is not the real grid coordinates, but internal coordinates used for algoim
421 double patch_pos = (key_g.get(d) - p_lo.get(d) + algoim_padding) * gd.spacing(d);
422 double cp_d = gd.template get<cp_field>(key)[d];
423 if(cp_d == -100.0)
424 std::cout<<"WARNING: Requesting closest point on nodes where it was not computed."<<std::endl;
425
426 distance += ((patch_pos - cp_d)*(patch_pos - cp_d));
427 }
428 distance = sqrt(distance);
429
430 gd.template get<phi_field>(key) = sign_fn*distance;
431 }
432 ++it;
433 }
434 }
435}
436
437#endif //__CLOSEST_POINT_HPP__
This is a distributed grid.
const openfpm::vector< GBoxes< device_grid::dims > > & getLocalGridsInfo() const
It return the informations about the local grids.
grid_key_dx< dim > getGKey(const grid_dist_key_dx< dim > &k) const
Convert a g_dist_key_dx into a global key.
grid_dist_iterator_sub< dim, device_grid > getSubDomainIterator(const grid_key_dx< dim > &start, const grid_key_dx< dim > &stop) const
It return an iterator that span the grid domain only in the specified part.
St spacing(size_t i) const
Get the spacing of the grid in direction i.
grid_dist_iterator< dim, device_grid, decltype(device_grid::type_of_subiterator()), FREE > getDomainIterator() const
It return an iterator that span the full grid domain (each processor span its local domain)
static const unsigned int dims
Number of dimensions.
Grid key for a distributed grid.
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
__device__ __host__ void set_d(index_type i, index_type id)
Set the i index.
Definition grid_key.hpp:516
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition grid_key.hpp:503
void reinitializeLS(grid_type &gd, const double nb_gamma)
Reinitializes the level set Phi field on a grid. The grid should have level set SDF and closest point...
void estimateClosestPoint(grid_type &gd, const double nb_gamma)
Computes the closest point coordinate for each grid point within nb_gamma from interface.
void extendLSField(grid_type &gd, const double nb_gamma)
Extends a (scalar) field to within nb_gamma from interface. The grid should have level set SDF and cl...
double operator()(const blitz::TinyVector< int, dim > idx) const
Call operator for the wrapper.
__device__ static __host__ void meta_copy_(const T &src, T &dst)
copy and object from src to dst
Definition meta_copy.hpp:60