From bd5ee1977eab791ffb94b56785e485a6aa96fc52 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 12 Mar 2026 10:14:10 -0400 Subject: [PATCH 1/7] Remove `fuzz` parameter from 2D and 3D grid point search methods and related call sites for simplification. --- src/pcms/point_search.cpp | 56 ++++++++++++++++----------------------- src/pcms/point_search.h | 16 +++++------ 2 files changed, 29 insertions(+), 43 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index eb135916..f00199b6 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -144,26 +144,26 @@ template } [[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( - const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords, Real fuzz) + const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords) { auto left = bbox.center[0] - bbox.half_width[0]; auto right = bbox.center[0] + bbox.half_width[0]; auto bot = bbox.center[1] - bbox.half_width[1]; auto top = bbox.center[1] + bbox.half_width[1]; auto xi = Omega_h::barycentric_from_global<2, 2>({left, bot}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { + if (Omega_h::is_barycentric_inside(xi)) { return true; } xi = Omega_h::barycentric_from_global<2, 2>({left, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { + if (Omega_h::is_barycentric_inside(xi)) { return true; } xi = Omega_h::barycentric_from_global<2, 2>({right, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { + if (Omega_h::is_barycentric_inside(xi)) { return true; } xi = Omega_h::barycentric_from_global<2, 2>({right, bot}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { + if (Omega_h::is_barycentric_inside(xi)) { return true; } return false; @@ -197,7 +197,7 @@ template 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, fuzz)) { + if (Omega_h::is_barycentric_inside(xi)) { return true; } } @@ -210,7 +210,7 @@ template */ [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( - const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox, Real fuzz) + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) { // triangle and grid cell bounding box intersect if (intersects(triangle_bbox(coords), bbox)) { @@ -220,7 +220,7 @@ KOKKOS_FUNCTION bool triangle_intersects_bbox( return true; } // if any of the bbox verts are within the triangle - if (bbox_verts_within_triangle(bbox, coords, fuzz)) { + if (bbox_verts_within_triangle(bbox, coords)) { return true; } // if any of the triangle's edges intersect with the bounding box @@ -257,12 +257,11 @@ namespace detail struct GridTriIntersectionFunctor2D { GridTriIntersectionFunctor2D(Omega_h::Mesh& mesh, - Kokkos::View grid, Real fuzz) + Kokkos::View grid) : mesh_(mesh), tris2verts_(mesh_.ask_elem_verts()), coords_(mesh_.coords()), grid_(grid), - fuzz_(fuzz), nelems_(mesh_.nelems()) { if (mesh_.dim() != 2) { @@ -286,7 +285,7 @@ struct GridTriIntersectionFunctor2D // 2d mesh with 2d coords, but 3 triangles const auto vertex_coords = Omega_h::gather_vectors<3, 2>(coords_, elem_tri2verts); - if (triangle_intersects_bbox(vertex_coords, grid_cell_bbox, fuzz_)) { + if (triangle_intersects_bbox(vertex_coords, grid_cell_bbox)) { if (fill) { fill[num_intersections] = elem_idx; } @@ -301,7 +300,6 @@ struct GridTriIntersectionFunctor2D Omega_h::LOs tris2verts_; Omega_h::Reals coords_; Kokkos::View grid_; - Real fuzz_; public: LO nelems_; @@ -362,10 +360,10 @@ struct GridTriIntersectionFunctor3D Kokkos::Crs construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, - int num_grid_cells, Real fuzz) + int num_grid_cells) { Kokkos::Crs intersection_map{}; - auto f = detail::GridTriIntersectionFunctor2D{mesh, grid, fuzz}; + auto f = detail::GridTriIntersectionFunctor2D{mesh, grid}; Kokkos::count_and_fill_crs(intersection_map, num_grid_cells, f); return intersection_map; } @@ -407,7 +405,6 @@ Kokkos::View GridPointSearch2D::operator()( auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; auto tolerances = tolerances_; - auto fuzz = fuzz_; Kokkos::parallel_for( points.extent(0), KOKKOS_LAMBDA(int p) { Omega_h::Vector<2> point( @@ -493,7 +490,7 @@ Kokkos::View GridPointSearch2D::operator()( if (edge_found) break; - if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { + if (Omega_h::is_barycentric_inside(parametric_coords)) { dimensionality = GridPointSearch2D::Result::Dimensionality::FACE; nearest_element_id = triangleID; parametric_coords_to_nearest = parametric_coords; @@ -511,17 +508,15 @@ Kokkos::View GridPointSearch2D::operator()( return results; } -GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, - Real fuzz) +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) : GridPointSearch2D(mesh, Nx, Ny, - PointSearchTolerances{"point search 2d tolerances"}, fuzz) + PointSearchTolerances{"point search 2d tolerances"}) { - Kokkos::deep_copy(tolerances_, 0); + Kokkos::deep_copy(tolerances_, 1E-12); } GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, - const PointSearchTolerances& tolerances, - Real fuzz) + const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); @@ -532,14 +527,13 @@ 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); - candidate_map_ = detail::construct_intersection_map_2d( - mesh, grid_, grid_h(0).GetNumCells(), fuzz_); + candidate_map_ = + detail::construct_intersection_map_2d(mesh, grid_, grid_h(0).GetNumCells()); coords_ = mesh.coords(); tris2verts_ = mesh.ask_elem_verts(); tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); - fuzz_ = fuzz; } Kokkos::View GridPointSearch3D::operator()( @@ -556,7 +550,6 @@ Kokkos::View GridPointSearch3D::operator()( auto tris2edges_adj = tris2edges_adj_; auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; - auto fuzz = fuzz_; Kokkos::parallel_for( points.extent(0), KOKKOS_LAMBDA(int p) { Omega_h::Vector point; @@ -585,7 +578,7 @@ Kokkos::View GridPointSearch3D::operator()( auto parametric_coords = Omega_h::barycentric_from_global(point, vertex_coords); - if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { + if (Omega_h::is_barycentric_inside(parametric_coords)) { results(p) = GridPointSearch3D::Result{ GridPointSearch3D::Result::Dimensionality::REGION, triangleID, parametric_coords}; @@ -605,17 +598,15 @@ Kokkos::View GridPointSearch3D::operator()( return results; } -GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, - Real fuzz) +GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) : GridPointSearch3D(mesh, Nx, Ny, Nz, - PointSearchTolerances{"point search 3d tolerances"}, fuzz) + PointSearchTolerances{"point search 3d tolerances"}) { Kokkos::deep_copy(tolerances_, 0); } GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, - const PointSearchTolerances& tolerances, - Real fuzz) + const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); @@ -641,6 +632,5 @@ GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); - fuzz_ = fuzz; } } // namespace pcms diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index d68f4432..20904ce3 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -22,12 +22,11 @@ namespace detail Kokkos::Crs construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, - int num_grid_cells, Real fuzz = 1E-12); + int num_grid_cells); } [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( - const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox, - Real fuzz = 1E-12); + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); template class PointLocalizationSearch @@ -77,9 +76,9 @@ class GridPointSearch2D : public PointLocalizationSearch2D public: using Result = PointLocalizationSearch2D::Result; - GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, Real fuzz = 1E-12); + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny); GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, - const PointSearchTolerances& tolerances, Real fuzz = 1E-12); + const PointSearchTolerances& tolerances); /** * given a point in global coordinates give the id of the triangle that the @@ -93,7 +92,6 @@ class GridPointSearch2D : public PointLocalizationSearch2D private: Omega_h::Mesh mesh_; - Real fuzz_; Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; Omega_h::Adj edges2verts_adj_; @@ -111,10 +109,9 @@ class GridPointSearch3D : public PointLocalizationSearch3D public: using Result = PointLocalizationSearch3D::Result; + GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz); GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, - Real fuzz = 1E-12); - GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, - const PointSearchTolerances& tolerances, Real fuzz = 1E-12); + const PointSearchTolerances& tolerances); /** * Given a point in global coordinates, returns the id of the tetrahedron (3D @@ -135,7 +132,6 @@ class GridPointSearch3D : public PointLocalizationSearch3D CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; Omega_h::Reals coords_; - Real fuzz_; }; } // namespace pcms From a57bbf8063843a201301a0e5326644ed43ff7e4c Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 1 Apr 2026 20:14:23 -0400 Subject: [PATCH 2/7] Add `distance_to_closest_edge_from_barycentric` function and corresponding tests for triangle edge distance calculation. --- src/pcms/point_search.h | 64 ++++++++++++++++++++++++++++++++++++++ test/test_point_search.cpp | 61 ++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+) diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 20904ce3..a6ea7813 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -6,6 +6,7 @@ #include #include #include +#include #include "pcms/utility/types.h" #include "pcms/uniform_grid.h" @@ -28,6 +29,69 @@ construct_intersection_map_2d(Omega_h::Mesh& mesh, [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); +/** + * Compute the Euclidean distance from a point to the closest edge of a triangle + * using only its barycentric coordinates with respect to that triangle. + * + * Given a triangle with vertex coordinates `coords` and a point represented by + * barycentric coordinates `xi` (with respect to those vertices), the distance + * from the point to the edge opposite vertex i is `xi[i] * h_i`, where `h_i` + * is the altitude from vertex i to its opposite edge. The altitude can be + * computed from the triangle geometry as `h_i = 2A / |e_i|`, where `A` is the + * triangle area and `|e_i|` is the length of the edge opposite vertex i. + * + * This function generalizes to N-dimensional embeddings of a 2-simplex + * (triangle) by computing the triangle area from two edge vectors using the + * parallelogram formula valid in any dimension: + * area = 0.5 * sqrt(||a||^2 ||b||^2 - (a·b)^2) + * + * Template parameter N is the embedding dimension (e.g., 2 for planar, 3 for + * spatial), while the simplex is always a triangle (3 vertices). + */ +template +[[nodiscard]] KOKKOS_FUNCTION inline Real +distance_to_closest_edge_from_barycentric( + const Omega_h::Matrix& coords, const Omega_h::Vector<3>& xi) +{ + // Edge vectors originating from vertex 0 + Omega_h::Vector a = coords[1] - coords[0]; + Omega_h::Vector b = coords[2] - coords[0]; + + // Parallelogram area magnitude squared = ||a||^2 ||b||^2 - (a·b)^2 + const Real aa = Omega_h::inner_product(a, a); + const Real bb = Omega_h::inner_product(b, b); + const Real ab = Omega_h::inner_product(a, b); + Real parallelogram_area_sq = aa * bb - ab * ab; + + // Guard against degeneracy and tiny negative due to FP roundoff + if (parallelogram_area_sq <= static_cast(0)) return 0; + const Real area = static_cast(0.5) * std::sqrt(parallelogram_area_sq); + const Real two_area = static_cast(2) * area; + + // Opposite edge lengths as a vector [|e(1,2)|, |e(0,2)|, |e(0,1)|] + Omega_h::Vector<3> ell; + for (int i = 0; i < 3; ++i) { + const int j = (i + 1) % 3; + const int k = (i + 2) % 3; + ell[i] = Omega_h::norm(coords[k] - coords[j]); + } + + // Altitudes h_i = 2A / |e_i| with degeneracy guard, stored as a vector + const Real eps = static_cast(1e-15); + Omega_h::Vector<3> h; + Omega_h::Vector<3> d; + for (int i = 0; i < 3; ++i) { + const Real hi = (ell[i] > eps) ? (two_area / ell[i]) : static_cast(0); + h[i] = hi; + d[i] = xi[i] * hi; // Distances to opposite edges using barycentric weights + } + + // Return the minimum component (avoid std::min for device-compatibility) + Real m = d[0]; + for (int i = 1; i < 3; ++i) m = (d[i] < m) ? d[i] : m; + return m; +} + template class PointLocalizationSearch { diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 58d5100e..f3d21b0c 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -100,6 +100,67 @@ TEST_CASE("Triangle BBox Intersection Regression") true); } +TEST_CASE("barycentric distance to closest edge") +{ + using pcms::distance_to_closest_edge_from_barycentric; + // Triangle embedded in 2D + Omega_h::Matrix<2, 3> coords{{0.0, 0.0}, {1.0, 0.0}, {0.5, 1.0}}; + + auto point_from_bary = [&](const Omega_h::Vector<3>& xi) { + Omega_h::Vector<2> p{0.0, 0.0}; + for (int i = 0; i < 3; ++i) { + p[0] += xi[i] * coords[i][0]; + p[1] += xi[i] * coords[i][1]; + } + return p; + }; + + auto dist_point_to_seg = [](const Omega_h::Vector<2>& p, + const Omega_h::Vector<2>& a, + const Omega_h::Vector<2>& b) { + auto ab = b - a; + auto ap = p - a; + auto ab2 = Omega_h::inner_product(ab, ab); + if (ab2 == 0.0) return Omega_h::norm(ap); + double t = Omega_h::inner_product(ap, ab) / ab2; + if (t < 0.0) t = 0.0; + if (t > 1.0) t = 1.0; + auto proj = a + t * ab; + return Omega_h::norm(p - proj); + }; + + auto check = [&](const Omega_h::Vector<3>& xi) { + auto p = point_from_bary(xi); + // Euclidean distance to each triangle edge + double d01 = dist_point_to_seg(p, coords[0], coords[1]); + double d12 = dist_point_to_seg(p, coords[1], coords[2]); + double d20 = dist_point_to_seg(p, coords[2], coords[0]); + double d_true = std::min(d01, std::min(d12, d20)); + + double d_bary = distance_to_closest_edge_from_barycentric<2>(coords, xi); + REQUIRE(d_bary == Catch::Approx(d_true).epsilon(1e-12)); + }; + + SECTION("interior points") + { + check(Omega_h::Vector<3>{1.0 / 3, 1.0 / 3, 1.0 / 3}); + check(Omega_h::Vector<3>{0.25, 0.25, 0.5}); + check(Omega_h::Vector<3>{0.1, 0.6, 0.3}); + check(Omega_h::Vector<3>{0.49, 0.49, 0.02}); + } + SECTION("on edges and at vertices") + { + // Edge midpoints + check(Omega_h::Vector<3>{0.0, 0.5, 0.5}); + check(Omega_h::Vector<3>{0.5, 0.0, 0.5}); + check(Omega_h::Vector<3>{0.5, 0.5, 0.0}); + // Vertices + check(Omega_h::Vector<3>{1.0, 0.0, 0.0}); + check(Omega_h::Vector<3>{0.0, 1.0, 0.0}); + check(Omega_h::Vector<3>{0.0, 0.0, 1.0}); + } +} + template bool num_candidates_within_range(const T& intersection_map, pcms::LO min, pcms::LO max) From 8afd7d5f32c0c92115b038799871325e31fcc6cb Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 1 Apr 2026 21:42:03 -0400 Subject: [PATCH 3/7] Remove inflation radius logic from 2D/3D intersection map construction and corresponding `construct_intersection_map` overloads for simplification. Refactor related grid search logic and adjust formatting for consistency. Replace `is_barycentric_inside_with_tolerance` with `Omega_h::is_barycentric_inside` in point search logic for simplified tolerance handling. Remove unused function implementation. Replace `element_id` with `face_id` across point search logic for improved entity identification consistency. Refactor tolerance checks and neighbor-cell fallback in `GridPointSearch2D`. Update related tests and adapt internal structures. Add compatibility mode to `GridPointSearch2D` for callers expecting face-based IDs. Update edge classification logic with min/max tracking and refine distance-based hierarchy. Add inflation radius to grid cell bounding boxes in 2D/3D intersection map construction for tolerance handling. Refactor candidate filtering logic in `GridPointSearch2D`. Add overloads for backward compatibility. --- src/pcms/adapter/omega_h/omega_h_field.h | 16 +- src/pcms/adapter/omega_h/omega_h_field2.cpp | 39 ++-- src/pcms/interpolator/adj_search.hpp | 10 +- src/pcms/point_search.cpp | 238 ++++++++++++++------ src/pcms/point_search.h | 16 +- test/test_point_search.cpp | 13 +- 6 files changed, 222 insertions(+), 110 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index 53af1256..1f361121 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -344,11 +344,11 @@ auto evaluate(const OmegaHField& field, Lagrange<1> /* method */, Kokkos::parallel_for( results.size(), KOKKOS_LAMBDA(LO i) { - auto [dim, elem_idx, coord] = results(i); - // TODO deal with case for elem_idx < 0 (point outside of mesh) - KOKKOS_ASSERT(elem_idx >= 0); + auto [dim, elem_idx, face_idx, coord] = results(i); + // TODO deal with case for face_idx < 0 (point outside of mesh) + KOKKOS_ASSERT(face_idx >= 0); const auto elem_tri2verts = - Omega_h::gather_verts<3>(tris2verts, elem_idx); + Omega_h::gather_verts<3>(tris2verts, face_idx); Real val = 0; for (int j = 0; j < 3; ++j) { val += field_values[elem_tri2verts[j]] * coord[j]; @@ -382,11 +382,11 @@ auto evaluate(const OmegaHField& field, NearestNeighbor /* method */, Kokkos::parallel_for( results.size(), KOKKOS_LAMBDA(LO i) { - auto [dim, elem_idx, coord] = results(i); - // TODO deal with case for elem_idx < 0 (point outside of mesh) - KOKKOS_ASSERT(elem_idx >= 0); + auto [dim, elem_idx, face_idx, coord] = results(i); + // TODO deal with case for face_idx < 0 (point outside of mesh) + KOKKOS_ASSERT(face_idx >= 0); const auto elem_tri2verts = - Omega_h::gather_verts<3>(tris2verts, elem_idx); + Omega_h::gather_verts<3>(tris2verts, face_idx); // value is closest to point has the largest coordinate int vert = 0; auto max_val = coord[0]; diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index 36db8729..f2370a4f 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -109,8 +109,8 @@ struct CountPointsPerElementFunctor KOKKOS_INLINE_FUNCTION void operator()(LO i) const { - auto [dim, elem_idx, coord] = search_results_(i); - Kokkos::atomic_add(&elem_counts_(elem_idx), 1); + auto [dim, elem_idx, face_idx, coord] = search_results_(i); + Kokkos::atomic_add(&elem_counts_(face_idx), 1); } }; @@ -142,14 +142,14 @@ struct FillCoordinatesAndIndicesFunctor KOKKOS_INLINE_FUNCTION void operator()(LO i) const { - auto [dim, elem_idx, coord] = search_results_(i); + auto [dim, elem_idx, face_idx, coord] = search_results_(i); // disable the host assertion macro for device code // currently don't handle case where point is on a boundary // PCMS_ALWAYS_ASSERT(static_cast(dim) == mesh_.dim()); - // element should be inside the domain (positive) - // PCMS_ALWAYS_ASSERT(elem_idx >= 0 && elem_idx < mesh_.nelems()); - LO count = Kokkos::atomic_sub_fetch(&elem_counts_(elem_idx), 1); - LO index = offsets_(elem_idx) + count - 1; + // face should be inside the domain (positive) + // PCMS_ALWAYS_ASSERT(face_idx >= 0 && face_idx < mesh_.nelems()); + LO count = Kokkos::atomic_sub_fetch(&elem_counts_(face_idx), 1); + LO index = offsets_(face_idx) + count - 1; for (int j = 0; j < (dim_ + 1); ++j) { coordinates_(index, j) = coord[j]; } @@ -172,18 +172,20 @@ struct OmegaHField2LocalizationHint if (mode_ == OutOfBoundsMode::ERROR) { // Error mode - throw error immediately if any point is out of bounds for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, coord] = search_results(i); - bool is_missing = - (static_cast(dim) != mesh.dim()) || (elem_idx < 0); + auto [dim, elem_idx, face_idx, coord] = search_results(i); + (void)dim; + (void)coord; + bool is_missing = (elem_idx < 0) || (face_idx < 0); PCMS_ALWAYS_ASSERT(!is_missing && "Points found outside mesh domain"); valid_point_indices.push_back(i); } } else { // Other modes - collect valid and missing points separately for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, coord] = search_results(i); - bool is_missing = - (static_cast(dim) != mesh.dim()) || (elem_idx < 0); + auto [dim, elem_idx, face_idx, coord] = search_results(i); + (void)dim; + (void)coord; + bool is_missing = (elem_idx < 0) || (face_idx < 0); if (is_missing) { missing_point_indices.push_back(i); } else { @@ -218,9 +220,10 @@ struct OmegaHField2LocalizationHint // Count points per element (valid points only) Kokkos::View elem_counts("elem_counts", mesh.nelems()); + Kokkos::deep_copy(elem_counts, 0); for (size_t i = 0; i < num_valid_; ++i) { - auto [dim, elem_idx, coord] = search_results(valid_point_indices[i]); - elem_counts[elem_idx] += 1; + auto [dim, elem_idx, face_idx, coord] = search_results(valid_point_indices[i]); + elem_counts[face_idx] += 1; } // Compute offsets @@ -235,9 +238,9 @@ struct OmegaHField2LocalizationHint // Fill coordinates and indices for valid points for (size_t i = 0; i < num_valid_; ++i) { size_t orig_idx = valid_point_indices[i]; - auto [dim, elem_idx, coord] = search_results(orig_idx); - elem_counts(elem_idx) -= 1; - LO index = offsets_(elem_idx) + elem_counts(elem_idx); + auto [dim, elem_idx, face_idx, coord] = search_results(orig_idx); + elem_counts(face_idx) -= 1; + LO index = offsets_(face_idx) + elem_counts(face_idx); for (int j = 0; j < (mesh.dim() + 1); ++j) { coordinates_(index, j) = coord[j]; } diff --git a/src/pcms/interpolator/adj_search.hpp b/src/pcms/interpolator/adj_search.hpp index 85baa236..2b56bc11 100644 --- a/src/pcms/interpolator/adj_search.hpp +++ b/src/pcms/interpolator/adj_search.hpp @@ -32,9 +32,9 @@ inline void checkTargetPoints( pcms::printInfo("INFO: Checking target points...\n"); auto check_target_points = OMEGA_H_LAMBDA(Omega_h::LO i) { - if (results(i).element_id < 0) { - OMEGA_H_CHECK_PRINTF(results(i).element_id >= 0, - "ERROR: Source cell id not found for target %d\n", + if (results(i).face_id < 0) { + OMEGA_H_CHECK_PRINTF(results(i).face_id >= 0, + "ERROR: Source face id not found for target %d\n", i); printf("%d, ", i); } @@ -130,10 +130,10 @@ inline void FindSupports::adjBasedSearch( Track visited; Omega_h::Real cutoffDistance = radii2[id]; - Omega_h::LO source_cell_id = results(id).element_id; + Omega_h::LO source_cell_id = results(id).face_id; OMEGA_H_CHECK_PRINTF( source_cell_id >= 0, - "ERROR: Source cell id not found for target %d (%f,%f)\n", id, + "ERROR: Source face id not found for target %d (%f,%f)\n", id, target_points(id, 0), target_points(id, 1)); const Omega_h::LO num_verts_in_dim = dim + 1; diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index f00199b6..77380219 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -1,6 +1,7 @@ #include "point_search.h" #include #include +#include // From // https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation @@ -276,7 +277,7 @@ struct GridTriIntersectionFunctor2D KOKKOS_INLINE_FUNCTION LO operator()(LO row, LO* fill) const { - const auto grid_cell_bbox = grid_(0).GetCellBBOX(row); + auto grid_cell_bbox = grid_(0).GetCellBBOX(row); LO num_intersections = 0; // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { @@ -327,7 +328,7 @@ struct GridTriIntersectionFunctor3D KOKKOS_INLINE_FUNCTION LO operator()(LO row, LO* fill) const { - const auto grid_cell_bbox = grid_(0).GetCellBBOX(row); + auto grid_cell_bbox = grid_(0).GetCellBBOX(row); LO num_intersections = 0; // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { @@ -414,95 +415,187 @@ Kokkos::View GridPointSearch2D::operator()( auto candidates_begin = candidate_map.row_map(cell_id); auto candidates_end = candidate_map.row_map(cell_id + 1); - bool vertex_found = false; - bool edge_found = false; - bool inside_cell = false; - - auto nearest_element_id = candidates_begin; - auto dimensionality = GridPointSearch2D::Result::Dimensionality::REGION; - Omega_h::Real distance_to_nearest{INFINITY}; - Omega_h::Vector<3> parametric_coords_to_nearest; - // create array that's size of number of candidates x num coords to store - // parametric inversion - for (auto i = candidates_begin; i < candidates_end; ++i) { - const int triangleID = candidate_map.entries(i); + // Track best entities across all candidates to ensure order invariance + Omega_h::Real best_vertex_dist = INFINITY; + LO best_vertex_id = -1; + int best_vertex_tid = -1; // triangle providing barycentric coords + Omega_h::Vector<3> best_vertex_bary{0.0, 0.0, 0.0}; + + Omega_h::Real best_edge_dist = INFINITY; + LO best_edge_id = -1; + int best_edge_tid_min = -1; + int best_edge_tid_max = -1; + Omega_h::Vector<3> best_edge_bary_min{0.0, 0.0, 0.0}; + Omega_h::Vector<3> best_edge_bary_max{0.0, 0.0, 0.0}; + + bool found_inside = false; + LO inside_face_id = -1; + Omega_h::Vector<3> inside_face_bary{0.0, 0.0, 0.0}; + + auto begin = candidate_map.row_map(cell_id); + auto end = candidate_map.row_map(cell_id + 1); + for (auto ii = begin; ii < end; ++ii) { + const int triangleID = candidate_map.entries(ii); const auto elem_tri2verts = Omega_h::gather_verts<3>(tris2verts, triangleID); - // 2d mesh with 2d coords, but 3 triangles auto vertex_coords = Omega_h::gather_vectors<3, 2>(coords, elem_tri2verts); auto parametric_coords = Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); - // Every triangle (face) is connected to 3 vertices + // Check vertices (hierarchy level 1): compute Euclidean distance for (int j = 0; j < 3; ++j) { - // Get the vertex ID from the connectivity array const int vertexID = tris2verts_adj.ab2b[triangleID * 3 + j]; - // Get the vertex coordinates from the mesh using vertexID - const Omega_h::Few vertex = - Omega_h::get_vector<2>(coords, vertexID); - - const auto distance = Omega_h::norm(point - vertex); - - if (distance < distance_to_nearest) { - dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; - nearest_element_id = vertexID; - distance_to_nearest = distance; - parametric_coords_to_nearest = parametric_coords; - - if (distance < tolerances(0)) { - vertex_found = true; - }; + const auto v = Omega_h::get_vector<2>(coords, vertexID); + const auto dv = Omega_h::norm(point - v); + if ((dv < best_vertex_dist) || + ((dv == best_vertex_dist) && (vertexID < best_vertex_id))) { + best_vertex_dist = dv; + best_vertex_id = vertexID; + best_vertex_tid = triangleID; + best_vertex_bary = parametric_coords; } } - if (vertex_found) - break; - + // Check edges (hierarchy level 2): only if projection falls within for (int j = 0; j < 3; ++j) { - // Every triangle (face) is connected to 3 edges const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; + const int va_id = edges2verts_adj.ab2b[edgeID * 2 + 0]; + const int vb_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; + const auto va = Omega_h::get_vector<2>(coords, va_id); + const auto vb = Omega_h::get_vector<2>(coords, vb_id); - auto vertex_a_id = edges2verts_adj.ab2b[edgeID * 2]; - auto vertex_b_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; - - auto vertex_a = Omega_h::get_vector<2>(coords, vertex_a_id); - auto vertex_b = Omega_h::get_vector<2>(coords, vertex_b_id); - - if (!normal_intersects_segment(vertex_a, vertex_b, point)) + if (!normal_intersects_segment(va, vb, point)) continue; - const auto distance_to_ab = - distance_from_line(vertex_a, vertex_b, point); - - if (distance_to_ab < distance_to_nearest) { - dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; - nearest_element_id = edgeID; - distance_to_nearest = distance_to_ab; - parametric_coords_to_nearest = parametric_coords; + const auto de = distance_from_line(va, vb, point); + if ((de < best_edge_dist) || + ((de == best_edge_dist) && (edgeID < best_edge_id)) || + ((de == best_edge_dist) && (edgeID == best_edge_id) && + (triangleID < best_edge_tid_min))) { + best_edge_dist = de; + best_edge_id = edgeID; + best_edge_tid_min = triangleID; + best_edge_tid_max = triangleID; + best_edge_bary_min = parametric_coords; + best_edge_bary_max = parametric_coords; + } else if ((de == best_edge_dist) && (edgeID == best_edge_id)) { + // Track the extents of face IDs sharing this nearest edge + if (triangleID < best_edge_tid_min) { + best_edge_tid_min = triangleID; + best_edge_bary_min = parametric_coords; + } + if (triangleID > best_edge_tid_max) { + best_edge_tid_max = triangleID; + best_edge_bary_max = parametric_coords; + } + } + } - if (distance_to_ab < tolerances(1)) { - edge_found = true; - }; + // Check face interior (hierarchy level 3) + if (Omega_h::is_barycentric_inside(parametric_coords)) { + if (!found_inside || (triangleID < inside_face_id)) { + found_inside = true; + inside_face_id = triangleID; + inside_face_bary = parametric_coords; } } + } - if (edge_found) - break; + const auto vtol = tolerances(0); + const auto etol = tolerances(1); + + // If we found an inside face, compute its edge distance using + // barycentric-only helper to check for edge classification while + // preserving the containing face ID. + Real inside_edge_dist = INFINITY; + int inside_edge_argmin = -1; + Omega_h::Matrix<2, 3> vcoords_in; + // Track Euclidean-nearest edge of the containing face (if any) + LO inside_edge_id = -1; + if (found_inside) { + const auto elem_tri2verts_in = + Omega_h::gather_verts<3>(tris2verts, inside_face_id); + vcoords_in = Omega_h::gather_vectors<3, 2>(coords, elem_tri2verts_in); + // Euclidean distance to the 3 edges of the containing triangle + inside_edge_argmin = -1; + inside_edge_dist = INFINITY; + inside_edge_id = -1; + for (int j = 0; j < 3; ++j) { + const int edgeID = tris2edges_adj.ab2b[inside_face_id * 3 + j]; + const int va_id = edges2verts_adj.ab2b[edgeID * 2 + 0]; + const int vb_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; + const auto va = Omega_h::get_vector<2>(coords, va_id); + const auto vb = Omega_h::get_vector<2>(coords, vb_id); + if (!normal_intersects_segment(va, vb, point)) + continue; + const auto de = distance_from_line(va, vb, point); + if (de < inside_edge_dist) { + inside_edge_dist = de; + inside_edge_argmin = j; + inside_edge_id = edgeID; + } + } + } - if (Omega_h::is_barycentric_inside(parametric_coords)) { - dimensionality = GridPointSearch2D::Result::Dimensionality::FACE; - nearest_element_id = triangleID; - parametric_coords_to_nearest = parametric_coords; - inside_cell = true; + GridPointSearch2D::Result::Dimensionality dim_out = + GridPointSearch2D::Result::Dimensionality::REGION; + LO element_id_out = -1; + LO face_id_out = -1; + Omega_h::Vector<3> bary_out{0.0, 0.0, 0.0}; + + // Apply hierarchy with tolerances + // Points within tolerance are considered "inside" with positive IDs + // Only points with no candidates at all get negative IDs + if (best_vertex_id >= 0 && best_vertex_dist <= vtol) { + // Point within vertex tolerance - still inside the mesh + dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; + element_id_out = best_vertex_id; + face_id_out = (best_vertex_tid >= 0) ? best_vertex_tid : -1; + bary_out = best_vertex_bary; + } else if ((best_edge_id >= 0 && best_edge_dist <= etol) || + (found_inside && inside_edge_dist <= etol)) { + // Point within edge tolerance - still inside the mesh + dim_out = GridPointSearch2D::Result::Dimensionality::EDGE; + if (found_inside && inside_edge_dist <= etol && + inside_edge_argmin >= 0) { + element_id_out = inside_edge_id; + face_id_out = inside_face_id; + bary_out = inside_face_bary; + } else { + element_id_out = best_edge_id; + face_id_out = (best_edge_tid_min >= 0) ? best_edge_tid_min : -1; + bary_out = best_edge_bary_min; + } + } else if (found_inside) { + // Point inside face + dim_out = GridPointSearch2D::Result::Dimensionality::FACE; + element_id_out = inside_face_id; + face_id_out = inside_face_id; + bary_out = inside_face_bary; + } else { + // Outside mesh - beyond tolerance of any entity + // Only negate if we truly have no candidates at all. + if (best_vertex_id >= 0 && + (best_vertex_dist <= best_edge_dist || best_edge_id < 0)) { + dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; + element_id_out = best_vertex_id; + face_id_out = (best_vertex_tid >= 0) ? best_vertex_tid : -1; + bary_out = best_vertex_bary; + element_id_out = -element_id_out; + } else if (best_edge_id >= 0) { + dim_out = GridPointSearch2D::Result::Dimensionality::EDGE; + element_id_out = best_edge_id; + face_id_out = (best_edge_tid_min >= 0) ? best_edge_tid_min : -1; + bary_out = best_edge_bary_min; + element_id_out = -element_id_out; + } else { + // No candidates at all: both IDs stay -1 } } - const int inside_mesh = - vertex_found || edge_found || inside_cell ? 1 : -1; - results(p) = GridPointSearch2D::Result{dimensionality, - inside_mesh * nearest_element_id, - parametric_coords_to_nearest}; + results(p) = GridPointSearch2D::Result{dim_out, element_id_out, + face_id_out, bary_out}; }); return results; @@ -527,6 +620,9 @@ 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()); coords_ = mesh.coords(); @@ -534,6 +630,7 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); + edges2faces_up_ = mesh.ask_up(Omega_h::EDGE, Omega_h::FACE); } Kokkos::View GridPointSearch3D::operator()( @@ -581,7 +678,7 @@ Kokkos::View GridPointSearch3D::operator()( if (Omega_h::is_barycentric_inside(parametric_coords)) { results(p) = GridPointSearch3D::Result{ GridPointSearch3D::Result::Dimensionality::REGION, triangleID, - parametric_coords}; + triangleID, parametric_coords}; found = true; break; } @@ -589,9 +686,10 @@ Kokkos::View GridPointSearch3D::operator()( // TODO: Get nearest element if no tetrahedron found } if (!found) { - results(p) = GridPointSearch3D::Result{ - dimensionality, -1 * candidate_map.entries(nearest_triangle), - parametric_coords_to_nearest}; + LO nearest_elem = candidate_map.entries(nearest_triangle); + results(p) = + GridPointSearch3D::Result{dimensionality, -nearest_elem, nearest_elem, + parametric_coords_to_nearest}; } }); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index a6ea7813..766b4e8d 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -24,7 +24,7 @@ Kokkos::Crs construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); -} +} // namespace detail [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); @@ -50,8 +50,8 @@ construct_intersection_map_2d(Omega_h::Mesh& mesh, */ template [[nodiscard]] KOKKOS_FUNCTION inline Real -distance_to_closest_edge_from_barycentric( - const Omega_h::Matrix& coords, const Omega_h::Vector<3>& xi) +distance_to_closest_edge_from_barycentric(const Omega_h::Matrix& coords, + const Omega_h::Vector<3>& xi) { // Edge vectors originating from vertex 0 Omega_h::Vector a = coords[1] - coords[0]; @@ -64,7 +64,8 @@ distance_to_closest_edge_from_barycentric( Real parallelogram_area_sq = aa * bb - ab * ab; // Guard against degeneracy and tiny negative due to FP roundoff - if (parallelogram_area_sq <= static_cast(0)) return 0; + if (parallelogram_area_sq <= static_cast(0)) + return 0; const Real area = static_cast(0.5) * std::sqrt(parallelogram_area_sq); const Real two_area = static_cast(2) * area; @@ -88,7 +89,8 @@ distance_to_closest_edge_from_barycentric( // Return the minimum component (avoid std::min for device-compatibility) Real m = d[0]; - for (int i = 1; i < 3; ++i) m = (d[i] < m) ? d[i] : m; + for (int i = 1; i < 3; ++i) + m = (d[i] < m) ? d[i] : m; return m; } @@ -107,7 +109,8 @@ class PointLocalizationSearch }; Dimensionality dimensionality; - LO element_id; + LO element_id; // ID of the actual entity (vertex/edge/face/region) + LO face_id; // ID of the containing face, or -1 only if no cell matched Omega_h::Vector parametric_coords; }; @@ -159,6 +162,7 @@ class GridPointSearch2D : public PointLocalizationSearch2D Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; Omega_h::Adj edges2verts_adj_; + Omega_h::Adj edges2faces_up_; Kokkos::View grid_{"uniform grid"}; CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index f3d21b0c..39c218c4 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -265,28 +265,31 @@ TEST_CASE("uniform grid search") SECTION("global coordinate within mesh") { { - auto [dim, idx, coords] = results_h(0); + auto [dim, idx, face_idx, coords] = results_h(0); CAPTURE(idx); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(idx == 0); + REQUIRE(face_idx >= 0); REQUIRE(coords[0] == Catch::Approx(1)); REQUIRE(coords[1] == Catch::Approx(0)); REQUIRE(coords[2] == Catch::Approx(0)); } { - auto [dim, idx, coords] = results_h(1); + auto [dim, idx, face_idx, coords] = results_h(1); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(idx == 156); + REQUIRE(face_idx >= 0); REQUIRE(coords[0] == Catch::Approx(0.5)); REQUIRE(coords[1] == Catch::Approx(0.1)); REQUIRE(coords[2] == Catch::Approx(0.4)); } { - auto [dim, idx, coords] = results_h(7); + auto [dim, idx, face_idx, coords] = results_h(7); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(idx == 0); + REQUIRE(face_idx >= 0); } } // feature needs to be added @@ -297,21 +300,25 @@ TEST_CASE("uniform grid search") REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(-1 * out_of_bounds.element_id == top_right.element_id); + REQUIRE(out_of_bounds.face_id >= 0); out_of_bounds = results_h(4); auto bot_left = results_h(0); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); + REQUIRE(out_of_bounds.face_id >= 0); out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(out_of_bounds.element_id == -219); + REQUIRE(out_of_bounds.face_id >= 0); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); + REQUIRE(out_of_bounds.face_id >= 0); } } From 6f71531905153f506df92d8072f22a6a58fae23c Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 2 Apr 2026 09:53:40 -0400 Subject: [PATCH 4/7] Remove unused `bbox_verts_within_triangle` function and refactor `bbox_verts_within_simplex` to replace it. Eliminate `fuzz` parameter in 2D intersection logic for simplification. --- src/pcms/localization/point_search.cpp | 33 +++----------------------- 1 file changed, 3 insertions(+), 30 deletions(-) diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index fd110e2f..7af091eb 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/point_search.cpp @@ -145,36 +145,9 @@ template return true; } -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( - const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords) -{ - auto left = bbox.center[0] - bbox.half_width[0]; - auto right = bbox.center[0] + bbox.half_width[0]; - auto bot = bbox.center[1] - bbox.half_width[1]; - auto top = bbox.center[1] + bbox.half_width[1]; - auto xi = Omega_h::barycentric_from_global<2, 2>({left, bot}, coords); - if (Omega_h::is_barycentric_inside(xi)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({left, top}, coords); - if (Omega_h::is_barycentric_inside(xi)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({right, top}, coords); - if (Omega_h::is_barycentric_inside(xi)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({right, bot}, coords); - if (Omega_h::is_barycentric_inside(xi)) { - return true; - } - return false; -} - template [[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_simplex( - const AABBox& bbox, const Omega_h::Matrix& coords, - Real fuzz) + 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 } @@ -211,7 +184,7 @@ template * intersects with a bounding box */ [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( - const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox, Real fuzz) + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) { // triangle and grid cell bounding box intersect if (intersects(triangle_bbox(coords), bbox)) { @@ -221,7 +194,7 @@ template return true; } // if any of the bbox verts are within the triangle - if (bbox_verts_within_triangle(bbox, coords)) { + if (bbox_verts_within_simplex(bbox, coords)) { return true; } // if any of the triangle's edges intersect with the bounding box From a92d0a8ecea6d9b56891dd3798a9d9d485b53591 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 2 Apr 2026 09:58:03 -0400 Subject: [PATCH 5/7] Apply consistent formatting to `test_point_search.cpp` by adjusting line splits in `distance_from_line` implementation. --- test/test_point_search.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index f2cf56a1..f2e6f02a 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -121,10 +121,13 @@ TEST_CASE("barycentric distance to closest edge") auto ab = b - a; auto ap = p - a; auto ab2 = Omega_h::inner_product(ab, ab); - if (ab2 == 0.0) return Omega_h::norm(ap); + if (ab2 == 0.0) + return Omega_h::norm(ap); double t = Omega_h::inner_product(ap, ab) / ab2; - if (t < 0.0) t = 0.0; - if (t > 1.0) t = 1.0; + if (t < 0.0) + t = 0.0; + if (t > 1.0) + t = 1.0; auto proj = a + t * ab; return Omega_h::norm(p - proj); }; From 7a0dd5146739817edcf5daa1fdaf4953bb1c0123 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 2 Apr 2026 10:01:32 -0400 Subject: [PATCH 6/7] Apply minor formatting fixes in `adj_search.hpp` and `mesh_fields_adapter2.h` for consistent line splitting. --- src/pcms/adapter/meshfields/mesh_fields_adapter2.h | 3 ++- src/pcms/transfer/adj_search.hpp | 5 +++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h index c894fb3c..fbd89b6e 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h @@ -246,7 +246,8 @@ struct MeshFieldsAdapter2LocalizationHint mesh.nelems()); Kokkos::deep_copy(elem_counts, 0); for (size_t i = 0; i < num_valid_; ++i) { - auto [dim, elem_idx, face_idx, coord] = search_results(valid_point_indices[i]); + auto [dim, elem_idx, face_idx, coord] = + search_results(valid_point_indices[i]); elem_counts[face_idx] += 1; } diff --git a/src/pcms/transfer/adj_search.hpp b/src/pcms/transfer/adj_search.hpp index 0f9b15d1..6d659159 100644 --- a/src/pcms/transfer/adj_search.hpp +++ b/src/pcms/transfer/adj_search.hpp @@ -22,9 +22,10 @@ class FindSupports public: FindSupports(Omega_h::Mesh& source_mesh_, Omega_h::Mesh& target_mesh_) - : source_mesh(source_mesh_), target_mesh(target_mesh_){}; + : source_mesh(source_mesh_), target_mesh(target_mesh_) {}; - FindSupports(Omega_h::Mesh& mesh_) : source_mesh(mesh_), target_mesh(mesh_){}; + FindSupports(Omega_h::Mesh& mesh_) + : source_mesh(mesh_), target_mesh(mesh_) {}; void adjBasedSearch(Omega_h::Write& supports_ptr, Omega_h::Write& nSupports, From 7cc1f656a7236cefd8813a2908985bfea0a26f33 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 15 Apr 2026 23:44:43 -0400 Subject: [PATCH 7/7] Refactor point search logic to use owning element IDs instead of face IDs for better clarity and consistency in results handling --- .../adapter/meshfields/mesh_fields_adapter.h | 45 +++- .../adapter/meshfields/mesh_fields_adapter2.h | 82 ++++++-- src/pcms/localization/point_search.cpp | 196 ++++++++++++++++-- src/pcms/localization/point_search.h | 20 +- src/pcms/transfer/adj_search.cpp | 17 +- src/pcms/transfer/mesh_intersection.cpp | 3 +- test/test_point_search.cpp | 17 +- 7 files changed, 319 insertions(+), 61 deletions(-) diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter.h b/src/pcms/adapter/meshfields/mesh_fields_adapter.h index d68fd56e..5c7e71ef 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter.h @@ -166,6 +166,15 @@ class MeshFieldsAdapter "search data structure must be constructed before use"); return (*search_)(points); } + [[nodiscard]] Kokkos::View GetOwningElementIds( + Kokkos::View results) + const + { + PCMS_FUNCTION_TIMER; + PCMS_ALWAYS_ASSERT(search_ != nullptr && + "search data structure must be constructed before use"); + return search_->GetOwningElementIds(results); + } [[nodiscard]] Omega_h::Read GetClassIDs() const { @@ -335,6 +344,7 @@ auto evaluate(const MeshFieldsAdapter& field, Lagrange<1> /* method */, PCMS_FUNCTION_TIMER; Omega_h::Write values(coordinates.size() / 2); auto tris2verts = field.GetMesh().ask_elem_verts(); + auto mesh_coords = field.GetMesh().coords(); auto field_values = field.GetMesh().template get_array(0, field.GetName()); Kokkos::View coords("coords", coordinates.size() / 2); @@ -344,17 +354,27 @@ auto evaluate(const MeshFieldsAdapter& field, Lagrange<1> /* method */, coords(i, 1) = coordinates(2 * i + 1); }); auto results = field.Search(coords); + auto owning_ids = field.GetOwningElementIds(results); Kokkos::parallel_for( results.size(), KOKKOS_LAMBDA(LO i) { - auto [dim, elem_idx, face_idx, coord] = results(i); - // TODO deal with case for face_idx < 0 (point outside of mesh) + auto [dim, elem_idx, coord] = results(i); + (void)dim; + (void)elem_idx; + (void)coord; + auto face_idx = owning_ids(i); + // TODO deal with case for elem_idx < 0 (point outside of mesh) KOKKOS_ASSERT(face_idx >= 0); const auto elem_tri2verts = Omega_h::gather_verts<3>(tris2verts, face_idx); + const auto vertex_coords = + Omega_h::gather_vectors<3, 2>(mesh_coords, elem_tri2verts); + const Omega_h::Vector<2> point{coords(i, 0), coords(i, 1)}; + const auto local = + Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); Real val = 0; for (int j = 0; j < 3; ++j) { - val += field_values[elem_tri2verts[j]] * coord[j]; + val += field_values[elem_tri2verts[j]] * local[j]; } if constexpr (std::is_integral_v) { val = std::round(val); @@ -373,6 +393,7 @@ auto evaluate(const MeshFieldsAdapter& field, NearestNeighbor /* method */, PCMS_FUNCTION_TIMER; Omega_h::Write values(coordinates.size() / 2); auto tris2verts = field.GetMesh().ask_elem_verts(); + auto mesh_coords = field.GetMesh().coords(); auto field_values = field.GetMesh().template get_array(0, field.GetName()); // TODO reuse coordinates_data if possible Kokkos::View coords("coords", coordinates.size() / 2); @@ -382,19 +403,29 @@ auto evaluate(const MeshFieldsAdapter& field, NearestNeighbor /* method */, coords(i, 1) = coordinates(2 * i + 1); }); auto results = field.Search(coords); + auto owning_ids = field.GetOwningElementIds(results); Kokkos::parallel_for( results.size(), KOKKOS_LAMBDA(LO i) { - auto [dim, elem_idx, face_idx, coord] = results(i); - // TODO deal with case for face_idx < 0 (point outside of mesh) + auto [dim, elem_idx, coord] = results(i); + (void)dim; + (void)elem_idx; + (void)coord; + auto face_idx = owning_ids(i); + // TODO deal with case for elem_idx < 0 (point outside of mesh) KOKKOS_ASSERT(face_idx >= 0); const auto elem_tri2verts = Omega_h::gather_verts<3>(tris2verts, face_idx); + const auto vertex_coords = + Omega_h::gather_vectors<3, 2>(mesh_coords, elem_tri2verts); + const Omega_h::Vector<2> point{coords(i, 0), coords(i, 1)}; + const auto local = + Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); // value is closest to point has the largest coordinate int vert = 0; - auto max_val = coord[0]; + auto max_val = local[0]; for (int j = 1; j <= 2; ++j) { - auto next_val = coord[j]; + auto next_val = local[j]; if (next_val > max_val) { max_val = next_val; vert = j; diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h index fbd89b6e..34dd8079 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h @@ -122,19 +122,27 @@ struct CountPointsPerElementFunctor { Kokkos::View elem_counts_; Kokkos::View search_results_; + Kokkos::View owning_elem_ids_; CountPointsPerElementFunctor( Kokkos::View elem_counts, - Kokkos::View search_results) - : elem_counts_(elem_counts), search_results_(search_results) + Kokkos::View search_results, + Kokkos::View owning_elem_ids) + : elem_counts_(elem_counts), + search_results_(search_results), + owning_elem_ids_(owning_elem_ids) { } KOKKOS_INLINE_FUNCTION void operator()(LO i) const { - auto [dim, elem_idx, face_idx, coord] = search_results_(i); - Kokkos::atomic_add(&elem_counts_(face_idx), 1); + auto [dim, elem_idx, coord] = search_results_(i); + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids_(i); + Kokkos::atomic_add(&elem_counts_(owner_idx), 1); } }; @@ -146,19 +154,22 @@ struct FillCoordinatesAndIndicesFunctor Kokkos::View coordinates_; Kokkos::View indices_; Kokkos::View search_results_; + Kokkos::View owning_elem_ids_; Omega_h::Int dim_; FillCoordinatesAndIndicesFunctor( Omega_h::Mesh& mesh, Kokkos::View elem_counts, Kokkos::View offsets, Kokkos::View coordinates, Kokkos::View indices, - Kokkos::View search_results) + Kokkos::View search_results, + Kokkos::View owning_elem_ids) : mesh_(mesh), elem_counts_(elem_counts), offsets_(offsets), coordinates_(coordinates), indices_(indices), search_results_(search_results), + owning_elem_ids_(owning_elem_ids), dim_(mesh.dim()) { } @@ -166,14 +177,15 @@ struct FillCoordinatesAndIndicesFunctor KOKKOS_INLINE_FUNCTION void operator()(LO i) const { - auto [dim, elem_idx, face_idx, coord] = search_results_(i); + auto [dim, elem_idx, coord] = search_results_(i); + const auto owner_idx = owning_elem_ids_(i); // disable the host assertion macro for device code // currently don't handle case where point is on a boundary // PCMS_ALWAYS_ASSERT(static_cast(dim) == mesh_.dim()); // face should be inside the domain (positive) // PCMS_ALWAYS_ASSERT(face_idx >= 0 && face_idx < mesh_.nelems()); - LO count = Kokkos::atomic_sub_fetch(&elem_counts_(face_idx), 1); - LO index = offsets_(face_idx) + count - 1; + LO count = Kokkos::atomic_sub_fetch(&elem_counts_(owner_idx), 1); + LO index = offsets_(owner_idx) + count - 1; for (int j = 0; j < (dim_ + 1); ++j) { coordinates_(index, j) = coord[j]; } @@ -186,7 +198,8 @@ struct MeshFieldsAdapter2LocalizationHint MeshFieldsAdapter2LocalizationHint( Omega_h::Mesh& mesh, Kokkos::View search_results, - OutOfBoundsMode mode) + Kokkos::View global_coords, + Kokkos::View owning_elem_ids, OutOfBoundsMode mode) : mode_(mode), num_valid_(0), num_missing_(0) { // First pass: count valid and invalid points @@ -196,20 +209,22 @@ struct MeshFieldsAdapter2LocalizationHint if (mode_ == OutOfBoundsMode::ERROR) { // Error mode - throw error immediately if any point is out of bounds for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, face_idx, coord] = search_results(i); + auto [dim, elem_idx, coord] = search_results(i); + const auto owner_idx = owning_elem_ids(i); (void)dim; (void)coord; - bool is_missing = (elem_idx < 0) || (face_idx < 0); + bool is_missing = (elem_idx < 0) || (owner_idx < 0); PCMS_ALWAYS_ASSERT(!is_missing && "Points found outside mesh domain"); valid_point_indices.push_back(i); } } else { // Other modes - collect valid and missing points separately for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, face_idx, coord] = search_results(i); + auto [dim, elem_idx, coord] = search_results(i); + const auto owner_idx = owning_elem_ids(i); (void)dim; (void)coord; - bool is_missing = (elem_idx < 0) || (face_idx < 0); + bool is_missing = (elem_idx < 0) || (owner_idx < 0); if (is_missing) { missing_point_indices.push_back(i); } else { @@ -246,9 +261,12 @@ struct MeshFieldsAdapter2LocalizationHint mesh.nelems()); Kokkos::deep_copy(elem_counts, 0); for (size_t i = 0; i < num_valid_; ++i) { - auto [dim, elem_idx, face_idx, coord] = - search_results(valid_point_indices[i]); - elem_counts[face_idx] += 1; + auto [dim, elem_idx, coord] = search_results(valid_point_indices[i]); + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids(valid_point_indices[i]); + elem_counts[owner_idx] += 1; } // Compute offsets @@ -261,14 +279,27 @@ struct MeshFieldsAdapter2LocalizationHint offsets_(mesh.nelems()) = total; // Fill coordinates and indices for valid points + const auto tris2verts = mesh.ask_elem_verts(); + const auto mesh_coords = mesh.coords(); for (size_t i = 0; i < num_valid_; ++i) { size_t orig_idx = valid_point_indices[i]; - auto [dim, elem_idx, face_idx, coord] = search_results(orig_idx); - elem_counts(face_idx) -= 1; - LO index = offsets_(face_idx) + elem_counts(face_idx); - for (int j = 0; j < (mesh.dim() + 1); ++j) { - coordinates_(index, j) = coord[j]; - } + auto [dim, elem_idx, coord] = search_results(orig_idx); + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids(orig_idx); + elem_counts(owner_idx) -= 1; + LO index = offsets_(owner_idx) + elem_counts(owner_idx); + const auto elem_tri2verts = + Omega_h::gather_verts<3>(tris2verts, owner_idx); + const auto vertex_coords = + Omega_h::gather_vectors<3, 2>(mesh_coords, elem_tri2verts); + const Omega_h::Vector<2> point{global_coords(orig_idx, 0), + global_coords(orig_idx, 1)}; + const auto local = + Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); + for (int j = 0; j < (mesh.dim() + 1); ++j) + coordinates_(index, j) = local[j]; indices_(index) = static_cast(orig_idx); } } @@ -432,11 +463,16 @@ inline LocalizationHint MeshFieldsAdapter2::GetLocalizationHint( coordinates.data_handle(), coordinates.extent(0), coordinates.extent(1)); deep_copy_mismatch_layouts(coords, coordinates_host); auto results = search_(coords); + auto owning_ids = search_.GetOwningElementIds(results); Kokkos::View results_h( "results_h", results.size()); + Kokkos::View owning_ids_h("owning_ids_h", + owning_ids.size()); Kokkos::deep_copy(results_h, results); + Kokkos::deep_copy(owning_ids_h, owning_ids); auto hint = std::make_shared( - mesh_, results_h, this->out_of_bounds_mode_); + mesh_, results_h, coordinates_host, owning_ids_h, + this->out_of_bounds_mode_); return LocalizationHint{hint}; } diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index 7af091eb..148b927b 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/point_search.cpp @@ -513,7 +513,6 @@ Kokkos::View GridPointSearch2D::operator()( GridPointSearch2D::Result::Dimensionality dim_out = GridPointSearch2D::Result::Dimensionality::REGION; LO element_id_out = -1; - LO face_id_out = -1; Omega_h::Vector<3> bary_out{0.0, 0.0, 0.0}; // Apply hierarchy with tolerances @@ -523,7 +522,6 @@ Kokkos::View GridPointSearch2D::operator()( // Point within vertex tolerance - still inside the mesh dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; element_id_out = best_vertex_id; - face_id_out = (best_vertex_tid >= 0) ? best_vertex_tid : -1; bary_out = best_vertex_bary; } else if ((best_edge_id >= 0 && best_edge_dist <= etol) || (found_inside && inside_edge_dist <= etol)) { @@ -532,18 +530,15 @@ Kokkos::View GridPointSearch2D::operator()( if (found_inside && inside_edge_dist <= etol && inside_edge_argmin >= 0) { element_id_out = inside_edge_id; - face_id_out = inside_face_id; bary_out = inside_face_bary; } else { element_id_out = best_edge_id; - face_id_out = (best_edge_tid_min >= 0) ? best_edge_tid_min : -1; bary_out = best_edge_bary_min; } } else if (found_inside) { // Point inside face dim_out = GridPointSearch2D::Result::Dimensionality::FACE; element_id_out = inside_face_id; - face_id_out = inside_face_id; bary_out = inside_face_bary; } else { // Outside mesh - beyond tolerance of any entity @@ -552,13 +547,11 @@ Kokkos::View GridPointSearch2D::operator()( (best_vertex_dist <= best_edge_dist || best_edge_id < 0)) { dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; element_id_out = best_vertex_id; - face_id_out = (best_vertex_tid >= 0) ? best_vertex_tid : -1; bary_out = best_vertex_bary; element_id_out = -element_id_out; } else if (best_edge_id >= 0) { dim_out = GridPointSearch2D::Result::Dimensionality::EDGE; element_id_out = best_edge_id; - face_id_out = (best_edge_tid_min >= 0) ? best_edge_tid_min : -1; bary_out = best_edge_bary_min; element_id_out = -element_id_out; } else { @@ -566,8 +559,7 @@ Kokkos::View GridPointSearch2D::operator()( } } - results(p) = GridPointSearch2D::Result{dim_out, element_id_out, - face_id_out, bary_out}; + results(p) = GridPointSearch2D::Result{dim_out, element_id_out, bary_out}; }); return results; @@ -603,6 +595,88 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); edges2faces_up_ = mesh.ask_up(Omega_h::EDGE, Omega_h::FACE); + verts2faces_up_ = mesh.ask_up(Omega_h::VERT, Omega_h::FACE); +} + +LO GridPointSearch2D::get_smallest_owner_face_id( + Result::Dimensionality dimensionality, LO element_id) const +{ + if (element_id < 0) + return -1; + if (dimensionality == Result::Dimensionality::FACE) + return element_id; + + const Omega_h::Adj* up_adj = nullptr; + if (dimensionality == Result::Dimensionality::EDGE) { + up_adj = &edges2faces_up_; + } else if (dimensionality == Result::Dimensionality::VERTEX) { + up_adj = &verts2faces_up_; + } else { + return -1; + } + + const auto begin = up_adj->a2ab[element_id]; + const auto end = up_adj->a2ab[element_id + 1]; + if (begin >= end) + return -1; + LO owner = up_adj->ab2b[begin]; + for (auto i = begin + 1; i < end; ++i) { + const LO candidate = up_adj->ab2b[i]; + if (candidate < owner) + owner = candidate; + } + return owner; +} + +LO GridPointSearch2D::GetOwningElementId(const Result& result) const +{ + const LO query_id = + (result.element_id < 0) ? -result.element_id : result.element_id; + return get_smallest_owner_face_id(result.dimensionality, query_id); +} + +Kokkos::View GridPointSearch2D::GetOwningElementIds( + Kokkos::View results) const +{ + Kokkos::View owners("point search owning face ids", results.extent(0)); + auto edges2faces_up = edges2faces_up_; + auto verts2faces_up = verts2faces_up_; + Kokkos::parallel_for( + results.extent(0), KOKKOS_LAMBDA(const LO i) { + const auto result = results(i); + LO element_id = result.element_id; + if (element_id < 0) + element_id = -element_id; + + LO owner = -1; + if (element_id >= 0) { + if (result.dimensionality == Result::Dimensionality::FACE) { + owner = element_id; + } else if (result.dimensionality == Result::Dimensionality::EDGE) { + const LO begin = edges2faces_up.a2ab[element_id]; + const LO end = edges2faces_up.a2ab[element_id + 1]; + if (begin < end) { + owner = edges2faces_up.ab2b[begin]; + for (LO j = begin + 1; j < end; ++j) { + const LO candidate = edges2faces_up.ab2b[j]; + owner = (candidate < owner) ? candidate : owner; + } + } + } else if (result.dimensionality == Result::Dimensionality::VERTEX) { + const LO begin = verts2faces_up.a2ab[element_id]; + const LO end = verts2faces_up.a2ab[element_id + 1]; + if (begin < end) { + owner = verts2faces_up.ab2b[begin]; + for (LO j = begin + 1; j < end; ++j) { + const LO candidate = verts2faces_up.ab2b[j]; + owner = (candidate < owner) ? candidate : owner; + } + } + } + } + owners(i) = owner; + }); + return owners; } Kokkos::View GridPointSearch3D::operator()( @@ -649,7 +723,7 @@ Kokkos::View GridPointSearch3D::operator()( if (Omega_h::is_barycentric_inside(parametric_coords)) { results(p) = GridPointSearch3D::Result{ GridPointSearch3D::Result::Dimensionality::REGION, triangleID, - triangleID, parametric_coords}; + parametric_coords}; found = true; break; } @@ -658,9 +732,8 @@ Kokkos::View GridPointSearch3D::operator()( } if (!found) { LO nearest_elem = candidate_map.entries(nearest_triangle); - results(p) = - GridPointSearch3D::Result{dimensionality, -nearest_elem, nearest_elem, - parametric_coords_to_nearest}; + results(p) = GridPointSearch3D::Result{dimensionality, -nearest_elem, + parametric_coords_to_nearest}; } }); @@ -701,5 +774,102 @@ GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, tris2edges_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::EDGE); tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); + verts2regions_up_ = mesh.ask_up(Omega_h::VERT, Omega_h::REGION); + edges2regions_up_ = mesh.ask_up(Omega_h::EDGE, Omega_h::REGION); + faces2regions_up_ = mesh.ask_up(Omega_h::FACE, Omega_h::REGION); +} + +LO GridPointSearch3D::get_smallest_owner_region_id( + Result::Dimensionality dimensionality, LO element_id) const +{ + if (element_id < 0) + return -1; + if (dimensionality == Result::Dimensionality::REGION) + return element_id; + + const Omega_h::Adj* up_adj = nullptr; + if (dimensionality == Result::Dimensionality::FACE) { + up_adj = &faces2regions_up_; + } else if (dimensionality == Result::Dimensionality::EDGE) { + up_adj = &edges2regions_up_; + } else if (dimensionality == Result::Dimensionality::VERTEX) { + up_adj = &verts2regions_up_; + } else { + return -1; + } + + const auto begin = up_adj->a2ab[element_id]; + const auto end = up_adj->a2ab[element_id + 1]; + if (begin >= end) + return -1; + LO owner = up_adj->ab2b[begin]; + for (auto i = begin + 1; i < end; ++i) { + const LO candidate = up_adj->ab2b[i]; + if (candidate < owner) + owner = candidate; + } + return owner; +} + +LO GridPointSearch3D::GetOwningElementId(const Result& result) const +{ + const LO query_id = + (result.element_id < 0) ? -result.element_id : result.element_id; + return get_smallest_owner_region_id(result.dimensionality, query_id); +} + +Kokkos::View GridPointSearch3D::GetOwningElementIds( + Kokkos::View results) const +{ + Kokkos::View owners("point search owning region ids", results.extent(0)); + auto verts2regions_up = verts2regions_up_; + auto edges2regions_up = edges2regions_up_; + auto faces2regions_up = faces2regions_up_; + Kokkos::parallel_for( + results.extent(0), KOKKOS_LAMBDA(const LO i) { + const auto result = results(i); + LO element_id = result.element_id; + if (element_id < 0) + element_id = -element_id; + + LO owner = -1; + if (element_id >= 0) { + if (result.dimensionality == Result::Dimensionality::REGION) { + owner = element_id; + } else if (result.dimensionality == Result::Dimensionality::FACE) { + const LO begin = faces2regions_up.a2ab[element_id]; + const LO end = faces2regions_up.a2ab[element_id + 1]; + if (begin < end) { + owner = faces2regions_up.ab2b[begin]; + for (LO j = begin + 1; j < end; ++j) { + const LO candidate = faces2regions_up.ab2b[j]; + owner = (candidate < owner) ? candidate : owner; + } + } + } else if (result.dimensionality == Result::Dimensionality::EDGE) { + const LO begin = edges2regions_up.a2ab[element_id]; + const LO end = edges2regions_up.a2ab[element_id + 1]; + if (begin < end) { + owner = edges2regions_up.ab2b[begin]; + for (LO j = begin + 1; j < end; ++j) { + const LO candidate = edges2regions_up.ab2b[j]; + owner = (candidate < owner) ? candidate : owner; + } + } + } else if (result.dimensionality == Result::Dimensionality::VERTEX) { + const LO begin = verts2regions_up.a2ab[element_id]; + const LO end = verts2regions_up.a2ab[element_id + 1]; + if (begin < end) { + owner = verts2regions_up.ab2b[begin]; + for (LO j = begin + 1; j < end; ++j) { + const LO candidate = verts2regions_up.ab2b[j]; + owner = (candidate < owner) ? candidate : owner; + } + } + } + } + owners(i) = owner; + }); + return owners; } } // namespace pcms diff --git a/src/pcms/localization/point_search.h b/src/pcms/localization/point_search.h index 752def22..ce9fa2b9 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -110,7 +110,6 @@ class PointLocalizationSearch Dimensionality dimensionality; LO element_id; // ID of the actual entity (vertex/edge/face/region) - LO face_id; // ID of the containing face, or -1 only if no cell matched Omega_h::Vector parametric_coords; }; @@ -126,6 +125,9 @@ class PointLocalizationSearch virtual Kokkos::View operator()( Kokkos::View point) const = 0; + [[nodiscard]] virtual LO GetOwningElementId(const Result& result) const = 0; + [[nodiscard]] virtual Kokkos::View GetOwningElementIds( + Kokkos::View results) const = 0; virtual ~PointLocalizationSearch() = default; protected: @@ -156,13 +158,20 @@ class GridPointSearch2D : public PointLocalizationSearch2D */ Kokkos::View operator()( Kokkos::View point) const override; + [[nodiscard]] LO GetOwningElementId(const Result& result) const override; + [[nodiscard]] Kokkos::View GetOwningElementIds( + Kokkos::View results) const override; private: + [[nodiscard]] LO get_smallest_owner_face_id( + Result::Dimensionality dimensionality, LO element_id) const; + Omega_h::Mesh mesh_; Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; Omega_h::Adj edges2verts_adj_; Omega_h::Adj edges2faces_up_; + Omega_h::Adj verts2faces_up_; Kokkos::View grid_{"uniform grid"}; CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; @@ -190,12 +199,21 @@ class GridPointSearch3D : public PointLocalizationSearch3D */ Kokkos::View operator()( Kokkos::View point) const override; + [[nodiscard]] LO GetOwningElementId(const Result& result) const override; + [[nodiscard]] Kokkos::View GetOwningElementIds( + Kokkos::View results) const override; private: + [[nodiscard]] LO get_smallest_owner_region_id( + Result::Dimensionality dimensionality, LO element_id) const; + Omega_h::Mesh mesh_; Omega_h::Adj tris2edges_adj_; Omega_h::Adj tris2verts_adj_; Omega_h::Adj edges2verts_adj_; + Omega_h::Adj verts2regions_up_; + Omega_h::Adj edges2regions_up_; + Omega_h::Adj faces2regions_up_; Kokkos::View[1]> grid_{"uniform grid"}; CandidateMapT candidate_map_; Omega_h::LOs tris2verts_; diff --git a/src/pcms/transfer/adj_search.cpp b/src/pcms/transfer/adj_search.cpp index 7060278a..e1ccfcdc 100644 --- a/src/pcms/transfer/adj_search.cpp +++ b/src/pcms/transfer/adj_search.cpp @@ -7,14 +7,16 @@ namespace pcms // Debug helper retained intentionally: useful for diagnosing point-localization // failures while developing support-search logic. [[maybe_unused]] static void checkTargetPoints( - const Kokkos::View& results) + const Kokkos::View& results, + const Kokkos::View& owning_cell_ids) { Kokkos::fence(); pcms::printInfo("INFO: Checking target points...\n"); auto check_target_points = OMEGA_H_LAMBDA(Omega_h::LO i) { - if (results(i).face_id < 0) { - OMEGA_H_CHECK_PRINTF(results(i).face_id >= 0, + (void)results; + if (owning_cell_ids(i) < 0) { + OMEGA_H_CHECK_PRINTF(owning_cell_ids(i) >= 0, "ERROR: Source face id not found for target %d\n", i); printf("%d, ", i); @@ -147,7 +149,8 @@ void FindSupports::adjBasedSearch(Omega_h::Write& supports_ptr, pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); auto results = search_cell(target_points); - // checkTargetPoints(results); + auto owning_cell_ids = search_cell.GetOwningElementIds(results); + // checkTargetPoints(results, owning_cell_ids); Omega_h::parallel_for( nvertices_target, @@ -155,11 +158,7 @@ void FindSupports::adjBasedSearch(Omega_h::Write& supports_ptr, Queue queue; Track visited; Omega_h::Real cutoffDistance = radii2[id]; - Omega_h::LO source_cell_id = results(id).face_id; - // if the point lies outside the mesh, it returns the negative cell id - // making the negative cell id to positive - if (source_cell_id < 0) - source_cell_id = Kokkos::abs(source_cell_id); + Omega_h::LO source_cell_id = owning_cell_ids(id); OMEGA_H_CHECK_PRINTF( source_cell_id >= 0, diff --git a/src/pcms/transfer/mesh_intersection.cpp b/src/pcms/transfer/mesh_intersection.cpp index 85e41967..0bad5cf5 100644 --- a/src/pcms/transfer/mesh_intersection.cpp +++ b/src/pcms/transfer/mesh_intersection.cpp @@ -29,6 +29,7 @@ void FindIntersections::adjBasedIntersectSearch( pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20); auto results = search_cell(centroids); + auto owning_cell_ids = search_cell.GetOwningElementIds(results); auto nfaces_target = target_mesh_.nfaces(); Omega_h::parallel_for( @@ -37,7 +38,7 @@ void FindIntersections::adjBasedIntersectSearch( Queue queue; Track visited; - auto current_cell_id = results(id).element_id; + auto current_cell_id = owning_cell_ids(id); auto current_tgt_elm_area = tgt_elem_areas[id]; OMEGA_H_CHECK_PRINTF(current_cell_id >= 0, diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index f2e6f02a..bfaaff78 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -268,7 +268,8 @@ TEST_CASE("uniform grid search") SECTION("global coordinate within mesh") { { - auto [dim, idx, face_idx, coords] = results_h(0); + auto [dim, idx, coords] = results_h(0); + const auto face_idx = search.GetOwningElementId(results_h(0)); CAPTURE(idx); @@ -280,7 +281,8 @@ TEST_CASE("uniform grid search") REQUIRE(coords[2] == Catch::Approx(0)); } { - auto [dim, idx, face_idx, coords] = results_h(1); + auto [dim, idx, coords] = results_h(1); + const auto face_idx = search.GetOwningElementId(results_h(1)); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(idx == 156); REQUIRE(face_idx >= 0); @@ -289,7 +291,8 @@ TEST_CASE("uniform grid search") REQUIRE(coords[2] == Catch::Approx(0.4)); } { - auto [dim, idx, face_idx, coords] = results_h(7); + auto [dim, idx, coords] = results_h(7); + const auto face_idx = search.GetOwningElementId(results_h(7)); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(idx == 0); REQUIRE(face_idx >= 0); @@ -303,25 +306,25 @@ TEST_CASE("uniform grid search") REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(-1 * out_of_bounds.element_id == top_right.element_id); - REQUIRE(out_of_bounds.face_id >= 0); + REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); out_of_bounds = results_h(4); auto bot_left = results_h(0); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); - REQUIRE(out_of_bounds.face_id >= 0); + REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(out_of_bounds.element_id == -219); - REQUIRE(out_of_bounds.face_id >= 0); + REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); - REQUIRE(out_of_bounds.face_id >= 0); + REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); } }