diff --git a/src/pcms/adapter/omega_h/omega_h_field.h b/src/pcms/adapter/omega_h/omega_h_field.h index c0aaa17e..82781a4a 100644 --- a/src/pcms/adapter/omega_h/omega_h_field.h +++ b/src/pcms/adapter/omega_h/omega_h_field.h @@ -154,13 +154,13 @@ class OmegaHField void ConstructSearch(int nx, int ny) { PCMS_FUNCTION_TIMER; - search_ = GridPointSearch(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() && + PCMS_ALWAYS_ASSERT(search_ != nullptr && "search data structure must be constructed before use"); return (*search_)(points); } @@ -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: 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_; LO size_; diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index ea45bf10..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), @@ -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()) @@ -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); 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_; }; diff --git a/src/pcms/interpolator/adj_search.hpp b/src/pcms/interpolator/adj_search.hpp index c39639fc..040d96a9 100644 --- a/src/pcms/interpolator/adj_search.hpp +++ b/src/pcms/interpolator/adj_search.hpp @@ -24,14 +24,14 @@ 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"); 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); @@ -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); @@ -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 e86cab36..01cbf745 100644 --- a/src/pcms/point_search.cpp +++ b/src/pcms/point_search.cpp @@ -57,6 +57,34 @@ 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; + 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)); + 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. @@ -117,15 +145,17 @@ 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 + const Omega_h::Vector coord, const AABBox& 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); + 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 bool bbox_verts_within_triangle( @@ -135,25 +165,59 @@ bool line_intersects_bbox(const Omega_h::Vector<2>& p0, 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; } +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 + constexpr unsigned num_verts = 1 << dim; + + // Each vertex is a just a unique combination of walls + // eg [left, bottom] (2D) or [right, top, front] (3) + for (unsigned i = 0; i < num_verts; ++i) { + // conveniently, i acts a bit field representing the current combination + Omega_h::Vector vert; + for (unsigned j = 0; j < dim; ++j) { + // eg 110 = 6 = [left, top, back] + // eg 01 = 1 = [right, top] + vert[j] = (i >> j) & 1 ? bbox_walls[j * 2] : bbox_walls[j * 2 + 1]; + } + auto xi = Omega_h::barycentric_from_global(vert, coords); + if (Omega_h::is_barycentric_inside(xi, fuzz)) { + return true; + } + } + return false; +} + /** * Check if a triangle element represented by 3 coordinates in two dimensions * intersects with a bounding box @@ -165,8 +229,8 @@ KOKKOS_FUNCTION bool triangle_intersects_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) || - 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 @@ -183,6 +247,15 @@ KOKKOS_FUNCTION bool triangle_intersects_bbox( 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 { /** @@ -195,10 +268,10 @@ 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()), @@ -206,7 +279,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(); } @@ -246,31 +319,80 @@ 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 " + "tetrahedral 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); + 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) +construct_intersection_map_2d(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; } -} // 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) +Kokkos::Crs +construct_intersection_map_3d(Omega_h::Mesh& mesh, + Kokkos::View grid, + int num_grid_cells) { - 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]}; + 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 OMEGA_H_INLINE double myreduce(const Omega_h::Vector& x, @@ -282,12 +404,11 @@ OMEGA_H_INLINE double myreduce(const Omega_h::Vector& x, 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_; @@ -297,6 +418,7 @@ Kokkos::View GridPointSearch::operator()( 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( @@ -305,10 +427,13 @@ Kokkos::View GridPointSearch::operator()( 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; + 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 @@ -320,16 +445,34 @@ Kokkos::View GridPointSearch::operator()( // 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) = GridPointSearch::Result{ - GridPointSearch::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 + 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; + for (int j = 0; j < 3; ++j) { // Every triangle (face) is connected to 3 edges const int edgeID = tris2edges_adj.ab2b[triangleID * 3 + j]; @@ -354,44 +497,53 @@ Kokkos::View GridPointSearch::operator()( 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) { + dimensionality = GridPointSearch2D::Result::Dimensionality::EDGE; + nearest_element_id = edgeID; + distance_to_nearest = distance_to_ab; + 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 (distance_to_ab < tolerances(1)) { + edge_found = true; + }; + } } - // 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 (edge_found) + break; - 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 (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; } } - if (!found) { - results(p) = GridPointSearch::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; } -GridPointSearch::GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny) +GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO 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, + const PointSearchTolerances& tolerances) + : PointLocalizationSearch(tolerances) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); @@ -402,7 +554,105 @@ GridPointSearch::GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny) .divisions = {Nx, Ny}}; Kokkos::deep_copy(grid_, grid_h); candidate_map_ = - detail::construct_intersection_map(mesh, grid_, grid_h(0).GetNumCells()); + 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::REGION, triangleID, + parametric_coords}; + found = true; + break; + } + + // TODO: Get nearest element if no tetrahedron found + } + if (!found) { + results(p) = GridPointSearch3D::Result{ + dimensionality, -1 * candidate_map.entries(nearest_triangle), + parametric_coords_to_nearest}; + } + }); + + return results; +} + +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_); + + 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 05aff657..4f231c16 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" @@ -18,39 +20,66 @@ namespace pcms 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 -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); -class GridPointSearch +template +class PointLocalizationSearch { - using CandidateMapT = - Kokkos::Crs; - public: - static constexpr auto dim = 2; struct Result { enum class Dimensionality { VERTEX = 0, EDGE = 1, - FACE = 2 + FACE = 2, + REGION = 3 }; Dimensionality dimensionality; - LO tri_id; + LO element_id; Omega_h::Vector parametric_coords; }; - GridPointSearch(Omega_h::Mesh& mesh, LO Nx, LO Ny); + static constexpr auto DIM = dim; + + using PointSearchTolerances = Kokkos::View; + + explicit PointLocalizationSearch(const PointSearchTolerances& tolerances) + : tolerances_(tolerances) + { + assert(tolerances_.is_allocated()); + } + + virtual Kokkos::View operator()( + Kokkos::View point) const = 0; + virtual ~PointLocalizationSearch() = default; + +protected: + PointSearchTolerances tolerances_; +}; + +using PointLocalizationSearch2D = PointLocalizationSearch<2>; +using PointLocalizationSearch3D = PointLocalizationSearch<3>; + +class GridPointSearch2D : public PointLocalizationSearch2D +{ + 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); + /** * 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 @@ -58,7 +87,8 @@ 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 +101,38 @@ class GridPointSearch 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); + GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, + 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. + */ + 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 pcms #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 diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 06bacf48..01093144 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]); @@ -142,7 +141,7 @@ TEST_CASE("construct intersection map") 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( + auto intersection_map = pcms::detail::construct_intersection_map_2d( mesh, grid_d, grid_h(0).GetNumCells()); // assert(cudaSuccess == cudaDeviceSynchronize()); REQUIRE(intersection_map.numRows() == 100); @@ -156,7 +155,7 @@ TEST_CASE("construct intersection map") .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( + auto intersection_map = pcms::detail::construct_intersection_map_2d( mesh, grid_d, grid_h(0).GetNumCells()); REQUIRE(intersection_map.numRows() == 3600); @@ -165,13 +164,19 @@ TEST_CASE("construct intersection map") } TEST_CASE("uniform grid search") { - using pcms::GridPointSearch; + 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); + 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", 8); // Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); points_h(0, 0) = 0; @@ -188,6 +193,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); @@ -196,7 +203,10 @@ TEST_CASE("uniform grid search") { { auto [dim, idx, coords] = results_h(0); - REQUIRE(dim == GridPointSearch::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)); @@ -204,12 +214,17 @@ TEST_CASE("uniform grid search") } { auto [dim, idx, coords] = results_h(1); - REQUIRE(dim == GridPointSearch::Result::Dimensionality::FACE); - REQUIRE(idx == 91); + REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); + 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]") @@ -217,23 +232,23 @@ TEST_CASE("uniform grid search") 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); + 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 == - GridPointSearch::Result::Dimensionality::VERTEX); - REQUIRE(-1 * out_of_bounds.tri_id == bot_left.tri_id); + 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 == - GridPointSearch::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.tri_id == top_right.tri_id); + GridPointSearch2D::Result::Dimensionality::EDGE); + REQUIRE(out_of_bounds.element_id == 219); out_of_bounds = results_h(6); REQUIRE(out_of_bounds.dimensionality == - GridPointSearch::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.tri_id == bot_left.tri_id); + GridPointSearch2D::Result::Dimensionality::EDGE); + REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); } }