diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter.h b/src/pcms/adapter/meshfields/mesh_fields_adapter.h index b1460b00..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, 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(elem_idx >= 0); + 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); + 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, 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(elem_idx >= 0); + 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); + 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 291b47de..34dd8079 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h @@ -122,11 +122,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) { } @@ -134,7 +138,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); } }; @@ -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()) { } @@ -167,13 +178,14 @@ struct FillCoordinatesAndIndicesFunctor void operator()(LO i) const { 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()); - // 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_(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 @@ -197,8 +210,10 @@ struct MeshFieldsAdapter2LocalizationHint // 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); + 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); } @@ -206,8 +221,10 @@ struct MeshFieldsAdapter2LocalizationHint // 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); + 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 { @@ -242,9 +259,14 @@ struct MeshFieldsAdapter2LocalizationHint // 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; + (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 @@ -257,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, 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); } } @@ -428,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 1b2bd45e..148b927b 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,9 @@ 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, - 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 +172,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 +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)) { @@ -220,7 +194,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 +230,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 +249,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 +258,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 +273,6 @@ struct GridTriIntersectionFunctor2D Omega_h::LOs tris2verts_; Omega_h::Reals coords_; Kokkos::View grid_; - Real fuzz_; public: LO nelems_; @@ -328,7 +300,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 +333,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 +378,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 +387,193 @@ 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, 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 +584,99 @@ 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; + + 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()( @@ -555,7 +693,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 +720,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 +731,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 +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); - 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; + + 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 cfcc2012..ce9fa2b9 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -6,6 +6,7 @@ #include #include #include +#include #include "pcms/utility/types.h" #include "pcms/utility/uniform_grid.h" @@ -22,12 +23,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); +} // namespace detail [[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 +109,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 +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: @@ -77,9 +145,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 +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_; - 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 +186,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,17 +199,25 @@ 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_; Omega_h::Reals coords_; - Real fuzz_; }; } // namespace pcms diff --git a/src/pcms/transfer/adj_search.cpp b/src/pcms/transfer/adj_search.cpp index 4ff86e02..e1ccfcdc 100644 --- a/src/pcms/transfer/adj_search.cpp +++ b/src/pcms/transfer/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); } @@ -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,15 +158,11 @@ 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).element_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, - "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/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, 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 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); } }