diff --git a/README.md b/README.md index be15e9bb..2a57876d 100644 --- a/README.md +++ b/README.md @@ -91,6 +91,14 @@ cmake -S redev -B build-redev \ -DADIOS2_DIR=$PWD/build-ADIOS2/install/lib64/cmake/adios2 \ -Dperfstubs_DIR=$PWD/build-perfstubs/install/lib/cmake cmake --build build-redev --target install + +git clone git@github.com:arborx/ArborX.git +cmake -S ArborX -B build-ArborX -D CMAKE_INSTALL_PREFIX="$PWD/build-ArborX/install" \ +-D ARBORX_ENABLE_MPI=ON \ +-D Kokkos_ROOT="$PWD/build-kokkos/install" \ +-D CMAKE_CXX_COMPILER=mpicxx \ +-D CMAKE_CXX_EXTENSIONS=OFF +cmake --build build-ArborX --target install ``` ### Build dependencies with CUDA ``` @@ -164,6 +172,14 @@ git clone git@github.com:catchorg/Catch2 cmake -S Catch2 -B build-Catch2 \ -DCMAKE_INSTALL_PREFIX=$PWD/build-Catch2/install cmake --build build-Catch2 --target install + +git clone git@github.com:arborx/ArborX.git +cmake -S ArborX -B build-ArborX -D CMAKE_INSTALL_PREFIX="$PWD/build-ArborX/install" \ +-D ARBORX_ENABLE_MPI=ON \ +-D Kokkos_ROOT="$PWD/build-kokkos/install" \ +-D CMAKE_CXX_COMPILER="$PWD/build-kokkos/install/bin/nvcc_wrapper" \ +-D CMAKE_CXX_EXTENSIONS=OFF +cmake --build build-ArborX --target install ``` ### Build, install, and test pcms ``` @@ -181,6 +197,7 @@ cmake -S pcms -B build-pcms \ -Dperfstubs_DIR=$PWD/build-perfstubs/install/lib/cmake \ -DCatch2_DIR=$PWD/build-Catch2/install/lib64/cmake/Catch2/ \ -DKokkosKernels_DIR=$PWD/build-kokkos-kernels/install/lib64/cmake/KokkosKernels/ \ + -DArborX_DIR=$PWD/build-ArborX/install/share/cmake/ArborX \ -DPCMS_TEST_DATA_DIR=$PWD/pcms_testcases cmake --build build-pcms -j 8 diff --git a/src/pcms/field/evaluator/omega_h_lagrange.h b/src/pcms/field/evaluator/omega_h_lagrange.h index 39ed0c94..37d1c99d 100644 --- a/src/pcms/field/evaluator/omega_h_lagrange.h +++ b/src/pcms/field/evaluator/omega_h_lagrange.h @@ -65,13 +65,11 @@ struct CopyCoordsFunctor template OmegaHLagrangeLocHint BuildLagrangeLocHint( Omega_h::Mesh& mesh, int mesh_dim, - Kokkos::View::Result*, - DeviceMemorySpace> - results, + PointSearch::Results& results, Kokkos::View coords_d, Kokkos::View owning_ids, OutOfBoundsMode mode) { - LO n = static_cast(results.size()); + LO n = static_cast(results.dimensionalities.size()); // First pass: count valid and missing Kokkos::View is_valid("is_valid", n); @@ -79,7 +77,7 @@ OmegaHLagrangeLocHint BuildLagrangeLocHint( "CheckValidity", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(LO i) { - bool out = (owning_ids(i) < 0) || (results(i).element_id < 0); + bool out = (owning_ids(i) < 0) || (results.element_ids(i) < 0); is_valid(i) = out ? 0 : 1; }); @@ -344,7 +342,7 @@ class OmegaHLagrangeEvaluatorFactory : public FieldEvaluatorFactory detail::CopyCoordsFunctor copy_functor(coords_d, raw_coords); Kokkos::parallel_for("copy_coords", n_pts, copy_functor); - auto results_d = search(coords_d); + auto results_d = search.apply(coords); auto owning_ids = search.GetOwningElementIds(results_d); return detail::BuildLagrangeLocHint( diff --git a/src/pcms/field/uniform_grid_binary_field.h b/src/pcms/field/uniform_grid_binary_field.h index 6a71de7a..61483053 100644 --- a/src/pcms/field/uniform_grid_binary_field.h +++ b/src/pcms/field/uniform_grid_binary_field.h @@ -47,7 +47,7 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid& grid) auto coords = coord_view.GetValues(); LO n = layout->GetNumOwnedDofHolder(); - Kokkos::View coords_d("coords_d", n); + Kokkos::View coords_d("coords_d", n, Dim); Kokkos::parallel_for( "CopyCoords", Kokkos::RangePolicy<>(0, n), KOKKOS_LAMBDA(int i) { for (int d = 0; d < Dim; ++d) { @@ -56,21 +56,21 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid& grid) }); // Run point-in-mesh search - Kokkos::View::Result*> results_d; + PointSearch::Results results_d; if constexpr (Dim == 2) { GridPointSearch2D search(mesh, grid.divisions[0], grid.divisions[1]); - results_d = search(coords_d); + results_d = search.apply(pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, pcms::MakeConstRank2View(coords_d))); } else { GridPointSearch3D search(mesh, grid.divisions[0], grid.divisions[1], grid.divisions[2]); - results_d = search(coords_d); + results_d = search.apply(pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, pcms::MakeConstRank2View(coords_d))); } - auto results_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, results_d); + auto result_ids_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, results_d.element_ids); Kokkos::View data("binary_mask", n); for (LO i = 0; i < n; ++i) - data(i) = (results_h(i).element_id >= 0) ? 1.0 : 0.0; + data(i) = (result_ids_h(i) >= 0) ? 1.0 : 0.0; field.SetDOFHolderDataHost( Rank2View(data.data(), n, 1)); diff --git a/src/pcms/localization/CMakeLists.txt b/src/pcms/localization/CMakeLists.txt index 269b1b77..a3ac9156 100644 --- a/src/pcms/localization/CMakeLists.txt +++ b/src/pcms/localization/CMakeLists.txt @@ -7,6 +7,7 @@ set( localization_factory.h point_cloud_localization.h mesh_localization.h + point_localization.h ) set( @@ -15,7 +16,9 @@ set( adj_search.cpp mls_support_helpers.cpp mesh_localization.cpp + point_localization.cpp ) +find_package(ArborX REQUIRED) add_library(pcms_localization ${PCMS_LOCALIZATION_SOURCES}) add_library(pcms::localization ALIAS pcms_localization) target_include_directories( @@ -32,6 +35,7 @@ target_link_libraries( pcms::utility pcms::discretization Omega_h::omega_h + ArborX::ArborX ) target_compile_features(pcms_localization PUBLIC cxx_std_20) diff --git a/src/pcms/localization/adj_search.cpp b/src/pcms/localization/adj_search.cpp index 7dd25364..8e400fdc 100644 --- a/src/pcms/localization/adj_search.cpp +++ b/src/pcms/localization/adj_search.cpp @@ -7,7 +7,7 @@ namespace pcms // Debug helper retained intentionally: useful for diagnosing point-localization // failures while developing support-search logic. [[maybe_unused]] static void checkTargetPoints( - const Kokkos::View& results, + const GridPointSearch2D::Results& results, const Kokkos::View& owning_cell_ids) { Kokkos::fence(); @@ -22,7 +22,7 @@ namespace pcms printf("%d, ", i); } }; - Omega_h::parallel_for(results.size(), check_target_points, + Omega_h::parallel_for(results.dimensionalities.size(), check_target_points, "check_target_points"); Kokkos::fence(); pcms::printInfo("\n"); @@ -81,8 +81,8 @@ static Omega_h::Write locate_target_cells( Omega_h::Write(nvertices_target, -1, "source cell ids"); if (dim == 2) { - Kokkos::View target_points("target_points", - nvertices_target); + Kokkos::View target_points("target_points", + nvertices_target, 2); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { target_points(i, 0) = target_coords[i * dim]; @@ -90,8 +90,10 @@ static Omega_h::Write locate_target_cells( }); Kokkos::fence(); - pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); - auto results = search_cell(target_points); + pcms::GridPointSearch2D search_cell(source_mesh, 10, 10); // + auto results = search_cell.apply(pcms::CoordinateView( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(target_points))); auto owning_cell_ids = search_cell.GetOwningElementIds(results); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { @@ -103,8 +105,8 @@ static Omega_h::Write locate_target_cells( source_cell_ids[i] = source_cell_id; }); } else if (dim == 3) { - Kokkos::View target_points("target_points", - nvertices_target); + Kokkos::View target_points("target_points", + nvertices_target, 3); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { target_points(i, 0) = target_coords[i * dim]; @@ -114,7 +116,9 @@ static Omega_h::Write locate_target_cells( Kokkos::fence(); pcms::GridPointSearch3D search_cell(source_mesh, 10, 10, 10); - auto results = search_cell(target_points); + auto results = search_cell.apply(pcms::CoordinateView( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(target_points))); auto owning_cell_ids = search_cell.GetOwningElementIds(results); Omega_h::parallel_for( nvertices_target, OMEGA_H_LAMBDA(const Omega_h::LO i) { diff --git a/src/pcms/localization/point_localization.cpp b/src/pcms/localization/point_localization.cpp new file mode 100644 index 00000000..6b11f6ce --- /dev/null +++ b/src/pcms/localization/point_localization.cpp @@ -0,0 +1,521 @@ +#include "point_localization.h" + +// Access traits for Omega_h +// constructs bounding boxes for elements +template +struct ArborX::AccessTraits> +{ + using memory_space = typename Omega_h::ExecSpace::memory_space; + + static KOKKOS_FUNCTION int size( + const pcms::detail::Omega_h_Mesh_Adapt& mesh) + { + return mesh.size_; + } + + static KOKKOS_FUNCTION auto get( + const pcms::detail::Omega_h_Mesh_Adapt& mesh, + int i) + { + ArborX::Point min = {INFINITY}; + ArborX::Point max = {-INFINITY}; + for (int j = 0; j < dim + 1; ++j) + { + auto cell_vert_id = mesh.adjacency[(dim+1)*i+j]; + for (int k = 0; k < dim; ++k) + { + Omega_h::Real curr_coord = mesh.coordinates[cell_vert_id*dim+k]; + if (min[k] > curr_coord) min[k] = curr_coord; + if (max[k] < curr_coord) max[k] = curr_coord; + } + } + return ArborX::Box(min, max); + } +}; + +template +struct ArborX::AccessTraits> +{ + using memory_space = MemorySpace; + static KOKKOS_FUNCTION int size( + pcms::detail::Coordinate_View_Adapt const &coords) + { + return coords.points.extent(0); + } + static KOKKOS_FUNCTION auto get( + pcms::detail::Coordinate_View_Adapt const &coords, + int i) + { + ArborX::Point ax_point; + for (int j = 0; j < dim; j++) ax_point[j] = coords.points(i,j); + return PredicateWithAttachment(intersects(ax_point), i); + } +}; + +namespace pcms +{ + +namespace detail +{ + +Mapping<2>::Mapping( + Omega_h::Matrix::DIM,Mapping<2>::DIM+1> const& triangle_, + const Kokkos::View& tolerances) + : triangle(triangle_) +{ + auto tol_h = Kokkos::create_mirror_view(tolerances); + Kokkos::deep_copy(tol_h, tolerances); + tolerances_[0] = tol_h(0); + tolerances_[1] = tol_h(1); + // printf("%d (%lf, %lf)\n", tol_h.size(), tolerances_[0], tolerances_[1]); + bary_transform = { triangle[0] - triangle[2], triangle[1] - triangle[2] }; + triangle_area = 0.5*fabs(Omega_h::determinant(bary_transform)); + bary_transform = Omega_h::invert(bary_transform); +} + +KOKKOS_FUNCTION +Omega_h::Vector::DIM + 1> Mapping<2>::get_bary(Omega_h::Vector::DIM> const& p) const +{ + Omega_h::Vector<2> coeffs = bary_transform*(p - triangle[2]); + return {coeffs[0], coeffs[1], 1 - coeffs[0] - coeffs[1]}; +} + +KOKKOS_FUNCTION +int Mapping<2>::which(int dim, + Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + if (dim == Omega_h::VERT) { + return which_vert(bary_coords); + } + if (dim == Omega_h::EDGE) { + return which_edge(bary_coords); + } + if (dim == Omega_h::FACE) { + return within_elem(bary_coords); + } + return -1; +} + +int Mapping<2>::which_vert(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + for (int i = 0; i < 3; i++) + { + // distance of the point to vertex i + Omega_h::Vector<2> error = (bary_coords[0] - kronecker(i,0)) * triangle[0] + + (bary_coords[1] - kronecker(i,1)) * triangle[1] + + (bary_coords[2] - kronecker(i,2)) * triangle[2]; + if (Omega_h::norm_squared(error) <= tolerances_[0]*tolerances_[0]) return i; + } + return -1; +} + +int Mapping<2>::which_edge(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + for (int i = 0; i <= 2; i++) + { + // Distance of the point to the edge opposite the ith vertex, + // see docs for derivation. + double dist = (2*bary_coords[i]*triangle_area)*(2*bary_coords[i]*triangle_area); + dist /= opposite_edge_len_sq(i); + if (dist <= tolerances_[1]*tolerances_[1] && + bary_coords[(i+1)%3] >= 0 && bary_coords[(i+2)%3] >= 0) return (i+1)%3; + } + return -1; +} + +int Mapping<2>::within_elem(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + return -1 * (int)!(bary_coords[0] > 0 && bary_coords[1] > 0 && bary_coords[2] > 0); +} + +double Mapping<2>::opposite_edge_len_sq(int i) const +{ + return Omega_h::norm_squared(triangle[(i+2)%3] - triangle[(i+1)%3]); +} + +Mapping<3>::Mapping( + Omega_h::Matrix::DIM,Mapping<3>::DIM+1> const& tetrahedron_, + const Kokkos::View& tolerances) + : tetrahedron(tetrahedron_) +{ + auto tol_h = Kokkos::create_mirror_view(tolerances); + Kokkos::deep_copy(tol_h, tolerances); + tolerances_[0] = tol_h(0); + tolerances_[1] = tol_h(1); + tolerances_[2] = tol_h(2); + set_triangle_areas(); + bary_transform = { tetrahedron[0] - tetrahedron[3], tetrahedron[1] - tetrahedron[3], tetrahedron[2] - tetrahedron[3] }; + tetrahedron_volume = fabs(Omega_h::determinant(bary_transform))/6.; + bary_transform = Omega_h::invert(bary_transform); +} + +KOKKOS_FUNCTION +Omega_h::Vector::DIM + 1> Mapping<3>::get_bary(Omega_h::Vector::DIM> const& p) const +{ + Omega_h::Vector<3> coeffs = bary_transform*(p - tetrahedron[3]); + return {coeffs[0], coeffs[1], coeffs[2], 1 - coeffs[0] - coeffs[1] - coeffs[2]}; +} + +KOKKOS_FUNCTION +int Mapping<3>::which(int dim, + Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + if (dim == Omega_h::VERT) { + return which_vert(bary_coords); + } + if (dim == Omega_h::EDGE) { + return which_edge(bary_coords); + } + if (dim == Omega_h::FACE) { + return which_face(bary_coords); + } + if (dim == Omega_h::REGION) { + return within_elem(bary_coords); + } + return -1; +} + +int Mapping<3>::which_vert(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + for (int i = 0; i < 4; i++) + { + // distance to the ith vertex + Omega_h::Vector<3> error = (bary_coords[0] - kronecker(i,0)) * tetrahedron[0] + (bary_coords[1] - kronecker(i,1)) * tetrahedron[1] + + (bary_coords[2] - kronecker(i,2)) * tetrahedron[2] + (bary_coords[3] - kronecker(i,3)) * tetrahedron[3]; + if (Omega_h::norm_squared(error) <= tolerances_[0]*tolerances_[0]) return i; + } + return -1; +} + + +int Mapping<3>::which_edge(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + int edge_ = 0; + for (int i = 0; i < 4; i++) + { + for (int j = i+1; j < 4; j++) + { + Omega_h::Vector<3> side1 = {0,0,0}, side2 = {0,0,0}, side3 = tetrahedron[i] - tetrahedron[j]; + for (int k = 0; k < 4; k++) + { + side1 += (bary_coords[k] - kronecker(i,k)) * tetrahedron[k]; + side2 += (bary_coords[k] - kronecker(j,k)) * tetrahedron[k]; + } + + // the norm of the cross product of two vectors is twice the area of the triangle + // those vectors form + double distnce_sq = Omega_h::norm_squared(Omega_h::cross(side1, side2)); + distnce_sq /= Omega_h::norm_squared(side3); + + if (distnce_sq <= tolerances_[1]*tolerances_[1] && bary_coords[i] > 0 && bary_coords[j] > 0) return edge(edge_); + edge_++; + } + } + return -1; +} + +int Mapping<3>::which_face(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + for (int i = 0; i < 4; i++) + { + if (fabs(3*tetrahedron_volume*bary_coords[i]/face_areas[i]) <= tolerances_[2] + && bary_coords[(i+1)%4] >= 0 && bary_coords[(i+2)%4] >= 0 && bary_coords[(i+3)%4] >= 0) + { + return face(i); + } + } + return -1; +} + +int Mapping<3>::within_elem(Omega_h::Vector::DIM + 1> const& bary_coords) const +{ + return -1 * (int)!(bary_coords[0] >= 0 && bary_coords[1] >= 0 && bary_coords[2] >= 0 && bary_coords[3] >= 0); +} + +// calculates the areas of each of the faces of the triangle +// the faces are "named" according to the vertex opposite +// (i.e., face 0 is defined by points 1, 2, 3) +void Mapping<3>::set_triangle_areas() +{ + for (int i = 0; i < 4; i++) + { + Omega_h::Vector<3> edge0 = tetrahedron[(i+1)%4] - tetrahedron[(i+3)%4], edge1 = tetrahedron[(i+2)%4] - tetrahedron[(i+3)%4]; + Omega_h::Vector<3> cross = Omega_h::cross(edge0, edge1); + face_areas[i] = 0.5*Omega_h::norm(cross); + } +} + +template +KOKKOS_FUNCTION void CallOnIntersect2D::operator()(Predicate const &predicate, Value const & val) const +{ + ArborX::Point const& ax = ArborX::getGeometry(predicate); + Omega_h::Vector point{ax[0], ax[1]}; + int point_ind = ArborX::getData(predicate); + + detail::Mapping<2> const& tm = mappings(val.index); + + // calculate the barycentric coefficients of the point + auto coeffs = tm.get_bary(point); + + for (int i = 0; i < DIM; i++) + { + int elem = tm.which(i, coeffs); + if (elem >= 0 && dimensionalities(point_ind) > (TreePointSearch::Dimensionality)i) + { + auto elem_ind = adjacencies[i][3*val.index + elem]; + dimensionalities(point_ind) = (TreePointSearch::Dimensionality)i; + element_ids(point_ind) = elem_ind; + for (int j = 0; j < DIM + 1; j++) + { + parametric_coords(point_ind, j) = coeffs(j); + } + return; + } + } + if (tm.which(DIM, coeffs) >= 0 + && dimensionalities(point_ind) > TreePointSearch::Dimensionality::FACE) + { + dimensionalities(point_ind) = TreePointSearch::Dimensionality::FACE; + element_ids(point_ind) = (LO)val.index; + for (int j = 0; j < DIM + 1; j++) + { + parametric_coords(point_ind, j) = coeffs(j); + } + } +} + +template +void CallOnIntersect3D::operator()(Predicate const &predicate, Value const & val) const +{ + ArborX::Point const& ax = ArborX::getGeometry(predicate); + Omega_h::Vector point{ax[0], ax[1], ax[2]}; + int point_ind = ArborX::getData(predicate); + + detail::Mapping<3> const& tm = mappings(val.index); + + // calculate the barycentric coefficients of the point + auto coeffs = tm.get_bary(point); + int offsets[DIM] = {4, 6, 4}; + for (int i = 0; i < DIM; i++) + { + int elem = tm.which(i, coeffs); + if (elem >= 0 && dimensionalities(point_ind) > (TreePointSearch::Dimensionality)i) + { + auto elem_ind = adjacencies[i][offsets[i]*val.index + elem]; + dimensionalities(point_ind) = (TreePointSearch::Dimensionality)i; + element_ids(point_ind) = elem_ind; + for (int j = 0; j < DIM + 1; j++) + { + parametric_coords(point_ind, j) = coeffs(j); + } + return; + } + } + + if (tm.which(DIM, coeffs) >= 0 + && dimensionalities(point_ind) > TreePointSearch::Dimensionality::REGION) + { + dimensionalities(point_ind) = TreePointSearch::Dimensionality::REGION; + element_ids(point_ind) = (LO)val.index; + for (int j = 0; j < DIM + 1; j++) + { + parametric_coords(point_ind, j) = coeffs(j); + } + } +} + +} //namespace detail + +TreePointSearch::Results TreePointSearch::apply( + const CoordinateView& coords) const +{ + if (coords.GetCoordinateSystem() != pcms::CoordinateSystem::Cartesian) + { + throw pcms_error("TreePointSearch::apply only implemented for" + " Cartesian coordinates"); + } + if (coords.GetValues().extent(1) != mesh_.dim()) + { + throw pcms_error("Input coordinate space dimension " + + std::to_string(coords.GetValues().extent(1)) + + " does not match query space dimension " + + std::to_string(mesh_.dim())); + } + + if (mesh_.dim() == 2) + { + static constexpr int DIM = 2; + Omega_h::ExecSpace execution_space; + + Kokkos::View dims("dimensionalities", coords.GetValues().extent(0)); + Kokkos::deep_copy(dims, TreePointSearch::Dimensionality::NO_INTERSECT); + + + Kokkos::View elem_ids("element IDs", coords.GetValues().extent(0)); + Kokkos::deep_copy(elem_ids, -1); + + Kokkos::View parametric_coords("parametric coordinates", coords.GetValues().extent(0), coords.GetValues().extent(1) + 1); + Kokkos::deep_copy(parametric_coords, -1.0); + + + tree->get_tree<2>().query( + execution_space, + detail::Coordinate_View_Adapt{coords.GetValues()}, + detail::CallOnIntersect2D( + tree->get_mappings<2>(), + mesh_.get_adj(Omega_h::FACE, 0).ab2b, + mesh_.get_adj(Omega_h::FACE, 1).ab2b, + dims, + elem_ids, + parametric_coords + )); + return Results{dims, elem_ids, parametric_coords}; + } + else + { + static constexpr int DIM = 3; + Omega_h::ExecSpace execution_space; + + Kokkos::View dims("dimensionalities", coords.GetValues().extent(0)); + Kokkos::deep_copy(dims, TreePointSearch::Dimensionality::NO_INTERSECT); + + + Kokkos::View elem_ids("element IDs", coords.GetValues().extent(0)); + Kokkos::deep_copy(elem_ids, -1); + + Kokkos::View parametric_coords("parametric coordinates", coords.GetValues().extent(0), coords.GetValues().extent(1) + 1); + Kokkos::deep_copy(parametric_coords, -1.0); + + tree->get_tree<3>().query( + execution_space, + detail::Coordinate_View_Adapt{coords.GetValues()}, + detail::CallOnIntersect3D( + tree->get_mappings<3>(), + mesh_.get_adj(Omega_h::REGION, 0).ab2b, + mesh_.get_adj(Omega_h::REGION, 1).ab2b, + mesh_.get_adj(Omega_h::REGION, 2).ab2b, + dims, + elem_ids, + parametric_coords + )); + return Results{dims, elem_ids, parametric_coords}; + } +} + +[[nodiscard]] LO TreePointSearch::GetOwningElementId( + const TreePointSearch::Results& results, int i) +{ + const Kokkos::View query_id{""}; + const Kokkos::View dim{""}; + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int){ + query_id(0) = (results.element_ids(i) < 0) ? -results.element_ids(i) : results.element_ids(i); + dim(0) = results.dimensionalities(i); + }); + auto query_id_h = Kokkos::create_mirror_view(query_id); + auto dim_h = Kokkos::create_mirror_view(dim); + Kokkos::deep_copy(query_id_h, query_id); + Kokkos::deep_copy(dim_h, dim); + return pcms::GetOwningElementId(mesh_, mesh_.dim(), static_cast(dim_h(0)), + query_id_h(0)); +} + +[[nodiscard]] Kokkos::View TreePointSearch::GetOwningElementIds( + const TreePointSearch::Results& results) +{ + Kokkos::View owning_ids("Owning element IDs", + results.dimensionalities.size()); + auto vert2elem = mesh_.ask_up(Omega_h::VERT, mesh_.dim()); + auto edge2elem = mesh_.ask_up(Omega_h::EDGE, mesh_.dim()); + auto face2elem = (mesh_.dim() == 3) ? mesh_.ask_up(Omega_h::FACE, mesh_.dim()) : Omega_h::Adj{}; + auto mesh_dim = mesh_.dim(); + Kokkos::parallel_for(owning_ids.size(), KOKKOS_LAMBDA(const int i){ + if (static_cast(results.dimensionalities(i)) == mesh_dim) { + owning_ids(i) = abs(results.element_ids(i)); + } else if (static_cast(results.dimensionalities(i)) == 2) { + owning_ids(i) = GetOwningElementIdFromAdj(face2elem, 2, mesh_dim, results.element_ids(i)); + } else if (static_cast(results.dimensionalities(i)) == 1) { + owning_ids(i) = GetOwningElementIdFromAdj(edge2elem, 1, mesh_dim, results.element_ids(i)); + } else if (static_cast(results.dimensionalities(i)) == 0) { + owning_ids(i) = GetOwningElementIdFromAdj(vert2elem, 0, mesh_dim, results.element_ids(i)); + } else { + owning_ids(i) = -1; + } + }); + return owning_ids; +} + +std::unique_ptr TreePointSearch::make_tree(const Omega_h::Mesh& mesh) const +{ + ExecSpace execution_space; + using DeviceType = Kokkos::Device; + + if (mesh.dim() == 2) + { + detail::Omega_h_Mesh_Adapt<2> tagged_mesh{ + mesh.nelems(), + mesh.get_adj(Omega_h::FACE, Omega_h::VERT).ab2b, + mesh.coords() + }; + + detail::TreeWrapper::Tree_t<2> tree = detail::TreeWrapper::Tree_t<2>( + execution_space, + ArborX::Experimental::attach_indices(tagged_mesh)); + + detail::TreeWrapper::Mappings_t<2> mappings = detail::TreeWrapper::Mappings_t<2>( + "mappings", + mesh.nelems()); + + auto face2vert = Omega_h::HostRead(mesh.get_adj(Omega_h::FACE, Omega_h::VERT).ab2b); + auto vert_coords = Omega_h::HostRead(mesh.coords()); + + auto mappings_h = Kokkos::create_mirror_view(mappings); + for (int i = 0; i < mesh.nelems(); i++) + { + Omega_h::Matrix<2,3> triangle = { + {vert_coords[face2vert[i*3]*2], vert_coords[face2vert[i*3]*2 + 1]}, + {vert_coords[face2vert[i*3 + 1]*2], vert_coords[face2vert[i*3 + 1]*2 + 1]}, + {vert_coords[face2vert[i*3 + 2]*2], vert_coords[face2vert[i*3 + 2]*2 + 1]}}; + mappings_h[i] = detail::Mapping<2>(triangle, tolerances_); + } + Kokkos::deep_copy(execution_space, mappings, mappings_h); + return std::make_unique(mappings, tree); + } + if (mesh.dim() == 3) + { + detail::Omega_h_Mesh_Adapt<3> tagged_mesh{ + mesh.nelems(), + mesh.get_adj(Omega_h::REGION, Omega_h::VERT).ab2b, + mesh.coords() + }; + detail::TreeWrapper::Tree_t<3> tree = detail::TreeWrapper::Tree_t<3>( + execution_space, + ArborX::Experimental::attach_indices(tagged_mesh)); + + detail::TreeWrapper::Mappings_t<3> mappings = detail::TreeWrapper::Mappings_t<3>( + "mappings", + mesh.nelems()); + + auto region2vert = Omega_h::HostRead(mesh.get_adj(Omega_h::REGION, Omega_h::VERT).ab2b); + auto vert_coords = Omega_h::HostRead(mesh.coords()); + + auto mappings_h = Kokkos::create_mirror_view(mappings); + for (int i = 0; i < mesh.nelems(); i++) + { + Omega_h::Matrix<3,4> tetrahedron = { + {vert_coords[region2vert[i*4]*3], vert_coords[region2vert[i*4]*3 + 1], vert_coords[region2vert[i*4]*3 + 2]}, + {vert_coords[region2vert[i*4 + 1]*3], vert_coords[region2vert[i*4 + 1]*3 + 1], vert_coords[region2vert[i*4 + 1]*3 + 2]}, + {vert_coords[region2vert[i*4 + 2]*3], vert_coords[region2vert[i*4 + 2]*3 + 1], vert_coords[region2vert[i*4 + 2]*3 + 2]}, + {vert_coords[region2vert[i*4 + 3]*3], vert_coords[region2vert[i*4 + 3]*3 + 1], vert_coords[region2vert[i*4 + 3]*3 + 2]}}; + + auto tol_h = Kokkos::create_mirror_view(tolerances_); + Kokkos::deep_copy(tol_h, tolerances_); + mappings_h[i] = detail::Mapping<3>(tetrahedron, tolerances_); + } + Kokkos::deep_copy(execution_space, mappings, mappings_h); + return std::make_unique(mappings, tree); + } + throw pcms_error("Invalid mesh dimension " + std::to_string(mesh.dim()) + ", TreePointSearch only implemented for 2D and 3D"); +} + +} // namespace pcms diff --git a/src/pcms/localization/point_localization.h b/src/pcms/localization/point_localization.h new file mode 100644 index 00000000..1125f396 --- /dev/null +++ b/src/pcms/localization/point_localization.h @@ -0,0 +1,479 @@ +#ifndef POINT_LOCALIZATION_H +#define POINT_LOCALIZATION_H + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "pcms/utility/assert.h" +#include "pcms/utility/types.h" +#include "pcms/utility/arrays.h" +#include "pcms/localization/point_search.h" + +namespace pcms +{ + +namespace detail +{ + +/** +* @brief Implementation of the Kronecker delta: +* https://en.wikipedia.org/wiki/Kronecker_delta +* @param i an integer +* @param j another integer +* @returns 1 if i==j, 0 if i i!= j +*/ +KOKKOS_INLINE_FUNCTION +LO kronecker(LO i, LO j) { return (LO)(i==j); } + +/** +* @brief Omega_h::Mesh tagged with the templated dimension required for +* compatability with ArborX::BVH +*/ +template +struct Omega_h_Mesh_Adapt +{ + const Omega_h::LO size_; + const Omega_h::LOs adjacency; + const Omega_h::Reals coordinates; +}; + +/** +* @brief pcms CoordinateView tagged with the templated dimension and +* memory space required for compatability with ArborX::BVH +*/ +template +struct Coordinate_View_Adapt +{ + const Rank2View> points; +}; + +template +class Mapping {}; + +/** +* @brief Mapping<2> represents a mapping from the global coordinate system +* to the barycentric coordinate system of an element in a 2D Omega_h::Mesh +* The class Mapping<2> calculates the mapping of any point in global coordinate +* space to the barycentric coordinate space belonging to a triangle (face) in +* an Omega_h::Mesh. +* An instance of Mapping<2> can, from the barycentric coordinates constructed by it, +* determine whether that point intersects the face, an edge, a vertex, or lies +* outside the face. It also determines which edge or vertex the point intersects +* by calculating the offset in the Omega_h::Adj::ab2b +*/ +template <> +class Mapping<2> +{ +public: + // the dimension of the space + static constexpr int DIM = 2; + // CONSTRUCTORS + /** + * @brief Default constructor + */ + KOKKOS_FUNCTION + Mapping() = default; + /** + * @brief Constructs a barycentric mapping of a triangle in an Omega_h::Mesh + * @param elem_index the index of the desired triangle or tetrahedron in the Omega_h::Mesh + * @param mesh the Omega_h::Mesh with the spatial information of the triangle + */ + + Mapping(Omega_h::Matrix const& triangle_, + Kokkos::View const& tolerances); + /** + * @brief Default destructor + */ + KOKKOS_FUNCTION + ~Mapping() = default; + /** + * @brief Computes the barycentric coordinates of a point in global space + * + * This function computes the barycentric coordinate of a point with respect + * to a triangle in an Omega_h::Mesh + * + * @param p the point to compute the barycentric coordinates of + * @returns an Omega_h::Vector<3> containing the barycentric coordinates + * of p + */ + KOKKOS_FUNCTION + Omega_h::Vector get_bary(Omega_h::Vector const& p) const; + /** + * @brief Determines which entity of a dimension `ent_dim` a point intersects + * from the barycentric coordinates + * + * If the point corresponding to the input barycentric coordinates intersects + * any entities with dimension `ent_dim` bordering the element this mapping is + * constructed from, this function determines the offset in the + * Omega_h::Adj::ab2b structure coresponding to FACE -> `ent_dim` adjacency. + * If the point does not intersect a border entity, this function returns -1 + * + * @param ent_dim the dimension of the entities we want to check intersection + * with + * @param bary_coords the barycentric coordinates computed by this mapping of + * a point in global space + * @returns -1 if the point does not intersect any entities, or the offset of + * the intersected entity in the Omega_h::Adj::ab2b structure + * coresponding to FACE -> `ent_dim` adjacency + */ + KOKKOS_FUNCTION + int which( + int ent_dim, + Omega_h::Vector const& bary_coords) const; +private: + /** + * @brief Computes the vertex offset of the point corresponding + * to the input barycentric coordinates + * @param bary_coords the input barycentric coordinates + * @returns the vertex offset if the point corresponding to the given bary. coordinates + * is within a certain (global) tolerance of a vertex + * -1 if the point does not lie within the tolerance of any vertex + */ + KOKKOS_FUNCTION + int which_vert(Omega_h::Vector const& bary_coords) const; + /** + * @brief Computes the edge offset of the point corresponding + * to the input barycentric coordinates + * @param bary_coords the input barycentric coordinates + * @returns the edge offset if the distnce from the point corresponding to the + * given bary. coordinates to the edge with the respective offset + * is within a certain (global) tolerance + * -1 if the point does not lie within the tolerance of any vertex + */ + KOKKOS_FUNCTION + int which_edge(Omega_h::Vector const& bary_coords) const; + /** + * @brief Computes whether a point with the given barycentric coordinates + * is within a triangle + * @param bary_coords the input barycentric coordinates + * @returns 0 if the point is within the highest order element and -1 otherwise + */ + KOKKOS_FUNCTION + int within_elem(Omega_h::Vector const& bary_coords) const; + /** + * @brief Constructs the triangle the mapping corresponds to + * @param index the element ID of the triangle + * @param mesh the Omega_h::Mesh with the spatial information for the triangle + * corresponding to index + */ + void set_mesh_triangle(int index, Omega_h::Mesh const& mesh); + /** + * @brief Calculates the length of the edge opposite vertex i + * @param i the vertex index opposite the desired edge + * @returns the length of the ith edge + */ + KOKKOS_FUNCTION + double opposite_edge_len_sq(int i) const; + + // REPRESENTATION + Omega_h::Vector tolerances_; + // representation of the barycentric coordinate mapping, source: + // https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Edge_approach + Omega_h::Matrix bary_transform; // Column-major order + // Points defining the triangle + Omega_h::Matrix triangle; + // area of the triangle + double triangle_area; +}; + +template <> +class Mapping<3> +{ +public: + static constexpr int DIM = 3; + /** + * @brief Default constructor + */ + KOKKOS_FUNCTION + Mapping() = default; + + Mapping(Omega_h::Matrix const& tetrahedron_, + const Kokkos::View& tolerances); + /** + * @brief Default destructor + */ + KOKKOS_FUNCTION + ~Mapping() = default; + /** + * @brief Computes the barycentric coordinates of a point in global space + * @param p the point to compute the barycentric coordinates of + * @returns an Omega_h::Vector containing the barycentric coordinates + * of p + */ + KOKKOS_FUNCTION + Omega_h::Vector get_bary(Omega_h::Vector const& p) const; + + /** + * @brief Determines which entity of dimension `ent_dim` a point intersects + * from the barycentric coordinates + * + * If the point corresponding to the input barycentric coordinates intersects + * any entities with dimension `ent_dim` bordering the element this mapping is + * constructed from, this function determines the offset in the + * Omega_h::Adj::ab2b structure coresponding to REGION -> `ent_dim` adjacency. + * If the point does not intersect a border entity, this function returns -1 + * + * @param ent_dim the dimension of the entities we want to check intersection + * with + * @param bary_coords the barycentric coordinates computed by this mapping of + * a point in global space + * @returns -1 if the point does not intersect any entities, or the offset of + * the intersected entity in the Omega_h::Adj::ab2b structure + * coresponding to REGION -> `ent_dim` adjacency + */ + KOKKOS_FUNCTION + int which(int ent_dim, + Omega_h::Vector const& bary_coords) const; +private: + /** + * @brief Computes the vertex offset of the point corresponding + * to the input barycentric coordinates + * @param bary_coords the input barycentric coordinates + * @returns the vertex offset if the point corresponding to the given bary. coordinates + * is within a certain (global) tolerance of a vertex + * -1 if the point does not lie within the tolerance of any vertex + */ + KOKKOS_FUNCTION + int which_vert(Omega_h::Vector const& bary_coords) const; + /** + * @brief Computes the edge offset of the point corresponding + * to the input barycentric coordinates + * @param bary_coords the input barycentric coordinates + * @returns the edge offset if the distnce from the point corresponding to the + * given bary. coordinates to the edge with the respective offset + * is within a certain (global) tolerance + * -1 if the point does not lie within the tolerance of any vertex + */ + KOKKOS_FUNCTION + int which_edge(Omega_h::Vector const& bary_coords) const; + /** + * @brief Computes the face offset of the point corresponding + * to the input barycentric coordinates + * @param bary_coords the input barycentric coordinates + * @returns if the point corresponding to the bary. coords is within a + * global tolerance of a face, returns the offset of that face + * -1 if the point does not lie within the tolerance of any face + * Algorithm source: + * C.E. Passerello, + * Interference detection using barycentric coordinates, + * Mechanics Research Communications, + * Volume 9, Issue 6, 1982, Pages 373-378, + * https://doi.org/10.1016/0093-6413(82)90034-9. + */ + KOKKOS_FUNCTION + int which_face(Omega_h::Vector const& bary_coords) const; + /** + * @brief Computes whether a point with the given barycentric coordinates + * is within a triangle + * @param bary_coords the input barycentric coordinates + * @returns 0 if the point is within the highet order element and -1 otherwise + */ + KOKKOS_FUNCTION + int within_elem(Omega_h::Vector const& bary_coords) const; + /** + * @brief Constructs the tetrahedron the mapping corresponds to + * @param index the element ID of the tetrahedron + * @param mesh the Omega_h::Mesh with the spatial information for the + * tetrahedron corresponding to `index` + */ + void set_mesh_tet(int index, Omega_h::Mesh const& mesh); + /** + * @brief Calculates the areas of each face in the tetrahedron the mapping + * corresponds to + * @param index the element ID of the tetrahedron + * @param mesh the Omega_h::Mesh with the spatial information for the + * tetrahedron corresponding to `index` + */ + void set_triangle_areas(); + // returns the actual face/edge offset from the barycentric offset + KOKKOS_INLINE_FUNCTION int face(int i) const + { return (OFFSETS & 3 << (i*2))>>(i*2); } + KOKKOS_INLINE_FUNCTION int edge(int i) const + { return (i < 4 && i > 0) ? (OFFSETS & 3 << ((i-1)*2))>>((i-1)*2) : i; } + // representation + Omega_h::Vector tolerances_; + Omega_h::Matrix bary_transform; // Column-major order + Omega_h::Matrix tetrahedron; + Omega_h::Vector face_areas; + double tetrahedron_volume; + // This is used to calculate the face and edge offsets in the cannonical ordering: + // https://user-images.githubusercontent.com/56453280/74203616-39d27a80-4c3e-11ea-885d-b0260490e184.png + static const char OFFSETS = 1 << 4 | 3 << 2 | 2; +}; + +/** +* Wrapper base class for the ArborX::BVH and Mappings classes, this is a *rough* +* workaround for ArborX::BVH being templated on dimension +*/ +struct TreeWrapper +{ + template + using Mappings_t = Kokkos::View*, Omega_h::ExecSpace::memory_space>; + template + using Tree_t = ArborX::BVH, unsigned>>; + + template + TreeWrapper(const Mappings_t& mappings, const Tree_t& tree) + { + dim_ = dim; + mappings_ = std::make_shared(mappings); + tree_ = std::make_shared(tree); + } + /** + * @brief returns a pointer to an ArborX tree + */ + template + inline Tree_t get_tree() + { + if (dim != dim_) + { + throw pcms_error("Requested get_tree return " + "type does not match internal Tree_t"); + } + return std::any_cast>(*tree_); + } + + /** + * @brief returns a pointer to a Kokkos::View of Mappings + */ + template + inline Mappings_t get_mappings() + { + if (dim != dim_) + { + throw pcms_error("Requested get_mappings return " + "type does not match internal Mappings_t"); + } + return std::any_cast>(*mappings_); + }; +private: + int dim_; + std::shared_ptr tree_; + std::shared_ptr mappings_; +}; +} //namespace detail + +class TreePointSearch : public PointSearch +{ +public: + using Results = PointSearch::Results; + using Dimensionality = Results::Dimensionality; + using ExecSpace = PointSearch::ExecSpace; + using MemorySpace = PointSearch::MemorySpace; + + TreePointSearch(Omega_h::Mesh& mesh) : + PointSearch(PointSearchTolerances("tree point search tolerances", mesh.dim())), + mesh_(mesh) + { + Kokkos::deep_copy(tolerances_, 1e-12); + tree = make_tree(mesh); + } + TreePointSearch(Omega_h::Mesh& mesh, const PointSearchTolerances& tolerances) : PointSearch(tolerances), mesh_(mesh), tree(make_tree(mesh)) {} + ~TreePointSearch() = default; + /** + * Given a set of points in global coordinates give the ids of the entities + * that the points lie within and the parametric coordinates of each point + * within a triangles or tetrahedra adjacent to the intersected entity. + * If the point does not lie within any triangle element, the id will + * be a negative number. + */ + Results apply(const CoordinateView& coords) const override; + [[nodiscard]] LO GetOwningElementId(const Results& results, int i) override; + [[nodiscard]] Kokkos::View GetOwningElementIds( + const Results& results) override; +private: + std::unique_ptr make_tree(const Omega_h::Mesh& mesh) const; + // Reference to the input mesh + Omega_h::Mesh &mesh_; + std::unique_ptr tree; +}; + +namespace detail +{ + + /** + * Functor to be called when a point intersects a triangle or tetrahedron + */ + class CallOnIntersect3D + { + public: + using MemorySpace = TreePointSearch::MemorySpace; + static constexpr int DIM = 3; + + CallOnIntersect3D( + const Kokkos::View*, MemorySpace>& mappings_, + const Omega_h::LOs& adjacencies0, + const Omega_h::LOs& adjacencies1, + const Omega_h::LOs& adjacencies2, + Kokkos::View& dimensionalities_, + Kokkos::View& element_ids_, + Kokkos::View& parametric_coords_ + ) : mappings(MakeRank1View(mappings_)), + adjacencies{adjacencies0, adjacencies1, adjacencies2}, + dimensionalities(MakeRank1View(dimensionalities_)), + element_ids(MakeRank1View(element_ids_)), + parametric_coords(MakeRank2View(parametric_coords_)) + { + }; + + /** + * Intersection callback, + */ + template + KOKKOS_FUNCTION void operator()(Predicate const &predicate, Value const & val) const; + private: + Rank1View, MemorySpace> mappings; + Omega_h::LOs adjacencies[3]; + Rank1View dimensionalities; + Rank1View element_ids; + Rank2View parametric_coords; + }; + + + class CallOnIntersect2D + { + public: + using MemorySpace = TreePointSearch::MemorySpace; + static constexpr int DIM = 2; + + CallOnIntersect2D( + const Kokkos::View*, MemorySpace>& mappings_, + const Omega_h::LOs& adjacencies0, + const Omega_h::LOs& adjacencies1, + Kokkos::View& dimensionalities_, + Kokkos::View& element_ids_, + Kokkos::View& parametric_coords_ + ) : mappings(MakeRank1View(mappings_)), + adjacencies{adjacencies0, adjacencies1}, + dimensionalities(MakeRank1View(dimensionalities_)), + element_ids(MakeRank1View(element_ids_)), + parametric_coords(MakeRank2View(parametric_coords_)) + { + }; + + template + KOKKOS_FUNCTION void operator()(Predicate const &predicate, Value const & val) const; + private: + Rank1View, MemorySpace> mappings; + Omega_h::LOs adjacencies[2]; + Rank1View dimensionalities; + Rank1View element_ids; + Rank2View parametric_coords; + }; + +} // namespace detail + +} // namespace pcms +#endif // POINT_LOCALIZATION_H diff --git a/src/pcms/localization/point_search.cpp b/src/pcms/localization/point_search.cpp index f7148107..844ae520 100644 --- a/src/pcms/localization/point_search.cpp +++ b/src/pcms/localization/point_search.cpp @@ -396,11 +396,16 @@ OMEGA_H_INLINE double myreduce(const Omega_h::Vector& x, return out; } -Kokkos::View GridPointSearch2D::operator()( - Kokkos::View points) const +GridPointSearch2D::Results GridPointSearch2D::apply( + const CoordinateView& points) const { - Kokkos::View results("point search result", - points.extent(0)); + Results results{ + Kokkos::View("dimensionalities", points.GetValues().extent(0)), + Kokkos::View("element IDs", points.GetValues().extent(0)), + Kokkos::View("parametric coordinates", points.GetValues().extent(0), + points.GetValues().extent(1) + 1) + }; + 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_; @@ -411,10 +416,14 @@ Kokkos::View GridPointSearch2D::operator()( auto edges2verts_adj = edges2verts_adj_; auto coords = coords_; auto tolerances = tolerances_; + auto point_coords = points.GetValues(); Kokkos::parallel_for( - points.extent(0), KOKKOS_LAMBDA(int p) { + point_coords.extent(0), KOKKOS_LAMBDA(int p) { + const auto vtol = tolerances(0); + const auto etol = tolerances(1); + Omega_h::Vector<2> point( - std::initializer_list{points(p, 0), points(p, 1)}); + std::initializer_list{point_coords(p, 0), point_coords(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); @@ -471,14 +480,14 @@ Kokkos::View GridPointSearch2D::operator()( const int vb_id = edges2verts_adj.ab2b[edgeID * 2 + 1]; const auto va = Omega_h::get_vector<2>(coords, va_id); const auto vb = Omega_h::get_vector<2>(coords, vb_id); - + if (!normal_intersects_segment(va, vb, point)) continue; const auto de = distance_from_line(va, vb, point); if ((de < best_edge_dist) || - ((de == best_edge_dist) && (edgeID < best_edge_id)) || - ((de == best_edge_dist) && (edgeID == best_edge_id) && + ((fabs(de - best_edge_dist) < etol) && (edgeID < best_edge_id)) || + ((fabs(de - best_edge_dist) < etol) && (edgeID == best_edge_id) && (triangleID < best_edge_tid_min))) { best_edge_dist = de; best_edge_id = edgeID; @@ -486,7 +495,7 @@ Kokkos::View GridPointSearch2D::operator()( best_edge_tid_max = triangleID; best_edge_bary_min = parametric_coords; best_edge_bary_max = parametric_coords; - } else if ((de == best_edge_dist) && (edgeID == best_edge_id)) { + } else if ((fabs(de - best_edge_dist) < etol) && (edgeID == best_edge_id)) { // Track the extents of face IDs sharing this nearest edge if (triangleID < best_edge_tid_min) { best_edge_tid_min = triangleID; @@ -508,10 +517,6 @@ Kokkos::View GridPointSearch2D::operator()( } } } - - const auto vtol = tolerances(0); - const auto etol = tolerances(1); - // If we found an inside face, compute its edge distance using // barycentric-only helper to check for edge classification while // preserving the containing face ID. @@ -545,8 +550,9 @@ Kokkos::View GridPointSearch2D::operator()( } } - GridPointSearch2D::Result::Dimensionality dim_out = - GridPointSearch2D::Result::Dimensionality::REGION; + // found_inside && inside_edge_dist <= etol) + GridPointSearch2D::Results::Dimensionality dim_out = + GridPointSearch2D::Results::Dimensionality::NO_INTERSECT; LO element_id_out = -1; Omega_h::Vector<3> bary_out{0.0, 0.0, 0.0}; @@ -555,13 +561,13 @@ Kokkos::View GridPointSearch2D::operator()( // Only points with no candidates at all get negative IDs if (best_vertex_id >= 0 && best_vertex_dist <= vtol) { // Point within vertex tolerance - still inside the mesh - dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; + dim_out = GridPointSearch2D::Results::Dimensionality::VERTEX; element_id_out = best_vertex_id; bary_out = best_vertex_bary; } else if ((best_edge_id >= 0 && best_edge_dist <= etol) || (found_inside && inside_edge_dist <= etol)) { // Point within edge tolerance - still inside the mesh - dim_out = GridPointSearch2D::Result::Dimensionality::EDGE; + dim_out = GridPointSearch2D::Results::Dimensionality::EDGE; if (found_inside && inside_edge_dist <= etol && inside_edge_argmin >= 0) { element_id_out = inside_edge_id; @@ -572,7 +578,7 @@ Kokkos::View GridPointSearch2D::operator()( } } else if (found_inside) { // Point inside face - dim_out = GridPointSearch2D::Result::Dimensionality::FACE; + dim_out = GridPointSearch2D::Results::Dimensionality::FACE; element_id_out = inside_face_id; bary_out = inside_face_bary; } else { @@ -580,12 +586,12 @@ Kokkos::View GridPointSearch2D::operator()( // Only negate if we truly have no candidates at all. if (best_vertex_id >= 0 && (best_vertex_dist <= best_edge_dist || best_edge_id < 0)) { - dim_out = GridPointSearch2D::Result::Dimensionality::VERTEX; + dim_out = GridPointSearch2D::Results::Dimensionality::VERTEX; element_id_out = best_vertex_id; bary_out = best_vertex_bary; element_id_out = -element_id_out; } else if (best_edge_id >= 0) { - dim_out = GridPointSearch2D::Result::Dimensionality::EDGE; + dim_out = GridPointSearch2D::Results::Dimensionality::EDGE; element_id_out = best_edge_id; bary_out = best_edge_bary_min; element_id_out = -element_id_out; @@ -594,7 +600,10 @@ Kokkos::View GridPointSearch2D::operator()( } } - results(p) = GridPointSearch2D::Result{dim_out, element_id_out, bary_out}; + results.dimensionalities(p) = dim_out; + results.element_ids(p) = element_id_out; + for (int j = 0; j < 3; j++) + results.parametric_coords(p, j) = bary_out(j); }); return results; @@ -602,14 +611,14 @@ Kokkos::View GridPointSearch2D::operator()( GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny) : GridPointSearch2D(mesh, Nx, Ny, - PointSearchTolerances{"point search 2d tolerances"}) + PointSearchTolerances("point search 2d tolerances", mesh.dim())) { Kokkos::deep_copy(tolerances_, 1E-12); } GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, const PointSearchTolerances& tolerances) - : PointLocalizationSearch(tolerances), mesh_(mesh) + : PointSearch(tolerances), mesh_(mesh) { auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); @@ -633,42 +642,50 @@ GridPointSearch2D::GridPointSearch2D(Omega_h::Mesh& mesh, LO Nx, LO Ny, verts2faces_up_ = mesh.ask_up(Omega_h::VERT, Omega_h::FACE); } -LO GridPointSearch2D::GetOwningElementId(const Result& result) +LO GridPointSearch2D::GetOwningElementId(const Results& results, int i) { - const LO query_id = - (result.element_id < 0) ? -result.element_id : result.element_id; + const Kokkos::View query_id{""}; + const Kokkos::View dim{""}; + Kokkos::parallel_for(1, KOKKOS_LAMBDA(const int){ + query_id(0) = (results.element_ids(i) < 0) ? -results.element_ids(i) : results.element_ids(i); + dim(0) = results.dimensionalities(i); + }); + auto query_id_h = Kokkos::create_mirror_view(query_id); + auto dim_h = Kokkos::create_mirror_view(dim); + Kokkos::deep_copy(query_id_h, query_id); + Kokkos::deep_copy(dim_h, dim); return pcms::GetOwningElementId( - mesh_, 2, static_cast(result.dimensionality), query_id); + mesh_, 2, static_cast(dim_h(0)), query_id_h(0)); } Kokkos::View GridPointSearch2D::GetOwningElementIds( - Kokkos::View results) + Results const& results) { - Kokkos::View owners("point search owning face ids", results.extent(0)); + Kokkos::View owners("point search owning face ids", + results.dimensionalities.extent(0)); auto edges2faces_up = edges2faces_up_; auto verts2faces_up = verts2faces_up_; constexpr int mesh_dim = 2; Kokkos::parallel_for( - results.extent(0), KOKKOS_LAMBDA(const LO i) { - const auto result = results(i); - LO element_id = result.element_id; + results.dimensionalities.extent(0), KOKKOS_LAMBDA(const LO i) { + LO element_id = results.element_ids(i); if (element_id < 0) element_id = -element_id; LO owner = -1; if (element_id >= 0) { - if (result.dimensionality == Result::Dimensionality::FACE) { + if (results.dimensionalities(i) == Results::Dimensionality::FACE) { owner = GetOwningElementIdFromAdj( - edges2faces_up, Result::Dimensionality::FACE, - Result::Dimensionality::FACE, element_id); - } else if (result.dimensionality == Result::Dimensionality::EDGE) { + edges2faces_up, Results::Dimensionality::FACE, + Results::Dimensionality::FACE, element_id); + } else if (results.dimensionalities(i) == Results::Dimensionality::EDGE) { owner = GetOwningElementIdFromAdj( - edges2faces_up, Result::Dimensionality::EDGE, - Result::Dimensionality::FACE, element_id); - } else if (result.dimensionality == Result::Dimensionality::VERTEX) { + edges2faces_up, Results::Dimensionality::EDGE, + Results::Dimensionality::FACE, element_id); + } else if (results.dimensionalities(i) == Results::Dimensionality::VERTEX) { owner = GetOwningElementIdFromAdj( - verts2faces_up, Result::Dimensionality::VERTEX, - Result::Dimensionality::FACE, element_id); + verts2faces_up, Results::Dimensionality::VERTEX, + Results::Dimensionality::FACE, element_id); } } owners(i) = owner; @@ -676,25 +693,45 @@ Kokkos::View GridPointSearch2D::GetOwningElementIds( return owners; } -Kokkos::View GridPointSearch3D::operator()( - Kokkos::View points) const +GridPointSearch3D::Results GridPointSearch3D::apply( + const CoordinateView& points) const { - Kokkos::View results("point search result", - points.extent(0)); + Results results{ + Kokkos::View("dimensionalities", points.GetValues().extent(0)), + Kokkos::View("element IDs", points.GetValues().extent(0)), + Kokkos::View("parametric coordinates", points.GetValues().extent(0), + points.GetValues().extent(1) + 1) + }; + Kokkos::deep_copy(results.dimensionalities, Dimensionality::NO_INTERSECT); + Kokkos::deep_copy(results.element_ids, -1); + Kokkos::deep_copy(results.parametric_coords, 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 tets2verts = tets2verts_; + auto tets2verts_adj = tets2verts_adj_; + auto tets2edges_adj = tets2edges_adj_; + auto tets2faces_adj = tets2faces_adj_; auto edges2verts_adj = edges2verts_adj_; + auto tolerances = tolerances_; auto coords = coords_; + auto point_coords = points.GetValues(); + auto get_face = KOKKOS_LAMBDA(int i) + { + // inverse face offsets as determined by the cannonical ordering: + // https://user-images.githubusercontent.com/56453280/74203616-39d27a80-4c3e-11ea-885d-b0260490e184.png + static const char OFFSETS = 1 << 6 | 2 << 2 | 3; + return (OFFSETS & 3 << (i*2))>>(i*2); + }; Kokkos::parallel_for( - points.extent(0), KOKKOS_LAMBDA(int p) { + point_coords.extent(0), KOKKOS_LAMBDA(int p) { + const auto vtol = tolerances(0); + const auto etol = tolerances(1); + const auto ftol = tolerances(2); Omega_h::Vector point; for (int i = 0; i < DIM; ++i) { - point[i] = points(p, i); + point[i] = point_coords(p, i); } auto cell_id = grid(0).ClosestCellID(point); @@ -703,34 +740,131 @@ Kokkos::View GridPointSearch3D::operator()( auto candidates_end = candidate_map.row_map(cell_id + 1); bool found = false; - auto nearest_triangle = candidates_begin; - auto dimensionality = GridPointSearch3D::Result::Dimensionality::EDGE; + auto nearest_tetrahedron = candidates_begin; + auto dimensionality = Dimensionality::EDGE; 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 int tetrahedronID = candidate_map.entries(i); const auto elem_tri2verts = - Omega_h::gather_verts(tris2verts, triangleID); - auto vertex_coords = + Omega_h::gather_verts(tets2verts, tetrahedronID); + const auto vertex_coords = Omega_h::gather_vectors(coords, elem_tri2verts); - auto parametric_coords = + const Omega_h::Vector parametric_coords = Omega_h::barycentric_from_global(point, vertex_coords); + LO vertexID = -1; + Real best_vdist_sq = INFINITY; + for (int k = 0; k < 4; k++) { + Real vdist_sq = Omega_h::norm_squared(vertex_coords[k] - point); + if (vdist_sq < vtol*vtol && vdist_sq < best_vdist_sq) + { + vertexID = elem_tri2verts[k]; + best_vdist_sq = vdist_sq; + } + } + if (vertexID > -1) + { + results.dimensionalities(p) = Dimensionality::VERTEX; + results.element_ids(p) = vertexID; + for (int j = 0; j < 4; j++) + results.parametric_coords(p, j) = parametric_coords[j]; + found = true; + } + + if (results.dimensionalities(p) == Dimensionality::VERTEX) continue; + + LO edgeID = -1; + Real best_edist = INFINITY; + for (int j = 0; j < 6; ++j) { + const int curr_edgeID = tets2edges_adj.ab2b[tetrahedronID * 6 + j]; + const int va_id = edges2verts_adj.ab2b[curr_edgeID * 2 + 0]; + const int vb_id = edges2verts_adj.ab2b[curr_edgeID * 2 + 1]; + const auto va = Omega_h::get_vector<3>(coords, va_id); + const auto vb = Omega_h::get_vector<3>(coords, vb_id); + if (!normal_intersects_segment(va, vb, point)) + continue; + const auto de = distance_from_line(va, vb, point); + if (de < etol && de < best_edist) { + best_edist = de; + edgeID = curr_edgeID; + } + } + if (edgeID > -1) + { + results.dimensionalities(p) = Dimensionality::EDGE; + results.element_ids(p) = edgeID; + for (int j = 0; j < 4; j++) + results.parametric_coords(p, j) = parametric_coords[j]; + found = true; + } + + if (results.dimensionalities(p) == Dimensionality::EDGE) continue; + + LO faceID = -1; + Real best_fdist_sq = INFINITY; + Omega_h::Matrix tet; + tet[0] = vertex_coords[0] - vertex_coords[3]; + tet[1] = vertex_coords[1] - vertex_coords[3]; + tet[2] = vertex_coords[2] - vertex_coords[3]; + Real tetrahedron_volume = fabs(Omega_h::determinant(tet))/6.; + + for (int k = 0; k < 4; k++) + { + LO bary_offset = get_face(k); + Real face_area_sq = Omega_h::norm_squared(Omega_h::cross( + vertex_coords[(bary_offset+1)%4] - vertex_coords[(bary_offset+3)%4], + vertex_coords[(bary_offset+2)%4] - vertex_coords[(bary_offset+3)%4]))/4.; + + Real fdist_sq = parametric_coords[bary_offset]*3*tetrahedron_volume; + fdist_sq *= fdist_sq; + fdist_sq /= face_area_sq; + + bool inside = true; + for (int j = 0; j < 4; j++) + { + if (j == bary_offset) continue; + if (parametric_coords[j] < 0) + { + inside = false; + break; + } + } + + if (fdist_sq < ftol*ftol && fdist_sq < best_fdist_sq && inside) + { + faceID = tets2faces_adj.ab2b[tetrahedronID * 4 + k]; + best_fdist_sq = fdist_sq; + } + } + if (faceID > -1) + { + results.dimensionalities(p) = Dimensionality::FACE; + results.element_ids(p) = faceID; + for (int j = 0; j < 4; j++) + results.parametric_coords(p, j) = parametric_coords[j]; + found = true; + } + + if (results.dimensionalities(p) == Dimensionality::FACE) continue; + if (Omega_h::is_barycentric_inside(parametric_coords)) { - results(p) = GridPointSearch3D::Result{ - GridPointSearch3D::Result::Dimensionality::REGION, triangleID, - parametric_coords}; + results.dimensionalities(p) = Dimensionality::REGION; + results.element_ids(p) = tetrahedronID; + for (int j = 0; j < 4; j++) + results.parametric_coords(p, j) = parametric_coords[j]; found = true; - break; } // TODO: Get nearest element if no tetrahedron found } if (!found) { - LO nearest_elem = candidate_map.entries(nearest_triangle); - results(p) = GridPointSearch3D::Result{dimensionality, -nearest_elem, - parametric_coords_to_nearest}; + LO nearest_elem = candidate_map.entries(nearest_tetrahedron); + results.dimensionalities(p) = dimensionality; + results.element_ids(p) = -nearest_elem; + for (int j = 0; j < 4; j++) + results.parametric_coords(p,j) = parametric_coords_to_nearest[j]; } }); @@ -739,14 +873,14 @@ Kokkos::View GridPointSearch3D::operator()( GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz) : GridPointSearch3D(mesh, Nx, Ny, Nz, - PointSearchTolerances{"point search 3d tolerances"}) + PointSearchTolerances("point search 3d tolerances", 3)) { - Kokkos::deep_copy(tolerances_, 0); + Kokkos::deep_copy(tolerances_, 1e-12); } GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, const PointSearchTolerances& tolerances) - : PointLocalizationSearch(tolerances), mesh_(mesh) + : PointSearch(tolerances), mesh_(mesh) { auto mesh_bbox = Omega_h::get_bounding_box<3>(&mesh); auto grid_h = Kokkos::create_mirror_view(grid_); @@ -767,56 +901,57 @@ GridPointSearch3D::GridPointSearch3D(Omega_h::Mesh& mesh, LO Nx, LO Ny, LO Nz, 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); - tris2verts_adj_ = mesh.ask_down(Omega_h::FACE, Omega_h::VERT); + tets2verts_ = mesh.ask_elem_verts(); + tets2faces_adj_ = mesh.ask_down(Omega_h::REGION, Omega_h::FACE); + tets2edges_adj_ = mesh.ask_down(Omega_h::REGION, Omega_h::EDGE); + tets2verts_adj_ = mesh.ask_down(Omega_h::REGION, Omega_h::VERT); edges2verts_adj_ = mesh.ask_down(Omega_h::EDGE, Omega_h::VERT); verts2regions_up_ = mesh.ask_up(Omega_h::VERT, Omega_h::REGION); edges2regions_up_ = mesh.ask_up(Omega_h::EDGE, Omega_h::REGION); faces2regions_up_ = mesh.ask_up(Omega_h::FACE, Omega_h::REGION); } -LO GridPointSearch3D::GetOwningElementId(const Result& result) +LO GridPointSearch3D::GetOwningElementId(const Results& results, int i) { const LO query_id = - (result.element_id < 0) ? -result.element_id : result.element_id; + (results.element_ids(i) < 0) ? -results.element_ids(i) : results.element_ids(i); return pcms::GetOwningElementId( - mesh_, 3, static_cast(result.dimensionality), query_id); + mesh_, 3, static_cast(results.dimensionalities(i)), query_id); } Kokkos::View GridPointSearch3D::GetOwningElementIds( - Kokkos::View results) + const Results& results) { - Kokkos::View owners("point search owning region ids", results.extent(0)); + Kokkos::View owners("point search owning region ids", + results.dimensionalities.extent(0)); auto verts2regions_up = verts2regions_up_; auto edges2regions_up = edges2regions_up_; auto faces2regions_up = faces2regions_up_; constexpr int mesh_dim = 3; Kokkos::parallel_for( - results.extent(0), KOKKOS_LAMBDA(const LO i) { - const auto result = results(i); - LO element_id = result.element_id; + owners.extent(0), KOKKOS_LAMBDA(const LO i) { + LO element_id = results.element_ids(i); if (element_id < 0) element_id = -element_id; LO owner = -1; if (element_id >= 0) { - if (result.dimensionality == Result::Dimensionality::REGION) { + if (results.dimensionalities(i) == Results::Dimensionality::REGION) { owner = GetOwningElementIdFromAdj( - faces2regions_up, Result::Dimensionality::REGION, - Result::Dimensionality::REGION, element_id); - } else if (result.dimensionality == Result::Dimensionality::FACE) { + faces2regions_up, Results::Dimensionality::REGION, + Results::Dimensionality::REGION, element_id); + } else if (results.dimensionalities(i) == Results::Dimensionality::FACE) { owner = GetOwningElementIdFromAdj( - faces2regions_up, Result::Dimensionality::FACE, - Result::Dimensionality::REGION, element_id); - } else if (result.dimensionality == Result::Dimensionality::EDGE) { + faces2regions_up, Results::Dimensionality::FACE, + Results::Dimensionality::REGION, element_id); + } else if (results.dimensionalities(i) == Results::Dimensionality::EDGE) { owner = GetOwningElementIdFromAdj( - edges2regions_up, Result::Dimensionality::EDGE, - Result::Dimensionality::REGION, element_id); - } else if (result.dimensionality == Result::Dimensionality::VERTEX) { + edges2regions_up, Results::Dimensionality::EDGE, + Results::Dimensionality::REGION, element_id); + } else if (results.dimensionalities(i) == Results::Dimensionality::VERTEX) { owner = GetOwningElementIdFromAdj( - verts2regions_up, Result::Dimensionality::VERTEX, - Result::Dimensionality::REGION, element_id); + verts2regions_up, Results::Dimensionality::VERTEX, + Results::Dimensionality::REGION, element_id); } } owners(i) = owner; diff --git a/src/pcms/localization/point_search.h b/src/pcms/localization/point_search.h index b1e0a9ef..e5054e8c 100644 --- a/src/pcms/localization/point_search.h +++ b/src/pcms/localization/point_search.h @@ -10,6 +10,7 @@ #include "pcms/utility/types.h" #include "pcms/utility/uniform_grid.h" #include "pcms/utility/bounding_box.h" +#include "pcms/field/coordinate_system.h" namespace pcms { @@ -22,33 +23,33 @@ namespace pcms * top-level dimension). */ LO GetOwningElementId(Omega_h::Mesh& mesh, int mesh_dim, int entity_dim, - LO element_id); + LO element_id); template KOKKOS_INLINE_FUNCTION LO GetOwningElementIdFromAdj(const Omega_h::Adj& upward_adj, DimEnum entity_dim, - DimEnum target_dim, LO element_id) + DimEnum target_dim, LO element_id) { - if (element_id < 0) - return -1; - - // If entity is already at target dimension, return it directly - if (static_cast(entity_dim) == static_cast(target_dim)) - return element_id; - - const auto begin = upward_adj.a2ab[element_id]; - const auto end = upward_adj.a2ab[element_id + 1]; - if (begin >= end) - return -1; - - // Find the smallest owning element ID - LO owner = upward_adj.ab2b[begin]; - for (auto i = begin + 1; i < end; ++i) { - const LO candidate = upward_adj.ab2b[i]; - if (candidate < owner) - owner = candidate; - } - return owner; + if (element_id < 0) + return -1; + + // If entity is already at target dimension, return it directly + if (static_cast(entity_dim) == static_cast(target_dim)) + return element_id; + + const auto begin = upward_adj.a2ab[element_id]; + const auto end = upward_adj.a2ab[element_id + 1]; + if (begin >= end) + return -1; + + // Find the smallest owning element ID + LO owner = upward_adj.ab2b[begin]; + for (auto i = begin + 1; i < end; ++i) { + const LO candidate = upward_adj.ab2b[i]; + if (candidate < owner) + owner = candidate; + } + return owner; } // @@ -59,142 +60,153 @@ namespace detail { Kokkos::Crs construct_intersection_map_2d(Omega_h::Mesh& mesh, - Kokkos::View grid, - int num_grid_cells); + Kokkos::View grid, + int num_grid_cells); } [[nodiscard]] KOKKOS_FUNCTION bool triangle_intersects_bbox( - const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); + const Omega_h::Matrix<2, 3>& coords, const AABBox<2>& bbox); -template -class PointLocalizationSearch +/** + * Point search base class + */ +class PointSearch { public: - struct Result - { - enum class Dimensionality - { - VERTEX = 0, - EDGE = 1, - FACE = 2, - REGION = 3 - }; - - Dimensionality dimensionality; - LO element_id; // ID of the actual entity (vertex/edge/face/region) - Omega_h::Vector parametric_coords; - }; - - 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; - - /** - * This function provides a temporary solution to retrieve the original - * behavior of the point search, which previously returned the face ID - * regardless of which entity the search result belonged to. Many parts of - * the codebase still rely on this legacy behavior. After updating the point - * search to return the exact entity, this function can be called to retrieve - * the old behavior and ensure correctness. Long term, the implementation - * should be updated to properly handle the new result. - */ - [[nodiscard]] virtual LO GetOwningElementId(const Result& result) = 0; - [[nodiscard]] virtual Kokkos::View GetOwningElementIds( - Kokkos::View results) = 0; - virtual ~PointLocalizationSearch() = default; - + using ExecSpace = Omega_h::ExecSpace; + using MemorySpace = ExecSpace::memory_space; + /** + * Results type gives dimensionalities of point intersection, the intersected + * element IDs, and the barycentric coordinate mappings of each point as follows: + * the `i`th input point has intersection dimensionality of + * `dimensionalities(i)`, intersected element ID of `element_ids(i)`, and + * likewise parametric coordinate mapping of `parametric_coords(i,0..d)` + * where `d` is the dimension of the space + */ + struct Results + { + enum class Dimensionality + { + VERTEX = 0, + EDGE = 1, + FACE = 2, + REGION = 3, + NO_INTERSECT = 4 + }; + + Kokkos::View dimensionalities; + Kokkos::View element_ids; + Kokkos::View parametric_coords; + }; + + using PointSearchTolerances = Kokkos::View; + + explicit PointSearch(const PointSearchTolerances& tolerances) + : tolerances_(tolerances) + { + assert(tolerances_.is_allocated()); + } + + ~PointSearch() = default; + + virtual Results apply(const CoordinateView& coords) const = 0; + /** + * This function provides a temporary solution to retrieve the original + * behavior of the point search, which previously returned the face ID + * regardless of which entity the search result belonged to. Many parts of + * the codebase still rely on this legacy behavior. After updating the point + * search to return the exact entity, this function can be called to retrieve + * the old behavior and ensure correctness. Long term, the implementation + * should be updated to properly handle the new result. + */ + [[nodiscard]] virtual LO GetOwningElementId(const Results& results, int i) = 0; + [[nodiscard]] virtual Kokkos::View GetOwningElementIds( + const Results& results) = 0; protected: - PointSearchTolerances tolerances_; + PointSearchTolerances tolerances_; }; -using PointLocalizationSearch2D = PointLocalizationSearch<2>; -using PointLocalizationSearch3D = PointLocalizationSearch<3>; - -class GridPointSearch2D : public PointLocalizationSearch2D +class GridPointSearch2D : public PointSearch { - 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); - - /** - * 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; - [[nodiscard]] LO GetOwningElementId(const Result& result) override; - [[nodiscard]] Kokkos::View GetOwningElementIds( - Kokkos::View results) override; + static constexpr int DIM = 2; + + using Results = PointSearch::Results; + using Dimensionality = Results::Dimensionality; + + 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 + * 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 + */ + Results apply( + const CoordinateView& coords) const override; + [[nodiscard]] LO GetOwningElementId(const Results& result, int i) override; + [[nodiscard]] Kokkos::View GetOwningElementIds(const Results& results) override; private: - Omega_h::Mesh mesh_; - Omega_h::Adj tris2edges_adj_; - Omega_h::Adj tris2verts_adj_; - Omega_h::Adj edges2verts_adj_; - Omega_h::Adj edges2faces_up_; - Omega_h::Adj verts2faces_up_; - Kokkos::View grid_{"uniform grid"}; - CandidateMapT candidate_map_; - Omega_h::LOs tris2verts_; - Omega_h::Reals coords_; + Omega_h::Mesh mesh_; + Omega_h::Adj tris2edges_adj_; + Omega_h::Adj tris2verts_adj_; + Omega_h::Adj edges2verts_adj_; + Omega_h::Adj edges2faces_up_; + Omega_h::Adj verts2faces_up_; + Kokkos::View grid_{"uniform grid"}; + CandidateMapT candidate_map_; + Omega_h::LOs tris2verts_; + Omega_h::Reals coords_; }; -class GridPointSearch3D : public PointLocalizationSearch3D +class GridPointSearch3D : public PointSearch { - 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); - - /** - * 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; - [[nodiscard]] LO GetOwningElementId(const Result& result) override; - [[nodiscard]] Kokkos::View GetOwningElementIds( - Kokkos::View results) override; + static constexpr int DIM = 3; + + using Results = PointSearch::Results; + using Dimensionality = Results::Dimensionality; + + 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. + */ + Results apply( + const CoordinateView& coords) const override; + [[nodiscard]] LO GetOwningElementId(const Results& result, int i) override; + [[nodiscard]] Kokkos::View GetOwningElementIds(const Results& results) override; private: - Omega_h::Mesh mesh_; - Omega_h::Adj tris2edges_adj_; - Omega_h::Adj tris2verts_adj_; - Omega_h::Adj edges2verts_adj_; - Omega_h::Adj verts2regions_up_; - Omega_h::Adj edges2regions_up_; - Omega_h::Adj faces2regions_up_; - Kokkos::View[1]> grid_{"uniform grid"}; - CandidateMapT candidate_map_; - Omega_h::LOs tris2verts_; - Omega_h::Reals coords_; - Real fuzz_; + Omega_h::Mesh mesh_; + Omega_h::Adj tets2faces_adj_; + Omega_h::Adj tets2edges_adj_; + Omega_h::Adj tets2verts_adj_; + Omega_h::Adj edges2verts_adj_; + Omega_h::Adj verts2regions_up_; + Omega_h::Adj edges2regions_up_; + Omega_h::Adj faces2regions_up_; + Kokkos::View[1]> grid_{"uniform grid"}; + CandidateMapT candidate_map_; + Omega_h::LOs tets2verts_; + Omega_h::Reals coords_; + Real fuzz_; }; } // namespace pcms diff --git a/src/pcms/transfer/mesh_intersection.cpp b/src/pcms/transfer/mesh_intersection.cpp index 564ae350..e96bfb3a 100644 --- a/src/pcms/transfer/mesh_intersection.cpp +++ b/src/pcms/transfer/mesh_intersection.cpp @@ -29,7 +29,9 @@ void FindIntersections::adjBasedIntersectSearch( auto centroids = ConvertCoordsTo2D(flat_centroids, target_mesh_.nfaces(), 2); pcms::GridPointSearch2D search_cell(source_mesh_, 20, 20); - auto results = search_cell(centroids); + auto results = search_cell.apply(CoordinateView( + CoordinateSystem::Cartesian, + MakeConstRank2View(centroids))); auto owning_cell_ids = search_cell.GetOwningElementIds(results); auto nfaces_target = target_mesh_.nfaces(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1c603c45..7a5293f1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -399,7 +399,9 @@ if(Catch2_FOUND) test_localization_factory.cpp test_omega_h_field2_outofbounds.cpp test_omega_h_lagrange_field.cpp - test_point_evaluator.cpp) + test_point_evaluator.cpp + test_tree_localization.cpp + ) endif() if(PCMS_ENABLE_MESHFIELDS) diff --git a/test/test_point_search.cpp b/test/test_point_search.cpp index 6bb9a1cc..2348deb8 100644 --- a/test/test_point_search.cpp +++ b/test/test_point_search.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -162,23 +163,23 @@ TEST_CASE("construct intersection map") REQUIRE(num_candidates_within_range(intersection_map, 1, 6)); } } -TEST_CASE("uniform grid search") +TEST_CASE("uniform grid search 2D") { - using pcms::GridPointSearch2D; + using pcms::PointSearch; 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"}; + pcms::PointSearch::PointSearchTolerances("point search 2d tolerances", 2); auto tolerances_h = Kokkos::create_mirror_view(tolerances); tolerances_h(0) = 0.01; tolerances_h(1) = 0.01; Kokkos::deep_copy(tolerances, tolerances_h); - GridPointSearch2D search{mesh, 10, 10, tolerances}; + std::unique_ptr search = std::make_unique(mesh, 10, 10, tolerances); - Kokkos::View points("test_points", 8); + Kokkos::View points("test_points", 8, 2); // Kokkos::View points("test_points", 1); auto points_h = Kokkos::create_mirror_view(points); points_h(0, 0) = 0; @@ -198,38 +199,49 @@ TEST_CASE("uniform grid search") 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); - Kokkos::deep_copy(results_h, results); + auto results = search->apply( + pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(points))); + auto result_dims_h = Kokkos::create_mirror_view(results.dimensionalities); + Kokkos::deep_copy(result_dims_h, results.dimensionalities); + auto result_ids_h = Kokkos::create_mirror_view(results.element_ids); + Kokkos::deep_copy(result_ids_h, results.element_ids); + auto result_coords_h = Kokkos::create_mirror_view(results.parametric_coords); + Kokkos::deep_copy(result_coords_h, results.parametric_coords); SECTION("global coordinate within mesh") { { - auto [dim, idx, coords] = results_h(0); - const auto face_idx = search.GetOwningElementId(results_h(0)); + auto dim = result_dims_h(0); + auto idx = result_ids_h(0); + auto coords = Omega_h::Vector<3>{result_coords_h(0,0), result_coords_h(0,1), result_coords_h(0,2)}; + const auto face_idx = search->GetOwningElementId(results, 0); CAPTURE(idx); - REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(dim == PointSearch::Results::Dimensionality::VERTEX); REQUIRE(idx == 0); REQUIRE(face_idx >= 0); - REQUIRE(coords[0] == Catch::Approx(1)); - REQUIRE(coords[1] == Catch::Approx(0)); - REQUIRE(coords[2] == Catch::Approx(0)); + REQUIRE(coords[0] == Catch::Approx(1).margin(0.00001)); + REQUIRE(coords[1] == Catch::Approx(0).margin(0.00001)); + REQUIRE(coords[2] == Catch::Approx(0).margin(0.00001)); } { - auto [dim, idx, coords] = results_h(1); - const auto face_idx = search.GetOwningElementId(results_h(1)); - REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::EDGE); + auto dim = result_dims_h(1); + auto idx = result_ids_h(1); + auto coords = Omega_h::Vector<3>{result_coords_h(1,0), result_coords_h(1,1), result_coords_h(1,2)}; + const auto face_idx = search->GetOwningElementId(results, 1); + REQUIRE(dim == PointSearch::Results::Dimensionality::EDGE); REQUIRE(idx == 156); REQUIRE(face_idx >= 0); - REQUIRE(coords[0] == Catch::Approx(0.5)); - REQUIRE(coords[1] == Catch::Approx(0.1)); - REQUIRE(coords[2] == Catch::Approx(0.4)); + REQUIRE(coords[0] == Catch::Approx(0.5).margin(0.00001)); + REQUIRE(coords[1] == Catch::Approx(0.1).margin(0.00001)); + REQUIRE(coords[2] == Catch::Approx(0.4).margin(0.00001)); } { - auto [dim, idx, coords] = results_h(7); - const auto face_idx = search.GetOwningElementId(results_h(7)); - REQUIRE(dim == GridPointSearch2D::Result::Dimensionality::FACE); + auto dim = result_dims_h(7); + auto idx = result_ids_h(7); + const auto face_idx = search->GetOwningElementId(results, 7); + REQUIRE(dim == PointSearch::Results::Dimensionality::FACE); REQUIRE(idx == 0); REQUIRE(face_idx >= 0); } @@ -237,51 +249,196 @@ TEST_CASE("uniform grid search") // feature needs to be added 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(search.GetOwningElementId(out_of_bounds) >= 0); + auto out_of_bounds_dim = result_dims_h(2); + auto out_of_bounds_id = result_ids_h(2); + auto top_right_id = result_ids_h(3); + REQUIRE(out_of_bounds_dim == + PointSearch::Results::Dimensionality::VERTEX); + REQUIRE(-1 * out_of_bounds_id == top_right_id); + REQUIRE(search->GetOwningElementId(results, 2) >= 0); - out_of_bounds = results_h(4); - auto bot_left = results_h(0); - REQUIRE(out_of_bounds.dimensionality == - GridPointSearch2D::Result::Dimensionality::VERTEX); - REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); - REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); + out_of_bounds_dim = result_dims_h(4); + out_of_bounds_id = result_ids_h(4); + auto bot_left_id = result_ids_h(0); + REQUIRE(out_of_bounds_dim == + PointSearch::Results::Dimensionality::VERTEX); + REQUIRE(-1 * out_of_bounds_id == bot_left_id); + REQUIRE(search->GetOwningElementId(results, 4) >= 0); - out_of_bounds = results_h(5); - REQUIRE(out_of_bounds.dimensionality == - GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(out_of_bounds.element_id == -219); - REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); + out_of_bounds_dim = result_dims_h(5); + out_of_bounds_id = result_ids_h(5); + REQUIRE(out_of_bounds_dim == + PointSearch::Results::Dimensionality::EDGE); + REQUIRE(out_of_bounds_id == -219); + REQUIRE(search->GetOwningElementId(results, 5) >= 0); - out_of_bounds = results_h(6); - REQUIRE(out_of_bounds.dimensionality == - GridPointSearch2D::Result::Dimensionality::EDGE); - REQUIRE(-1 * out_of_bounds.element_id == bot_left.element_id); - REQUIRE(search.GetOwningElementId(out_of_bounds) >= 0); + out_of_bounds_dim = result_dims_h(6); + out_of_bounds_id = result_ids_h(6); + REQUIRE(out_of_bounds_dim == + PointSearch::Results::Dimensionality::EDGE); + REQUIRE(-1 * out_of_bounds_id == bot_left_id); + REQUIRE(search->GetOwningElementId(results, 6) >= 0); } SECTION("point on extension of an edge") { - Kokkos::View ext_points("ext_test_points", 1); + Kokkos::View ext_points("ext_test_points", 1, 2); auto ext_h = Kokkos::create_mirror_view(ext_points); ext_h(0, 0) = 1.5; ext_h(0, 1) = 0.0; Kokkos::deep_copy(ext_points, ext_h); - auto ext_results = search(ext_points); - auto ext_results_h = Kokkos::create_mirror_view(ext_results); - Kokkos::deep_copy(ext_results_h, ext_results); + auto ext_results = search->apply(pcms::CoordinateView(pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(ext_points))); - auto res = ext_results_h(0); + auto result_dims_h = Kokkos::create_mirror_view(ext_results.dimensionalities); + Kokkos::deep_copy(result_dims_h, ext_results.dimensionalities); + auto result_ids_h = Kokkos::create_mirror_view(ext_results.element_ids); + Kokkos::deep_copy(result_ids_h, ext_results.element_ids); // Must be out-of-bounds (negative element id) - REQUIRE(res.element_id < 0); + REQUIRE(result_ids_h(0) < 0); // Dimensionality must be VERTEX — the nearest entity is the // rightmost bottom vertex (1.0, 0). - REQUIRE(res.dimensionality == - GridPointSearch2D::Result::Dimensionality::VERTEX); + REQUIRE(result_dims_h(0) == + PointSearch::Results::Dimensionality::VERTEX); // Owning element should still resolve to a valid face - REQUIRE(search.GetOwningElementId(res) >= 0); + REQUIRE(search->GetOwningElementId(results, 0) >= 0); + } +} +TEST_CASE("uniform grid search 3D") +{ + using pcms::GridPointSearch3D; + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 1, 1, 1, false); + auto tolerances = + GridPointSearch3D::PointSearchTolerances("point search 3d tolerances", 3); + auto tolerances_h = Kokkos::create_mirror_view(tolerances); + tolerances_h(0) = 0.01; + tolerances_h(1) = 0.01; + tolerances_h(2) = 0.01; + Kokkos::deep_copy(tolerances, tolerances_h); + + GridPointSearch3D search{mesh, 5, 5, 5, tolerances}; + + auto check_res = [] (auto const& results, int dim) + { + auto dims = Kokkos::create_mirror_view(results.dimensionalities); + auto ids = Kokkos::create_mirror_view(results.element_ids); + Kokkos::deep_copy(dims, results.dimensionalities); + Kokkos::deep_copy(ids, results.element_ids); + for (int i = 0; i < dims.size(); i++) + { + CHECK(dims(i) == (pcms::PointSearch::Results::Dimensionality)dim); + CHECK(ids(i) == i); + } + }; + SECTION ("Vertex intersection") { + Kokkos::View coords("vertex coordinates", mesh.nverts(), 3); + auto coords_h = Kokkos::create_mirror_view(coords); + auto meshcoords = Omega_h::HostRead(mesh.coords()); + for (int i = 0; i < mesh.nverts(); i++) + { + coords_h(i, 0) = meshcoords[i*3]; + coords_h(i, 1) = meshcoords[i*3 + 1]; + coords_h(i, 2) = meshcoords[i*3 + 2]; + } + Kokkos::deep_copy(coords, coords_h); + pcms::CoordinateView> vert_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::PointSearch::Results results = search.apply(vert_cv); + + check_res(results, 0); + } + SECTION ("Edge intersection") { + Kokkos::View coords("intersection coords", mesh.nedges(), 3); + auto edge2vert = Omega_h::HostRead(mesh.get_adj(Omega_h::EDGE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (int i = 0; i < mesh.nedges(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[edge2vert[2*i]*3], meshcoords[edge2vert[2*i]*3 + 1], meshcoords[edge2vert[2*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[edge2vert[2*i + 1]*3], meshcoords[edge2vert[2*i + 1]*3 + 1], meshcoords[edge2vert[2*i + 1]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1)/2.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + + Kokkos::deep_copy(coords, coords_h); + + pcms::CoordinateView> edge_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::PointSearch::Results results = search.apply(edge_cv); + check_res(results, 1); + } + SECTION ("Face intersection") { + Kokkos::View coords("intersection coords", mesh.nfaces(), 3); + auto face2vert = Omega_h::HostRead(mesh.get_adj(Omega_h::FACE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (int i = 0; i < mesh.nfaces(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[face2vert[3*i]*3], meshcoords[face2vert[3*i]*3 + 1], meshcoords[face2vert[3*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[face2vert[3*i + 1]*3], meshcoords[face2vert[3*i + 1]*3 + 1], meshcoords[face2vert[3*i + 1]*3 + 2] + }; + Omega_h::Vector<3> v2{ + meshcoords[face2vert[3*i + 2]*3], meshcoords[face2vert[3*i + 2]*3 + 1], meshcoords[face2vert[3*i + 2]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1 + v2)/3.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + Kokkos::deep_copy(coords, coords_h); + + pcms::CoordinateView> face_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::GridPointSearch3D::Results results = search.apply(face_cv); + check_res(results, 2); + } + SECTION ("Region intersection") { + Kokkos::View coords("intersection coords", mesh.nelems(), 3); + auto region2vert = Omega_h::HostRead(mesh.get_adj(Omega_h::REGION, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (unsigned i = 0; i < mesh.nelems(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[region2vert[4*i]*3], meshcoords[region2vert[4*i]*3 + 1], meshcoords[region2vert[4*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[region2vert[4*i + 1]*3], meshcoords[region2vert[4*i + 1]*3 + 1], meshcoords[region2vert[4*i + 1]*3 + 2] + }; + Omega_h::Vector<3> v2{ + meshcoords[region2vert[4*i + 2]*3], meshcoords[region2vert[4*i + 2]*3 + 1], meshcoords[region2vert[4*i + 2]*3 + 2] + }; + Omega_h::Vector<3> v3{ + meshcoords[region2vert[4*i + 3]*3], meshcoords[region2vert[4*i + 3]*3 + 1], meshcoords[region2vert[4*i + 3]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1 + v2 + v3)/4.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + Kokkos::deep_copy(coords, coords_h); + + pcms::CoordinateView> region_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::PointSearch::Results results = search.apply(region_cv); + check_res(results, 3); } } diff --git a/test/test_tree_localization.cpp b/test/test_tree_localization.cpp new file mode 100644 index 00000000..dd854dfb --- /dev/null +++ b/test/test_tree_localization.cpp @@ -0,0 +1,291 @@ +#include +#include +#include +#include + +#include +#include + +#include + +TEST_CASE ("Test 10x10 2D tree point classification") { + // Setup for 10x10 tree test cases + auto lib = Omega_h::Library{}; + auto world = lib.world(); + Omega_h::Mesh mesh_1x1 = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false); + REQUIRE(mesh_1x1.dim() == 2); + Omega_h::ExecSpace execs; + + pcms::TreePointSearch tree_search(mesh_1x1); + // End setup for 10x10 tree test cases + + auto check_res = [] (auto const& results, int dim) + { + auto dims = Kokkos::create_mirror_view(results.dimensionalities); + auto ids = Kokkos::create_mirror_view(results.element_ids); + Kokkos::deep_copy(dims, results.dimensionalities); + Kokkos::deep_copy(ids, results.element_ids); + for (int i = 0; i < dims.size(); i++) + { + CHECK(dims(i) == (pcms::TreePointSearch::Dimensionality)dim); + CHECK(ids(i) == i); + // CHECK_THAT(results(i).parametric_coords, + // Catch::Matchers::Contains(Catch::Matchers::WithinAbs(1.0, 10e-7))); + } + }; + + SECTION ("Vertex intersection") { + Kokkos::View coords("vertex coordinates", mesh_1x1.nverts(), 2); + auto coords_h = Kokkos::create_mirror_view(coords); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + for (int i = 0; i < mesh_1x1.nverts(); i++) + { + coords_h(i,0) = meshcoords[i*2]; + coords_h(i,1) = meshcoords[i*2 + 1]; + } + Kokkos::deep_copy(execs, coords, coords_h); + pcms::CoordinateView> vert_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(vert_cv); + check_res(results, 0); + } + SECTION ("Edge intersection") { + Kokkos::View coords("edge coordinates", mesh_1x1.nedges(), 2); + auto coords_h = Kokkos::create_mirror_view(coords); + auto edge2vert = Omega_h::HostRead(mesh_1x1.get_adj(Omega_h::EDGE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + for (int i = 0; i < mesh_1x1.nedges(); i++) + { + Omega_h::Vector<2> v0{ + meshcoords[edge2vert[2*i]*2], meshcoords[edge2vert[2*i]*2 + 1] + }; + Omega_h::Vector<2> v1{ + meshcoords[edge2vert[2*i + 1]*2], meshcoords[edge2vert[2*i + 1]*2 + 1] + }; + Omega_h::Vector<2> midpoint = (v0 + v1)/2.; + coords_h(i,0) = midpoint[0]; + coords_h(i,1) = midpoint[1]; + } + Kokkos::deep_copy(execs, coords, coords_h); + + pcms::CoordinateView> edge_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(edge_cv); + check_res(results, 1); + } + SECTION ("Face intersection") { + Kokkos::View coords("face coordinates", mesh_1x1.nfaces(), 2); + auto face2vert = Omega_h::HostRead(mesh_1x1.get_adj(Omega_h::FACE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + for (int i = 0; i < mesh_1x1.nfaces(); i++) + { + Omega_h::Vector<2> v0{ + meshcoords[face2vert[3*i]*2], meshcoords[face2vert[3*i]*2 + 1] + }; + Omega_h::Vector<2> v1{ + meshcoords[face2vert[3*i + 1]*2], meshcoords[face2vert[3*i + 1]*2 + 1] + }; + Omega_h::Vector<2> v2{ + meshcoords[face2vert[3*i + 2]*2], meshcoords[face2vert[3*i + 2]*2 + 1] + }; + Omega_h::Vector<2> midpoint = (v0 + v1 + v2)/3.; + coords_h(i,0) = midpoint[0]; + coords_h(i,1) = midpoint[1]; + } + Kokkos::deep_copy(execs, coords, coords_h); + + pcms::CoordinateView> face_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(face_cv); + check_res(results, 2); + } + // SECTION("Outside Mesh") { + // Omega_h::Write coords + // { + // -0.5, -0.5, + // -0.5, 0.0, + // -0.5, 0.5, + // -0.5, 1.0, + // -0.5, 1.5, + // 0.0, -0.5, + // 0.0, 1.5, + // 0.5, -0.5, + // 0.5, 1.5, + // 1.0, -0.5, + // 1.0, 1.5, + // 1.5, -0.5, + // 1.5, 0.0, + // 1.5, 0.5, + // 1.5, 1.0, + // 1.5, 1.5 + // }; + // pcms::CoordinateView outside_cv( + // pcms::CoordinateSystem::Cartesian, + // pcms::MakeConstRank2View(Omega_h::Read(coords), 2)); + // pcms::TreePointSearch::Results results = tree_search.apply(outside_cv); + // check_res(results, 4); + // } +} + +TEST_CASE ("Test 5x5x5 3D tree point classification") { + // Setup for 5x5x5 tree test cases + auto lib = Omega_h::Library{}; + auto world = lib.world(); + Omega_h::Mesh mesh_1x1 = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 5, 5, 5, false); + REQUIRE(mesh_1x1.dim() == 3); + Omega_h::ExecSpace execs; + + pcms::TreePointSearch tree_search(mesh_1x1); + + auto check_res = [] (auto const& results, int dim) + { + auto dims = Kokkos::create_mirror_view(results.dimensionalities); + auto ids = Kokkos::create_mirror_view(results.element_ids); + Kokkos::deep_copy(dims, results.dimensionalities); + Kokkos::deep_copy(ids, results.element_ids); + for (int i = 0; i < dims.size(); i++) + { + CHECK(dims(i) == (pcms::TreePointSearch::Dimensionality)dim); + CHECK(ids(i) == i); + // CHECK_THAT(results(i).parametric_coords, + // Catch::Matchers::Contains(Catch::Matchers::WithinAbs(1.0, 10e-7))); + } + }; + + SECTION ("Vertex intersection") { + Kokkos::View coords("vertex coordinates", mesh_1x1.nverts(), 3); + auto coords_h = Kokkos::create_mirror_view(coords); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + for (int i = 0; i < mesh_1x1.nverts(); i++) + { + coords_h(i, 0) = meshcoords[i*3]; + coords_h(i, 1) = meshcoords[i*3 + 1]; + coords_h(i, 2) = meshcoords[i*3 + 2]; + } + Kokkos::deep_copy(execs, coords, coords_h); + pcms::CoordinateView> vert_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(vert_cv); + + check_res(results, 0); + } + SECTION ("Edge intersection") { + Kokkos::View coords("intersection coords", mesh_1x1.nedges(), 3); + auto edge2vert = Omega_h::HostRead(mesh_1x1.get_adj(Omega_h::EDGE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (int i = 0; i < mesh_1x1.nedges(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[edge2vert[2*i]*3], meshcoords[edge2vert[2*i]*3 + 1], meshcoords[edge2vert[2*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[edge2vert[2*i + 1]*3], meshcoords[edge2vert[2*i + 1]*3 + 1], meshcoords[edge2vert[2*i + 1]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1)/2.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + + Kokkos::deep_copy(execs, coords, coords_h); + + pcms::CoordinateView> edge_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(edge_cv); + check_res(results, 1); + } + SECTION ("Face intersection") { + Kokkos::View coords("intersection coords", mesh_1x1.nfaces(), 3); + auto face2vert = Omega_h::HostRead(mesh_1x1.get_adj(Omega_h::FACE, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (int i = 0; i < mesh_1x1.nfaces(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[face2vert[3*i]*3], meshcoords[face2vert[3*i]*3 + 1], meshcoords[face2vert[3*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[face2vert[3*i + 1]*3], meshcoords[face2vert[3*i + 1]*3 + 1], meshcoords[face2vert[3*i + 1]*3 + 2] + }; + Omega_h::Vector<3> v2{ + meshcoords[face2vert[3*i + 2]*3], meshcoords[face2vert[3*i + 2]*3 + 1], meshcoords[face2vert[3*i + 2]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1 + v2)/3.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + Kokkos::deep_copy(execs, coords, coords_h); + + pcms::CoordinateView> face_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(face_cv); + check_res(results, 2); + } + SECTION ("Region intersection") { + Kokkos::View coords("intersection coords", mesh_1x1.nelems(), 3); + auto region2vert = Omega_h::HostRead(mesh_1x1.get_adj(Omega_h::REGION, Omega_h::VERT).ab2b); + auto meshcoords = Omega_h::HostRead(mesh_1x1.coords()); + auto coords_h = Kokkos::create_mirror_view(coords); + + for (int i = 0; i < mesh_1x1.nelems(); i++) + { + Omega_h::Vector<3> v0{ + meshcoords[region2vert[4*i]*3], meshcoords[region2vert[4*i]*3 + 1], meshcoords[region2vert[4*i]*3 + 2] + }; + Omega_h::Vector<3> v1{ + meshcoords[region2vert[4*i + 1]*3], meshcoords[region2vert[4*i + 1]*3 + 1], meshcoords[region2vert[4*i + 1]*3 + 2] + }; + Omega_h::Vector<3> v2{ + meshcoords[region2vert[4*i + 2]*3], meshcoords[region2vert[4*i + 2]*3 + 1], meshcoords[region2vert[4*i + 2]*3 + 2] + }; + Omega_h::Vector<3> v3{ + meshcoords[region2vert[4*i + 3]*3], meshcoords[region2vert[4*i + 3]*3 + 1], meshcoords[region2vert[4*i + 3]*3 + 2] + }; + Omega_h::Vector<3> midpoint = (v0 + v1 + v2 + v3)/4.; + coords_h(i, 0) = midpoint[0]; + coords_h(i, 1) = midpoint[1]; + coords_h(i, 2) = midpoint[2]; + } + Kokkos::deep_copy(execs, coords, coords_h); + + pcms::CoordinateView> region_cv( + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(coords)); + pcms::TreePointSearch::Results results = tree_search.apply(region_cv); + check_res(results, 3); + } + + // SECTION ("Region intersection") { + // Kokkos::View coords + // { + // -0.5, -0.5, -0.5, + // -0.5, -0.5, 1.5, + // -0.5, 1.5, -0.5, + // -0.5, 1.5, 1.5, + // 1.5, -0.5, -0.5, + // 1.5, -0.5, 1.5, + // 1.5, 1.5, -0.5, + // 1.5, 1.5, 1.5 + // }; + + // pcms::CoordinateView> outside_cv( + // pcms::CoordinateSystem::Cartesian, + // pcms::MakeConstRank2View(coords)); + // Kokkos::View results + // = tree_search.apply(outside_cv); +// auto results_h = Kokkos::create_mirror_view(results); +// Kokkos::deep_copy(execs, results_h, results); + // check_outside(results_h); + // } +}