#ifndef ALGOIM_HOCP_HPP #define ALGOIM_HOCP_HPP /* Contains the main drivers for the high-order closest point algorithm, mainly for rectangular Cartesian grids, as described in the paper R. I. Saye, High-order methods for computing distances to implicitly defined surfaces, Communications in Applied Mathematics and Computational Science, 9(1), 107-141 (2014), http://dx.doi.org/10.2140/camcos.2014.9.107 */ #include #include "multiloop.hpp" #include "kdtree.hpp" #include "newtoncp.hpp" #include "stencilpoly.hpp" #include "utility.hpp" namespace algoim { namespace detail { // A CellPoly struct contains the interpolated polynomial for a particular grid cell template struct CellPoly { Poly poly; // The polynomial uvector i; // Coordinates of the cell CellPoly(const Poly& poly, const uvector& i) : poly(poly), i(i) {} }; // Given a Cartesian grid-defined level set function phi, interpolate phi on each relevant grid cell to // create a collection of polynomials. If interpAllCells is true, then every grid cell in the domain // is interpolated, otherwise only those cells for which phi changes sign on the vertices are considered. template void createCellPolynomials(const uvector& ext, const Phi& phi, const uvector& dx, bool interpAllCells, std::vector>& cells) { if (interpAllCells) { for (MultiLoop i(0, ext - 1); ~i; ++i) cells.push_back(CellPoly(Poly(i(), phi, dx), i())); } else { for (MultiLoop i(0, ext - 1); ~i; ++i) { int s = 0; for (MultiLoop j(i(), i() + 2); ~j; ++j) s += (phi(j()) >= 0.0) ? 1 : -1; if (std::abs(s) < 1 << N) cells.push_back(CellPoly(Poly(i(), phi, dx), i())); } } } // Given a collection of cell-generated polynomials, sample their zero level set to generate the // seed points in the cloud, where subcellExt is the number of subcells to use per cell. The output // is a collection of points and the corresponding index of the polynomial from which they were generated. template void samplePolynomials(const std::vector>& cells, int subcellExt, const uvector& dx, const uvector& xmin, std::vector>& points, std::vector& pointcells) { points.clear(); pointcells.clear(); // subcellsqrradius is the induced squared radius of each subcell (used to discard points if they go too far) double subcellsqrradius = util::sqr(1.25*norm(dx)*0.5/subcellExt); // tolerance 1% of subcell size double isotolsqr = std::max(1.0e-25, 1e-4*subcellsqrradius); for (size_t i = 0, len = cells.size(); i < len; ++i) { for (MultiLoop j(0, subcellExt); ~j; ++j) { uvector x = (j() + 0.5)*dx/subcellExt; uvector y = x; if (newtonIso(x, cells[i].poly, isotolsqr, 10) > 0 && sqrnorm(x - y) < subcellsqrradius) { points.push_back(cells[i].i*dx + xmin + x); pointcells.push_back(static_cast(i)); } } } } } /* The ComputeHighOrderCP struct is the main engine for computing high-order accurate closest points. It accepts as parameters: - bandrsqr: squared radius of the narrow band - overlapr: radius of the bounding ball used in Newton's method - cptolsqr: tolerance used to determine convergence in Newton's method - cells: collection of polynomials - kdtree: tree containing the seed points - points: collection of seed points - pointcells: index of the polynomial corresponding to each seed point - dx: grid cell size - xmin: coordinate of lower-left vertex of cell 0 */ template struct ComputeHighOrderCP { const double bandrsqr, overlapr, cptolsqr; const std::vector>& cells; const KDTree& kdtree; const std::vector>& points; const std::vector& pointcells; const uvector dx; const uvector xmin; ComputeHighOrderCP(double bandrsqr, double overlapr, double cptolsqr, const std::vector>& cells, const KDTree& kdtree, const std::vector>& points, const std::vector& pointcells, const uvector& dx, const uvector& xmin) : bandrsqr(bandrsqr), overlapr(overlapr), cptolsqr(cptolsqr), cells(cells), kdtree(kdtree), points(points), pointcells(pointcells), dx(dx), xmin(xmin) {} // Given x, compute the closest point to x; returns false if the closest seed point is outside // the narrow band bool compute(const uvector& x, uvector& cp, double* signedDist = 0) const { // Find the closest point in the cloud int index = kdtree.nearest(x, bandrsqr); // If there is no point within the defined band radius, return false if (index < 0) return false; // Fetch the polynomial corresponding to the closest point in the cloud and calculate the reference point // and initial guess of the closest point in the local coordinate system const detail::CellPoly& cell = cells[pointcells[index]]; uvector ref = x - (cell.i*dx + xmin); cp = points[index] - (cell.i*dx + xmin); // Execute Newton's method newtonCP(cp, ref, cell.poly, overlapr, cptolsqr, 20); // If the signed distance is also requested, the sign can (usually) be determined // according to alignment of the normal and polynomial gradient at the closest point if (signedDist) { uvector grad = cell.poly.grad(cp); if (dot(ref - cp, grad) >= 0.0) *signedDist = norm(ref - x); else *signedDist = -norm(ref - x); } // Adjust the closest point such that it is in the global coordinate system cp += cell.i*dx + xmin; return true; } }; // Reinitialiase a Cartesian grid-defined level set function phi in a narrow band of radius r using the // high-order closest point algorithm. The template parameter Degree specifies the type of // polynomial to use. The given phi should implement operator()(uvector i) accepting a grid // coordinate 'i', where 0 <= i < extent template void reinit(Phi& phi, const uvector& extent, const uvector& dx, double r) { // Determine the type of polynomial to use based on given Degree and dimension N typedef typename StencilPoly::T_Poly Poly; // In this simple implementation of a reinitialisation algorithm, boundaries of phi // are treated with a simple type of constant extrapolation, similar to how Neumann // boundary conditions are specified with ghost layers. In other words, if phi(i) // is requested and i is outside the domain, then the nearest grid point to i is // used to define phi(i). auto phiextrap = [&](uvector i) { for (int dim = 0; dim < N; ++dim) { if (i(dim) < 0) i(dim) = 0; else if (i(dim) >= extent(dim)) i(dim) = extent(dim) - 1; } return phi(i); }; // Find all cells containing the interface and construct the high-order polynomials std::vector> cells; detail::createCellPolynomials(extent, phiextrap, dx, false, cells); // Using the polynomials, sample the zero level set in each cell to create a cloud of seed points std::vector> points; std::vector pointcells; detail::samplePolynomials(cells, 2, dx, 0.0, points, pointcells); // Construct a k-d tree from the seed points KDTree kdtree(points); // Pass everything to the closest point computation engine ComputeHighOrderCP hocp(r < std::numeric_limits::max() ? r*r : std::numeric_limits::max(), // squared bandradius 0.5*max(dx), // amount that each polynomial overlaps / size of the bounding ball in Newton's method util::sqr(std::max(1.0e-14, std::pow(max(dx), Poly::order))), // tolerance to determine convergence cells, kdtree, points, pointcells, dx, 0.0); // For each grid point in the domain, compute the closest point and use this to re-define phi as the signed distance function for (MultiLoop i(0, extent); ~i; ++i) { uvector x = i()*dx, cp; if (hocp.compute(x, cp)) phi(i()) = (phi(i()) >= 0.0)? norm(x - cp) : -norm(x - cp); else phi(i()) = (phi(i()) >= 0.0)? std::numeric_limits::max() : -std::numeric_limits::max(); } } } // namespace algoim #endif