From aaf9e25e7bb2e7eeeb0b23e440c2939e3940e1cd Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 15 May 2026 01:32:24 -0400 Subject: [PATCH 1/3] remove fuzz and return face id alongside element id --- src/pcms/field/evaluator/mesh_fields.h | 4 +- .../field/evaluator/mesh_fields_backend.h | 157 ++++-- src/pcms/field/evaluator/omega_h_lagrange.h | 63 ++- src/pcms/localization/adj_search.cpp | 24 +- src/pcms/localization/point_search.cpp | 520 ++++++++++++++---- src/pcms/localization/point_search.h | 100 +++- src/pcms/transfer/mesh_intersection.cpp | 3 +- test/test_point_search.cpp | 74 +++ 8 files changed, 750 insertions(+), 195 deletions(-) diff --git a/src/pcms/field/evaluator/mesh_fields.h b/src/pcms/field/evaluator/mesh_fields.h index cbe69f5b..54dc5c35 100644 --- a/src/pcms/field/evaluator/mesh_fields.h +++ b/src/pcms/field/evaluator/mesh_fields.h @@ -140,8 +140,10 @@ class MeshFieldsEvaluatorFactory : public FieldEvaluatorFactory coords_search(i, 1) = coordinates(i, 1); }); auto results = search_(coords_search); + auto owning_ids = search_.GetOwningElementIds(results); - MeshFieldsAdapter2LocalizationHint hint(mesh_, results, policy.mode); + MeshFieldsAdapter2LocalizationHint hint(mesh_, results, coords_search, + owning_ids, policy.mode); return std::make_unique>( layout_, std::move(hint), policy.fill_value); diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index fdce6d6c..a895381e 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/src/pcms/field/evaluator/mesh_fields_backend.h @@ -183,11 +183,15 @@ 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) { } @@ -195,7 +199,11 @@ struct CountPointsPerElementFunctor void operator()(LO i) const { auto [dim, elem_idx, coord] = search_results_(i); - Kokkos::atomic_add(&elem_counts_(elem_idx), 1); + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids_(i); + Kokkos::atomic_add(&elem_counts_(owner_idx), 1); } }; @@ -207,19 +215,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()) { } @@ -228,8 +239,9 @@ struct FillCoordinatesAndIndicesFunctor void operator()(LO i) const { auto [dim, elem_idx, coord] = search_results_(i); - LO count = Kokkos::atomic_sub_fetch(&elem_counts_(elem_idx), 1); - LO index = offsets_(elem_idx) + count - 1; + const auto owner_idx = owning_elem_ids_(i); + 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]; } @@ -242,14 +254,17 @@ struct FilterValidPointsFunctor { Kokkos::View search_results_; Kokkos::View valid_flags_; + Kokkos::View owning_elem_ids_; Omega_h::Int mesh_dim_; OutOfBoundsMode mode_; FilterValidPointsFunctor(Kokkos::View results, - Kokkos::View flags, Omega_h::Int mesh_dim, - OutOfBoundsMode mode) + Kokkos::View flags, + Kokkos::View owning_elem_ids, + Omega_h::Int mesh_dim, OutOfBoundsMode mode) : search_results_(results), valid_flags_(flags), + owning_elem_ids_(owning_elem_ids), mesh_dim_(mesh_dim), mode_(mode) { @@ -259,14 +274,17 @@ struct FilterValidPointsFunctor void operator()(LO i, LO& num_valid, LO& num_missing) const { auto [dim, elem_idx, coord] = search_results_(i); - bool is_valid = (static_cast(dim) == mesh_dim_) && (elem_idx >= 0); + const auto owner_idx = owning_elem_ids_(i); + (void)dim; + (void)coord; + bool is_missing = (elem_idx < 0) || (owner_idx < 0); - if (mode_ == OutOfBoundsMode::ERROR && !is_valid) { + if (mode_ == OutOfBoundsMode::ERROR && is_missing) { Kokkos::abort("Points found outside mesh domain"); } - valid_flags_(i) = is_valid ? 1 : 0; - if (is_valid) { + valid_flags_(i) = !is_missing ? 1 : 0; + if (!is_missing) { num_valid++; } else { num_missing++; @@ -311,13 +329,16 @@ struct CountPerElementFunctor { Kokkos::View search_results_; Kokkos::View valid_indices_; + Kokkos::View owning_elem_ids_; Kokkos::View elem_counts_; CountPerElementFunctor(Kokkos::View results, Kokkos::View valid_indices, + Kokkos::View owning_elem_ids, Kokkos::View elem_counts) : search_results_(results), valid_indices_(valid_indices), + owning_elem_ids_(owning_elem_ids), elem_counts_(elem_counts) { } @@ -327,7 +348,11 @@ struct CountPerElementFunctor { LO orig_idx = valid_indices_(i); auto [dim, elem_idx, coord] = search_results_(orig_idx); - Kokkos::atomic_add(&elem_counts_(elem_idx), 1); + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids_(orig_idx); + Kokkos::atomic_add(&elem_counts_(owner_idx), 1); } }; @@ -335,24 +360,36 @@ struct FillCoordinatesDeviceFunctor { Kokkos::View search_results_; Kokkos::View valid_indices_; + Kokkos::View owning_elem_ids_; Kokkos::View elem_counts_; Kokkos::View offsets_; Kokkos::View coordinates_; Kokkos::View indices_; + Omega_h::LOs tris2verts_; + Omega_h::Reals mesh_coords_; + Kokkos::View global_coords_; + Omega_h::Mesh& mesh_; Omega_h::Int dim_; - FillCoordinatesDeviceFunctor(Kokkos::View results, - Kokkos::View valid_indices, - Kokkos::View elem_counts, - Kokkos::View offsets, - Kokkos::View coordinates, - Kokkos::View indices, Omega_h::Int dim) + FillCoordinatesDeviceFunctor( + Kokkos::View results, + Kokkos::View valid_indices, Kokkos::View owning_elem_ids, + Kokkos::View elem_counts, Kokkos::View offsets, + Kokkos::View coordinates, Kokkos::View indices, + Omega_h::LOs tris2verts, Omega_h::Reals mesh_coords, + Kokkos::View global_coords, + Omega_h::Mesh& mesh, Omega_h::Int dim) : search_results_(results), valid_indices_(valid_indices), + owning_elem_ids_(owning_elem_ids), elem_counts_(elem_counts), offsets_(offsets), coordinates_(coordinates), indices_(indices), + tris2verts_(tris2verts), + mesh_coords_(mesh_coords), + global_coords_(global_coords), + mesh_(mesh), dim_(dim) { } @@ -362,12 +399,23 @@ struct FillCoordinatesDeviceFunctor { LO orig_idx = valid_indices_(i); auto [dim, elem_idx, coord] = search_results_(orig_idx); - LO count = Kokkos::atomic_fetch_sub(&elem_counts_(elem_idx), 1); - LO index = offsets_(elem_idx) + count - 1; - - for (int j = 0; j < (dim_ + 1); ++j) { - coordinates_(index, j) = coord[j]; - } + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids_(orig_idx); + LO count = + Kokkos::atomic_fetch_add(&elem_counts_(owner_idx), static_cast(-1)); + LO index = offsets_(owner_idx) + count - 1; + 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 < (dim_ + 1); ++j) + coordinates_(index, j) = local[j]; indices_(index) = orig_idx; } }; @@ -416,7 +464,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) { std::vector valid_point_indices; @@ -425,16 +474,20 @@ struct MeshFieldsAdapter2LocalizationHint if (mode_ == OutOfBoundsMode::ERROR) { 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); + const auto owner_idx = owning_elem_ids(i); + (void)dim; + (void)coord; + 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 { 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); + const auto owner_idx = owning_elem_ids(i); + (void)dim; + (void)coord; + bool is_missing = (elem_idx < 0) || (owner_idx < 0); if (is_missing) { missing_point_indices.push_back(i); } else { @@ -465,9 +518,14 @@ struct MeshFieldsAdapter2LocalizationHint 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; + (void)dim; + (void)elem_idx; + (void)coord; + const auto owner_idx = owning_elem_ids(valid_point_indices[i]); + elem_counts[owner_idx] += 1; } LO total; @@ -478,14 +536,27 @@ struct MeshFieldsAdapter2LocalizationHint functor, total); offsets_(mesh.nelems()) = total; + 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, coord] = search_results(orig_idx); - elem_counts(elem_idx) -= 1; - LO index = offsets_(elem_idx) + elem_counts(elem_idx); - for (int j = 0; j < (mesh.dim() + 1); ++j) { - coordinates_(index, j) = coord[j]; - } + (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); } @@ -511,7 +582,8 @@ struct MeshFieldsAdapter2LocalizationHint MeshFieldsAdapter2LocalizationHint( Omega_h::Mesh& mesh, Kokkos::View search_results_d, - OutOfBoundsMode mode) + Kokkos::View global_coords, + Kokkos::View owning_elem_ids, OutOfBoundsMode mode) : mode_(mode), num_valid_(0), num_missing_(0) { const LO n_points = search_results_d.size(); @@ -523,7 +595,7 @@ struct MeshFieldsAdapter2LocalizationHint // Step 1: Filter valid/missing points on device Kokkos::View valid_flags("valid_flags", n_points); FilterValidPointsFunctor filter_functor(search_results_d, valid_flags, - mesh.dim(), mode_); + owning_elem_ids, mesh.dim(), mode_); LO num_valid = 0, num_missing = 0; Kokkos::parallel_reduce( "FilterValidPoints", @@ -555,7 +627,7 @@ struct MeshFieldsAdapter2LocalizationHint Kokkos::deep_copy(elem_counts_d, 0); CountPerElementFunctor count_functor(search_results_d, valid_indices_d, - elem_counts_d); + owning_elem_ids, elem_counts_d); Kokkos::parallel_for( "CountPerElement", Kokkos::RangePolicy(0, num_valid_), @@ -580,10 +652,13 @@ struct MeshFieldsAdapter2LocalizationHint coordinates_d_ = Kokkos::View("coordinates_d", num_valid_, mesh.dim() + 1); indices_d_ = Kokkos::View("indices_d", num_valid_); + const auto tris2verts = mesh.ask_elem_verts(); + const auto mesh_coords = mesh.coords(); FillCoordinatesDeviceFunctor fill_functor( - search_results_d, valid_indices_d, elem_counts_d, offsets_d_, - coordinates_d_, indices_d_, mesh.dim()); + search_results_d, valid_indices_d, owning_elem_ids, elem_counts_d, + offsets_d_, coordinates_d_, indices_d_, tris2verts, mesh_coords, + global_coords, mesh, mesh.dim()); Kokkos::parallel_for( "FillCoordinatesAndIndices", Kokkos::RangePolicy(0, num_valid_), diff --git a/src/pcms/field/evaluator/omega_h_lagrange.h b/src/pcms/field/evaluator/omega_h_lagrange.h index 827a99d3..781eff35 100644 --- a/src/pcms/field/evaluator/omega_h_lagrange.h +++ b/src/pcms/field/evaluator/omega_h_lagrange.h @@ -64,11 +64,12 @@ struct CopyCoordsFunctor template OmegaHLagrangeLocHint BuildLagrangeLocHint( - int mesh_dim, + Omega_h::Mesh& mesh, int mesh_dim, Kokkos::View::Result*, DeviceMemorySpace> results, - OutOfBoundsMode mode) + Kokkos::View coords_d, + Kokkos::View owning_ids, OutOfBoundsMode mode) { LO n = static_cast(results.size()); @@ -78,8 +79,7 @@ OmegaHLagrangeLocHint BuildLagrangeLocHint( "CheckValidity", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(LO i) { - bool out = (static_cast(results(i).dimensionality) != mesh_dim) || - (results(i).element_id < 0); + bool out = (owning_ids(i) < 0) || (results(i).element_id < 0); is_valid(i) = out ? 0 : 1; }); @@ -117,16 +117,58 @@ OmegaHLagrangeLocHint BuildLagrangeLocHint( update += val; }); + auto elem_verts = mesh.ask_elem_verts(); + auto mesh_coords = mesh.coords(); + Kokkos::parallel_for( "CompactData", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(LO i) { if (is_valid(i)) { LO k = valid_offsets(i); - elem_ids(k) = results(i).element_id; + LO owner_elem = owning_ids(i); + elem_ids(k) = owner_elem; orig_indices(k) = i; - for (int d = 0; d <= mesh_dim; ++d) - bary(k, d) = results(i).parametric_coords[d]; + + // Recompute barycentric coordinates in the owning element + // using global coordinates + int nvpe = mesh_dim + 1; + // Get vertices of the owning element + Omega_h::Few verts; + for (int v = 0; v < nvpe; ++v) { + verts[v] = elem_verts[owner_elem * nvpe + v]; + } + + // Get vertex coordinates + if (mesh_dim == 2) { + Omega_h::Few, 3> vertex_coords; + for (int v = 0; v < 3; ++v) { + vertex_coords[v] = Omega_h::get_vector<2>(mesh_coords, verts[v]); + } + Omega_h::Vector<2> point; + for (int d = 0; d < 2; ++d) { + point[d] = coords_d(i, d); + } + auto local = + Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); + for (int d = 0; d < 3; ++d) { + bary(k, d) = local[d]; + } + } else if (mesh_dim == 3) { + Omega_h::Few, 4> vertex_coords; + for (int v = 0; v < 4; ++v) { + vertex_coords[v] = Omega_h::get_vector<3>(mesh_coords, verts[v]); + } + Omega_h::Vector<3> point; + for (int d = 0; d < 3; ++d) { + point[d] = coords_d(i, d); + } + auto local = + Omega_h::barycentric_from_global<3, 3>(point, vertex_coords); + for (int d = 0; d < 4; ++d) { + bary(k, d) = local[d]; + } + } } else { LO k = missing_offsets(i); missing_indices(k) = i; @@ -291,6 +333,7 @@ class OmegaHLagrangeEvaluatorFactory : public FieldEvaluatorFactory auto raw_coords = coords.GetCoordinates(); LO n_pts = static_cast(raw_coords.extent(0)); int mesh_dim = layout_->GetMesh().dim(); + Omega_h::Mesh& mesh = const_cast(layout_->GetMesh()); OmegaHLagrangeLocHint hint = std::visit( [&](auto& search) { @@ -303,8 +346,10 @@ class OmegaHLagrangeEvaluatorFactory : public FieldEvaluatorFactory Kokkos::parallel_for("copy_coords", n_pts, copy_functor); auto results_d = search(coords_d); - return detail::BuildLagrangeLocHint(mesh_dim, results_d, - policy.mode); + auto owning_ids = search.GetOwningElementIds(results_d); + + return detail::BuildLagrangeLocHint( + mesh, mesh_dim, results_d, coords_d, owning_ids, policy.mode); }, search_); diff --git a/src/pcms/localization/adj_search.cpp b/src/pcms/localization/adj_search.cpp index e1b937d7..7dd25364 100644 --- a/src/pcms/localization/adj_search.cpp +++ b/src/pcms/localization/adj_search.cpp @@ -7,15 +7,17 @@ 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).element_id < 0) { - OMEGA_H_CHECK_PRINTF(results(i).element_id >= 0, - "ERROR: Source cell id not found for target %d\n", + (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); } @@ -90,14 +92,13 @@ static Omega_h::Write locate_target_cells( pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); auto results = search_cell(target_points); + auto owning_cell_ids = search_cell.GetOwningElementIds(results); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { - auto source_cell_id = results(i).element_id; - if (source_cell_id < 0) - source_cell_id = Kokkos::abs(source_cell_id); + auto source_cell_id = owning_cell_ids(i); OMEGA_H_CHECK_PRINTF( source_cell_id >= 0, - "ERROR: Source cell id not found for target %d (%f,%f)\n", i, + "ERROR: Source face id not found for target %d (%f,%f)\n", i, target_points(i, 0), target_points(i, 1)); source_cell_ids[i] = source_cell_id; }); @@ -114,13 +115,12 @@ static Omega_h::Write locate_target_cells( pcms::GridPointSearch3D search_cell(source_mesh, 10, 10, 10); auto results = search_cell(target_points); + auto owning_cell_ids = search_cell.GetOwningElementIds(results); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { - auto source_cell_id = results(i).element_id; - if (source_cell_id < 0) - source_cell_id = Kokkos::abs(source_cell_id); + auto source_cell_id = owning_cell_ids(i); OMEGA_H_CHECK_PRINTF(source_cell_id >= 0, - "ERROR: Source cell id not found for target %d\n", + "ERROR: Source face id not found for target %d\n", i); source_cell_ids[i] = source_cell_id; }); diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index 1b2bd45e..84a4dc04 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/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 @@ -144,36 +145,35 @@ template return true; } -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( - const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords, Real fuzz) -{ - 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)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({left, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({right, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - xi = Omega_h::barycentric_from_global<2, 2>({right, bot}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - return false; -} +// [[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( +// const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords, Real fuzz) +// { +// 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)) { +// return true; +// } +// xi = Omega_h::barycentric_from_global<2, 2>({left, top}, coords); +// if (Omega_h::is_barycentric_inside(xi, fuzz)) { +// return true; +// } +// xi = Omega_h::barycentric_from_global<2, 2>({right, top}, coords); +// if (Omega_h::is_barycentric_inside(xi, fuzz)) { +// return true; +// } +// xi = Omega_h::barycentric_from_global<2, 2>({right, bot}, coords); +// if (Omega_h::is_barycentric_inside(xi, fuzz)) { +// 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 } @@ -198,7 +198,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 * 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)) { @@ -220,7 +220,7 @@ template return true; } // if any of the bbox verts are within the triangle - if (bbox_verts_within_triangle(bbox, coords, fuzz)) { + if (bbox_verts_within_simplex(bbox, coords)) { return true; } // if any of the triangle's edges intersect with the bounding box @@ -256,12 +256,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) { @@ -276,7 +275,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) { @@ -285,7 +284,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; } @@ -300,7 +299,6 @@ struct GridTriIntersectionFunctor2D Omega_h::LOs tris2verts_; Omega_h::Reals coords_; Kokkos::View grid_; - Real fuzz_; public: LO nelems_; @@ -328,7 +326,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) { @@ -361,10 +359,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; } @@ -406,7 +404,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( @@ -416,111 +413,195 @@ 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]; - 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); + 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(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, fuzz)) { - 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; + 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; + 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; + bary_out = inside_face_bary; + } else { + element_id_out = best_edge_id; + 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; + 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; + 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; + 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, bary_out}; }); 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); @@ -531,14 +612,107 @@ 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_); + // 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(); 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; + 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; + + auto edges2faces_a2ab_h_ = Omega_h::HostRead(edges2faces_up_.a2ab); + auto edges2faces_ab2b_h_ = Omega_h::HostRead(edges2faces_up_.ab2b); + auto verts2faces_a2ab_h_ = Omega_h::HostRead(verts2faces_up_.a2ab); + auto verts2faces_ab2b_h_ = Omega_h::HostRead(verts2faces_up_.ab2b); + + const Omega_h::HostRead* a2ab_h = nullptr; + const Omega_h::HostRead* ab2b_h = nullptr; + if (dimensionality == Result::Dimensionality::EDGE) { + a2ab_h = &edges2faces_a2ab_h_; + ab2b_h = &edges2faces_ab2b_h_; + } else if (dimensionality == Result::Dimensionality::VERTEX) { + a2ab_h = &verts2faces_a2ab_h_; + ab2b_h = &verts2faces_ab2b_h_; + } else { + return -1; + } + + const auto begin = (*a2ab_h)[element_id]; + const auto end = (*a2ab_h)[element_id + 1]; + if (begin >= end) + return -1; + LO owner = (*ab2b_h)[begin]; + for (auto i = begin + 1; i < end; ++i) { + const LO candidate = (*ab2b_h)[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()( @@ -555,7 +729,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; @@ -583,7 +756,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}; @@ -594,26 +767,24 @@ 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, + parametric_coords_to_nearest}; } }); 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); @@ -639,6 +810,113 @@ 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; + 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; + + auto verts2regions_a2ab_h_ = Omega_h::HostRead(verts2regions_up_.a2ab); + auto verts2regions_ab2b_h_ = Omega_h::HostRead(verts2regions_up_.ab2b); + auto edges2regions_a2ab_h_ = Omega_h::HostRead(edges2regions_up_.a2ab); + auto edges2regions_ab2b_h_ = Omega_h::HostRead(edges2regions_up_.ab2b); + auto faces2regions_a2ab_h_ = Omega_h::HostRead(faces2regions_up_.a2ab); + auto faces2regions_ab2b_h_ = Omega_h::HostRead(faces2regions_up_.ab2b); + + const Omega_h::HostRead* a2ab_h = nullptr; + const Omega_h::HostRead* ab2b_h = nullptr; + if (dimensionality == Result::Dimensionality::FACE) { + a2ab_h = &faces2regions_a2ab_h_; + ab2b_h = &faces2regions_ab2b_h_; + } else if (dimensionality == Result::Dimensionality::EDGE) { + a2ab_h = &edges2regions_a2ab_h_; + ab2b_h = &edges2regions_ab2b_h_; + } else if (dimensionality == Result::Dimensionality::VERTEX) { + a2ab_h = &verts2regions_a2ab_h_; + ab2b_h = &verts2regions_ab2b_h_; + } else { + return -1; + } + + const auto begin = (*a2ab_h)[element_id]; + const auto end = (*a2ab_h)[element_id + 1]; + if (begin >= end) + return -1; + LO owner = (*ab2b_h)[begin]; + for (auto i = begin + 1; i < end; ++i) { + const LO candidate = (*ab2b_h)[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 cfcc2012..835d6c46 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -22,12 +22,76 @@ 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); + +/** + * 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 @@ -44,7 +108,7 @@ class PointLocalizationSearch }; Dimensionality dimensionality; - LO element_id; + LO element_id; // ID of the actual entity (vertex/edge/face/region) Omega_h::Vector parametric_coords; }; @@ -60,6 +124,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: @@ -77,9 +144,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 @@ -90,13 +157,19 @@ 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_; - Real fuzz_; 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_; @@ -111,10 +184,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 @@ -125,12 +197,20 @@ 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/mesh_intersection.cpp b/src/pcms/transfer/mesh_intersection.cpp index 9cb37de3..564ae350 100644 --- a/src/pcms/transfer/mesh_intersection.cpp +++ b/src/pcms/transfer/mesh_intersection.cpp @@ -30,6 +30,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( @@ -38,7 +39,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 6d8dc6a1..bfaaff78 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -100,6 +100,70 @@ 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) @@ -205,27 +269,33 @@ TEST_CASE("uniform grid search") { { auto [dim, idx, coords] = results_h(0); + const auto face_idx = search.GetOwningElementId(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); + const auto face_idx = search.GetOwningElementId(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); + const auto face_idx = search.GetOwningElementId(results_h(7)); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(idx == 0); + REQUIRE(face_idx >= 0); } } // feature needs to be added @@ -236,21 +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(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(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(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(search.GetOwningElementId(out_of_bounds) >= 0); } } From f4bedd1975e6a9f466cd2a33c35dddba94f70576 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 15 May 2026 17:10:58 -0400 Subject: [PATCH 2/3] add a new out of bound test --- .../field/evaluator/mesh_fields_backend.h | 24 +------ test/test_omega_h_field2_outofbounds.cpp | 71 +++++++++++++++++++ 2 files changed, 74 insertions(+), 21 deletions(-) diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index a895381e..84eb206d 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/src/pcms/field/evaluator/mesh_fields_backend.h @@ -198,10 +198,6 @@ struct CountPointsPerElementFunctor KOKKOS_INLINE_FUNCTION void operator()(LO i) const { - 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); } @@ -239,6 +235,8 @@ struct FillCoordinatesAndIndicesFunctor void operator()(LO i) const { auto [dim, elem_idx, coord] = search_results_(i); + (void)dim; + (void)elem_idx; const auto owner_idx = owning_elem_ids_(i); LO count = Kokkos::atomic_sub_fetch(&elem_counts_(owner_idx), 1); LO index = offsets_(owner_idx) + count - 1; @@ -274,9 +272,9 @@ struct FilterValidPointsFunctor void operator()(LO i, LO& num_valid, LO& num_missing) const { auto [dim, elem_idx, coord] = search_results_(i); - const auto owner_idx = owning_elem_ids_(i); (void)dim; (void)coord; + const auto owner_idx = owning_elem_ids_(i); bool is_missing = (elem_idx < 0) || (owner_idx < 0); if (mode_ == OutOfBoundsMode::ERROR && is_missing) { @@ -347,10 +345,6 @@ struct CountPerElementFunctor void operator()(LO i) const { LO orig_idx = valid_indices_(i); - 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); Kokkos::atomic_add(&elem_counts_(owner_idx), 1); } @@ -398,10 +392,6 @@ struct FillCoordinatesDeviceFunctor void operator()(LO i) const { LO orig_idx = valid_indices_(i); - 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); LO count = Kokkos::atomic_fetch_add(&elem_counts_(owner_idx), static_cast(-1)); @@ -520,10 +510,6 @@ struct MeshFieldsAdapter2LocalizationHint 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]); - (void)dim; - (void)elem_idx; - (void)coord; const auto owner_idx = owning_elem_ids(valid_point_indices[i]); elem_counts[owner_idx] += 1; } @@ -540,10 +526,6 @@ struct MeshFieldsAdapter2LocalizationHint 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, 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); diff --git a/test/test_omega_h_field2_outofbounds.cpp b/test/test_omega_h_field2_outofbounds.cpp index da0cd338..5e4c3813 100644 --- a/test/test_omega_h_field2_outofbounds.cpp +++ b/test/test_omega_h_field2_outofbounds.cpp @@ -4,7 +4,9 @@ #include #include #include "pcms/field/function_space/lagrange.h" +#include "pcms/field/function_space/spline.h" #include "pcms/field/field_metadata.h" +#include "pcms/utility/uniform_grid.h" #include "field_test_utils.h" #include #include @@ -40,3 +42,72 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") factory, field, coords, is_inside, OMEGA_H_LAMBDA(Real x, Real y) { return x + y; }, fill_value, 1.0e-10); } + +TEST_CASE("uniform_grid_field out of bounds FILL mode") +{ + const int N = 10; + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {1.0, 1.0}; + grid.divisions = {N, N}; + + auto factory = pcms::LagrangeFunctionSpace::FromUniformGrid( + grid, 1, pcms::CoordinateSystem::Cartesian, 1); + auto field = factory.CreateField(pcms::FieldMetadata{}); + pcms::test::SetField( + field.GetData(), *factory.GetLayout(), + OMEGA_H_LAMBDA(Real x, Real y) { return x + y; }); + + Real fill_value = -999.0; + + // Test points - mix of inside and outside the [0,1]^2 grid + std::vector coords = { + 0.5, 0.5, // inside - should evaluate normally + 1.5, 0.5, // outside (x > 1) - should return fill_value + 0.5, -0.1, // outside (y < 0) - should return fill_value + 0.3, 0.7, // inside - should evaluate normally + -0.1, 0.5, // outside (x < 0) - should return fill_value + }; + + std::vector is_inside = {true, false, false, true, false}; + pcms::test::CheckEvaluationWithFill( + factory, field, coords, is_inside, + OMEGA_H_LAMBDA(Real x, Real y) { return x + y; }, fill_value, 1.0e-8); +} + +TEST_CASE("spline_field out of bounds FILL mode") +{ + const int N = 10; + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {1.0, 1.0}; + grid.divisions = {N, N}; + + auto factory = pcms::SplineFunctionSpace::FromUniformGrid( + grid, pcms::CoordinateSystem::Cartesian); + auto field = factory.CreateField(pcms::FieldMetadata{}); + pcms::test::SetField( + field.GetData(), *factory.GetLayout(), + OMEGA_H_LAMBDA(Real x, Real y) { return x + y; }); + + Real fill_value = -999.0; + + // Test points - mix of inside and outside the [0,1]^2 grid + std::vector coords = { + 0.5, 0.5, // inside - should evaluate normally + 1.5, 0.5, // outside (x > 1) - should return fill_value + 0.5, -0.1, // outside (y < 0) - should return fill_value + 0.3, 0.7, // inside - should evaluate normally + -0.1, 0.5, // outside (x < 0) - should return fill_value + }; + + std::vector is_inside = {true, false, false, true, false}; + pcms::test::CheckEvaluationWithFill( + factory, field, coords, is_inside, + OMEGA_H_LAMBDA(Real x, Real y) { return x + y; }, fill_value, 1.0e-8); +} + +// Note: polynomial_reconstruction_field does not support out-of-bounds FILL +// mode. The MLS (Moving Least Squares) implementation always attempts +// interpolation/ extrapolation for all points and does not check the +// OutOfBoundsPolicy. From ce0289e925fb3cc862c43d08683b34e7ac6a5d1e Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Mon, 18 May 2026 15:15:06 -0400 Subject: [PATCH 3/3] comment GetOwningElementId --- src/pcms/localization/point_search.cpp | 26 -------------------------- src/pcms/localization/point_search.h | 10 ++++++++++ 2 files changed, 10 insertions(+), 26 deletions(-) diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index 84a4dc04..9c8ab150 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/point_search.cpp @@ -145,32 +145,6 @@ template return true; } -// [[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( -// const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords, Real fuzz) -// { -// 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)) { -// return true; -// } -// xi = Omega_h::barycentric_from_global<2, 2>({left, top}, coords); -// if (Omega_h::is_barycentric_inside(xi, fuzz)) { -// return true; -// } -// xi = Omega_h::barycentric_from_global<2, 2>({right, top}, coords); -// if (Omega_h::is_barycentric_inside(xi, fuzz)) { -// return true; -// } -// xi = Omega_h::barycentric_from_global<2, 2>({right, bot}, coords); -// if (Omega_h::is_barycentric_inside(xi, fuzz)) { -// return true; -// } -// return false; -// } - template [[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_simplex( const AABBox& bbox, const Omega_h::Matrix& coords) diff --git a/src/pcms/localization/point_search.h b/src/pcms/localization/point_search.h index 835d6c46..f66cc89a 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -124,6 +124,16 @@ class PointLocalizationSearch virtual Kokkos::View operator()( Kokkos::View point) const = 0; + + /** + * This function provides a temporary solution to retrieve the original + * behavior of the point search, which previously returned the face ID + * regardless of which entity the search result belonged to. Many parts of + * the codebase still rely on this legacy behavior. After updating the point + * search to return the exact entity, this function can be called to retrieve + * the old behavior and ensure correctness. Long term, the implementation + * should be updated to properly handle the new result. + */ [[nodiscard]] virtual LO GetOwningElementId(const Result& result) const = 0; [[nodiscard]] virtual Kokkos::View GetOwningElementIds( Kokkos::View results) const = 0;