From a7a929723b7c474b18946925e90c38a3c2a2c23a Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Wed, 4 Nov 2026 22:39:51 -0500 Subject: [PATCH] Factor overlap predicates out of the grid point search --- src/pcms/localization/CMakeLists.txt | 1 + src/pcms/localization/overlap.h | 427 ++++++++++++++++++++++++ src/pcms/localization/point_search.cpp | 400 ++++++++-------------- src/pcms/localization/point_search.h | 98 ++++++ src/pcms/transfer/mesh_intersection.hpp | 62 +--- src/pcms/utility/uniform_grid.h | 33 +- test/CMakeLists.txt | 1 + test/test_cell_range.cpp | 270 +++++++++++++++ test/test_point_search.cpp | 285 ++++++++++++++++ test/test_uniform_grid.cpp | 59 ++++ 10 files changed, 1304 insertions(+), 332 deletions(-) create mode 100644 src/pcms/localization/overlap.h create mode 100644 test/test_cell_range.cpp diff --git a/src/pcms/localization/CMakeLists.txt b/src/pcms/localization/CMakeLists.txt index 269b1b776..f415e054f 100644 --- a/src/pcms/localization/CMakeLists.txt +++ b/src/pcms/localization/CMakeLists.txt @@ -1,6 +1,7 @@ set( PCMS_LOCALIZATION_HEADERS point_search.h + overlap.h queue_visited.hpp adj_search.hpp mls_support_helpers.h diff --git a/src/pcms/localization/overlap.h b/src/pcms/localization/overlap.h new file mode 100644 index 000000000..929c44817 --- /dev/null +++ b/src/pcms/localization/overlap.h @@ -0,0 +1,427 @@ +#ifndef PCMS_LOCALIZATION_OVERLAP_H +#define PCMS_LOCALIZATION_OVERLAP_H +#include +#include +#include +#include +#include +#include + +#include + +#include "pcms/utility/bounding_box.h" +#include "pcms/utility/types.h" +#include "pcms/utility/uniform_grid.h" + +namespace pcms +{ + +/// Predicate tags. Each tag names the operand pair with the exact type words +/// from the original function name (Triangle/Bbox, Simplex/Bbox, +/// Simplex/Simplex). A `_Method` suffix is added only when two tags would +/// otherwise collide. +/// +/// Adding a new predicate is a three-step, local change: define a tag here, +/// declare its padding if it bins by bounding box (see below), and add an +/// `Overlap` specialization. +struct TriangleVsBbox +{ + /// Slack in grid cells that `candidate_cells` must add around the box. + /// + /// The predicate accepts a shared boundary plane as an intersection, so a + /// cell that merely *touches* the box is a candidate. One cell is enough: a + /// cell two or more away cannot touch the box at all. + static constexpr LO boundary_pad = 1; +}; +struct SimplexVsBbox +{ + /// \see TriangleVsBbox::boundary_pad + static constexpr LO boundary_pad = 1; +}; +struct SimplexVsSimplex +{ + // No `boundary_pad` +}; + +/// Inclusive per-axis range of grid cells overlapped by a bounding box. +/// +/// The rules that define what a range *means* live here with the type rather +/// than in a consumer. Note that `CellRange{}` is not the empty range: its +/// zero initialization denotes the single cell at axis index 0, so use +/// `Empty()` when no cell should belong to it. +template +struct CellRange +{ + Kokkos::Array lo; + Kokkos::Array hi; + + /// The range no cell belongs to, as published by a rank with no geometry. + [[nodiscard]] static CellRange Empty() noexcept + { + CellRange range{}; + for (int d = 0; d < Dim; ++d) { + range.lo[d] = 0; + range.hi[d] = -1; + } + return range; + } + + /// True when no cell is in the range, so no query can match it. + [[nodiscard]] KOKKOS_INLINE_FUNCTION bool IsEmpty() const noexcept + { + for (int d = 0; d < Dim; ++d) { + if (hi[d] < lo[d]) { + return true; + } + } + return false; + } + + /// True when the cell at per-axis indices `axis_indices` falls in the range. + [[nodiscard]] KOKKOS_INLINE_FUNCTION bool Contains( + const Kokkos::Array& axis_indices) const noexcept + { + // An empty range has hi < lo on every axis, so this is false for it too. + for (int d = 0; d < Dim; ++d) { + if (axis_indices[d] < lo[d] || axis_indices[d] > hi[d]) { + return false; + } + } + return true; + } +}; + +namespace detail +{ + +KOKKOS_INLINE_FUNCTION +AABBox<2> triangle_bbox(const Omega_h::Matrix<2, 3>& coords) +{ + Kokkos::Array max{coords(0, 0), coords(1, 0)}; + Kokkos::Array min{coords(0, 0), coords(1, 0)}; + for (int i = 1; i < 3; ++i) { + max[0] = std::fmax(max[0], coords(0, i)); + max[1] = std::fmax(max[1], coords(1, i)); + min[0] = std::fmin(min[0], coords(0, i)); + min[1] = std::fmin(min[1], coords(1, i)); + } + return {.center = {(max[0] + min[0]) / 2.0, (max[1] + min[1]) / 2.0}, + .half_width = {(max[0] - min[0]) / 2.0, (max[1] - min[1]) / 2.0}}; +} + +template +KOKKOS_INLINE_FUNCTION AABBox simplex_bbox( + const Omega_h::Matrix& coords) +{ + Kokkos::Array max; + Kokkos::Array min; + for (int j = 0; j < dim; ++j) { + max[j] = coords(j, 0); + min[j] = coords(j, 0); + } + for (int i = 1; i < dim + 1; ++i) { + for (int j = 0; j < dim; ++j) { + max[j] = std::fmax(max[j], coords(j, i)); + min[j] = std::fmin(min[j], coords(j, i)); + } + } + + Kokkos::Array center; + Kokkos::Array half_width; + + for (int j = 0; j < dim; ++j) { + center[j] = (max[j] + min[j]) / 2.0; + half_width[j] = (max[j] - min[j]) / 2.0; + } + + return {.center = center, .half_width = half_width}; +} + +// Liang, You-Dong, and B. A. Barsky. "A New Concept and Method for Line +// Clipping." ACM Transactions on Graphics 3, no. 1 (January 1984): 1-22. +// https://doi.org/10.1145/357332.357333. +KOKKOS_INLINE_FUNCTION +bool clipt(Real p, Real q, Real& t0, Real& t1) +{ + if (p < 0) { + auto r = q / p; + if (r > t1) { + return false; + } + if (r > t0) { + t0 = r; + } + } else if (p > 0) { + auto r = q / p; + if (r < t0) { + return false; + } + if (r < t1) { + t1 = r; + } + } + // p==0 + else { + if (q < 0) { + return false; + } + } + return true; +} + +// Liang, You-Dong, and B. A. Barsky. "A New Concept and Method for Line +// Clipping." ACM Transactions on Graphics 3, no. 1 (January 1984): 1-22. +// https://doi.org/10.1145/357332.357333. +KOKKOS_INLINE_FUNCTION +bool line_intersects_bbox(const Omega_h::Vector<2>& p0, + const Omega_h::Vector<2>& p1, const AABBox<2>& bbox) +{ + Real t0 = 0; + Real t1 = 1; + auto xleft = bbox.center[0] - bbox.half_width[0]; + auto deltax = p1[0] - p0[0]; + if (clipt(-deltax, p0[0] - xleft, t0, t1)) { + auto xright = bbox.center[0] + bbox.half_width[0]; + if (clipt(deltax, xright - p0[0], t0, t1)) { + auto deltay = p1[1] - p0[1]; + auto ybottom = bbox.center[1] - bbox.half_width[1]; + if (clipt(-deltay, p0[1] - ybottom, t0, t1)) { + auto ytop = bbox.center[1] + bbox.half_width[1]; + if (clipt(deltay, ytop - p0[1], t0, t1)) { + // full liang-barksy algorithm computes new x-y coordinates here, + // but we just need to check intersections + return true; + } + } + } + } + return false; +} + +template +[[nodiscard]] KOKKOS_INLINE_FUNCTION bool within_bbox( + const Omega_h::Vector coord, const AABBox& bbox) noexcept +{ + for (int i = 0; i < dim; ++i) { + if (coord[i] < bbox.center[i] - bbox.half_width[i]) + return false; + if (coord[i] > bbox.center[i] + bbox.half_width[i]) + return false; + } + return true; +} + +template +[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_simplex( + const AABBox& bbox, const Omega_h::Matrix& coords) +{ + // each dimension has a pair of opposing "walls" + // 2D: { [left, right], [top, bottom] } -> { left, right, top, bottom } + // 3D: { [left, right], [top, bottom], [front, back] } -> { left, ..., back } + Kokkos::Array bbox_walls{}; + for (int i = 0; i < dim; i++) { + bbox_walls[i * 2] = bbox.center[i] - bbox.half_width[i]; + bbox_walls[i * 2 + 1] = bbox.center[i] + bbox.half_width[i]; + } + + // 1 << dim == 2 ** dim == num vertices in ndim bounding box / hypercube + constexpr unsigned num_verts = 1 << dim; + + // Each vertex is a just a unique combination of walls + // eg [left, bottom] (2D) or [right, top, front] (3) + for (unsigned i = 0; i < num_verts; ++i) { + // conveniently, i acts a bit field representing the current combination + Omega_h::Vector vert; + for (unsigned j = 0; j < dim; ++j) { + // eg 110 = 6 = [left, top, back] + // eg 01 = 1 = [right, top] + vert[j] = (i >> j) & 1 ? bbox_walls[j * 2] : bbox_walls[j * 2 + 1]; + } + auto xi = Omega_h::barycentric_from_global(vert, coords); + if (Omega_h::is_barycentric_inside(xi)) { + return true; + } + } + return false; +} + +} // namespace detail + +/// An overlap predicate, selected as `Overlap` and invoked through +/// its static `apply()`. +/// +/// This is the single entry point for overlap tests: the primary template is +/// intentionally undefined, so each `Dim x Tag` combination that exists has an +/// explicit specialization whose `apply()` holds the original function body +/// verbatim, and asking for a combination that does not exist is a compile +/// error rather than a silent wrong-dimension call. +/// +/// `Tag` names the operand pair (see the tags above); `Dim` fixes the +/// dimension. Because the tag is a template argument rather than a value, a tag +/// that carries no data is never passed as a dummy parameter. +template +struct Overlap; + +template <> +struct Overlap<2, TriangleVsBbox> +{ + /** + * Check if a triangle element represented by 3 coordinates in two dimensions + * intersects with a bounding box + */ + [[nodiscard]] KOKKOS_FUNCTION static bool apply( + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) + { + // triangle and grid cell bounding box intersect + if (intersects(detail::triangle_bbox(coords), bbox)) { + // if any of the triangle verts inside of bbox + if (detail::within_bbox<2>(coords[0], bbox) || + detail::within_bbox<2>(coords[1], bbox) || + detail::within_bbox<2>(coords[2], bbox)) { + return true; + } + // if any of the bbox verts are within the triangle + if (detail::bbox_verts_within_simplex(bbox, coords)) { + return true; + } + // if any of the triangle's edges intersect with the bounding box + if (detail::line_intersects_bbox(coords[0], coords[1], bbox) || + detail::line_intersects_bbox(coords[1], coords[2], bbox) || + detail::line_intersects_bbox(coords[2], coords[0], bbox)) { + return true; + } + } + return false; + } +}; + +template +struct Overlap +{ + [[nodiscard]] KOKKOS_FUNCTION static bool apply( + const Omega_h::Matrix& coords, const AABBox& bbox) + { + return intersects(detail::simplex_bbox(coords), bbox); + // TODO: Add refined cases from triangle_intersects_bbox + } +}; + +template +struct Overlap +{ + /** + * @brief Conservative test for a zero-volume overlap between two simplices. + * + * If verticies of one simplex lies in the closed outer half-space of some + * face plane of the other, their intersection is contained in that plane and + * so has no volume. + * + * @return true when the overlap is provably degenerate and the clip can be + * skipped + */ + [[nodiscard]] OMEGA_H_INLINE static bool apply( + const r3d::Few, Dim + 1>& a, + const r3d::Few, Dim + 1>& b) + { + // The face normal below is a 2D perpendicular or a 3D cross product, so + // this predicate is 2D/3D only even though everything else about it reads + // as dimension-generic. Without this the higher-dimension case is silent: + // at Dim >= 4 `normal[Dim - 1]` is never written and is read as garbage in + // the loop below, and at Dim == 1 writing normal[1] and normal[2] runs off + // the end of the vector. + static_assert( + Dim == 2 || Dim == 3, + "simplex-vs-simplex overlap is implemented for 2D and 3D only"); + + // Each iteration treats one simplex as the half-space owner. Both are + // needed: a face of `a` may separate the pair while no face of `b` does. + for (int side = 0; side < 2; ++side) { + const auto& owner = (side == 0) ? a : b; + const auto& other = (side == 0) ? b : a; + + // Face `f` is the face opposite vertex `f`. + for (int f = 0; f < Dim + 1; ++f) { + // Vertices spanning the face, in index order with `f` skipped. + int face_verts[Dim]; + int n_face_verts = 0; + for (int v = 0; v < Dim + 1; ++v) { + if (v != f) { + face_verts[n_face_verts++] = v; + } + } + const auto& origin = owner[face_verts[0]]; + + // Face normal, oriented to point away from the opposite vertex so the + // owner's interior lies on the negative side. + r3d::Vector normal; + if constexpr (Dim == 2) { + const auto edge = owner[face_verts[1]] - origin; + normal[0] = edge[1]; + normal[1] = -edge[0]; + } else { + const auto e1 = owner[face_verts[1]] - origin; + const auto e2 = owner[face_verts[2]] - origin; + normal[0] = e1[1] * e2[2] - e1[2] * e2[1]; + normal[1] = e1[2] * e2[0] - e1[0] * e2[2]; + normal[2] = e1[0] * e2[1] - e1[1] * e2[0]; + } + + Omega_h::Real opposite_side = 0.0; + Omega_h::Real normal_sq = 0.0; + for (int d = 0; d < Dim; ++d) { + opposite_side += normal[d] * (owner[f][d] - origin[d]); + normal_sq += normal[d] * normal[d]; + } + // A degenerate face gives no usable axis. + if (normal_sq == 0.0) { + continue; + } + const Omega_h::Real orientation = (opposite_side > 0.0) ? -1.0 : 1.0; + + bool all_outside = true; + for (int v = 0; v < Dim + 1 && all_outside; ++v) { + Omega_h::Real signed_distance = 0.0; + for (int d = 0; d < Dim; ++d) { + signed_distance += normal[d] * (other[v][d] - origin[d]); + } + all_outside = (orientation * signed_distance) >= 0.0; + } + if (all_outside) { + return true; + } + } + } + return false; + } +}; + +/// Low-level cell range: cells overlapped by `box` on `grid`, widened by `pad` +/// cells per side and clamped to the grid. +template +[[nodiscard]] KOKKOS_INLINE_FUNCTION CellRange overlapping_cells( + const AABBox& box, const UniformGrid& grid, LO pad = 1) +{ + CellRange range; + for (int d = 0; d < Dim; ++d) { + range.lo[d] = Kokkos::max( + grid.AxisCellIndex(d, box.center[d] - box.half_width[d]) - pad, 0); + range.hi[d] = Kokkos::min( + grid.AxisCellIndex(d, box.center[d] + box.half_width[d]) + pad, + grid.divisions[d] - 1); + } + return range; +} + +/// Cells on `grid` that could satisfy a `Pred` test within `bound`: the cells +/// `bound` spans, widened by the predicate's own `boundary_pad` so that a cell +/// which merely touches the bound is included. +template +[[nodiscard]] KOKKOS_INLINE_FUNCTION CellRange candidate_cells( + const AABBox& bound, const UniformGrid& grid) +{ + return overlapping_cells(bound, grid, Pred::boundary_pad); +} + +} // namespace pcms + +#endif // PCMS_LOCALIZATION_OVERLAP_H diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index da74f055a..cd5a140d4 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/point_search.cpp @@ -1,7 +1,10 @@ #include "point_search.h" +#include "pcms/localization/overlap.h" #include #include #include +#include +#include // From // https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation @@ -112,217 +115,63 @@ LO GetOwningElementId(Omega_h::Mesh& mesh, int mesh_dim, int entity_dim, return owner; } -KOKKOS_INLINE_FUNCTION -AABBox<2> triangle_bbox(const Omega_h::Matrix<2, 3>& coords) -{ - Kokkos::Array max{coords(0, 0), coords(1, 0)}; - Kokkos::Array min{coords(0, 0), coords(1, 0)}; - for (int i = 1; i < 3; ++i) { - max[0] = std::fmax(max[0], coords(0, i)); - max[1] = std::fmax(max[1], coords(1, i)); - min[0] = std::fmin(min[0], coords(0, i)); - min[1] = std::fmin(min[1], coords(1, i)); - } - return {.center = {(max[0] + min[0]) / 2.0, (max[1] + min[1]) / 2.0}, - .half_width = {(max[0] - min[0]) / 2.0, (max[1] - min[1]) / 2.0}}; -} - -template -KOKKOS_INLINE_FUNCTION AABBox simplex_bbox( - const Omega_h::Matrix& coords) -{ - Kokkos::Array max; - Kokkos::Array min; - for (int j = 0; j < dim; ++j) { - max[j] = coords(j, 0); - min[j] = coords(j, 0); - } - for (int i = 1; i < dim + 1; ++i) { - for (int j = 0; j < dim; ++j) { - max[j] = std::fmax(max[j], coords(j, i)); - min[j] = std::fmin(min[j], coords(j, i)); - } - } - - Kokkos::Array center; - Kokkos::Array half_width; - - for (int j = 0; j < dim; ++j) { - center[j] = (max[j] + min[j]) / 2.0; - half_width[j] = (max[j] - min[j]) / 2.0; - } - - return {.center = center, .half_width = half_width}; -} - -// Liang, You-Dong, and B. A. Barsky. “A New Concept and Method for Line -// Clipping.” ACM Transactions on Graphics 3, no. 1 (January 1984): 1–22. -// https://doi.org/10.1145/357332.357333. -KOKKOS_INLINE_FUNCTION -bool clipt(Real p, Real q, Real& t0, Real& t1) -{ - if (p < 0) { - auto r = q / p; - if (r > t1) { - return false; - } - if (r > t0) { - t0 = r; - } - } else if (p > 0) { - auto r = q / p; - if (r < t0) { - return false; - } - if (r < t1) { - t1 = r; - } - } - // p==0 - else { - if (q < 0) { - return false; - } - } - return true; -} -// Liang, You-Dong, and B. A. Barsky. “A New Concept and Method for Line -// Clipping.” ACM Transactions on Graphics 3, no. 1 (January 1984): 1–22. -// https://doi.org/10.1145/357332.357333. -KOKKOS_INLINE_FUNCTION -bool line_intersects_bbox(const Omega_h::Vector<2>& p0, - const Omega_h::Vector<2>& p1, const AABBox<2>& bbox) -{ - Real t0 = 0; - Real t1 = 1; - auto xleft = bbox.center[0] - bbox.half_width[0]; - auto deltax = p1[0] - p0[0]; - if (clipt(-deltax, p0[0] - xleft, t0, t1)) { - auto xright = bbox.center[0] + bbox.half_width[0]; - if (clipt(deltax, xright - p0[0], t0, t1)) { - auto deltay = p1[1] - p0[1]; - auto ybottom = bbox.center[1] - bbox.half_width[1]; - if (clipt(-deltay, p0[1] - ybottom, t0, t1)) { - auto ytop = bbox.center[1] + bbox.half_width[1]; - if (clipt(deltay, ytop - p0[1], t0, t1)) { - // full liang-barksy algorithm computes new x-y coordinates here, but - // we just need to check intersections - return true; - } - } - } - } - return false; -} - -template -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool within_bbox( - const Omega_h::Vector coord, const AABBox& bbox) noexcept -{ - for (int i = 0; i < dim; ++i) { - if (coord[i] < bbox.center[i] - bbox.half_width[i]) - return false; - if (coord[i] > bbox.center[i] + bbox.half_width[i]) - return false; - } - return true; -} - -template -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_simplex( - const AABBox& bbox, const Omega_h::Matrix& coords) -{ - // each dimension has a pair of opposing "walls" - // 2D: { [left, right], [top, bottom] } -> { left, right, top, bottom } - // 3D: { [left, right], [top, bottom], [front, back] } -> { left, ..., back } - Kokkos::Array bbox_walls{}; - for (int i = 0; i < dim; i++) { - bbox_walls[i * 2] = bbox.center[i] - bbox.half_width[i]; - bbox_walls[i * 2 + 1] = bbox.center[i] + bbox.half_width[i]; - } - - // 1 << dim == 2 ** dim == num vertices in ndim bounding box / hypercube - constexpr unsigned num_verts = 1 << dim; - - // Each vertex is a just a unique combination of walls - // eg [left, bottom] (2D) or [right, top, front] (3) - for (unsigned i = 0; i < num_verts; ++i) { - // conveniently, i acts a bit field representing the current combination - Omega_h::Vector vert; - for (unsigned j = 0; j < dim; ++j) { - // eg 110 = 6 = [left, top, back] - // eg 01 = 1 = [right, top] - vert[j] = (i >> j) & 1 ? bbox_walls[j * 2] : bbox_walls[j * 2 + 1]; - } - auto xi = Omega_h::barycentric_from_global(vert, coords); - if (Omega_h::is_barycentric_inside(xi)) { - return true; - } - } - return false; -} - /** * Check if a triangle element represented by 3 coordinates in two dimensions - * intersects with a bounding box + * intersects with a bounding box. */ [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) { - // triangle and grid cell bounding box intersect - if (intersects(triangle_bbox(coords), bbox)) { - // if any of the triangle verts inside of bbox - if (within_bbox<2>(coords[0], bbox) || within_bbox<2>(coords[1], bbox) || - within_bbox<2>(coords[2], bbox)) { - return true; - } - // if any of the bbox verts are within the triangle - if (bbox_verts_within_simplex(bbox, coords)) { - return true; - } - // if any of the triangle's edges intersect with the bounding box - if (line_intersects_bbox(coords[0], coords[1], bbox) || - line_intersects_bbox(coords[1], coords[2], bbox) || - line_intersects_bbox(coords[2], coords[0], bbox)) { - return true; - } - } - return false; -} - -template -[[nodiscard]] KOKKOS_FUNCTION bool simplex_intersects_bbox( - const Omega_h::Matrix& coords, const AABBox& bbox) -{ - return intersects(simplex_bbox(coords), bbox); - // TODO: Add refined cases from triangle_intersects_bbox + return Overlap<2, TriangleVsBbox>::apply(coords, bbox); } namespace detail { +// The candidate-map machinery. Declared in the header, defined here, and +// explicitly instantiated for the supported predicates at the end of this +// file -- see the note there. + /** - * Functor for constructing the mapping from grid cells to intersecting - * triangles. Each row of the resulting CSR structure represents a grid cell and - * each row entry corresponds to an ID of an element that intersects that grid - * cell. - * \Warning since this works on Omega_h meshes, we currently assume each element - * is a 2D simplex (triangle) - * \Warning since this uses Omega_h data which is only available in the - * "Default" Execution space, the should not be used in an alternative EXE space + * The bridge between "an element id" and "the geometry a predicate tests". + * + * This is the only shape-specific part of building a candidate map. Everything + * downstream -- the broad-phase binning and the predicate test -- works on + * whatever `operands()` returns, so supporting a different element shape is a + * new policy with these three members and nothing else changes: + * + * struct MyElements { + * static MyElements FromMesh(Omega_h::Mesh&); + * Operands operands(LO elem) const; // what Pred::apply + * takes AABBox bound(const Operands&) const; // a box containing it + * }; + * + * `bound` is allowed to be conservative: the predicate rejects what it + * over-includes, whereas a bound that is too tight silently drops candidates. */ -/// The bbox-overlap predicate for a simplex against one grid cell. 2D applies -/// the refined edge tests, 3D is a plain bbox overlap; both match what the -/// previous cell-major functors used, so candidate lists are unchanged. template -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool simplex_overlaps_cell( - const Omega_h::Matrix& coords, const AABBox& bbox) +struct MeshSimplexElements { - if constexpr (Dim == 2) { - return triangle_intersects_bbox(coords, bbox); - } else { - return simplex_intersects_bbox(coords, bbox); + Omega_h::LOs elems2verts; + Omega_h::Reals coords; + + [[nodiscard]] static MeshSimplexElements FromMesh(Omega_h::Mesh& mesh) + { + return {mesh.ask_elem_verts(), mesh.coords()}; } -} + + [[nodiscard]] KOKKOS_INLINE_FUNCTION Omega_h::Matrix operands( + LO elem) const + { + const auto verts = Omega_h::gather_verts(elems2verts, elem); + return Omega_h::gather_vectors(coords, verts); + } + + [[nodiscard]] KOKKOS_INLINE_FUNCTION AABBox bound( + const Omega_h::Matrix& ops) const + { + return simplex_bbox(ops); + } +}; /** * Functor for the element -> grid cell map, the transpose of what the search @@ -332,19 +181,18 @@ template * testing every element against a cell, each element only visits the cells its * own bounding box spans. That makes construction O(nelems * cells_per_elem) * instead of O(num_cells * nelems). + * + * Operand-agnostic: `Elements` decides what an element *is* and how to bound + * it, `Pred` decides when it overlaps a cell. Neither is fixed here, so the + * same functor bins any shape a caller can supply a policy for. */ -template +template > struct ElemCellOverlapFunctor { - ElemCellOverlapFunctor(Omega_h::Mesh& mesh, + ElemCellOverlapFunctor(Elements elements, Kokkos::View[1]> grid) - : elems2verts_(mesh.ask_elem_verts()), coords_(mesh.coords()), grid_(grid) + : elements_(elements), grid_(grid) { - if (mesh.dim() != Dim) { - std::cerr << "ElemCellOverlapFunctor: mesh dimension does not match the " - "functor dimension\n"; - std::terminate(); - } } /// Two-pass functor. On the first pass we count the cells this element @@ -352,36 +200,16 @@ struct ElemCellOverlapFunctor KOKKOS_INLINE_FUNCTION LO operator()(LO elem, LO* fill) const { - const auto elem_verts = Omega_h::gather_verts(elems2verts_, elem); - const auto vertex_coords = - Omega_h::gather_vectors(coords_, elem_verts); const auto& grid = grid_(0); + const auto ops = elements_.operands(elem); - // Per-axis inclusive range of cells the element's bbox spans, widened by - // one cell on each side. The overlap predicate treats a shared boundary - // plane as an intersection, so a cell that only touches the bbox counts as - // a candidate and would be missed by the unwidened range -- which happens - // constantly on a structured mesh whose vertices land on grid lines. One - // cell is enough: a cell two or more away cannot touch the bbox at all. - // The predicate below rejects whatever the widening over-includes. - Kokkos::Array lo, hi, idx; - for (int d = 0; d < Dim; ++d) { - Real min = vertex_coords[0][d]; - Real max = vertex_coords[0][d]; - for (int v = 1; v < Dim + 1; ++v) { - min = Kokkos::min(min, vertex_coords[v][d]); - max = Kokkos::max(max, vertex_coords[v][d]); - } - lo[d] = Kokkos::max(grid.AxisCellIndex(d, min) - 1, 0); - hi[d] = - Kokkos::min(grid.AxisCellIndex(d, max) + 1, grid.divisions[d] - 1); - idx[d] = lo[d]; - } + const auto range = candidate_cells(elements_.bound(ops), grid); + Kokkos::Array idx = range.lo; LO num_overlaps = 0; while (true) { const LO cell = grid.CellIndexFromAxisIndices(idx); - if (simplex_overlaps_cell(vertex_coords, grid.GetCellBBOX(cell))) { + if (Overlap::apply(ops, grid.GetCellBBOX(cell))) { if (fill) { fill[num_overlaps] = cell; } @@ -390,10 +218,10 @@ struct ElemCellOverlapFunctor // odometer over the inclusive range [lo, hi] int d = 0; for (; d < Dim; ++d) { - if (++idx[d] <= hi[d]) { + if (++idx[d] <= range.hi[d]) { break; } - idx[d] = lo[d]; + idx[d] = range.lo[d]; } if (d == Dim) { break; @@ -403,8 +231,7 @@ struct ElemCellOverlapFunctor } private: - Omega_h::LOs elems2verts_; - Omega_h::Reals coords_; + Elements elements_; Kokkos::View[1]> grid_; }; @@ -416,18 +243,26 @@ struct ElemCellOverlapFunctor * in.numRows() (the element count), which both over-allocates the row map and * writes out of bounds once num_grid_cells exceeds the element count. */ -template +template Kokkos::Crs construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View[1]> grid, int num_grid_cells) { + if (mesh.dim() != Dim) { + std::cerr + << "construct_intersection_map: mesh dimension does not match the " + "map dimension\n"; + std::terminate(); + } + using CrsT = Kokkos::Crs; const auto nelems = mesh.nelems(); CrsT elem2cells{}; Kokkos::count_and_fill_crs(elem2cells, nelems, - ElemCellOverlapFunctor{mesh, grid}); + ElemCellOverlapFunctor{ + Elements::FromMesh(mesh), grid}); Kokkos::View counts("candidate_counts", num_grid_cells); Kokkos::parallel_for( @@ -483,15 +318,12 @@ construct_intersection_map(Omega_h::Mesh& mesh, * ground truth for what the candidate map should contain, so it is kept and * exposed as an oracle for the equivalence test. */ -template +template > struct ReferenceCellMajorFunctor { - ReferenceCellMajorFunctor(Omega_h::Mesh& mesh, - Kokkos::View[1]> grid) - : elems2verts_(mesh.ask_elem_verts()), - coords_(mesh.coords()), - grid_(grid), - nelems_(mesh.nelems()) + ReferenceCellMajorFunctor(Elements elements, + Kokkos::View[1]> grid, LO nelems) + : elements_(elements), grid_(grid), nelems_(nelems) { } @@ -501,11 +333,8 @@ struct ReferenceCellMajorFunctor const auto cell_bbox = grid_(0).GetCellBBOX(row); LO num_intersections = 0; for (LO elem = 0; elem < nelems_; ++elem) { - const auto elem_verts = - Omega_h::gather_verts(elems2verts_, elem); - const auto vertex_coords = - Omega_h::gather_vectors(coords_, elem_verts); - if (simplex_overlaps_cell(vertex_coords, cell_bbox)) { + const auto ops = elements_.operands(elem); + if (Overlap::apply(ops, cell_bbox)) { if (fill) { fill[num_intersections] = elem; } @@ -516,13 +345,12 @@ struct ReferenceCellMajorFunctor } private: - Omega_h::LOs elems2verts_; - Omega_h::Reals coords_; + Elements elements_; Kokkos::View[1]> grid_; LO nelems_; }; -template +template Kokkos::Crs construct_intersection_map_reference(Omega_h::Mesh& mesh, Kokkos::View[1]> grid, @@ -530,16 +358,17 @@ construct_intersection_map_reference(Omega_h::Mesh& mesh, { Kokkos::Crs intersection_map{}; Kokkos::count_and_fill_crs(intersection_map, num_grid_cells, - ReferenceCellMajorFunctor{mesh, grid}); + ReferenceCellMajorFunctor{ + Elements::FromMesh(mesh), grid, mesh.nelems()}); return intersection_map; } - Kokkos::Crs construct_intersection_map_reference_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { - return construct_intersection_map_reference<2>(mesh, grid, num_grid_cells); + return construct_intersection_map_reference<2, TriangleVsBbox>( + mesh, grid, num_grid_cells); } Kokkos::Crs @@ -547,17 +376,24 @@ construct_intersection_map_reference_3d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { - return construct_intersection_map_reference<3>(mesh, grid, num_grid_cells); + return construct_intersection_map_reference<3, SimplexVsBbox>(mesh, grid, + num_grid_cells); } // num_grid_cells should be result of grid.GetNumCells(), take as argument to // avoid extra copy of grid from gpu to cpu +// +// The 2D and 3D overloads name their predicate here rather than reading it from +// a shared mapping: a "2D refined candidate map" is what this function *is*, so +// the choice belongs to its definition and cannot leak into another call site. +// Callers needing a different predicate use the template directly. Kokkos::Crs construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { - return construct_intersection_map<2>(mesh, grid, num_grid_cells); + return construct_intersection_map<2, TriangleVsBbox>(mesh, grid, + num_grid_cells); } Kokkos::Crs @@ -565,8 +401,38 @@ construct_intersection_map_3d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { - return construct_intersection_map<3>(mesh, grid, num_grid_cells); + return construct_intersection_map<3, SimplexVsBbox>(mesh, grid, + num_grid_cells); } +// The supported predicates, instantiated here so that the search classes' +// predicate-taking constructors and the tests can use them from other +// translation units. This is the complete set: naming a predicate that is not +// listed below is a link error, which is the accepted cost of keeping the +// machinery out of the header. Add a line here when adding a predicate that +// callers outside this file need. +template Kokkos::Crs +construct_intersection_map<2, TriangleVsBbox>(Omega_h::Mesh&, + Kokkos::View[1]>, + int); +template Kokkos::Crs +construct_intersection_map<2, SimplexVsBbox>(Omega_h::Mesh&, + Kokkos::View[1]>, + int); +template Kokkos::Crs +construct_intersection_map<3, SimplexVsBbox>(Omega_h::Mesh&, + Kokkos::View[1]>, + int); + +template Kokkos::Crs +construct_intersection_map_reference<2, TriangleVsBbox>( + Omega_h::Mesh&, Kokkos::View[1]>, int); +template Kokkos::Crs +construct_intersection_map_reference<2, SimplexVsBbox>( + Omega_h::Mesh&, Kokkos::View[1]>, int); +template Kokkos::Crs +construct_intersection_map_reference<3, SimplexVsBbox>( + Omega_h::Mesh&, Kokkos::View[1]>, int); + } // namespace detail template @@ -791,6 +657,13 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances), mesh_(mesh) +{ + InitializeGeometry(mesh, Nx, Ny); + // This search's default predicate: the refined triangle-vs-box test. + InitializeMap(); +} + +void GridPointSearch2D::InitializeGeometry(Omega_h::Mesh& mesh, LO Nx, LO Ny) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); @@ -800,11 +673,7 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, .bot_left = {mesh_bbox.min[0], mesh_bbox.min[1]}, .divisions = {Nx, Ny}}; Kokkos::deep_copy(grid_, grid_h); - // Determine inflation radius from tolerances (max of vertex/edge tol) - auto tol_h = Kokkos::create_mirror_view(tolerances_); - Kokkos::deep_copy(tol_h, tolerances_); - candidate_map_ = - detail::construct_intersection_map_2d(mesh, grid_, grid_h(0).GetNumCells()); + num_grid_cells_ = grid_h(0).GetNumCells(); coords_ = mesh.coords(); tris2verts_ = mesh.ask_elem_verts(); tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); @@ -932,6 +801,14 @@ GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances), mesh_(mesh) +{ + InitializeGeometry(mesh, Nx, Ny, Nz); + // This search's default predicate: the conservative bounding-box test. + InitializeMap(); +} + +void GridPointSearch3D::InitializeGeometry(Omega_h::Mesh& mesh, LO Nx, LO Ny, + LO Nz) { auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); @@ -949,8 +826,7 @@ GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, .divisions = {Nx, Ny, Nz}}; Kokkos::deep_copy(grid_, grid_h); - candidate_map_ = - detail::construct_intersection_map_3d(mesh, grid_, grid_h(0).GetNumCells()); + num_grid_cells_ = grid_h(0).GetNumCells(); coords_ = mesh.coords(); tris2verts_ = mesh.ask_elem_verts(); tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); diff --git a/src/pcms/localization/point_search.h b/src/pcms/localization/point_search.h index 602c99fb3..220d14020 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -5,10 +5,12 @@ #include #include +#include #include #include #include +#include "pcms/localization/overlap.h" #include "pcms/utility/types.h" #include "pcms/utility/uniform_grid.h" #include "pcms/utility/bounding_box.h" @@ -81,6 +83,33 @@ Kokkos::Crs construct_intersection_map_reference_3d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); + +/// How the candidate map learns an element's geometry, and a conservative bound +/// on it. The default policy for an Omega_h simplex mesh; defined in the source +/// file alongside the rest of the map machinery. A different element shape +/// means a different policy with the same three members, which is the only +/// thing that has to change to bin that shape against the same grid. +template +struct MeshSimplexElements; + +/** + * Build the grid cell -> element candidate map with `Pred`: the element-major + * build, and the point search's hot path at construction time. + */ +template > +Kokkos::Crs +construct_intersection_map(Omega_h::Mesh& mesh, + Kokkos::View[1]> grid, + int num_grid_cells); + +/// O(num_cells * nelems) brute-force counterpart, for tests. \see +/// construct_intersection_map +template > +Kokkos::Crs +construct_intersection_map_reference(Omega_h::Mesh& mesh, + Kokkos::View[1]> grid, + int num_grid_cells); + } // namespace detail [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( @@ -175,6 +204,20 @@ class GridPointSearch2D : public PointLocalizationSearch2D GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances); + /** + * Build the candidate map with a predicate of the caller's choosing: + * + * GridPointSearch2D search{mesh, nx, ny, tolerances, MyPredicate{}}; + * + * The predicate decides which elements are offered as candidates for a grid + * cell. The query itself is predicate-independent -- it tests containment of + * the point in each candidate -- so a narrower predicate costs a pass over + * more candidates but cannot change which answer is reported. + */ + template + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, + const PointSearchTolerances& tolerances, Pred); + /** * Given points in global coordinates, returns for each the entity (vertex, * edge or face, chosen with the vertex and edge tolerances) containing it @@ -187,6 +230,13 @@ class GridPointSearch2D : public PointLocalizationSearch2D Kokkos::View results) const override; private: + /// Grid and adjacency: everything the map needs except the predicate. + void InitializeGeometry(Omega_h::Mesh& mesh, LO Nx, LO Ny); + + /// The candidate map, which is the only predicate-dependent part. + template + void InitializeMap(); + Omega_h::Mesh mesh_; Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; @@ -194,6 +244,7 @@ class GridPointSearch2D : public PointLocalizationSearch2D Omega_h::Adj edges2faces_up_; Omega_h::Adj verts2faces_up_; Kokkos::View grid_{"uniform grid"}; + LO num_grid_cells_ = 0; CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; Omega_h::Reals coords_; @@ -211,6 +262,11 @@ class GridPointSearch3D : public PointLocalizationSearch3D GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, const PointSearchTolerances& tolerances); + /// \see GridPointSearch2D's predicate-taking constructor. + template + GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, + const PointSearchTolerances& tolerances, Pred); + /** * Given points in global coordinates, returns for each the id of a * tetrahedron containing it and the point's barycentric coordinates in that @@ -234,6 +290,13 @@ class GridPointSearch3D : public PointLocalizationSearch3D Kokkos::View results) const override; private: + /// \see GridPointSearch2D::InitializeGeometry + void InitializeGeometry(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz); + + /// \see GridPointSearch2D::InitializeMap + template + void InitializeMap(); + Omega_h::Mesh mesh_; Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; @@ -242,11 +305,46 @@ class GridPointSearch3D : public PointLocalizationSearch3D Omega_h::Adj edges2regions_up_; Omega_h::Adj faces2regions_up_; Kokkos::View[1]> grid_{"uniform grid"}; + LO num_grid_cells_ = 0; CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; Omega_h::Reals coords_; }; +template +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, + const PointSearchTolerances& tolerances, + Pred) + : PointLocalizationSearch2D(tolerances), mesh_(mesh) +{ + InitializeGeometry(mesh, Nx, Ny); + InitializeMap(); +} + +template +void GridPointSearch2D::InitializeMap() +{ + candidate_map_ = + detail::construct_intersection_map<2, Pred>(mesh_, grid_, num_grid_cells_); +} + +template +GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, + const PointSearchTolerances& tolerances, + Pred) + : PointLocalizationSearch3D(tolerances), mesh_(mesh) +{ + InitializeGeometry(mesh, Nx, Ny, Nz); + InitializeMap(); +} + +template +void GridPointSearch3D::InitializeMap() +{ + candidate_map_ = + detail::construct_intersection_map<3, Pred>(mesh_, grid_, num_grid_cells_); +} + /// A grid search of either spatial dimension, as owned by components that are /// dimension-agnostic at compile time (e.g. an evaluator factory over an /// Omega_h mesh) and handed to consumers that need to locate points on the same diff --git a/src/pcms/transfer/mesh_intersection.hpp b/src/pcms/transfer/mesh_intersection.hpp index 9d5ead08d..18bd5e882 100644 --- a/src/pcms/transfer/mesh_intersection.hpp +++ b/src/pcms/transfer/mesh_intersection.hpp @@ -2,11 +2,11 @@ #define PCMS_TRANSFER_MESH_INTERSECTION_HPP #include +#include #include #include #include #include -#include #include #include #include @@ -51,65 +51,7 @@ template const r3d::Few, Dim + 1>& a, const r3d::Few, Dim + 1>& b) { - // Each iteration treats one simplex as the half-space owner. Both are - // needed: a face of `a` may separate the pair while no face of `b` does. - for (int side = 0; side < 2; ++side) { - const auto& owner = (side == 0) ? a : b; - const auto& other = (side == 0) ? b : a; - - // Face `f` is the face opposite vertex `f`. - for (int f = 0; f < Dim + 1; ++f) { - // Vertices spanning the face, in index order with `f` skipped. - int face_verts[Dim]; - int n_face_verts = 0; - for (int v = 0; v < Dim + 1; ++v) { - if (v != f) { - face_verts[n_face_verts++] = v; - } - } - const auto& origin = owner[face_verts[0]]; - - // Face normal, oriented to point away from the opposite vertex so the - // owner's interior lies on the negative side. - r3d::Vector normal; - if constexpr (Dim == 2) { - const auto edge = owner[face_verts[1]] - origin; - normal[0] = edge[1]; - normal[1] = -edge[0]; - } else { - const auto e1 = owner[face_verts[1]] - origin; - const auto e2 = owner[face_verts[2]] - origin; - normal[0] = e1[1] * e2[2] - e1[2] * e2[1]; - normal[1] = e1[2] * e2[0] - e1[0] * e2[2]; - normal[2] = e1[0] * e2[1] - e1[1] * e2[0]; - } - - Omega_h::Real opposite_side = 0.0; - Omega_h::Real normal_sq = 0.0; - for (int d = 0; d < Dim; ++d) { - opposite_side += normal[d] * (owner[f][d] - origin[d]); - normal_sq += normal[d] * normal[d]; - } - // A degenerate face gives no usable axis. - if (normal_sq == 0.0) { - continue; - } - const Omega_h::Real orientation = (opposite_side > 0.0) ? -1.0 : 1.0; - - bool all_outside = true; - for (int v = 0; v < Dim + 1 && all_outside; ++v) { - Omega_h::Real signed_distance = 0.0; - for (int d = 0; d < Dim; ++d) { - signed_distance += normal[d] * (other[v][d] - origin[d]); - } - all_outside = (orientation * signed_distance) >= 0.0; - } - if (all_outside) { - return true; - } - } - } - return false; + return Overlap::apply(a, b); } /** diff --git a/src/pcms/utility/uniform_grid.h b/src/pcms/utility/uniform_grid.h index e6730fcec..7f5731310 100644 --- a/src/pcms/utility/uniform_grid.h +++ b/src/pcms/utility/uniform_grid.h @@ -65,6 +65,20 @@ struct UniformGrid return GetCellIndex(indexes); } + /// Per-axis cell indices in coordinate order (x,y,z) for a linear cell ID. + /// + /// The exact inverse of `CellIndexFromAxisIndices`. Provided so callers that + /// work with linear cell IDs -- a routing table keyed by cell, for example -- + /// do not have to reimplement the grid's internal index ordering, which is + /// reversed relative to coordinates and easy to get wrong silently. + [[nodiscard]] KOKKOS_INLINE_FUNCTION Kokkos::Array + AxisIndicesFromCellIndex(LO idx) const + { + auto indices = GetDimensionedIndex(idx); + reverse(indices); + return indices; + } + [[nodiscard]] KOKKOS_INLINE_FUNCTION AABBox GetCellBBOX(LO idx) const { auto index = GetDimensionedIndex(idx); @@ -98,21 +112,20 @@ struct UniformGrid return true; } + /// Per-axis indices of a linear cell ID, in reversed order (z, y, x), the + /// order `GetCellIndex` takes. The exact inverse of `GetCellIndex`, which + /// numbers cells as `idx = x + Nx * (y + Ny * z)`. + /// + /// This is a numbering of cells, not a memory layout: it is the same on host + /// and device, whatever layout a Kokkos view uses. [[nodiscard]] KOKKOS_INLINE_FUNCTION Kokkos::Array GetDimensionedIndex(LO idx) const { - LO stride = 1; - for (std::size_t i = 0; i < dim - 1; ++i) { - stride *= divisions[i]; - } Kokkos::Array result; - - for (size_t i = 0; i < dim; ++i) { - result[i] = idx / stride; - idx -= result[i] * stride; - stride /= divisions[i]; + for (size_t axis = 0; axis < dim; ++axis) { + result[dim - 1 - axis] = idx % divisions[axis]; + idx /= divisions[axis]; } - return result; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 12ec4f11b..2d270abaf 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -429,6 +429,7 @@ if(Catch2_FOUND) test_error_handling.cpp test_eqdsk.cpp test_uniform_grid.cpp + test_cell_range.cpp test_field_evaluation.cpp test_field_interpolation.cpp test_field_copy.cpp diff --git a/test/test_cell_range.cpp b/test/test_cell_range.cpp new file mode 100644 index 000000000..24e75c1ed --- /dev/null +++ b/test/test_cell_range.cpp @@ -0,0 +1,270 @@ +#include +#include +#include +#include + +// Serial tests of the cell-range layer the routing table is built on: which +// grid cells a bounding box covers (`candidate_cells`, `overlapping_cells`) and +// the invariants of `CellRange`. + +using pcms::AABBox; +using pcms::candidate_cells; +using pcms::CellRange; +using pcms::LO; +using pcms::Overlap; +using pcms::overlapping_cells; +using pcms::Real; +using pcms::SimplexVsBbox; +using pcms::TriangleVsBbox; +using pcms::Uniform2DGrid; +using pcms::UniformGrid; + +TEST_CASE("bbox vs grid cell range") +{ + Uniform2DGrid grid{ + .edge_length = {10, 10}, .bot_left = {0, 0}, .divisions = {10, 10}}; + + SECTION("box inside one cell still covers padded neighbors") + { + AABBox<2> box{.center = {0.5, 0.5}, .half_width = {0.1, 0.1}}; + const auto range = overlapping_cells<2>(box, grid); + // default pad = 1 clamps at the grid border + REQUIRE(range.lo[0] == 0); + REQUIRE(range.lo[1] == 0); + REQUIRE(range.hi[0] == 1); + REQUIRE(range.hi[1] == 1); + } + + SECTION("no padding gives the exact span") + { + AABBox<2> box{.center = {2.5, 2.5}, .half_width = {0.4, 0.4}}; + const auto range = overlapping_cells<2>(box, grid, 0); + REQUIRE(range.lo[0] == 2); + REQUIRE(range.lo[1] == 2); + REQUIRE(range.hi[0] == 2); + REQUIRE(range.hi[1] == 2); + } + + SECTION("box spanning several cells") + { + AABBox<2> box{.center = {2.5, 2.5}, .half_width = {1.0, 1.0}}; + const auto range = overlapping_cells<2>(box, grid, 0); + REQUIRE(range.lo[0] == 1); + REQUIRE(range.lo[1] == 1); + REQUIRE(range.hi[0] == 3); + REQUIRE(range.hi[1] == 3); + } +} + +namespace +{ +/// True when grid cell `(cx, cy)` and `box` meet, each treated as closed. +/// +/// A shared boundary counts, which is the semantics the overlap predicate has +/// and therefore the set `candidate_cells` has to cover. +[[nodiscard]] bool cell_touches_box(const Uniform2DGrid& grid, LO cx, LO cy, + const AABBox<2>& box) +{ + const Kokkos::Array cell{cx, cy}; + for (int d = 0; d < 2; ++d) { + const Real width = grid.edge_length[d] / grid.divisions[d]; + const Real cell_lo = grid.bot_left[d] + cell[d] * width; + const Real cell_hi = cell_lo + width; + const Real box_lo = box.center[d] - box.half_width[d]; + const Real box_hi = box.center[d] + box.half_width[d]; + if (cell_hi < box_lo || cell_lo > box_hi) { + return false; + } + } + return true; +} +} // namespace + +/// The property `candidate_cells` exists to guarantee: its range contains every +/// cell a closed-boundary predicate can accept. A range that is too narrow +/// drops candidates silently -- the element is never offered to the query and +/// nothing downstream can tell -- so this is checked as a coverage property +/// rather than against the exact bounds, which are allowed to over-include. +TEST_CASE( + "candidate_cells covers every cell a closed-boundary predicate accepts") +{ + const Uniform2DGrid grid{ + .edge_length = {10, 10}, .bot_left = {0, 0}, .divisions = {10, 10}}; + + // Centres at, near and between grid lines, and extents that land exactly on a + // grid line: that is when an unwidened range loses the touching cell. + for (const Real center : {0.5, 2.5, 5.0, 9.5}) { + for (const Real half : {0.0, 0.25, 0.5, 1.0, 3.0}) { + const AABBox<2> box{.center = {center, center}, + .half_width = {half, half}}; + const auto range = candidate_cells<2, TriangleVsBbox>(box, grid); + for (LO cx = 0; cx < grid.divisions[0]; ++cx) { + for (LO cy = 0; cy < grid.divisions[1]; ++cy) { + if (cell_touches_box(grid, cx, cy, box)) { + REQUIRE(range.Contains(Kokkos::Array{cx, cy})); + } + } + } + } + } +} + +TEST_CASE("candidate_cells bins a supplied bound") +{ + const Uniform2DGrid grid{ + .edge_length = {10, 10}, .bot_left = {0, 0}, .divisions = {10, 10}}; + + // A degenerate triangle whose bounding box is the single grid line x = y = 3: + // the cells on either side of that line only touch it, so the padding must + // come from the predicate rather than from the caller. + const Omega_h::Matrix<2, 3> degenerate{{3.0, 3.0}, {3.0, 3.0}, {3.0, 3.0}}; + // `candidate_cells` bins a bound it is given; the element side supplies it. + const auto range = candidate_cells<2, TriangleVsBbox>( + pcms::detail::simplex_bbox<2>(degenerate), grid); + REQUIRE(range.Contains(Kokkos::Array{3, 3})); + REQUIRE(range.Contains(Kokkos::Array{2, 2})); +} + +TEST_CASE("CellRange states its own invariants") +{ + const auto empty = pcms::CellRange<2>::Empty(); + REQUIRE(empty.IsEmpty()); + REQUIRE(!empty.Contains(Kokkos::Array{0, 0})); + + const pcms::CellRange<2> unit{.lo = {2, 2}, .hi = {2, 2}}; + REQUIRE(!unit.IsEmpty()); + REQUIRE(unit.Contains(Kokkos::Array{2, 2})); + REQUIRE(!unit.Contains(Kokkos::Array{3, 2})); + + // The empty range must not be what value initialization produces: a + // zero-initialized CellRange denotes the single cell at axis index 0. + const pcms::CellRange<2> zeroed{}; + REQUIRE(!zeroed.IsEmpty()); + REQUIRE(zeroed.Contains(Kokkos::Array{0, 0})); +} + +// What dimension does the overlap + grid-cell-range layer actually support? +// +// `Dim` is a template parameter throughout and every loop runs to `Dim`, so +// nothing in these functions is written for 2 or 3 specifically. That is easy +// to lose by accident -- one hardcoded `[2]` or `<3, 3>` and the layer silently +// becomes 2D/3D only -- so it is pinned here at dimensions no mesh can +// represent, which is exactly where a hidden assumption would show up. +// +// Findings, and the boundary they draw: +// +// generic, and covered here AABBox, intersects, simplex_bbox, CellRange, +// overlapping_cells, candidate_cells, +// UniformGrid, Overlap +// generic over Dim, but the MeshSimplexElements, ElemCellOverlapFunctor, +// map builder needs a mesh construct_intersection_map -- all take `Dim`, +// but are fed an Omega_h mesh, which is <= 3D, so +// dimensions above 3 cannot be reached through +// them +// NOT generic GridPointSearch2D / GridPointSearch3D, which +// are +// concrete classes with hand-written +// per-dimension query logic (VERTEX/EDGE/FACE for +// 2D, REGION for 3D) +// NOT generic, was silent Overlap: nominally +// Dim-generic but its face normal is a 2D +// perpendicular or a 3D cross product. Now +// guarded by a static_assert, so asking for +// another dimension is a compile error instead of +// reading uninitialized memory. +TEST_CASE("overlap and cell-range layers are dimension-generic") +{ + SECTION("1D") + { + const UniformGrid<1> grid{ + .edge_length = {4.0}, .bot_left = {0.0}, .divisions = {4}}; + const AABBox<1> box{.center = {1.5}, .half_width = {0.25}}; + // Cell 1, padded by 1, clamped at the low edge. + const auto range = candidate_cells<1, SimplexVsBbox>(box, grid); + REQUIRE(range.lo[0] == 0); + REQUIRE(range.hi[0] == 2); + REQUIRE(range.Contains(Kokkos::Array{1})); + } + + SECTION("4D") + { + const UniformGrid<4> grid{.edge_length = {4.0, 4.0, 4.0, 4.0}, + .bot_left = {0.0, 0.0, 0.0, 0.0}, + .divisions = {4, 4, 4, 4}}; + const AABBox<4> box{.center = {1.5, 1.5, 1.5, 1.5}, + .half_width = {0.25, 0.25, 0.25, 0.25}}; + const auto range = candidate_cells<4, SimplexVsBbox>(box, grid); + for (int d = 0; d < 4; ++d) { + REQUIRE(range.lo[d] == 0); + REQUIRE(range.hi[d] == 2); + } + REQUIRE(range.Contains(Kokkos::Array{1, 1, 1, 1})); + // Same axis order as 2D: CellIndexFromAxisIndices round-trips. + const auto indices = Kokkos::Array{1, 1, 1, 1}; + const auto linear = grid.CellIndexFromAxisIndices(indices); + const auto back = grid.AxisIndicesFromCellIndex(linear); + for (int d = 0; d < 4; ++d) { + REQUIRE(back[d] == indices[d]); + } + } + + SECTION("5D range invariants") + { + const auto empty = CellRange<5>::Empty(); + REQUIRE(empty.IsEmpty()); + REQUIRE(!empty.Contains(Kokkos::Array{0, 0, 0, 0, 0})); + + const CellRange<5> unit{.lo = {3, 3, 3, 3, 3}, .hi = {3, 3, 3, 3, 3}}; + REQUIRE(!unit.IsEmpty()); + REQUIRE(unit.Contains(Kokkos::Array{3, 3, 3, 3, 3})); + REQUIRE(!unit.Contains(Kokkos::Array{3, 3, 3, 3, 4})); + } + + SECTION("SimplexVsBbox predicate in 4D") + { + const AABBox<4> box{.center = {1.0, 1.0, 1.0, 1.0}, + .half_width = {0.5, 0.5, 0.5, 0.5}}; + + // A degenerate simplex: all 5 vertices at a point inside the box. + Omega_h::Matrix<4, 5> inside{}; + for (int v = 0; v < 5; ++v) { + for (int d = 0; d < 4; ++d) { + inside(d, v) = 1.0; + } + } + REQUIRE(Overlap<4, SimplexVsBbox>::apply(inside, box)); + + // Far outside on every axis. + Omega_h::Matrix<4, 5> outside{}; + for (int v = 0; v < 5; ++v) { + for (int d = 0; d < 4; ++d) { + outside(d, v) = 100.0; + } + } + REQUIRE_FALSE(Overlap<4, SimplexVsBbox>::apply(outside, box)); + + // Straddling the box on one axis only: the conservative predicate is true. + Omega_h::Matrix<4, 5> straddling{}; + for (int v = 0; v < 5; ++v) { + for (int d = 0; d < 4; ++d) { + straddling(d, v) = 1.0; + } + } + straddling(3, 0) = 0.5; + straddling(3, 1) = 1.9; + REQUIRE(Overlap<4, SimplexVsBbox>::apply(straddling, box)); + } + + SECTION("SimplexVsBbox predicate in 5D") + { + const AABBox<5> box{.center = {2.0, 2.0, 2.0, 2.0, 2.0}, + .half_width = {1.0, 1.0, 1.0, 1.0, 1.0}}; + Omega_h::Matrix<5, 6> inside{}; + for (int v = 0; v < 6; ++v) { + for (int d = 0; d < 5; ++d) { + inside(d, v) = 2.0; + } + } + REQUIRE(Overlap<5, SimplexVsBbox>::apply(inside, box)); + } +} diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index f19fa3141..1c6d09317 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -527,3 +527,288 @@ TEST_CASE("point search encodes outside points unambiguously", REQUIRE(owners_h(2) == 0); } } + +// Baseline for the reused local search: on a mesh with no ghost layer every +// element is owned and every element centroid must localize to its own element. +// The distributed search treats a local answer as authoritative only when this +// holds, so it is the contract the distributed layer builds on. +TEST_CASE("ghost-free search localizes every element centroid") +{ + using pcms::GridPointSearch2D; + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 4, 4, 0, false); + REQUIRE(mesh.nghost_layers() == 0); + + const auto n_elems = mesh.nelems(); + const auto owned = Omega_h::HostRead(mesh.owned(2)); + REQUIRE(static_cast(owned.size()) == n_elems); + for (Omega_h::LO i = 0; i < n_elems; ++i) { + REQUIRE(owned[i] != 0); + } + + const auto n = pcms::DivisionsPerAxisForMesh<2>(mesh.nelems()); + GridPointSearch2D search{mesh, n, n}; + const auto centroids = + Omega_h::HostRead(pcms::get_entity_centroids(mesh, 2)); + Kokkos::View points("query points", n_elems); + auto points_h = Kokkos::create_mirror_view(points); + for (Omega_h::LO i = 0; i < n_elems; ++i) { + points_h(i, 0) = centroids[i * 2 + 0]; + points_h(i, 1) = centroids[i * 2 + 1]; + } + Kokkos::deep_copy(points, points_h); + + const auto owners_h = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace{}, search.GetOwningElementIds(search(points))); + for (Omega_h::LO i = 0; i < n_elems; ++i) { + REQUIRE(owners_h(i) == i); + } +} + +// A predicate chosen at the map-building call must govern *both* halves: the +// range of cells visited and the test applied to each. Checking the map against +// a brute-force build made with that same predicate is what catches the case +// where one half ignores the argument -- the fast map would then disagree. +// +// Only predicates instantiated in point_search.cpp can be named from here; the +// set is small and internal, which is why the map machinery no longer lives in +// the header. See the instantiations at the end of that file. +TEST_CASE("candidate map accepts a predicate other than the default") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false); + + Kokkos::View grid_d("uniform grid"); + auto grid_h = Kokkos::create_mirror_view(grid_d); + grid_h(0) = Uniform2DGrid{ + .edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; + Kokkos::deep_copy(grid_d, grid_h); + const auto ncells = grid_h(0).GetNumCells(); + + // SimplexVsBbox is the 3D default; in 2D it is the conservative bounding-box + // test, which is a legitimate thing to want. No other caller asks for it in + // 2D. + auto conservative = + pcms::detail::construct_intersection_map<2, pcms::SimplexVsBbox>( + mesh, grid_d, ncells); + require_maps_identical( + conservative, + pcms::detail::construct_intersection_map_reference<2, pcms::SimplexVsBbox>( + mesh, grid_d, ncells)); + + // The conservative test accepts everything the refined default accepts and + // more, wherever an element's box reaches a cell its triangle does not. If + // the predicate were ignored, the two maps would be identical. + auto refined = + pcms::detail::construct_intersection_map_2d(mesh, grid_d, ncells); + REQUIRE(conservative.entries.size() > refined.entries.size()); +} + +// The search itself must accept a predicate the same way, chosen once when it +// is built. The predicate must be one instantiated in point_search.cpp. +TEST_CASE("point search accepts a predicate other than its default") +{ + using pcms::GridPointSearch2D; + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false); + + auto tolerances = + GridPointSearch2D::PointSearchTolerances{"point search 2d tolerances"}; + auto tolerances_h = Kokkos::create_mirror_view(tolerances); + tolerances_h(0) = 0.01; + tolerances_h(1) = 0.01; + Kokkos::deep_copy(tolerances, tolerances_h); + + GridPointSearch2D default_search{mesh, 10, 10, tolerances}; + // The conservative 2D predicate instead of this search's default. + GridPointSearch2D custom_search{mesh, 10, 10, tolerances, + pcms::SimplexVsBbox{}}; + + Kokkos::View points("test_points", 4); + auto points_h = Kokkos::create_mirror_view(points); + points_h(0, 0) = 0.05; + points_h(0, 1) = 0.02; + points_h(1, 0) = 0.55; + points_h(1, 1) = 0.54; + points_h(2, 0) = 0.9; + points_h(2, 1) = 0.9; + points_h(3, 0) = 1.0; + points_h(3, 1) = 1.0; + Kokkos::deep_copy(points, points_h); + + auto default_results = default_search(points); + auto custom_results = custom_search(points); + auto default_h = Kokkos::create_mirror_view(default_results); + auto custom_h = Kokkos::create_mirror_view(custom_results); + Kokkos::deep_copy(default_h, default_results); + Kokkos::deep_copy(custom_h, custom_results); + + // The query tests point containment in each candidate, so the predicate only + // decides which elements are offered. For points inside the mesh both + // searches must therefore answer identically: the value-constructor path is a + // candidate-set choice, not a change of semantics. + for (size_t i = 0; i < 4; ++i) { + CAPTURE(i); + REQUIRE(custom_h(i).dimensionality == default_h(i).dimensionality); + REQUIRE(custom_h(i).element_id == default_h(i).element_id); + } +} + +// The 3D search takes a division count per axis. Unequal counts must locate +// points exactly as equal counts do: elements are binned by each cell's box, so +// a cell box computed for the wrong axis would drop candidates silently. +TEST_CASE("3D point search with unequal divisions per axis", + "[point_search][3d]") +{ + using pcms::GridPointSearch3D; + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 4, 4, 4, false); + GridPointSearch3D search{mesh, 3, 5, 7}; + + const auto n_elems = mesh.nelems(); + const auto centroids = + Omega_h::HostRead(pcms::get_entity_centroids(mesh, 3)); + Kokkos::View points("query points", n_elems); + auto points_h = Kokkos::create_mirror_view(points); + for (Omega_h::LO i = 0; i < n_elems; ++i) { + for (int d = 0; d < 3; ++d) { + points_h(i, d) = centroids[i * 3 + d]; + } + } + Kokkos::deep_copy(points, points_h); + + const auto owners_h = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace{}, search.GetOwningElementIds(search(points))); + Omega_h::LO n_wrong = 0; + for (Omega_h::LO i = 0; i < n_elems; ++i) { + n_wrong += owners_h(i) != i; + } + CAPTURE(n_wrong, n_elems); + REQUIRE(n_wrong == 0); +} + +// The grid and overlap layer under the point search is written for any +// dimension, though the search classes themselves are 2D and 3D only. These +// checks pin that down above 3D. Every grid has a different division count on +// each axis, and every probe uses different indices per axis, so an axis mix-up +// cannot hide behind a symmetric case. +namespace +{ + +/// Cells of `grid` whose id, per-axis indices and box disagree, counted on the +/// device: the same functions the search runs in its kernels. +template +pcms::LO CountCellRoundTripFailuresOnDevice(const pcms::UniformGrid& grid) +{ + pcms::LO failures = 0; + Kokkos::parallel_reduce( + "cell round trip", grid.GetNumCells(), + KOKKOS_LAMBDA(pcms::LO cell, pcms::LO & n) { + const auto axis = grid.AxisIndicesFromCellIndex(cell); + bool ok = grid.CellIndexFromAxisIndices(axis) == cell; + Omega_h::Vector center; + const auto box = grid.GetCellBBOX(cell); + for (unsigned d = 0; d < Dim; ++d) { + ok = ok && axis[d] >= 0 && axis[d] < grid.divisions[d]; + center[d] = box.center[d]; + } + ok = ok && grid.ClosestCellID(center) == cell; + n += ok ? 0 : 1; + }, + failures); + return failures; +} + +} // namespace + +TEST_CASE("grid and overlap layer above 3D", "[point_search][nd]") +{ + using pcms::LO; + + SECTION("4D grid: cell ids, per-axis indices and cell boxes agree") + { + // Cell size 1 on every axis, so a cell's box is easy to state. + const pcms::UniformGrid<4> grid{.edge_length = {2.0, 3.0, 4.0, 5.0}, + .bot_left = {0.0, 0.0, 0.0, 0.0}, + .divisions = {2, 3, 4, 5}}; + REQUIRE(grid.GetNumCells() == 120); + + // Ids number cells x-fastest: x + Nx * (y + Ny * (z + Nz * w)). + const Kokkos::Array last{1, 2, 3, 4}; + REQUIRE(grid.CellIndexFromAxisIndices(last) == 119); + const Kokkos::Array probe{1, 0, 2, 3}; + const LO probe_id = 1 + 2 * (0 + 3 * (2 + 4 * 3)); + REQUIRE(grid.CellIndexFromAxisIndices(probe) == probe_id); + const auto back = grid.AxisIndicesFromCellIndex(probe_id); + for (int d = 0; d < 4; ++d) { + REQUIRE(back[d] == probe[d]); + } + const auto box = grid.GetCellBBOX(probe_id); + for (int d = 0; d < 4; ++d) { + REQUIRE(box.center[d] == Catch::Approx(probe[d] + 0.5)); + REQUIRE(box.half_width[d] == Catch::Approx(0.5)); + } + + REQUIRE(CountCellRoundTripFailuresOnDevice(grid) == 0); + } + + SECTION("5D grid: every cell round-trips on the device") + { + const pcms::UniformGrid<5> grid{.edge_length = {1.0, 2.0, 3.0, 2.0, 1.0}, + .bot_left = {0.0, -1.0, 2.0, 0.5, 0.0}, + .divisions = {2, 3, 4, 3, 2}}; + REQUIRE(grid.GetNumCells() == 144); + REQUIRE(CountCellRoundTripFailuresOnDevice(grid) == 0); + } + + SECTION("4D candidate cells keep each axis separate") + { + const pcms::UniformGrid<4> grid{.edge_length = {2.0, 3.0, 4.0, 5.0}, + .bot_left = {0.0, 0.0, 0.0, 0.0}, + .divisions = {2, 3, 4, 5}}; + // Inside cell (0, 1, 2, 3); padded by one cell, clamped to the grid. + const AABBox<4> box{.center = {0.5, 1.5, 2.5, 3.5}, + .half_width = {0.25, 0.25, 0.25, 0.25}}; + const auto range = pcms::candidate_cells<4, pcms::SimplexVsBbox>(box, grid); + const Kokkos::Array lo{0, 0, 1, 2}; + const Kokkos::Array hi{1, 2, 3, 4}; + for (int d = 0; d < 4; ++d) { + CAPTURE(d); + REQUIRE(range.lo[d] == lo[d]); + REQUIRE(range.hi[d] == hi[d]); + } + REQUIRE(range.Contains(Kokkos::Array{0, 1, 2, 3})); + REQUIRE_FALSE(range.Contains(Kokkos::Array{0, 1, 2, 1})); + } + + SECTION("4D simplex-vs-box overlap needs every axis to overlap") + { + const AABBox<4> box{.center = {1.0, 1.0, 1.0, 1.0}, + .half_width = {0.5, 0.5, 0.5, 0.5}}; + // A simplex spanning [0.8, 1.2] on axes 0-2. On axis 3 it either reaches + // into the box or stops short of it. + const auto simplex = [](pcms::Real axis3_lo, pcms::Real axis3_hi) { + Omega_h::Matrix<4, 5> coords{}; + for (int v = 0; v < 5; ++v) { + for (int d = 0; d < 3; ++d) { + coords(d, v) = (v % 2 == 0) ? 0.8 : 1.2; + } + coords(3, v) = (v == 0) ? axis3_lo : axis3_hi; + } + return coords; + }; + REQUIRE( + pcms::Overlap<4, pcms::SimplexVsBbox>::apply(simplex(1.4, 2.0), box)); + // Overlapping on axes 0-2 but entirely past the box on axis 3. + REQUIRE_FALSE( + pcms::Overlap<4, pcms::SimplexVsBbox>::apply(simplex(1.6, 2.0), box)); + } +} diff --git a/test/test_uniform_grid.cpp b/test/test_uniform_grid.cpp index d580c9f68..7fcdc4c28 100644 --- a/test/test_uniform_grid.cpp +++ b/test/test_uniform_grid.cpp @@ -53,6 +53,40 @@ TEST_CASE("uniform grid") REQUIRE(0 == uniform_grid.GetCellIndex({0, 0})); REQUIRE(119 == uniform_grid.GetCellIndex({11, 9})); } + SECTION("AxisIndicesFromCellIndex") + { + // The inverse of CellIndexFromAxisIndices, in coordinate order. Cell 1 is + // the second cell along x, which the cell-BBOX test below also pins down + // (its center is (1.5, 0.5)): the grid stores its indices reversed relative + // to coordinates, so getting this wrong is silent. + { + const auto [i, j] = uniform_grid.AxisIndicesFromCellIndex(0); + REQUIRE(i == 0); + REQUIRE(j == 0); + } + { + const auto [i, j] = uniform_grid.AxisIndicesFromCellIndex(1); + REQUIRE(i == 1); + REQUIRE(j == 0); + } + { + const auto [i, j] = uniform_grid.AxisIndicesFromCellIndex(119); + REQUIRE(i == 9); + REQUIRE(j == 11); + } + + // Round-trip over every cell of the grid, which is the property a caller + // keyed by linear cell id depends on. + for (pcms::LO cell = 0; cell < uniform_grid.GetNumCells(); ++cell) { + const auto axis_indices = uniform_grid.AxisIndicesFromCellIndex(cell); + REQUIRE(uniform_grid.CellIndexFromAxisIndices(axis_indices) == cell); + // ...and it must agree with the binning rule used to find a point's cell. + for (int d = 0; d < 2; ++d) { + const auto center = uniform_grid.GetCellBBOX(cell).center[d]; + REQUIRE(uniform_grid.AxisCellIndex(d, center) == axis_indices[d]); + } + } + } } TEST_CASE("Create uniform grid from omega_h mesh") @@ -170,3 +204,28 @@ TEST_CASE("Create uniform grid from 3D mesh") REQUIRE(grid.divisions[2] == 12); } } + +// Grid resolution has to be a property of the mesh as a whole, not of one +// partition: search cells are named by a linear id, and an id identifies the +// same region of space on two ranks only if both resolved the grid the same +// way. A resolution derived from `mesh.nelems()` breaks that, because +// `nelems()` counts ghost copies and therefore changes with the partition. +TEST_CASE("cell index round trip on a grid with unequal divisions") +{ + const UniformGrid<3> grid{.edge_length = {3.0, 4.0, 5.0}, + .bot_left = {0.0, 0.0, 0.0}, + .divisions = {3, 4, 5}}; + for (pcms::LO cell = 0; cell < grid.GetNumCells(); ++cell) { + CAPTURE(cell); + const auto axis_indices = grid.AxisIndicesFromCellIndex(cell); + for (int d = 0; d < 3; ++d) { + REQUIRE(axis_indices[d] >= 0); + REQUIRE(axis_indices[d] < grid.divisions[d]); + } + REQUIRE(grid.CellIndexFromAxisIndices(axis_indices) == cell); + // The cell's own box must bin back to the same cell. + const auto box = grid.GetCellBBOX(cell); + REQUIRE(grid.ClosestCellID(Omega_h::Vector<3>{box.center[0], box.center[1], + box.center[2]}) == cell); + } +}