17#ifndef __CLOSEST_POINT_HPP__
18#define __CLOSEST_POINT_HPP__
20#include "Grid/grid_dist_key.hpp"
21#include "algoim_hocp.hpp"
24constexpr int algoim_padding = 4;
34template<
size_t wrapping_field,
typename gr
id_type,
typename wrapping_field_type =
typename boost::mpl::at<
typename gr
id_type::value_type::type,boost::mpl::
int_<wrapping_field>>::type>
43 double operator() (
const blitz::TinyVector<int, dim> idx)
const
45 long int local_key[dim];
49 for (
int d = 0; d < dim; ++d)
50 local_key[d] = idx(d) - algoim_padding;
55 return gd.template get<wrapping_field>(
grid_key);
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) {
61 using Poly =
typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
63 Poly field_poly = Poly(coord, *
this, dx);
65 gd.template get<extend_field_temp>(key) = field_poly(pos);
70template<
size_t wrapping_field,
typename gr
id_type,
typename wrapping_field_type,
size_t N1>
80 double operator() (
const blitz::TinyVector<int, dim> idx)
const
82 long int local_key[dim];
86 for (
int d = 0; d < dim; ++d)
87 local_key[d] = idx(d) - algoim_padding;
92 return gd.template get<wrapping_field>(
grid_key)[comp_i];
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) {
98 using Poly =
typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
100 for (
int i = 0; i < N1; ++i) {
102 Poly field_poly = Poly(coord, *
this, dx);
104 gd.template get<extend_field_temp>(key)[i] = field_poly(pos);
109template<
size_t wrapping_field,
typename gr
id_type,
typename wrapping_field_type,
size_t N1,
size_t N2>
115 size_t comp_i, comp_j;
119 double operator() (
const blitz::TinyVector<int, dim> idx)
const
121 long int local_key[dim];
125 for (
int d = 0; d < dim; ++d)
126 local_key[d] = idx(d) - algoim_padding;
131 return gd.template get<wrapping_field>(
grid_key)[comp_i][comp_j];
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) {
137 using Poly =
typename Algoim::StencilPoly<grid_type::dims, poly_order>::T_Poly;
139 for (
int i = 0; i < N1; ++i) {
140 for (
int j = 0; j < N2; ++j) {
143 Poly field_poly = Poly(coord, *
this, dx);
145 gd.template get<extend_field_temp>(key)[i][j] = field_poly(pos);
151template<
size_t wrapping_field,
typename gr
id_type,
typename wrapping_field_type,
size_t N1,
size_t N2,
size_t N3>
157 size_t comp_i, comp_j, comp_k;
161 double operator() (
const blitz::TinyVector<int, dim> idx)
const
163 long int local_key[dim];
167 for (
int d = 0; d < dim; ++d)
168 local_key[d] = idx(d) - algoim_padding;
173 return gd.template get<wrapping_field>(
grid_key)[comp_i][comp_j][comp_k];
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) {
179 using Poly =
typename Algoim::StencilPoly<grid_type::dims, poly_order>::T_Poly;
181 for (
int i = 0; i < N1; ++i) {
182 for (
int j = 0; j < N2; ++j) {
183 for (
int k = 0; k < N3; ++k) {
187 Poly field_poly = Poly(coord, *
this, dx);
189 gd.template get<extend_field_temp>(key)[i][j][k] = field_poly(pos);
208template<
size_t phi_field,
size_t cp_field,
int poly_order,
typename gr
id_type>
213 gd.template ghost_get<phi_field>(KEEP_PROPERTIES);
216 using Poly =
typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
219 blitz::TinyVector<double,dim> dx;
220 for(
int d = 0; d < dim; ++d)
228 for(
int i = 0; i < patches.size();i++)
230 for(
int d = 0; d < dim; ++d)
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]);
239 std::vector<Algoim::detail::CellPoly<dim,Poly>> cells;
241 blitz::TinyVector<int,dim> ext;
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);
246 Algoim::detail::createCellPolynomials(ext, phiwrap, dx,
false, cells);
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);
252 Algoim::KDTree<double,dim> kdtree(points);
255 double nb_gamma_plus_dx = nb_gamma + gd.
spacing(0);
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(),
259 Algoim::sqr(std::max(1.0e-14, std::pow(blitz::max(dx), Poly::order))),
260 cells, kdtree, points, pointcells, dx, 0.0);
266 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
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);
274 if (hocp.compute(patch_pos, cp))
276 for(
int d = 0; d < dim; ++d)
277 gd.template get<cp_field>(key)[d] = cp(d);
281 std::cout<<
"WARN: Closest point computation fails at : ";
282 for(
int d = 0; d < dim; ++d)
284 std::cout<<key_g.get(d)<<
" ";
285 gd.template get<cp_field>(key)[d] = -100.0;
308template<
size_t phi_field,
size_t cp_field,
size_t extend_field,
size_t extend_field_temp,
int poly_order,
typename gr
id_type>
313 gd.template ghost_get<phi_field, cp_field, extend_field>(KEEP_PROPERTIES);
316 using Poly =
typename Algoim::StencilPoly<dim, poly_order>::T_Poly;
318 blitz::TinyVector<double,dim> dx;
319 for(
int d = 0; d < dim; ++d)
325 for(
int i = 0; i < patches.size();i++)
327 for(
int d = 0; d < dim; ++d)
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]);
338 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
340 blitz::TinyVector<int,dim> coord;
341 blitz::TinyVector<double,dim> pos;
343 for(
int d = 0; d < dim; ++d)
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);
351 fieldwrap.template extend<extend_field_temp,poly_order>(coord,dx,pos,key);
361 typedef typename boost::mpl::at<typename grid_type::value_type::type,boost::mpl::int_<extend_field>>::type type_to_copy;
366 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
381template<
size_t phi_field,
size_t cp_field,
typename gr
id_type>
387 gd.template ghost_get<cp_field>(KEEP_PROPERTIES);
390 blitz::TinyVector<double,dim> dx;
391 for(
int d = 0; d < dim; ++d)
397 for(
int i = 0; i < patches.size();i++)
399 for(
int d = 0; d < dim; ++d)
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]);
410 if(std::abs(gd.template get<phi_field>(key)) < nb_gamma)
413 double sign_fn = (gd.template get<phi_field>(key) >= 0.0)?1.0:-1.0;
417 double distance = 0.0;
418 for(
int d = 0; d < dim; ++d)
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];
424 std::cout<<
"WARNING: Requesting closest point on nodes where it was not computed."<<std::endl;
426 distance += ((patch_pos - cp_d)*(patch_pos - cp_d));
428 distance = sqrt(distance);
430 gd.template get<phi_field>(key) = sign_fn*distance;
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
__device__ __host__ void set_d(index_type i, index_type id)
Set the i index.
__device__ __host__ index_type get(index_type i) const
Get the i index.
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.