From 42bf41cab0d35ce0c82a39098dc35802c64d0dfc Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 21 Mar 2025 13:06:42 -0400 Subject: [PATCH 01/22] Refactor GridPointSearch to support dimensional template specialization. Updated `GridPointSearch` to use a templated design (`GridPointSearch`) and introduced `GridPointSearch2D` as a specific case for 2D point localization. This change ensures clearer extensibility for other dimensions and improves code maintainability by centralizing shared logic in a new base class (`PointLocalizationSearch`). Adjusted test cases, headers, and related implementations accordingly. --- src/pcms/adapter/omega_h/omega_h_field.h | 7 +- src/pcms/interpolator/adj_search.hpp | 4 +- src/pcms/point_search.cpp | 265 ++++++++++------------- src/pcms/point_search.h | 42 ++-- test/test_point_search.cpp | 75 +++---- 5 files changed, 178 insertions(+), 215 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index c0aaa17e..124f43e5 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -154,7 +154,7 @@ class OmegaHField void ConstructSearch(int nx, int ny) { PCMS_FUNCTION_TIMER; - search_ = GridPointSearch(mesh_, nx, ny); + search_ = GridPointSearch2D(mesh_, nx, ny); } // pass through to search function [[nodiscard]] auto Search(Kokkos::View points) const @@ -220,9 +220,8 @@ class OmegaHField private: std::string name_; Omega_h::Mesh& mesh_; - // TODO make this a pointer and introduce base class to Search for alternative - // search methods - std::optional search_; + // TODO make this a pointer and introduce base class to Search for alternative search methods + std::optional search_; // bitmask array that specifies a filter on the field Omega_h::Read mask_; LO size_; diff --git a/src/pcms/interpolator/adj_search.hpp b/src/pcms/interpolator/adj_search.hpp index c39639fc..7f5cc2e2 100644 --- a/src/pcms/interpolator/adj_search.hpp +++ b/src/pcms/interpolator/adj_search.hpp @@ -24,7 +24,7 @@ Omega_h::Real calculateDistance(const Omega_h::Real* p1, } inline void checkTargetPoints( - const Kokkos::View& results) + const Kokkos::View& results) { Kokkos::fence(); printf("INFO: Checking target points...\n"); @@ -117,7 +117,7 @@ inline void FindSupports::adjBasedSearch( }); Kokkos::fence(); - pcms::GridPointSearch search_cell(source_mesh, 10, 10); + pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); auto results = search_cell(target_points); checkTargetPoints(results); diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index e86cab36..8f5b9b96 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -2,35 +2,29 @@ #include #include -// From -// https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points +// From https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points KOKKOS_INLINE_FUNCTION -double distance_from_line(const double x0, const double y0, const double x1, - const double y1, const double x2, const double y2) +double distance_from_line(const double x0, const double y0, const double x1, const double y1, const double x2, const double y2) { - const Omega_h::Vector<2> p1 = {x1, y1}; - const Omega_h::Vector<2> p2 = {x2, y2}; + const Omega_h::Vector<2> p1 = { x1, y1 }; + const Omega_h::Vector<2> p2 = { x2, y2 }; auto disp = p2 - p1; - return std::abs(disp[1] * x0 - disp[0] * y0 + x2 * y1 - y2 * x1) / - Omega_h::norm(disp); + return std::abs(disp[1]*x0 - disp[0]*y0 + x2*y1 - y2*x1) / Omega_h::norm(disp); } -// Law of Cosines, where a, b, c and gamma are defined here: -// https://en.wikipedia.org/wiki/Law_of_cosines#Use_in_solving_triangles +// Law of Cosines, where a, b, c and gamma are defined here: https://en.wikipedia.org/wiki/Law_of_cosines#Use_in_solving_triangles KOKKOS_INLINE_FUNCTION double angle_from_side_lengths(const double a, const double b, const double c) { - return std::acos((a * a + b * b - c * c) / 2 * a * b); + return std::acos((a*a + b*b - c*c) / 2*a*b); } KOKKOS_INLINE_FUNCTION -bool normal_intersects_segment(const Omega_h::Few a, - const Omega_h::Few b, - const Omega_h::Few c) +bool normal_intersects_segment(const Omega_h::Few a, const Omega_h::Few b, const Omega_h::Few c) { const auto ab_len = Omega_h::norm(a - b); - const auto bc_len = Omega_h::norm(b - c); + const auto bc_len = Omega_h::norm(b -c); const auto ac_len = Omega_h::norm(a - c); const double angle1 = angle_from_side_lengths(bc_len, ac_len, ab_len); @@ -117,40 +111,29 @@ bool line_intersects_bbox(const Omega_h::Vector<2>& p0, return false; } -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool within_bbox( - const Omega_h::Vector<2> coord, const AABBox<2>& bbox) noexcept -{ +[[nodiscard]] KOKKOS_INLINE_FUNCTION +bool within_bbox(const Omega_h::Vector<2> coord, const AABBox<2> & bbox) noexcept { 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]; - return (coord[0] >= left) && (coord[0] <= right) && (coord[1] >= bot) && - (coord[1] <= top); + return (coord[0]>=left) && (coord[0]<=right) && (coord[1] >= bot) && (coord[1]<=top); } -[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( - const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords) -{ +[[nodiscard]] KOKKOS_INLINE_FUNCTION +bool bbox_verts_within_triangle(const AABBox<2>& bbox, const Omega_h::Matrix<2,3>& coords) { auto left = bbox.center[0] - bbox.half_width[0]; auto right = bbox.center[0] + bbox.half_width[0]; auto bot = bbox.center[1] - bbox.half_width[1]; auto top = bbox.center[1] + bbox.half_width[1]; - auto xi = barycentric_from_global({left, bot}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - xi = barycentric_from_global({left, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - xi = barycentric_from_global({right, top}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } - xi = barycentric_from_global({right, bot}, coords); - if (Omega_h::is_barycentric_inside(xi, fuzz)) { - return true; - } + auto xi = barycentric_from_global({left,bot}, coords); + if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} + xi = barycentric_from_global({left,top}, coords); + if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} + xi = barycentric_from_global({right,top}, coords); + if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} + xi = barycentric_from_global({right,bot}, coords); + if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} return false; } @@ -159,13 +142,15 @@ bool line_intersects_bbox(const Omega_h::Vector<2>& p0, * intersects with a bounding box */ [[nodiscard]] -KOKKOS_FUNCTION bool triangle_intersects_bbox( - const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) +KOKKOS_FUNCTION +bool triangle_intersects_bbox(const Omega_h::Matrix<2, 3>& coords, + const AABBox<2>& bbox) { // triangle and grid cell bounding box intersect if (intersects(triangle_bbox(coords), bbox)) { // if any of the triangle verts inside of bbox - if (within_bbox(coords[0], bbox) || within_bbox(coords[1], bbox) || + if (within_bbox(coords[0], bbox) || + within_bbox(coords[1], bbox) || within_bbox(coords[2], bbox)) { return true; } @@ -197,8 +182,7 @@ namespace detail */ struct GridTriIntersectionFunctor { - GridTriIntersectionFunctor(Omega_h::Mesh& mesh, - Kokkos::View grid) + GridTriIntersectionFunctor(Omega_h::Mesh& mesh, Kokkos::View grid) : mesh_(mesh), tris2verts_(mesh_.ask_elem_verts()), coords_(mesh_.coords()), @@ -221,11 +205,9 @@ struct GridTriIntersectionFunctor LO num_intersections = 0; // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { - const auto elem_tri2verts = - Omega_h::gather_verts<3>(tris2verts_, elem_idx); + const auto elem_tri2verts = Omega_h::gather_verts<3>(tris2verts_, elem_idx); // 2d mesh with 2d coords, but 3 triangles - const auto vertex_coords = - Omega_h::gather_vectors<3, 2>(coords_, elem_tri2verts); + const auto vertex_coords = Omega_h::gather_vectors<3, 2>(coords_, elem_tri2verts); if (triangle_intersects_bbox(vertex_coords, grid_cell_bbox)) { if (fill) { fill[num_intersections] = elem_idx; @@ -241,17 +223,14 @@ struct GridTriIntersectionFunctor Omega_h::LOs tris2verts_; Omega_h::Reals coords_; Kokkos::View grid_; - public: LO nelems_; }; -// num_grid_cells should be result of grid.GetNumCells(), take as argument to -// avoid extra copy of grid from gpu to cpu +// num_grid_cells should be result of grid.GetNumCells(), take as argument to avoid extra copy +// of grid from gpu to cpu Kokkos::Crs -construct_intersection_map(Omega_h::Mesh& mesh, - Kokkos::View grid, - int num_grid_cells) +construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { Kokkos::Crs intersection_map{}; auto f = detail::GridTriIntersectionFunctor{mesh, grid}; @@ -272,137 +251,117 @@ Omega_h::Vector<3> barycentric_from_global( return {1 - xi[0] - xi[1], xi[0], xi[1]}; } -template -OMEGA_H_INLINE double myreduce(const Omega_h::Vector& x, - Op op) OMEGA_H_NOEXCEPT -{ +template +OMEGA_H_INLINE double myreduce(const Omega_h::Vector & x, Op op) OMEGA_H_NOEXCEPT { auto out = x[0]; - for (int i = 1; i < n; ++i) - out = op(out, x[i]); + for (int i = 1; i < n; ++i) out = op(out, x[i]); return out; -} +} -Kokkos::View GridPointSearch::operator()( - Kokkos::View points) const +Kokkos::View GridPointSearch2D::operator()(Kokkos::View points) const { - static_assert(dim == 2, "point search assumes dim==2"); - Kokkos::View results("point search result", - points.extent(0)); + Kokkos::View results("point search result", points.extent(0)); auto num_rows = candidate_map_.numRows(); // needed so that we don't capture this ptr which will be memory error on cuda - auto grid = grid_; + auto grid = grid_; auto candidate_map = candidate_map_; auto tris2verts = tris2verts_; auto tris2verts_adj = tris2verts_adj_; auto tris2edges_adj = tris2edges_adj_; auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; - Kokkos::parallel_for( - points.extent(0), KOKKOS_LAMBDA(int p) { - Omega_h::Vector<2> point( - std::initializer_list{points(p, 0), points(p, 1)}); - auto cell_id = grid(0).ClosestCellID(point); - assert(cell_id < num_rows && cell_id >= 0); - auto candidates_begin = candidate_map.row_map(cell_id); - auto candidates_end = candidate_map.row_map(cell_id + 1); - bool found = false; - - auto nearest_triangle = candidates_begin; - auto dimensionality = GridPointSearch::Result::Dimensionality::EDGE; - 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); - 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 = barycentric_from_global(point, vertex_coords); - - if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { - results(p) = GridPointSearch::Result{ - GridPointSearch::Result::Dimensionality::FACE, triangleID, - parametric_coords}; - found = true; - break; - } + Kokkos::parallel_for(points.extent(0), KOKKOS_LAMBDA(int p) { + Omega_h::Vector<2> point(std::initializer_list{points(p,0), points(p,1)}); + auto cell_id = grid(0).ClosestCellID(point); + assert(cell_id < num_rows && cell_id >= 0); + auto candidates_begin = candidate_map.row_map(cell_id); + auto candidates_end = candidate_map.row_map(cell_id + 1); + bool found = false; + + auto nearest_triangle = candidates_begin; + auto dimensionality = GridPointSearch::Result::Dimensionality::EDGE; + 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); + 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 = barycentric_from_global(point, vertex_coords); - for (int j = 0; j < 3; ++j) { - // Every triangle (face) is connected to 3 edges - const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; + if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { + results(p) = GridPointSearch::Result{GridPointSearch::Result::Dimensionality::FACE, triangleID, parametric_coords}; + found = true; + break; + } - auto vertex_a_id = edges2verts_adj.ab2b[edgeID * 2]; - auto vertex_b_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; + 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 = Omega_h::get_vector<2>(coords, vertex_a_id); - auto vertex_b = Omega_h::get_vector<2>(coords, vertex_b_id); + auto vertex_a_id = edges2verts_adj.ab2b[edgeID * 2]; + auto vertex_b_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; - if (!normal_intersects_segment(point, vertex_a, vertex_b)) - continue; + 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 auto xa = vertex_a[0]; - const auto ya = vertex_a[1]; - const auto xb = vertex_b[0]; - const auto yb = vertex_b[1]; + if (!normal_intersects_segment(point, vertex_a, vertex_b)) continue; - const auto xp = point[0]; - const auto yp = point[1]; + const auto xa = vertex_a[0]; + const auto ya = vertex_a[1]; + const auto xb = vertex_b[0]; + const auto yb = vertex_b[1]; - const auto distance_to_ab = - distance_from_line(xp, yp, xa, ya, xb, yb); + const auto xp = point[0]; + const auto yp = point[1]; - if (distance_to_ab >= distance_to_nearest) { - continue; - } + const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); - dimensionality = GridPointSearch::Result::Dimensionality::EDGE; - nearest_triangle = i; - distance_to_nearest = distance_to_ab; - parametric_coords_to_nearest = parametric_coords; - } + if (distance_to_ab >= distance_to_nearest) { continue; } - // Every triangle (face) is connected to 3 vertices - 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); - - if (const auto distance = Omega_h::norm(point - vertex); - distance < distance_to_nearest) { - dimensionality = GridPointSearch::Result::Dimensionality::VERTEX; - nearest_triangle = i; - distance_to_nearest = distance; - parametric_coords_to_nearest = parametric_coords; - } - } + dimensionality = GridPointSearch::Result::Dimensionality::EDGE; + nearest_triangle = i; + distance_to_nearest = distance_to_ab; + parametric_coords_to_nearest = parametric_coords; } - if (!found) { - results(p) = GridPointSearch::Result{ - dimensionality, -1 * candidate_map.entries(nearest_triangle), - parametric_coords_to_nearest}; + + // Every triangle (face) is connected to 3 vertices + 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); + + if (const auto distance = Omega_h::norm(point - vertex);distance < distance_to_nearest) { + dimensionality = GridPointSearch::Result::Dimensionality::VERTEX; + nearest_triangle = i; + distance_to_nearest = distance; + parametric_coords_to_nearest = parametric_coords; + } } - }); + } + if(!found) + { + results(p) = GridPointSearch::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; + } + }); return results; } -GridPointSearch::GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny) +GridPointSearch<2>::GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); - grid_h(0) = - Uniform2DGrid{.edge_length = {mesh_bbox.max[0] - mesh_bbox.min[0], - mesh_bbox.max[1] - mesh_bbox.min[1]}, - .bot_left = {mesh_bbox.min[0], mesh_bbox.min[1]}, - .divisions = {Nx, Ny}}; + grid_h(0) = Uniform2DGrid{.edge_length = {mesh_bbox.max[0] - mesh_bbox.min[0], + mesh_bbox.max[1] - mesh_bbox.min[1]}, + .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(mesh, grid_, grid_h(0).GetNumCells()); + candidate_map_ = detail::construct_intersection_map(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); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 05aff657..33b7ab47 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -15,12 +15,9 @@ namespace pcms // TODO take a bounding box as we may want a bbox that's bigger than the mesh! // this function is in the public header for testing, but should not be directly // used -namespace detail -{ +namespace detail { Kokkos::Crs -construct_intersection_map(Omega_h::Mesh& mesh, - Kokkos::View grid, - int num_grid_cells); +construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); } KOKKOS_FUNCTION Omega_h::Vector<3> barycentric_from_global( @@ -29,15 +26,11 @@ Omega_h::Vector<3> barycentric_from_global( [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); -class GridPointSearch +template +class PointLocalizationSearch { - using CandidateMapT = - Kokkos::Crs; - public: - static constexpr auto dim = 2; - struct Result - { + struct Result { enum class Dimensionality { VERTEX = 0, @@ -50,6 +43,25 @@ class GridPointSearch Omega_h::Vector parametric_coords; }; + static constexpr auto DIM = 2; + + virtual Kokkos::View operator()(Kokkos::View point) const = 0; +}; + +template +class GridPointSearch : public PointLocalizationSearch +{ + static_assert(false, "Not implemented"); +}; + +template <> +class GridPointSearch<2> : public PointLocalizationSearch<2> +{ + using CandidateMapT = Kokkos::Crs; + +public: + using Result = PointLocalizationSearch::Result; + GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny); /** * given a point in global coordinates give the id of the triangle that the @@ -58,7 +70,7 @@ class GridPointSearch * id will be a negative number and (TODO) will return a negative id of the * closest element */ - Kokkos::View operator()(Kokkos::View point) const; + Kokkos::View operator()(Kokkos::View point) const override; private: Omega_h::Mesh mesh_; @@ -71,5 +83,7 @@ class GridPointSearch Omega_h::Reals coords_; }; -} // namespace pcms +using GridPointSearch2D = GridPointSearch<2>; + +} // namespace detail #endif // PCMS_COUPLING_POINT_SEARCH_H diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 06bacf48..7590efe3 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -112,8 +112,7 @@ bool num_candidates_within_range(const T& intersection_map, pcms::LO min, Kokkos::parallel_reduce( intersection_map.numRows(), KOKKOS_LAMBDA(const int i, result_type& update) { - auto num_candidates = - intersection_map.row_map(i + 1) - intersection_map.row_map(i); + auto num_candidates = intersection_map.row_map(i + 1) - intersection_map.row_map(i); if (num_candidates > update.max_val) update.max_val = num_candidates; if (num_candidates < update.min_val) @@ -125,11 +124,11 @@ bool num_candidates_within_range(const T& intersection_map, pcms::LO min, std::cerr << result.min_val << ' ' << result.max_val << '\n'; return within_range; } -// extern Omega_h::Library omega_h_library; +//extern Omega_h::Library omega_h_library; TEST_CASE("construct intersection map") { - // auto world = omega_h_library.world(); + //auto world = omega_h_library.world(); auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = @@ -139,11 +138,9 @@ TEST_CASE("construct intersection map") { Kokkos::View grid_d("uniform grid"); auto grid_h = Kokkos::create_mirror_view(grid_d); - grid_h(0) = Uniform2DGrid{ - .edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; + grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; Kokkos::deep_copy(grid_d, grid_h); - auto intersection_map = pcms::detail::construct_intersection_map( - mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map(mesh, grid_d, grid_h(0).GetNumCells()); // assert(cudaSuccess == cudaDeviceSynchronize()); REQUIRE(intersection_map.numRows() == 100); REQUIRE(num_candidates_within_range(intersection_map, 2, 16)); @@ -152,38 +149,35 @@ TEST_CASE("construct intersection map") { Kokkos::View grid_d("uniform grid"); auto grid_h = Kokkos::create_mirror_view(grid_d); - grid_h(0) = Uniform2DGrid{ - .edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {60, 60}}; + grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {60, 60}}; Kokkos::deep_copy(grid_d, grid_h); // require number of candidates is >=1 and <=6 - auto intersection_map = pcms::detail::construct_intersection_map( - mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map(mesh, grid_d, grid_h(0).GetNumCells()); REQUIRE(intersection_map.numRows() == 3600); REQUIRE(num_candidates_within_range(intersection_map, 1, 6)); } } -TEST_CASE("uniform grid search") -{ - using pcms::GridPointSearch; +TEST_CASE("uniform grid search") { + using pcms::GridPointSearch2D; auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); - GridPointSearch search{mesh, 10, 10}; - Kokkos::View points("test_points", 7); - // Kokkos::View points("test_points", 1); + GridPointSearch2D search{mesh,10,10}; + Kokkos::View points("test_points", 7); + //Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); - points_h(0, 0) = 0; - points_h(0, 1) = 0; - points_h(1, 0) = 0.55; - points_h(1, 1) = 0.54; - points_h(2, 0) = 100; - points_h(2, 1) = 100; - points_h(3, 0) = 1; - points_h(3, 1) = 1; - points_h(4, 0) = -1; - points_h(4, 1) = -1; + points_h(0,0) = 0; + points_h(0,1) = 0; + points_h(1,0) = 0.55; + points_h(1,1) = 0.54; + points_h(2,0) = 100; + points_h(2,1) = 100; + points_h(3,0) = 1; + points_h(3,1) = 1; + points_h(4,0) = -1; + points_h(4,1) = -1; points_h(5, 0) = 1.01; points_h(5, 1) = 0.95; points_h(6, 0) = 0.05; @@ -195,16 +189,16 @@ TEST_CASE("uniform grid search") SECTION("global coordinate within mesh") { { - auto [dim, idx, coords] = results_h(0); - REQUIRE(dim == GridPointSearch::Result::Dimensionality::FACE); + auto [dim, idx,coords] = results_h(0); + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(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); - REQUIRE(dim == GridPointSearch::Result::Dimensionality::FACE); + auto [dim, idx,coords] = results_h(1); + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(idx == 91); REQUIRE(coords[0] == Catch::Approx(0.5)); REQUIRE(coords[1] == Catch::Approx(0.1)); @@ -212,28 +206,25 @@ TEST_CASE("uniform grid search") } } // feature needs to be added - SECTION("Global coordinate outside mesh", "[!mayfail]") - { + SECTION("Global coordinate outside mesh", "[!mayfail]") { auto out_of_bounds = results_h(2); auto top_right = results_h(3); - REQUIRE(out_of_bounds.dimensionality == - GridPointSearch::Result::Dimensionality::VERTEX); - REQUIRE(-1 * out_of_bounds.tri_id == top_right.tri_id); + REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(-1*out_of_bounds.tri_id == top_right.tri_id); out_of_bounds = results_h(4); auto bot_left = results_h(0); - REQUIRE(out_of_bounds.dimensionality == - GridPointSearch::Result::Dimensionality::VERTEX); - REQUIRE(-1 * out_of_bounds.tri_id == bot_left.tri_id); + REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(-1*out_of_bounds.tri_id == bot_left.tri_id); out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == - GridPointSearch::Result::Dimensionality::EDGE); + GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(-1 * out_of_bounds.tri_id == top_right.tri_id); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == - GridPointSearch::Result::Dimensionality::EDGE); + GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(-1 * out_of_bounds.tri_id == bot_left.tri_id); } } From 2e6c364d4c38192dfb5e82e353ff028b5629fd7d Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 21 Mar 2025 13:14:53 -0400 Subject: [PATCH 02/22] Switch `optional` to `unique_ptr` Replaced `std::optional` with `std::unique_ptr` for `search_` to better handle polymorphism and align with future extensions. Updated type aliasing and constructors accordingly to enhance code clarity and maintainability. --- src/pcms/adapter/omega_h/omega_h_field.h | 7 +++---- src/pcms/point_search.h | 6 ++++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index 124f43e5..52784e89 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -154,14 +154,13 @@ class OmegaHField void ConstructSearch(int nx, int ny) { PCMS_FUNCTION_TIMER; - search_ = GridPointSearch2D(mesh_, nx, ny); + search_ = std::make_unique(mesh_, nx, ny); } // pass through to search function [[nodiscard]] auto Search(Kokkos::View points) const { PCMS_FUNCTION_TIMER; - PCMS_ALWAYS_ASSERT(search_.has_value() && - "search data structure must be constructed before use"); + PCMS_ALWAYS_ASSERT(search_ != nullptr && "search data structure must be constructed before use"); return (*search_)(points); } @@ -221,7 +220,7 @@ class OmegaHField std::string name_; Omega_h::Mesh& mesh_; // TODO make this a pointer and introduce base class to Search for alternative search methods - std::optional search_; + std::unique_ptr search_; // bitmask array that specifies a filter on the field Omega_h::Read mask_; LO size_; diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 33b7ab47..ba785ca2 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -48,6 +48,8 @@ class PointLocalizationSearch virtual Kokkos::View operator()(Kokkos::View point) const = 0; }; +using PointLocalizationSearch2D = PointLocalizationSearch<2>; + template class GridPointSearch : public PointLocalizationSearch { @@ -55,12 +57,12 @@ class GridPointSearch : public PointLocalizationSearch }; template <> -class GridPointSearch<2> : public PointLocalizationSearch<2> +class GridPointSearch<2> : public PointLocalizationSearch2D { using CandidateMapT = Kokkos::Crs; public: - using Result = PointLocalizationSearch::Result; + using Result = PointLocalizationSearch2D::Result; GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny); /** From b848d2e1b3ddd7a43bcc3fd5b379312178423a46 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Sat, 5 Apr 2025 13:14:14 -0400 Subject: [PATCH 03/22] Add 3D grid point search functionality and refactor 2D implementation Introduced support for 3D point localization and grid intersection with a new `GridTriIntersectionFunctor3D`. Refactored 2D implementations for consistency, generalizing methods to handle dimensions dynamically. Added templates and utilities to support both 2D and 3D use cases effectively. --- src/pcms/point_search.cpp | 122 +++++++++++++++++++++++++++++++------- src/pcms/point_search.h | 15 +---- src/pcms/uniform_grid.h | 5 +- 3 files changed, 108 insertions(+), 34 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 8f5b9b96..c38aaf84 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -51,6 +51,31 @@ AABBox<2> triangle_bbox(const Omega_h::Matrix<2, 3>& coords) return {.center = {(max[0] + min[0]) / 2.0, (max[1] + min[1]) / 2.0}, .half_width = {(max[0] - min[0]) / 2.0, (max[1] - min[1]) / 2.0}}; } + +template +AABBox simplex_bbox(const Omega_h::Matrix& coords) +{ + std::array max{coords(0, 0), coords(1, 0)}; + std::array min{coords(0, 0), coords(1, 0)}; + + for (int i = 1; i < dim + 1; ++i) { + for (int j = 0; j < dim; ++j) { + max[j] = std::fmax(max[j], coords(j, i)); + min[j] = std::fmin(min[j], coords(j, i)); + } + } + + std::array center; + std::array half_width; + + for (int j = 0; j < dim; ++j) { + center[j] = (max[j] + min[j]) / 2.0; + half_width[j] = (max[j] - min[j]) / 2.0; + } + + return {.center = center, .half_width = half_width}; +} + // Liang, You-Dong, and B. A. Barsky. “A New Concept and Method for Line // Clipping.” ACM Transactions on Graphics 3, no. 1 (January 1984): 1–22. // https://doi.org/10.1145/357332.357333. @@ -111,13 +136,14 @@ bool line_intersects_bbox(const Omega_h::Vector<2>& p0, return false; } +template [[nodiscard]] KOKKOS_INLINE_FUNCTION -bool within_bbox(const Omega_h::Vector<2> coord, const AABBox<2> & bbox) noexcept { - 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]; - return (coord[0]>=left) && (coord[0]<=right) && (coord[1] >= bot) && (coord[1]<=top); +bool within_bbox(const Omega_h::Vector coord, const AABBox & bbox) noexcept { + for (int i = 0; i < dim; ++i) { + if (coord[i] < bbox.center[i] - bbox.half_width[i]) return false; + if (coord[i] > bbox.center[i] + bbox.half_width[i]) return false; + } + return true; } [[nodiscard]] KOKKOS_INLINE_FUNCTION @@ -149,9 +175,9 @@ bool triangle_intersects_bbox(const Omega_h::Matrix<2, 3>& coords, // triangle and grid cell bounding box intersect if (intersects(triangle_bbox(coords), bbox)) { // if any of the triangle verts inside of bbox - if (within_bbox(coords[0], bbox) || - within_bbox(coords[1], bbox) || - within_bbox(coords[2], bbox)) { + if (within_bbox<2>(coords[0], bbox) || + within_bbox<2>(coords[1], bbox) || + within_bbox<2>(coords[2], bbox)) { return true; } // if any of the bbox verts are within the triangle @@ -168,6 +194,15 @@ bool triangle_intersects_bbox(const Omega_h::Matrix<2, 3>& coords, return false; } +template +[[nodiscard]] +KOKKOS_FUNCTION +bool simplex_intersects_bbox(const Omega_h::Matrix& coords, const AABBox& bbox) +{ + return intersects(simplex_bbox(coords), bbox); + // TODO: Add refined cases from triangle_intersects_bbox +} + namespace detail { /** @@ -180,9 +215,9 @@ namespace detail * \Warning since this uses Omega_h data which is only available in the * "Default" Execution space, the should not be used in an alternative EXE space */ -struct GridTriIntersectionFunctor +struct GridTriIntersectionFunctor2D { - GridTriIntersectionFunctor(Omega_h::Mesh& mesh, Kokkos::View grid) + GridTriIntersectionFunctor2D(Omega_h::Mesh& mesh, Kokkos::View grid) : mesh_(mesh), tris2verts_(mesh_.ask_elem_verts()), coords_(mesh_.coords()), @@ -190,7 +225,7 @@ struct GridTriIntersectionFunctor nelems_(mesh_.nelems()) { if (mesh_.dim() != 2) { - std::cerr << "GridTriIntersection currently only developed for 2D " + std::cerr << "GridTriIntersection2D currently only developed for 2D " "triangular meshes\n"; std::terminate(); } @@ -227,13 +262,60 @@ struct GridTriIntersectionFunctor LO nelems_; }; +struct GridTriIntersectionFunctor3D +{ + GridTriIntersectionFunctor3D(Omega_h::Mesh& mesh, Kokkos::View grid) + : mesh_(mesh), + tets2verts_(mesh_.ask_elem_verts()), + coords_(mesh_.coords()), + grid_(grid), + nelems_(mesh_.nelems()) + { + if (mesh_.dim() != 3) { + std::cerr << "GridTriIntersection3D currently only developed for 3D " + "triangular meshes\n"; + std::terminate(); + } + } + /// Two-pass functor. On the first pass we set the number of grid/triangle + /// intersections. On the second pass we fill the CSR array with the indexes + /// to the triangles intersecting with the current grid cell (row) + KOKKOS_INLINE_FUNCTION + LO operator()(LO row, LO* fill) const + { + const 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) { + const auto elem_tet2verts = Omega_h::gather_verts<4>(tets2verts_, elem_idx); + // 2d mesh with 2d coords, but 3 triangles + const auto vertex_coords = Omega_h::gather_vectors<4, 3>(coords_, elem_tet2verts); + if (simplex_intersects_bbox<3>(vertex_coords, grid_cell_bbox)) { + if (fill) { + fill[num_intersections] = elem_idx; + } + ++num_intersections; + } + } + return num_intersections; + } + +private: + Omega_h::Mesh& mesh_; + Omega_h::LOs tets2verts_; + Omega_h::Reals coords_; + Kokkos::View grid_; +public: + LO nelems_; +}; + // num_grid_cells should be result of grid.GetNumCells(), take as argument to avoid extra copy // of grid from gpu to cpu Kokkos::Crs construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { Kokkos::Crs intersection_map{}; - auto f = detail::GridTriIntersectionFunctor{mesh, grid}; + auto f = detail::GridTriIntersectionFunctor2D{mesh, grid}; Kokkos::count_and_fill_crs(intersection_map, num_grid_cells, f); return intersection_map; } @@ -260,7 +342,7 @@ OMEGA_H_INLINE double myreduce(const Omega_h::Vector & x, Op op) OMEGA_H_NOEX Kokkos::View GridPointSearch2D::operator()(Kokkos::View points) const { - Kokkos::View results("point search result", points.extent(0)); + Kokkos::View results("point search result", points.extent(0)); auto num_rows = candidate_map_.numRows(); // needed so that we don't capture this ptr which will be memory error on cuda auto grid = grid_; @@ -279,7 +361,7 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V bool found = false; auto nearest_triangle = candidates_begin; - auto dimensionality = GridPointSearch::Result::Dimensionality::EDGE; + auto dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; 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 @@ -292,7 +374,7 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V auto parametric_coords = barycentric_from_global(point, vertex_coords); if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { - results(p) = GridPointSearch::Result{GridPointSearch::Result::Dimensionality::FACE, triangleID, parametric_coords}; + results(p) = GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, triangleID, parametric_coords}; found = true; break; } @@ -321,7 +403,7 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V if (distance_to_ab >= distance_to_nearest) { continue; } - dimensionality = GridPointSearch::Result::Dimensionality::EDGE; + dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; nearest_triangle = i; distance_to_nearest = distance_to_ab; parametric_coords_to_nearest = parametric_coords; @@ -336,7 +418,7 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V Omega_h::get_vector<2>(coords, vertexID); if (const auto distance = Omega_h::norm(point - vertex);distance < distance_to_nearest) { - dimensionality = GridPointSearch::Result::Dimensionality::VERTEX; + dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; nearest_triangle = i; distance_to_nearest = distance; parametric_coords_to_nearest = parametric_coords; @@ -345,14 +427,14 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V } if(!found) { - results(p) = GridPointSearch::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; + results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; } }); return results; } -GridPointSearch<2>::GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny) +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index ba785ca2..f6dd5d9e 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -43,28 +43,21 @@ class PointLocalizationSearch Omega_h::Vector parametric_coords; }; - static constexpr auto DIM = 2; + static constexpr auto DIM = dim; virtual Kokkos::View operator()(Kokkos::View point) const = 0; }; using PointLocalizationSearch2D = PointLocalizationSearch<2>; -template -class GridPointSearch : public PointLocalizationSearch -{ - static_assert(false, "Not implemented"); -}; - -template <> -class GridPointSearch<2> : public PointLocalizationSearch2D +class GridPointSearch2D : public PointLocalizationSearch2D { using CandidateMapT = Kokkos::Crs; public: using Result = PointLocalizationSearch2D::Result; - GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny); + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny); /** * given a point in global coordinates give the id of the triangle that the * point lies within and the parametric coordinate of the point within the @@ -85,7 +78,5 @@ class GridPointSearch<2> : public PointLocalizationSearch2D Omega_h::Reals coords_; }; -using GridPointSearch2D = GridPointSearch<2>; - } // namespace detail #endif // PCMS_COUPLING_POINT_SEARCH_H diff --git a/src/pcms/uniform_grid.h b/src/pcms/uniform_grid.h index 5ec13d19..52680254 100644 --- a/src/pcms/uniform_grid.h +++ b/src/pcms/uniform_grid.h @@ -6,7 +6,7 @@ namespace pcms { -template +template struct UniformGrid { // Make private? @@ -117,7 +117,8 @@ struct UniformGrid } }; -using Uniform2DGrid = UniformGrid<>; +using Uniform2DGrid = UniformGrid<2>; +using Uniform3DGrid = UniformGrid<3>; } // namespace pcms From 8f1dfa96d6d4549fa314d2871cd0f1b6599a25f7 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 10 Apr 2025 22:26:19 -0400 Subject: [PATCH 04/22] Replace custom barycentric_from_global with Omega_h implementation Replaced the local implementation of barycentric_from_global with Omega_h's standard implementation. Removed the custom function from both source and header files, and updated all calls to use Omega_h::barycentric_from_global<2, 2>. This change enhances code maintainability and aligns functionality with the Omega_h library. --- src/pcms/point_search.cpp | 22 +++++----------------- src/pcms/point_search.h | 2 -- test/test_point_search.cpp | 5 ++--- 3 files changed, 7 insertions(+), 22 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index c38aaf84..3188ae04 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -152,13 +152,13 @@ bool bbox_verts_within_triangle(const AABBox<2>& bbox, const Omega_h::Matrix<2,3 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 = barycentric_from_global({left,bot}, coords); + auto xi = Omega_h::barycentric_from_global<2, 2>({left,bot}, coords); if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} - xi = barycentric_from_global({left,top}, coords); + xi = Omega_h::barycentric_from_global<2, 2>({left,top}, coords); if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} - xi = barycentric_from_global({right,top}, coords); + xi = Omega_h::barycentric_from_global<2, 2>({right,top}, coords); if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} - xi = barycentric_from_global({right,bot}, coords); + xi = Omega_h::barycentric_from_global<2, 2>({right,bot}, coords); if(Omega_h::is_barycentric_inside(xi, fuzz)){ return true;} return false; } @@ -321,18 +321,6 @@ construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View g } } // namespace detail -KOKKOS_FUNCTION -Omega_h::Vector<3> barycentric_from_global( - const Omega_h::Vector<2>& point, const Omega_h::Matrix<2, 3>& vertex_coords) -{ - const auto inverse_basis = - Omega_h::pseudo_invert(Omega_h::simplex_basis<2, 2>(vertex_coords)); - auto xi = inverse_basis * (point - vertex_coords[0]); - // note omega_h form_barycentric is currently broken. - // see https://github.com/sandialabs/omega_h/issues/389 - return {1 - xi[0] - xi[1], xi[0], xi[1]}; -} - template OMEGA_H_INLINE double myreduce(const Omega_h::Vector & x, Op op) OMEGA_H_NOEXCEPT { auto out = x[0]; @@ -371,7 +359,7 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V 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 = barycentric_from_global(point, vertex_coords); + auto parametric_coords = Omega_h::barycentric_from_global<2, 2>(point, vertex_coords); if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { results(p) = GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, triangleID, parametric_coords}; diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index f6dd5d9e..4659bd79 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -20,8 +20,6 @@ Kokkos::Crs construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); } KOKKOS_FUNCTION -Omega_h::Vector<3> barycentric_from_global( - const Omega_h::Vector<2>& point, const Omega_h::Matrix<2, 3>& vertex_coords); [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 7590efe3..6652e71a 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -5,7 +5,6 @@ #include using pcms::AABBox; -using pcms::barycentric_from_global; using pcms::Uniform2DGrid; TEST_CASE("global to local") @@ -14,7 +13,7 @@ TEST_CASE("global to local") SECTION("check verts") { for (int i = 0; i < coords.size(); ++i) { - auto xi = barycentric_from_global(coords[i], coords); + auto xi = Omega_h::barycentric_from_global<2, 2>(coords[i], coords); Omega_h::Vector<3> hand_xi{0, 0, 0}; hand_xi[i] = 1; printf("[%f,%f,%f] == [%f,%f,%f]\n", xi[0], xi[1], xi[2], hand_xi[0], @@ -25,7 +24,7 @@ TEST_CASE("global to local") SECTION("check point") { Omega_h::Vector<2> point{0.5, 0.5}; - auto xi = barycentric_from_global(point, coords); + auto xi = Omega_h::barycentric_from_global<2, 2>(point, coords); Omega_h::Vector<3> hand_xi{0.25, 0.25, 0.5}; printf("[%f,%f,%f] == [%f,%f,%f]\n", xi[0], xi[1], xi[2], hand_xi[0], hand_xi[1], hand_xi[2]); From 19fd48254d3855df0259074cea2ed49923249585 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 11 Apr 2025 12:01:38 -0400 Subject: [PATCH 05/22] Add bbox_verts_within_simplex function and virtual destructor Introduce bbox_verts_within_simplex to check if an AABBox's vertices are within a simplex using barycentric coordinates. Also, add a virtual destructor to the PointLocalizationSearch class for proper cleanup in derived classes. --- src/pcms/point_search.cpp | 32 ++++++++++++++++++++++++++++++++ src/pcms/point_search.h | 1 + 2 files changed, 33 insertions(+) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 3188ae04..63b750b0 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -163,6 +163,38 @@ bool bbox_verts_within_triangle(const AABBox<2>& bbox, const Omega_h::Matrix<2,3 return false; } +template +[[nodiscard]] KOKKOS_INLINE_FUNCTION +bool bbox_verts_within_simplex(const AABBox& bbox, const Omega_h::Matrix& coords) +{ + // each dimension has a pair of opposing "walls" + // 2D: { [left, right], [top, bottom] } -> { left, right, top, bottom } + // 3D: { [left, right], [top, bottom], [front, back] } -> { left, ..., back } + std::array bbox_walls {}; + for (int i = 0; i < dim; i++) { + bbox_walls[i * 2] = bbox.center[i] - bbox.half_width[i]; + bbox_walls[i * 2 + 1] = bbox.center[i] + bbox.half_width[i]; + } + + // 1 << dim == 2 ** dim == num vertices in ndim bounding box / hypercube + // Each vertex is a just a unique combination of walls + // eg [left, bottom] (2D) or [right, top, front] (3) + for (unsigned i = 0; i < 1 << dim; ++i) { + // conveniently, i acts a bit field representing the current combination + Omega_h::Vector vert; + for (unsigned j = 0; j < dim; ++j) { + // eg 110 = 6 = [left, top, back] + // eg 01 = 1 = [right, top] + vert[j] = (i >> j) & 1 ? bbox_walls[j * 2] : bbox_walls[j * 2 + 1]; + } + auto xi = Omega_h::barycentric_from_global(vert, coords); + if (Omega_h::is_barycentric_inside(xi, fuzz)) { + return true; + } + } + return false; +} + /** * Check if a triangle element represented by 3 coordinates in two dimensions * intersects with a bounding box diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 4659bd79..23c58f8c 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -44,6 +44,7 @@ class PointLocalizationSearch static constexpr auto DIM = dim; virtual Kokkos::View operator()(Kokkos::View point) const = 0; + virtual ~PointLocalizationSearch() = default; }; using PointLocalizationSearch2D = PointLocalizationSearch<2>; From 1ae14727a69df901e521228144ef02c607d6ce6c Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 11 Apr 2025 13:05:09 -0400 Subject: [PATCH 06/22] Refactor point search to support 3D grids. Renamed `construct_intersection_map` to `construct_intersection_map_2d` and added `construct_intersection_map_3d` for 3D support. Introduced `GridPointSearch3D` class, extending point localization for 3D grids with associated functionality and data structures. Updated tests to reflect these changes. --- src/pcms/point_search.cpp | 96 ++++++++++++++++++++++++++++++++++++-- src/pcms/point_search.h | 31 +++++++++++- test/test_point_search.cpp | 4 +- 3 files changed, 125 insertions(+), 6 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 63b750b0..f029fee6 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -320,7 +320,6 @@ struct GridTriIntersectionFunctor3D // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { const auto elem_tet2verts = Omega_h::gather_verts<4>(tets2verts_, elem_idx); - // 2d mesh with 2d coords, but 3 triangles const auto vertex_coords = Omega_h::gather_vectors<4, 3>(coords_, elem_tet2verts); if (simplex_intersects_bbox<3>(vertex_coords, grid_cell_bbox)) { if (fill) { @@ -344,13 +343,22 @@ struct GridTriIntersectionFunctor3D // num_grid_cells should be result of grid.GetNumCells(), take as argument to avoid extra copy // of grid from gpu to cpu Kokkos::Crs -construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) +construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) { Kokkos::Crs intersection_map{}; auto f = detail::GridTriIntersectionFunctor2D{mesh, grid}; Kokkos::count_and_fill_crs(intersection_map, num_grid_cells, f); return intersection_map; } + +Kokkos::Crs +construct_intersection_map_3d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) +{ + Kokkos::Crs intersection_map{}; + auto f = detail::GridTriIntersectionFunctor3D{mesh, grid}; + Kokkos::count_and_fill_crs(intersection_map, num_grid_cells, f); + return intersection_map; +} } // namespace detail template @@ -463,7 +471,89 @@ 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(mesh, grid_, grid_h(0).GetNumCells()); + 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); +} + +Kokkos::View GridPointSearch3D::operator()(Kokkos::View points) const +{ + Kokkos::View results("point search result", points.extent(0)); + auto num_rows = candidate_map_.numRows(); + // needed so that we don't capture this ptr which will be memory error on cuda + auto grid = grid_; + auto candidate_map = candidate_map_; + auto tris2verts = tris2verts_; + auto tris2verts_adj = tris2verts_adj_; + auto tris2edges_adj = tris2edges_adj_; + auto edges2verts_adj = edges2verts_adj_; + auto coords = coords_; + Kokkos::parallel_for(points.extent(0), KOKKOS_LAMBDA(int p) { + + Omega_h::Vector point; + for (int i = 0; i < DIM; ++i) { + point[i] = points(p, i); + } + + auto cell_id = grid(0).ClosestCellID(point); + assert(cell_id < num_rows && cell_id >= 0); + auto candidates_begin = candidate_map.row_map(cell_id); + auto candidates_end = candidate_map.row_map(cell_id + 1); + bool found = false; + + auto nearest_triangle = candidates_begin; + auto dimensionality = GridPointSearch3D::Result::Dimensionality::EDGE; + Omega_h::Real distance_to_nearest { INFINITY }; + Omega_h::Vector 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); + const auto elem_tri2verts = Omega_h::gather_verts(tris2verts, triangleID); + auto vertex_coords = Omega_h::gather_vectors(coords, elem_tri2verts); + auto parametric_coords = Omega_h::barycentric_from_global(point, vertex_coords); + + if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { + results(p) = GridPointSearch3D::Result{GridPointSearch3D::Result::Dimensionality::FACE, triangleID, parametric_coords}; + found = true; + break; + } + + // TODO: Get nearest element if no triangle found + } + if(!found) + { + results(p) = GridPointSearch3D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; + } + }); + + return results; +} + +GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) +{ + auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); + auto grid_h = Kokkos::create_mirror_view(grid_); + + std::array edge_lengths {}; + std::array bot_left {}; + + for (int i = 0; i < DIM; ++i) { + edge_lengths[i] = mesh_bbox.max[i] - mesh_bbox.min[i]; + bot_left[i] = mesh_bbox.min[i]; + } + + grid_h(0) = Uniform3DGrid { + .edge_length = edge_lengths, + .bot_left = bot_left, + .divisions = {Nx, Ny, Nz} + }; + + Kokkos::deep_copy(grid_, grid_h); + candidate_map_ = detail::construct_intersection_map_3d(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); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 23c58f8c..e046bf38 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -17,7 +17,7 @@ namespace pcms // used namespace detail { Kokkos::Crs -construct_intersection_map(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); +construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); } KOKKOS_FUNCTION @@ -48,6 +48,7 @@ class PointLocalizationSearch }; using PointLocalizationSearch2D = PointLocalizationSearch<2>; +using PointLocalizationSearch3D = PointLocalizationSearch<3>; class GridPointSearch2D : public PointLocalizationSearch2D { @@ -77,5 +78,33 @@ class GridPointSearch2D : public PointLocalizationSearch2D Omega_h::Reals coords_; }; +class GridPointSearch3D : public PointLocalizationSearch3D +{ + using CandidateMapT = Kokkos::Crs; + +public: + using Result = PointLocalizationSearch3D::Result; + + GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz); + /** + * given a point in global coordinates give the id of the triangle that the + * point lies within and the parametric coordinate of the point within the + * triangle. If the point does not lie within any triangle element. Then the + * id will be a negative number and (TODO) will return a negative id of the + * closest element + */ + Kokkos::View operator()(Kokkos::View point) const override; + +private: + Omega_h::Mesh mesh_; + Omega_h::Adj tris2edges_adj_; + Omega_h::Adj tris2verts_adj_; + Omega_h::Adj edges2verts_adj_; + Kokkos::View[1]> grid_{"uniform grid"}; + CandidateMapT candidate_map_; + Omega_h::LOs tris2verts_; + Omega_h::Reals coords_; +}; + } // namespace detail #endif // PCMS_COUPLING_POINT_SEARCH_H diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 6652e71a..20b16eb1 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -139,7 +139,7 @@ TEST_CASE("construct intersection map") auto grid_h = Kokkos::create_mirror_view(grid_d); grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; Kokkos::deep_copy(grid_d, grid_h); - auto intersection_map = pcms::detail::construct_intersection_map(mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map_2d(mesh, grid_d, grid_h(0).GetNumCells()); // assert(cudaSuccess == cudaDeviceSynchronize()); REQUIRE(intersection_map.numRows() == 100); REQUIRE(num_candidates_within_range(intersection_map, 2, 16)); @@ -151,7 +151,7 @@ TEST_CASE("construct intersection map") grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {60, 60}}; Kokkos::deep_copy(grid_d, grid_h); // require number of candidates is >=1 and <=6 - auto intersection_map = pcms::detail::construct_intersection_map(mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map_2d(mesh, grid_d, grid_h(0).GetNumCells()); REQUIRE(intersection_map.numRows() == 3600); REQUIRE(num_candidates_within_range(intersection_map, 1, 6)); From ab88b05c2bcf36849e518a8f4146a77c489c0bf4 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 18 Apr 2025 12:23:54 -0400 Subject: [PATCH 07/22] Add REGION dimensionality to point search Introduce a new REGION dimensionality to represent 3D regions in point search operations. Adjust loops and result assignment logic to incorporate this new dimension for improved handling of 3D cases. --- src/pcms/point_search.cpp | 6 ++++-- src/pcms/point_search.h | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index f029fee6..cf6cb00b 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -177,9 +177,11 @@ bool bbox_verts_within_simplex(const AABBox& bbox, const Omega_h::Matrix vert; for (unsigned j = 0; j < dim; ++j) { @@ -517,7 +519,7 @@ Kokkos::View GridPointSearch3D::operator()(Kokkos::V auto parametric_coords = Omega_h::barycentric_from_global(point, vertex_coords); if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { - results(p) = GridPointSearch3D::Result{GridPointSearch3D::Result::Dimensionality::FACE, triangleID, parametric_coords}; + results(p) = GridPointSearch3D::Result{GridPointSearch3D::Result::Dimensionality::REGION, triangleID, parametric_coords}; found = true; break; } diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index e046bf38..fac0a50d 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -33,7 +33,8 @@ class PointLocalizationSearch { VERTEX = 0, EDGE = 1, - FACE = 2 + FACE = 2, + REGION = 3 }; Dimensionality dimensionality; From ca5caa02a0b9aa9457bb063a17e9f9164e6ee4af Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 8 Oct 2025 12:17:33 -0400 Subject: [PATCH 08/22] Refactor edge and vertex search loops in 2D grid point localization Simplified indentation and adjusted line breaks in edge and vertex search loops for improved readability. Added constructor overloads to `GridPointSearch2D` for handling tolerances via a new `PointSearchTolerances` abstraction. Updated `PointLocalizationSearch` class with support for tolerances. --- src/pcms/point_search.cpp | 61 ++++++++++++++++++++++----------------- src/pcms/point_search.h | 9 ++++++ 2 files changed, 44 insertions(+), 26 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index cf6cb00b..6cce3627 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -409,34 +409,35 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V break; } - for (int j = 0; j < 3; ++j) { - // Every triangle (face) is connected to 3 edges - const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; + 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_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); + 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(point, vertex_a, vertex_b)) continue; + if (!normal_intersects_segment(point, vertex_a, vertex_b)) + continue; - const auto xa = vertex_a[0]; - const auto ya = vertex_a[1]; - const auto xb = vertex_b[0]; - const auto yb = vertex_b[1]; + const auto xa = vertex_a[0]; + const auto ya = vertex_a[1]; + const auto xb = vertex_b[0]; + const auto yb = vertex_b[1]; - const auto xp = point[0]; - const auto yp = point[1]; + const auto xp = point[0]; + const auto yp = point[1]; - const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); + const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); - if (distance_to_ab >= distance_to_nearest) { continue; } + if (distance_to_ab >= distance_to_nearest) continue; - dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; - nearest_triangle = i; - distance_to_nearest = distance_to_ab; - parametric_coords_to_nearest = parametric_coords; + dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; + nearest_triangle = i; + distance_to_nearest = distance_to_ab; + parametric_coords_to_nearest = parametric_coords; } // Every triangle (face) is connected to 3 vertices @@ -447,12 +448,13 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V const Omega_h::Few vertex = Omega_h::get_vector<2>(coords, vertexID); - if (const auto distance = Omega_h::norm(point - vertex);distance < distance_to_nearest) { - dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; - nearest_triangle = i; - distance_to_nearest = distance; - parametric_coords_to_nearest = parametric_coords; - } + if (const auto distance = Omega_h::norm(point - vertex); + distance < distance_to_nearest) { + dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; + nearest_triangle = i; + distance_to_nearest = distance; + parametric_coords_to_nearest = parametric_coords; + } } } if(!found) @@ -465,6 +467,13 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V } GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) + : GridPointSearch2D(mesh, Nx, Ny, { }) +{ + Kokkos::deep_copy(tolerances_, 0); +} + +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, PointSearchTolerances tolerances) + : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index fac0a50d..71af9a97 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -44,8 +44,15 @@ class PointLocalizationSearch static constexpr auto DIM = dim; + using PointSearchTolerances = Kokkos::View; + + explicit PointLocalizationSearch(PointSearchTolerances tolerances) : tolerances_(tolerances){ } + virtual Kokkos::View operator()(Kokkos::View point) const = 0; virtual ~PointLocalizationSearch() = default; + +protected: + PointSearchTolerances tolerances_; }; using PointLocalizationSearch2D = PointLocalizationSearch<2>; @@ -59,6 +66,8 @@ class GridPointSearch2D : public PointLocalizationSearch2D using Result = PointLocalizationSearch2D::Result; GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny); + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, PointSearchTolerances tolerances); + /** * given a point in global coordinates give the id of the triangle that the * point lies within and the parametric coordinate of the point within the From 26bd542dc254d3f3fc5509d44edb745e5cbc6049 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 8 Oct 2025 12:42:15 -0400 Subject: [PATCH 09/22] Reversed algorithm logic and implemented tolerances Previous algorithm checked edges then vertices. Now we check vertices, then edges. --- src/pcms/point_search.cpp | 55 +++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 23 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 6cce3627..17790d68 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -409,6 +409,32 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V break; } + // Every triangle (face) is connected to 3 vertices + 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 > tolerances_[0]) continue; + + if (distance < distance_to_nearest) { + dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; + nearest_triangle = i; + distance_to_nearest = distance; + parametric_coords_to_nearest = parametric_coords; + } + } + + if (dimensionality == Result::Dimensionality::VERTEX) { + results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; + continue; + } + + dimensionality = Result::Dimensionality::FACE; + for (int j = 0; j < 3; ++j) { // Every triangle (face) is connected to 3 edges const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; @@ -432,35 +458,18 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); - if (distance_to_ab >= distance_to_nearest) continue; - - dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; - nearest_triangle = i; - distance_to_nearest = distance_to_ab; - parametric_coords_to_nearest = parametric_coords; - } + if (distance_to_ab > tolerances_[1]) continue; - // Every triangle (face) is connected to 3 vertices - 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); - - if (const auto distance = Omega_h::norm(point - vertex); - distance < distance_to_nearest) { - dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; + if (distance_to_ab < distance_to_nearest) { + dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; nearest_triangle = i; - distance_to_nearest = distance; + distance_to_nearest = distance_to_ab; parametric_coords_to_nearest = parametric_coords; } } } - if(!found) - { - results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; - } + + results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; }); return results; From dc655c44d06c3385d8de44b9d72de4ee0cee06d4 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 30 Oct 2025 15:38:29 -0400 Subject: [PATCH 10/22] Refactor point search tolerances handling and rename Result fields Converted `PointSearchTolerances` to use `const` references for safer and more efficient handling. Added default tolerances initialization in 2D and 3D grid search constructors. Renamed `tri_id` to `element_id` in `Result` for improved clarity and consistency across dimensionalities. --- src/pcms/interpolator/adj_search.hpp | 6 +++--- src/pcms/point_search.cpp | 12 ++++++++++-- src/pcms/point_search.h | 17 ++++++++++++----- test/test_point_search.cpp | 8 ++++---- 4 files changed, 29 insertions(+), 14 deletions(-) diff --git a/src/pcms/interpolator/adj_search.hpp b/src/pcms/interpolator/adj_search.hpp index 7f5cc2e2..040d96a9 100644 --- a/src/pcms/interpolator/adj_search.hpp +++ b/src/pcms/interpolator/adj_search.hpp @@ -30,8 +30,8 @@ inline void checkTargetPoints( printf("INFO: Checking target points...\n"); auto check_target_points = OMEGA_H_LAMBDA(Omega_h::LO i) { - if (results(i).tri_id < 0) { - OMEGA_H_CHECK_PRINTF(results(i).tri_id >= 0, + 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", i); printf("%d, ", i); @@ -128,7 +128,7 @@ inline void FindSupports::adjBasedSearch( Track visited; Omega_h::Real cutoffDistance = radii2[id]; - Omega_h::LO source_cell_id = results(id).tri_id; + Omega_h::LO source_cell_id = results(id).element_id; OMEGA_H_CHECK_PRINTF( source_cell_id >= 0, "ERROR: Source cell id not found for target %d (%f,%f)\n", id, diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 17790d68..f75f38f0 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -476,12 +476,12 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V } GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) - : GridPointSearch2D(mesh, Nx, Ny, { }) + : GridPointSearch2D(mesh, Nx, Ny, PointSearchTolerances { "point search 2d tolerances" }) { Kokkos::deep_copy(tolerances_, 0); } -GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, PointSearchTolerances tolerances) +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); @@ -554,6 +554,14 @@ Kokkos::View GridPointSearch3D::operator()(Kokkos::V } GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) + : GridPointSearch3D(mesh, Nx, Ny, Nz, 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) + : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 71af9a97..0910b771 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -1,11 +1,13 @@ #ifndef PCMS_COUPLING_POINT_SEARCH_H #define PCMS_COUPLING_POINT_SEARCH_H -#include +#include + #include #include -#include "types.h" #include #include + +#include "types.h" #include "pcms/uniform_grid.h" #include "pcms/bounding_box.h" @@ -38,7 +40,7 @@ class PointLocalizationSearch }; Dimensionality dimensionality; - LO tri_id; + LO element_id; Omega_h::Vector parametric_coords; }; @@ -46,7 +48,10 @@ class PointLocalizationSearch using PointSearchTolerances = Kokkos::View; - explicit PointLocalizationSearch(PointSearchTolerances tolerances) : tolerances_(tolerances){ } + explicit PointLocalizationSearch(const PointSearchTolerances& tolerances) : tolerances_(tolerances) + { + assert(tolerances_.is_allocated()); + } virtual Kokkos::View operator()(Kokkos::View point) const = 0; virtual ~PointLocalizationSearch() = default; @@ -66,7 +71,7 @@ class GridPointSearch2D : public PointLocalizationSearch2D using Result = PointLocalizationSearch2D::Result; GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny); - GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, PointSearchTolerances tolerances); + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances); /** * given a point in global coordinates give the id of the triangle that the @@ -96,6 +101,8 @@ class GridPointSearch3D : public PointLocalizationSearch3D 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, const PointSearchTolerances& tolerances); + /** * given a point in global coordinates give the id of the triangle that the * point lies within and the parametric coordinate of the point within the diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 20b16eb1..54d1bb45 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -209,21 +209,21 @@ TEST_CASE("uniform grid search") { auto out_of_bounds = results_h(2); auto top_right = results_h(3); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); - REQUIRE(-1*out_of_bounds.tri_id == top_right.tri_id); + REQUIRE(-1*out_of_bounds.element_id == top_right.element_id); 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.tri_id == bot_left.tri_id); + REQUIRE(-1*out_of_bounds.element_id == bot_left.element_id); out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.tri_id == top_right.tri_id); + REQUIRE(-1 * out_of_bounds.element_id == top_right.element_id); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.tri_id == bot_left.tri_id); + REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); } } From 18e3fab4081cb8ece511ad7cb73c2a49435a4513 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Thu, 30 Oct 2025 15:39:21 -0400 Subject: [PATCH 11/22] Refactor GridPointSearch2D result assignment logic and update tolerances handling Revised result assignment logic to incorporate detailed checks for vertex, edge, and face proximity, using new tolerances-based conditions. Updated tests to reflect changes in dimensionality and result handling for improved accuracy verification. --- src/pcms/point_search.cpp | 56 ++++++++++++++++++++++++-------------- test/test_point_search.cpp | 14 ++++++++-- 2 files changed, 46 insertions(+), 24 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index f75f38f0..79799421 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -382,16 +382,20 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V auto tris2edges_adj = tris2edges_adj_; auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; + auto tolerances = tolerances_; Kokkos::parallel_for(points.extent(0), KOKKOS_LAMBDA(int p) { Omega_h::Vector<2> point(std::initializer_list{points(p,0), points(p,1)}); auto cell_id = grid(0).ClosestCellID(point); assert(cell_id < num_rows && cell_id >= 0); auto candidates_begin = candidate_map.row_map(cell_id); auto candidates_end = candidate_map.row_map(cell_id + 1); - bool found = false; - auto nearest_triangle = candidates_begin; - auto dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; + 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 @@ -403,12 +407,6 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V 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); - if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { - results(p) = GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, triangleID, parametric_coords}; - found = true; - break; - } - // Every triangle (face) is connected to 3 vertices for (int j = 0; j < 3; ++j) { // Get the vertex ID from the connectivity array @@ -418,22 +416,21 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V Omega_h::get_vector<2>(coords, vertexID); const auto distance = Omega_h::norm(point - vertex); - if (distance > tolerances_[0]) continue; if (distance < distance_to_nearest) { dimensionality = GridPointSearch2D::Result::Dimensionality::VERTEX; - nearest_triangle = i; + nearest_element_id = vertexID; distance_to_nearest = distance; parametric_coords_to_nearest = parametric_coords; - } - } - if (dimensionality == Result::Dimensionality::VERTEX) { - results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; - continue; + if (distance < tolerances(0)) { + vertex_found = true; + }; + } } - dimensionality = Result::Dimensionality::FACE; + if (vertex_found) + break; for (int j = 0; j < 3; ++j) { // Every triangle (face) is connected to 3 edges @@ -458,18 +455,35 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); - if (distance_to_ab > tolerances_[1]) continue; - if (distance_to_ab < distance_to_nearest) { + edge_found = true; dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; - nearest_triangle = i; + nearest_element_id = edgeID; distance_to_nearest = distance_to_ab; parametric_coords_to_nearest = parametric_coords; + + if (distance_to_ab < tolerances(1)) { + edge_found = true; + }; } } + + if (edge_found) + break; + + 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; + + // results(p) = GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, triangleID, parametric_coords}; + // return; + } } - results(p) = GridPointSearch2D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; + 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}; }); return results; diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 54d1bb45..80814ca4 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -163,7 +163,12 @@ TEST_CASE("uniform grid search") { auto world = lib.world(); auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); - GridPointSearch2D search{mesh,10,10}; + auto tolerances = GridPointSearch2D::PointSearchTolerances { "point search 2d tolerances" }; + tolerances(0) = 0.01; + tolerances(1) = 0.01; + + GridPointSearch2D search{mesh,10,10, tolerances}; + Kokkos::View points("test_points", 7); //Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); @@ -189,7 +194,10 @@ TEST_CASE("uniform grid search") { { { auto [dim, idx,coords] = results_h(0); - REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); + + CAPTURE(idx); + + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::VERTEX); REQUIRE(idx == 0); REQUIRE(coords[0] == Catch::Approx(1)); REQUIRE(coords[1] == Catch::Approx(0)); @@ -197,7 +205,7 @@ TEST_CASE("uniform grid search") { } { auto [dim, idx,coords] = results_h(1); - REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(idx == 91); REQUIRE(coords[0] == Catch::Approx(0.5)); REQUIRE(coords[1] == Catch::Approx(0.1)); From bccb55e557d3c30cadb93fbe495b2a52b39f2cba Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 5 Nov 2025 14:54:50 -0500 Subject: [PATCH 12/22] Refactor edge handling in GridPointSearch2D and extend test cases Removed unused `edge_found` variable for cleaner logic in result assignment. Updated test cases to add new points and verify dimensionality, element IDs, and coordinates for various scenarios, enhancing accuracy and coverage. --- src/pcms/point_search.cpp | 1 - test/test_point_search.cpp | 13 ++++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 79799421..a54aca44 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -456,7 +456,6 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); if (distance_to_ab < distance_to_nearest) { - edge_found = true; dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; nearest_element_id = edgeID; distance_to_nearest = distance_to_ab; diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 80814ca4..22f4fcff 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -169,7 +169,7 @@ TEST_CASE("uniform grid search") { GridPointSearch2D search{mesh,10,10, tolerances}; - Kokkos::View points("test_points", 7); + Kokkos::View points("test_points", 8); //Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); points_h(0,0) = 0; @@ -186,6 +186,8 @@ TEST_CASE("uniform grid search") { points_h(5, 1) = 0.95; points_h(6, 0) = 0.05; points_h(6, 1) = -0.01; + points_h(7, 0) = 0.05; + points_h(7, 1) = 0.02; Kokkos::deep_copy(points, points_h); auto results = search(points); auto results_h = Kokkos::create_mirror_view(results); @@ -206,11 +208,16 @@ TEST_CASE("uniform grid search") { { auto [dim, idx,coords] = results_h(1); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(idx == 91); + REQUIRE(idx == 156); 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); + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); + REQUIRE(idx == 0); + } } // feature needs to be added SECTION("Global coordinate outside mesh", "[!mayfail]") { @@ -227,7 +234,7 @@ TEST_CASE("uniform grid search") { out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.element_id == top_right.element_id); + REQUIRE(out_of_bounds.element_id == 219); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == From c1f9e4f92b95b9bee5959779db7f43a4bc926148 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 5 Nov 2025 15:44:04 -0500 Subject: [PATCH 13/22] Ran clang format --- src/pcms/point_search.cpp | 431 +++++++++++++++++++++---------------- src/pcms/point_search.h | 36 ++-- test/test_point_search.cpp | 70 +++--- 3 files changed, 307 insertions(+), 230 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index a54aca44..d484ba4f 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -2,29 +2,35 @@ #include #include -// From https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points +// From +// https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points KOKKOS_INLINE_FUNCTION -double distance_from_line(const double x0, const double y0, const double x1, const double y1, const double x2, const double y2) +double distance_from_line(const double x0, const double y0, const double x1, + const double y1, const double x2, const double y2) { - const Omega_h::Vector<2> p1 = { x1, y1 }; - const Omega_h::Vector<2> p2 = { x2, y2 }; + const Omega_h::Vector<2> p1 = {x1, y1}; + const Omega_h::Vector<2> p2 = {x2, y2}; auto disp = p2 - p1; - return std::abs(disp[1]*x0 - disp[0]*y0 + x2*y1 - y2*x1) / Omega_h::norm(disp); + return std::abs(disp[1] * x0 - disp[0] * y0 + x2 * y1 - y2 * x1) / + Omega_h::norm(disp); } -// Law of Cosines, where a, b, c and gamma are defined here: https://en.wikipedia.org/wiki/Law_of_cosines#Use_in_solving_triangles +// Law of Cosines, where a, b, c and gamma are defined here: +// https://en.wikipedia.org/wiki/Law_of_cosines#Use_in_solving_triangles KOKKOS_INLINE_FUNCTION double angle_from_side_lengths(const double a, const double b, const double c) { - return std::acos((a*a + b*b - c*c) / 2*a*b); + return std::acos((a * a + b * b - c * c) / 2 * a * b); } KOKKOS_INLINE_FUNCTION -bool normal_intersects_segment(const Omega_h::Few a, const Omega_h::Few b, const Omega_h::Few c) +bool normal_intersects_segment(const Omega_h::Few a, + const Omega_h::Few b, + const Omega_h::Few c) { const auto ab_len = Omega_h::norm(a - b); - const auto bc_len = Omega_h::norm(b -c); + const auto bc_len = Omega_h::norm(b - c); const auto ac_len = Omega_h::norm(a - c); const double angle1 = angle_from_side_lengths(bc_len, ac_len, ab_len); @@ -137,40 +143,52 @@ bool line_intersects_bbox(const Omega_h::Vector<2>& p0, } template -[[nodiscard]] KOKKOS_INLINE_FUNCTION -bool within_bbox(const Omega_h::Vector coord, const AABBox & bbox) noexcept { +[[nodiscard]] KOKKOS_INLINE_FUNCTION bool within_bbox( + const Omega_h::Vector coord, const AABBox& bbox) noexcept +{ for (int i = 0; i < dim; ++i) { - if (coord[i] < bbox.center[i] - bbox.half_width[i]) return false; - if (coord[i] > bbox.center[i] + bbox.half_width[i]) return false; + if (coord[i] < bbox.center[i] - bbox.half_width[i]) + return false; + if (coord[i] > bbox.center[i] + bbox.half_width[i]) + return false; } return true; } -[[nodiscard]] KOKKOS_INLINE_FUNCTION -bool bbox_verts_within_triangle(const AABBox<2>& bbox, const Omega_h::Matrix<2,3>& coords) { +[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_triangle( + const AABBox<2>& bbox, const Omega_h::Matrix<2, 3>& coords) +{ auto left = bbox.center[0] - bbox.half_width[0]; auto right = bbox.center[0] + bbox.half_width[0]; auto bot = bbox.center[1] - bbox.half_width[1]; auto top = bbox.center[1] + bbox.half_width[1]; - auto xi = Omega_h::barycentric_from_global<2, 2>({left,bot}, coords); - if(Omega_h::is_barycentric_inside(xi, 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;} + 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) +[[nodiscard]] KOKKOS_INLINE_FUNCTION bool bbox_verts_within_simplex( + const AABBox& bbox, const Omega_h::Matrix& coords) { // each dimension has a pair of opposing "walls" // 2D: { [left, right], [top, bottom] } -> { left, right, top, bottom } // 3D: { [left, right], [top, bottom], [front, back] } -> { left, ..., back } - std::array bbox_walls {}; + std::array bbox_walls{}; for (int i = 0; i < dim; i++) { bbox_walls[i * 2] = bbox.center[i] - bbox.half_width[i]; bbox_walls[i * 2 + 1] = bbox.center[i] + bbox.half_width[i]; @@ -202,15 +220,13 @@ bool bbox_verts_within_simplex(const AABBox& bbox, const Omega_h::Matrix& coords, - const AABBox<2>& bbox) +KOKKOS_FUNCTION bool triangle_intersects_bbox( + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox) { // triangle and grid cell bounding box intersect if (intersects(triangle_bbox(coords), bbox)) { // if any of the triangle verts inside of bbox - if (within_bbox<2>(coords[0], bbox) || - within_bbox<2>(coords[1], bbox) || + if (within_bbox<2>(coords[0], bbox) || within_bbox<2>(coords[1], bbox) || within_bbox<2>(coords[2], bbox)) { return true; } @@ -230,8 +246,8 @@ bool triangle_intersects_bbox(const Omega_h::Matrix<2, 3>& coords, template [[nodiscard]] -KOKKOS_FUNCTION -bool simplex_intersects_bbox(const Omega_h::Matrix& coords, const AABBox& bbox) +KOKKOS_FUNCTION bool simplex_intersects_bbox( + const Omega_h::Matrix& coords, const AABBox& bbox) { return intersects(simplex_bbox(coords), bbox); // TODO: Add refined cases from triangle_intersects_bbox @@ -251,7 +267,8 @@ namespace detail */ struct GridTriIntersectionFunctor2D { - GridTriIntersectionFunctor2D(Omega_h::Mesh& mesh, Kokkos::View grid) + GridTriIntersectionFunctor2D(Omega_h::Mesh& mesh, + Kokkos::View grid) : mesh_(mesh), tris2verts_(mesh_.ask_elem_verts()), coords_(mesh_.coords()), @@ -274,9 +291,11 @@ struct GridTriIntersectionFunctor2D LO num_intersections = 0; // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { - const auto elem_tri2verts = Omega_h::gather_verts<3>(tris2verts_, elem_idx); + const auto elem_tri2verts = + Omega_h::gather_verts<3>(tris2verts_, elem_idx); // 2d mesh with 2d coords, but 3 triangles - const auto vertex_coords = Omega_h::gather_vectors<3, 2>(coords_, elem_tri2verts); + const auto vertex_coords = + Omega_h::gather_vectors<3, 2>(coords_, elem_tri2verts); if (triangle_intersects_bbox(vertex_coords, grid_cell_bbox)) { if (fill) { fill[num_intersections] = elem_idx; @@ -292,13 +311,15 @@ struct GridTriIntersectionFunctor2D Omega_h::LOs tris2verts_; Omega_h::Reals coords_; Kokkos::View grid_; + public: LO nelems_; }; struct GridTriIntersectionFunctor3D { - GridTriIntersectionFunctor3D(Omega_h::Mesh& mesh, Kokkos::View grid) + GridTriIntersectionFunctor3D(Omega_h::Mesh& mesh, + Kokkos::View grid) : mesh_(mesh), tets2verts_(mesh_.ask_elem_verts()), coords_(mesh_.coords()), @@ -321,8 +342,10 @@ struct GridTriIntersectionFunctor3D LO num_intersections = 0; // hierarchical parallel may make be very beneficial here... for (LO elem_idx = 0; elem_idx < nelems_; ++elem_idx) { - const auto elem_tet2verts = Omega_h::gather_verts<4>(tets2verts_, elem_idx); - const auto vertex_coords = Omega_h::gather_vectors<4, 3>(coords_, elem_tet2verts); + const auto elem_tet2verts = + Omega_h::gather_verts<4>(tets2verts_, elem_idx); + const auto vertex_coords = + Omega_h::gather_vectors<4, 3>(coords_, elem_tet2verts); if (simplex_intersects_bbox<3>(vertex_coords, grid_cell_bbox)) { if (fill) { fill[num_intersections] = elem_idx; @@ -338,14 +361,17 @@ struct GridTriIntersectionFunctor3D Omega_h::LOs tets2verts_; Omega_h::Reals coords_; Kokkos::View grid_; + public: LO nelems_; }; -// num_grid_cells should be result of grid.GetNumCells(), take as argument to avoid extra copy -// of grid from gpu to cpu +// num_grid_cells should be result of grid.GetNumCells(), take as argument to +// avoid extra copy of grid from gpu to cpu Kokkos::Crs -construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) +construct_intersection_map_2d(Omega_h::Mesh& mesh, + Kokkos::View grid, + int num_grid_cells) { Kokkos::Crs intersection_map{}; auto f = detail::GridTriIntersectionFunctor2D{mesh, grid}; @@ -354,7 +380,9 @@ construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View -construct_intersection_map_3d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells) +construct_intersection_map_3d(Omega_h::Mesh& mesh, + Kokkos::View grid, + int num_grid_cells) { Kokkos::Crs intersection_map{}; auto f = detail::GridTriIntersectionFunctor3D{mesh, grid}; @@ -363,19 +391,24 @@ construct_intersection_map_3d(Omega_h::Mesh& mesh, Kokkos::View -OMEGA_H_INLINE double myreduce(const Omega_h::Vector & x, Op op) OMEGA_H_NOEXCEPT { +template +OMEGA_H_INLINE double myreduce(const Omega_h::Vector& x, + Op op) OMEGA_H_NOEXCEPT +{ auto out = x[0]; - for (int i = 1; i < n; ++i) out = op(out, x[i]); + for (int i = 1; i < n; ++i) + out = op(out, x[i]); return out; -} +} -Kokkos::View GridPointSearch2D::operator()(Kokkos::View points) const +Kokkos::View GridPointSearch2D::operator()( + Kokkos::View points) const { - Kokkos::View results("point search result", points.extent(0)); + Kokkos::View results("point search result", + points.extent(0)); auto num_rows = candidate_map_.numRows(); // needed so that we don't capture this ptr which will be memory error on cuda - auto grid = grid_; + auto grid = grid_; auto candidate_map = candidate_map_; auto tris2verts = tris2verts_; auto tris2verts_adj = tris2verts_adj_; @@ -383,128 +416,142 @@ Kokkos::View GridPointSearch2D::operator()(Kokkos::V auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; auto tolerances = tolerances_; - Kokkos::parallel_for(points.extent(0), KOKKOS_LAMBDA(int p) { - Omega_h::Vector<2> point(std::initializer_list{points(p,0), points(p,1)}); - auto cell_id = grid(0).ClosestCellID(point); - assert(cell_id < num_rows && cell_id >= 0); - 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); - 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 - 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; - }; + Kokkos::parallel_for( + points.extent(0), KOKKOS_LAMBDA(int p) { + Omega_h::Vector<2> point( + std::initializer_list{points(p, 0), points(p, 1)}); + auto cell_id = grid(0).ClosestCellID(point); + assert(cell_id < num_rows && cell_id >= 0); + 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); + 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 + 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; + }; + } } - } - if (vertex_found) - break; + if (vertex_found) + break; - for (int j = 0; j < 3; ++j) { - // Every triangle (face) is connected to 3 edges - const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; + 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_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); + 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(point, vertex_a, vertex_b)) - continue; + if (!normal_intersects_segment(point, vertex_a, vertex_b)) + continue; - const auto xa = vertex_a[0]; - const auto ya = vertex_a[1]; - const auto xb = vertex_b[0]; - const auto yb = vertex_b[1]; + const auto xa = vertex_a[0]; + const auto ya = vertex_a[1]; + const auto xb = vertex_b[0]; + const auto yb = vertex_b[1]; - const auto xp = point[0]; - const auto yp = point[1]; + const auto xp = point[0]; + const auto yp = point[1]; - const auto distance_to_ab = distance_from_line(xp, yp, xa, ya, xb, yb); + const auto distance_to_ab = + distance_from_line(xp, yp, xa, ya, xb, yb); - 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; + 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; - if (distance_to_ab < tolerances(1)) { - edge_found = true; - }; + if (distance_to_ab < tolerances(1)) { + edge_found = true; + }; + } } - } - if (edge_found) - break; + if (edge_found) + break; - 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; + 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; - // results(p) = GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, triangleID, parametric_coords}; - // return; + // results(p) = + // GridPointSearch2D::Result{GridPointSearch2D::Result::Dimensionality::FACE, + // triangleID, parametric_coords}; return; + } } - } - 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}; - }); + 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}; + }); return results; } GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) - : GridPointSearch2D(mesh, Nx, Ny, PointSearchTolerances { "point search 2d tolerances" }) + : GridPointSearch2D(mesh, Nx, Ny, + PointSearchTolerances{"point search 2d tolerances"}) { Kokkos::deep_copy(tolerances_, 0); } -GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances) +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, + const PointSearchTolerances& tolerances) : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); - grid_h(0) = Uniform2DGrid{.edge_length = {mesh_bbox.max[0] - mesh_bbox.min[0], - mesh_bbox.max[1] - mesh_bbox.min[1]}, - .bot_left = {mesh_bbox.min[0], mesh_bbox.min[1]}, - .divisions = {Nx, Ny}}; + grid_h(0) = + Uniform2DGrid{.edge_length = {mesh_bbox.max[0] - mesh_bbox.min[0], + mesh_bbox.max[1] - mesh_bbox.min[1]}, + .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()); + 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); @@ -512,9 +559,11 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const Po edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); } -Kokkos::View GridPointSearch3D::operator()(Kokkos::View points) const +Kokkos::View GridPointSearch3D::operator()( + Kokkos::View points) const { - Kokkos::View results("point search result", points.extent(0)); + Kokkos::View results("point search result", + points.extent(0)); auto num_rows = candidate_map_.numRows(); // needed so that we don't capture this ptr which will be memory error on cuda auto grid = grid_; @@ -524,77 +573,83 @@ Kokkos::View GridPointSearch3D::operator()(Kokkos::V auto tris2edges_adj = tris2edges_adj_; auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; - Kokkos::parallel_for(points.extent(0), KOKKOS_LAMBDA(int p) { + Kokkos::parallel_for( + points.extent(0), KOKKOS_LAMBDA(int p) { + Omega_h::Vector point; + for (int i = 0; i < DIM; ++i) { + point[i] = points(p, i); + } - Omega_h::Vector point; - for (int i = 0; i < DIM; ++i) { - point[i] = points(p, i); - } + auto cell_id = grid(0).ClosestCellID(point); + assert(cell_id < num_rows && cell_id >= 0); + auto candidates_begin = candidate_map.row_map(cell_id); + auto candidates_end = candidate_map.row_map(cell_id + 1); + bool found = false; + + auto nearest_triangle = candidates_begin; + auto dimensionality = GridPointSearch3D::Result::Dimensionality::EDGE; + Omega_h::Real distance_to_nearest{INFINITY}; + Omega_h::Vector 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); + const auto elem_tri2verts = + Omega_h::gather_verts(tris2verts, triangleID); + auto vertex_coords = + Omega_h::gather_vectors(coords, elem_tri2verts); + auto parametric_coords = + Omega_h::barycentric_from_global(point, vertex_coords); + + if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { + results(p) = GridPointSearch3D::Result{ + GridPointSearch3D::Result::Dimensionality::REGION, triangleID, + parametric_coords}; + found = true; + break; + } - auto cell_id = grid(0).ClosestCellID(point); - assert(cell_id < num_rows && cell_id >= 0); - auto candidates_begin = candidate_map.row_map(cell_id); - auto candidates_end = candidate_map.row_map(cell_id + 1); - bool found = false; - - auto nearest_triangle = candidates_begin; - auto dimensionality = GridPointSearch3D::Result::Dimensionality::EDGE; - Omega_h::Real distance_to_nearest { INFINITY }; - Omega_h::Vector 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); - const auto elem_tri2verts = Omega_h::gather_verts(tris2verts, triangleID); - auto vertex_coords = Omega_h::gather_vectors(coords, elem_tri2verts); - auto parametric_coords = Omega_h::barycentric_from_global(point, vertex_coords); - - if (Omega_h::is_barycentric_inside(parametric_coords, fuzz)) { - results(p) = GridPointSearch3D::Result{GridPointSearch3D::Result::Dimensionality::REGION, triangleID, parametric_coords}; - found = true; - break; + // TODO: Get nearest element if no triangle found } - - // TODO: Get nearest element if no triangle found - } - if(!found) - { - results(p) = GridPointSearch3D::Result{dimensionality, -1 * candidate_map.entries(nearest_triangle), parametric_coords_to_nearest}; - } - }); + if (!found) { + results(p) = GridPointSearch3D::Result{ + dimensionality, -1 * candidate_map.entries(nearest_triangle), + parametric_coords_to_nearest}; + } + }); return results; } GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) - : GridPointSearch3D(mesh, Nx, Ny, Nz, PointSearchTolerances { "point search 3d tolerances" }) + : GridPointSearch3D(mesh, Nx, Ny, Nz, + 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) - : PointLocalizationSearch(tolerances) + : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); - std::array edge_lengths {}; - std::array bot_left {}; + std::array edge_lengths{}; + std::array bot_left{}; for (int i = 0; i < DIM; ++i) { edge_lengths[i] = mesh_bbox.max[i] - mesh_bbox.min[i]; bot_left[i] = mesh_bbox.min[i]; } - grid_h(0) = Uniform3DGrid { - .edge_length = edge_lengths, - .bot_left = bot_left, - .divisions = {Nx, Ny, Nz} - }; + grid_h(0) = Uniform3DGrid{.edge_length = edge_lengths, + .bot_left = bot_left, + .divisions = {Nx, Ny, Nz}}; Kokkos::deep_copy(grid_, grid_h); - candidate_map_ = detail::construct_intersection_map_3d(mesh, grid_, grid_h(0).GetNumCells()); + candidate_map_ = + detail::construct_intersection_map_3d(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); diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 0910b771..6925d77f 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -17,9 +17,12 @@ namespace pcms // TODO take a bounding box as we may want a bbox that's bigger than the mesh! // this function is in the public header for testing, but should not be directly // used -namespace detail { +namespace detail +{ Kokkos::Crs -construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); +construct_intersection_map_2d(Omega_h::Mesh& mesh, + Kokkos::View grid, + int num_grid_cells); } KOKKOS_FUNCTION @@ -30,7 +33,8 @@ template class PointLocalizationSearch { public: - struct Result { + struct Result + { enum class Dimensionality { VERTEX = 0, @@ -48,12 +52,14 @@ class PointLocalizationSearch using PointSearchTolerances = Kokkos::View; - explicit PointLocalizationSearch(const PointSearchTolerances& tolerances) : tolerances_(tolerances) + explicit PointLocalizationSearch(const PointSearchTolerances& tolerances) + : tolerances_(tolerances) { assert(tolerances_.is_allocated()); } - virtual Kokkos::View operator()(Kokkos::View point) const = 0; + virtual Kokkos::View operator()( + Kokkos::View point) const = 0; virtual ~PointLocalizationSearch() = default; protected: @@ -65,13 +71,15 @@ using PointLocalizationSearch3D = PointLocalizationSearch<3>; class GridPointSearch2D : public PointLocalizationSearch2D { - using CandidateMapT = Kokkos::Crs; + using CandidateMapT = + Kokkos::Crs; public: using Result = PointLocalizationSearch2D::Result; GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny); - GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances); + GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, + const PointSearchTolerances& tolerances); /** * given a point in global coordinates give the id of the triangle that the @@ -80,7 +88,8 @@ class GridPointSearch2D : public PointLocalizationSearch2D * id will be a negative number and (TODO) will return a negative id of the * closest element */ - Kokkos::View operator()(Kokkos::View point) const override; + Kokkos::View operator()( + Kokkos::View point) const override; private: Omega_h::Mesh mesh_; @@ -95,13 +104,15 @@ class GridPointSearch2D : public PointLocalizationSearch2D class GridPointSearch3D : public PointLocalizationSearch3D { - using CandidateMapT = Kokkos::Crs; + using CandidateMapT = + Kokkos::Crs; 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, const PointSearchTolerances& tolerances); + GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, + const PointSearchTolerances& tolerances); /** * given a point in global coordinates give the id of the triangle that the @@ -110,7 +121,8 @@ class GridPointSearch3D : public PointLocalizationSearch3D * id will be a negative number and (TODO) will return a negative id of the * closest element */ - Kokkos::View operator()(Kokkos::View point) const override; + Kokkos::View operator()( + Kokkos::View point) const override; private: Omega_h::Mesh mesh_; @@ -123,5 +135,5 @@ class GridPointSearch3D : public PointLocalizationSearch3D Omega_h::Reals coords_; }; -} // namespace detail +} // namespace pcms #endif // PCMS_COUPLING_POINT_SEARCH_H diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 22f4fcff..01093144 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -111,7 +111,8 @@ bool num_candidates_within_range(const T& intersection_map, pcms::LO min, Kokkos::parallel_reduce( intersection_map.numRows(), KOKKOS_LAMBDA(const int i, result_type& update) { - auto num_candidates = intersection_map.row_map(i + 1) - intersection_map.row_map(i); + auto num_candidates = + intersection_map.row_map(i + 1) - intersection_map.row_map(i); if (num_candidates > update.max_val) update.max_val = num_candidates; if (num_candidates < update.min_val) @@ -123,11 +124,11 @@ bool num_candidates_within_range(const T& intersection_map, pcms::LO min, std::cerr << result.min_val << ' ' << result.max_val << '\n'; return within_range; } -//extern Omega_h::Library omega_h_library; +// extern Omega_h::Library omega_h_library; TEST_CASE("construct intersection map") { - //auto world = omega_h_library.world(); + // auto world = omega_h_library.world(); auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = @@ -137,9 +138,11 @@ TEST_CASE("construct intersection map") { Kokkos::View grid_d("uniform grid"); auto grid_h = Kokkos::create_mirror_view(grid_d); - grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; + grid_h(0) = Uniform2DGrid{ + .edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {10, 10}}; Kokkos::deep_copy(grid_d, grid_h); - auto intersection_map = pcms::detail::construct_intersection_map_2d(mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map_2d( + mesh, grid_d, grid_h(0).GetNumCells()); // assert(cudaSuccess == cudaDeviceSynchronize()); REQUIRE(intersection_map.numRows() == 100); REQUIRE(num_candidates_within_range(intersection_map, 2, 16)); @@ -148,40 +151,44 @@ TEST_CASE("construct intersection map") { Kokkos::View grid_d("uniform grid"); auto grid_h = Kokkos::create_mirror_view(grid_d); - grid_h(0) = Uniform2DGrid{.edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {60, 60}}; + grid_h(0) = Uniform2DGrid{ + .edge_length{1, 1}, .bot_left = {0, 0}, .divisions = {60, 60}}; Kokkos::deep_copy(grid_d, grid_h); // require number of candidates is >=1 and <=6 - auto intersection_map = pcms::detail::construct_intersection_map_2d(mesh, grid_d, grid_h(0).GetNumCells()); + auto intersection_map = pcms::detail::construct_intersection_map_2d( + mesh, grid_d, grid_h(0).GetNumCells()); REQUIRE(intersection_map.numRows() == 3600); REQUIRE(num_candidates_within_range(intersection_map, 1, 6)); } } -TEST_CASE("uniform grid search") { +TEST_CASE("uniform grid search") +{ using pcms::GridPointSearch2D; auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); - auto tolerances = GridPointSearch2D::PointSearchTolerances { "point search 2d tolerances" }; + auto tolerances = + GridPointSearch2D::PointSearchTolerances{"point search 2d tolerances"}; tolerances(0) = 0.01; tolerances(1) = 0.01; - GridPointSearch2D search{mesh,10,10, tolerances}; + GridPointSearch2D search{mesh, 10, 10, tolerances}; - Kokkos::View points("test_points", 8); - //Kokkos::View points("test_points", 1); + Kokkos::View points("test_points", 8); + // Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); - points_h(0,0) = 0; - points_h(0,1) = 0; - points_h(1,0) = 0.55; - points_h(1,1) = 0.54; - points_h(2,0) = 100; - points_h(2,1) = 100; - points_h(3,0) = 1; - points_h(3,1) = 1; - points_h(4,0) = -1; - points_h(4,1) = -1; + points_h(0, 0) = 0; + points_h(0, 1) = 0; + points_h(1, 0) = 0.55; + points_h(1, 1) = 0.54; + points_h(2, 0) = 100; + points_h(2, 1) = 100; + points_h(3, 0) = 1; + points_h(3, 1) = 1; + points_h(4, 0) = -1; + points_h(4, 1) = -1; points_h(5, 0) = 1.01; points_h(5, 1) = 0.95; points_h(6, 0) = 0.05; @@ -195,7 +202,7 @@ TEST_CASE("uniform grid search") { SECTION("global coordinate within mesh") { { - auto [dim, idx,coords] = results_h(0); + auto [dim, idx, coords] = results_h(0); CAPTURE(idx); @@ -206,7 +213,7 @@ TEST_CASE("uniform grid search") { REQUIRE(coords[2] == Catch::Approx(0)); } { - auto [dim, idx,coords] = results_h(1); + auto [dim, idx, coords] = results_h(1); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); REQUIRE(idx == 156); REQUIRE(coords[0] == Catch::Approx(0.5)); @@ -214,22 +221,25 @@ TEST_CASE("uniform grid search") { REQUIRE(coords[2] == Catch::Approx(0.4)); } { - auto [dim, idx,coords] = results_h(7); + auto [dim, idx, coords] = results_h(7); REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); REQUIRE(idx == 0); } } // feature needs to be added - SECTION("Global coordinate outside mesh", "[!mayfail]") { + SECTION("Global coordinate outside mesh", "[!mayfail]") + { auto out_of_bounds = results_h(2); auto top_right = results_h(3); - REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); - REQUIRE(-1*out_of_bounds.element_id == top_right.element_id); + REQUIRE(out_of_bounds.dimensionality == + GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(-1 * out_of_bounds.element_id == top_right.element_id); out_of_bounds = results_h(4); auto bot_left = results_h(0); - REQUIRE(out_of_bounds.dimensionality == GridPointSearch2D::Result::Dimensionality::VERTEX); - REQUIRE(-1*out_of_bounds.element_id == bot_left.element_id); + REQUIRE(out_of_bounds.dimensionality == + GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); out_of_bounds = results_h(5); REQUIRE(out_of_bounds.dimensionality == From 23678a5c6f01a5a30256a4c008b5b178a72fc0b6 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 5 Nov 2025 15:49:16 -0500 Subject: [PATCH 14/22] fix narrowing conversion --- src/pcms/point_search.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index d484ba4f..aac27341 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -188,7 +188,7 @@ template // each dimension has a pair of opposing "walls" // 2D: { [left, right], [top, bottom] } -> { left, right, top, bottom } // 3D: { [left, right], [top, bottom], [front, back] } -> { left, ..., back } - std::array bbox_walls{}; + std::array bbox_walls{}; for (int i = 0; i < dim; i++) { bbox_walls[i * 2] = bbox.center[i] - bbox.half_width[i]; bbox_walls[i * 2 + 1] = bbox.center[i] + bbox.half_width[i]; From 2e7e4e15dc6633d2a0771e99a89f4168ea58eaa3 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Mon, 10 Nov 2025 00:16:48 -0500 Subject: [PATCH 15/22] Correctly initialize initial min and max for n dimensions Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/pcms/point_search.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index aac27341..49fe8999 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -61,9 +61,12 @@ AABBox<2> triangle_bbox(const Omega_h::Matrix<2, 3>& coords) template AABBox simplex_bbox(const Omega_h::Matrix& coords) { - std::array max{coords(0, 0), coords(1, 0)}; - std::array min{coords(0, 0), coords(1, 0)}; - + std::array max; + std::array min; + for (int j = 0; j < dim; ++j) { + max[j] = coords(j, 0); + min[j] = coords(j, 0); + } for (int i = 1; i < dim + 1; ++i) { for (int j = 0; j < dim; ++j) { max[j] = std::fmax(max[j], coords(j, i)); From 879586f5b31a51aed97df1dae5b6eb8c3699f4b9 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Mon, 10 Nov 2025 00:17:34 -0500 Subject: [PATCH 16/22] correct wording for 3D meshes Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/pcms/point_search.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 49fe8999..447ae128 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -331,7 +331,7 @@ struct GridTriIntersectionFunctor3D { if (mesh_.dim() != 3) { std::cerr << "GridTriIntersection3D currently only developed for 3D " - "triangular meshes\n"; + "tetrahedral meshes\n"; std::terminate(); } } From 30b7e54ffe6fdb319845a8dbf62932ba4c08f861 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Mon, 10 Nov 2025 00:18:12 -0500 Subject: [PATCH 17/22] correct working for 3D meshes Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/pcms/point_search.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 6925d77f..de745374 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -115,11 +115,11 @@ class GridPointSearch3D : public PointLocalizationSearch3D const PointSearchTolerances& tolerances); /** - * given a point in global coordinates give the id of the triangle that the - * point lies within and the parametric coordinate of the point within the - * triangle. If the point does not lie within any triangle element. Then the + * Given a point in global coordinates, returns the id of the tetrahedron (3D element) + * that the point lies within and the parametric coordinate of the point within the + * tetrahedron. If the point does not lie within any tetrahedron element, then the * id will be a negative number and (TODO) will return a negative id of the - * closest element + * closest element. */ Kokkos::View operator()( Kokkos::View point) const override; From 32e886f5a8e05e94f70a4301730b362b020d58c1 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Mon, 10 Nov 2025 00:24:38 -0500 Subject: [PATCH 18/22] update TODO comment in omega_h_field.h --- src/pcms/adapter/omega_h/omega_h_field.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index 52784e89..745c723b 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -219,7 +219,7 @@ class OmegaHField private: std::string name_; Omega_h::Mesh& mesh_; - // TODO make this a pointer and introduce base class to Search for alternative search methods + // TODO: introduce base class to Search for alternative search methods std::unique_ptr search_; // bitmask array that specifies a filter on the field Omega_h::Read mask_; From a88e058f9fb837fdb8a71295dcddc4aaaefeb6aa Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 12 Nov 2025 12:12:15 -0500 Subject: [PATCH 19/22] Update OmegaHField2 to use GridPointSearch2D Replaced `GridPointSearch` with `GridPointSearch2D` in `omega_h_field2` for improved dimensionality handling. --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 2 +- src/pcms/adapter/omega_h/omega_h_field2.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index ea45bf10..ed3606c7 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -158,7 +158,7 @@ struct OmegaHField2LocalizationHint { OmegaHField2LocalizationHint( Omega_h::Mesh& mesh, - Kokkos::View search_results) + Kokkos::View search_results) : offsets_("", mesh.nelems() + 1), coordinates_("", search_results.size(), mesh.dim() + 1), indices_("", search_results.size()) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.h b/src/pcms/adapter/omega_h/omega_h_field2.h index 2c3365c8..bc6248e4 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.h +++ b/src/pcms/adapter/omega_h/omega_h_field2.h @@ -60,7 +60,7 @@ class OmegaHField2 : public FieldT const OmegaHFieldLayout& layout_; Omega_h::Mesh& mesh_; std::unique_ptr mesh_field_; - GridPointSearch search_; + GridPointSearch2D search_; Kokkos::View dof_holder_data_; }; From c0a3180d7b8dc4a9553b4f2403b9b4c5aa69b363 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 19 Nov 2025 14:14:05 -0500 Subject: [PATCH 20/22] Apply minor formatting adjustments and line splits in point_search and omega_h files --- src/pcms/adapter/omega_h/omega_h_field.h | 3 ++- src/pcms/point_search.h | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index 745c723b..82781a4a 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -160,7 +160,8 @@ class OmegaHField [[nodiscard]] auto Search(Kokkos::View points) const { PCMS_FUNCTION_TIMER; - PCMS_ALWAYS_ASSERT(search_ != nullptr && "search data structure must be constructed before use"); + PCMS_ALWAYS_ASSERT(search_ != nullptr && + "search data structure must be constructed before use"); return (*search_)(points); } diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index de745374..854a46a2 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -115,11 +115,11 @@ class GridPointSearch3D : public PointLocalizationSearch3D const PointSearchTolerances& tolerances); /** - * Given a point in global coordinates, returns the id of the tetrahedron (3D element) - * that the point lies within and the parametric coordinate of the point within the - * tetrahedron. If the point does not lie within any tetrahedron element, then the - * id will be a negative number and (TODO) will return a negative id of the - * closest element. + * Given a point in global coordinates, returns the id of the tetrahedron (3D + * element) that the point lies within and the parametric coordinate of the + * point within the tetrahedron. If the point does not lie within any + * tetrahedron element, then the id will be a negative number and (TODO) will + * return a negative id of the closest element. */ Kokkos::View operator()( Kokkos::View point) const override; From 9c458a76b166687f27e9060036b50ae1d51be91e Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 19 Nov 2025 14:36:15 -0500 Subject: [PATCH 21/22] Remove unused `KOKKOS_FUNCTION` annotation and correct TODO comment in point_search files --- src/pcms/point_search.cpp | 2 +- src/pcms/point_search.h | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pcms/point_search.cpp b/src/pcms/point_search.cpp index 447ae128..01cbf745 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -612,7 +612,7 @@ Kokkos::View GridPointSearch3D::operator()( break; } - // TODO: Get nearest element if no triangle found + // TODO: Get nearest element if no tetrahedron found } if (!found) { results(p) = GridPointSearch3D::Result{ diff --git a/src/pcms/point_search.h b/src/pcms/point_search.h index 854a46a2..4f231c16 100644 --- a/src/pcms/point_search.h +++ b/src/pcms/point_search.h @@ -24,7 +24,6 @@ construct_intersection_map_2d(Omega_h::Mesh& mesh, Kokkos::View grid, int num_grid_cells); } -KOKKOS_FUNCTION [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); From 2926a7994598d3392a4078bc8812af2ead6f3a58 Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Wed, 10 Dec 2025 11:15:24 -0500 Subject: [PATCH 22/22] Replace `GridPointSearch` with `GridPointSearch2D` in `omega_h_field2` for 2D-specific compatibility improvements and clarity --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index ed3606c7..d000c314 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -97,11 +97,11 @@ struct ComputeOffsetsFunctor struct CountPointsPerElementFunctor { Kokkos::View elem_counts_; - Kokkos::View search_results_; + Kokkos::View search_results_; CountPointsPerElementFunctor( Kokkos::View elem_counts, - Kokkos::View search_results) + Kokkos::View search_results) : elem_counts_(elem_counts), search_results_(search_results) { } @@ -121,13 +121,13 @@ struct FillCoordinatesAndIndicesFunctor Kokkos::View offsets_; Kokkos::View coordinates_; Kokkos::View indices_; - Kokkos::View search_results_; + Kokkos::View search_results_; 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) : mesh_(mesh), elem_counts_(elem_counts), offsets_(offsets), @@ -299,7 +299,7 @@ LocalizationHint OmegaHField2::GetLocalizationHint( coordinates.data_handle(), coordinates.extent(0), coordinates.extent(1)); deep_copy_mismatch_layouts(coords, coordinates_host); auto results = search_(coords); - Kokkos::View results_h( + Kokkos::View results_h( "results_h", results.size()); Kokkos::deep_copy(results_h, results); auto hint = std::make_shared(mesh_, results_h);