diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index 8a45cf32c..37c4edb53 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -40,6 +40,7 @@ if(NOT WIN32) add_subdirectory(develop) add_subdirectory(union_find) endif() +add_subdirectory(triangulated_surface_distance) if (ARBORX_ENABLE_MPI) add_subdirectory(distributed_tree_driver) diff --git a/benchmarks/triangulated_surface_distance/CMakeLists.txt b/benchmarks/triangulated_surface_distance/CMakeLists.txt new file mode 100644 index 000000000..1a32bb8b6 --- /dev/null +++ b/benchmarks/triangulated_surface_distance/CMakeLists.txt @@ -0,0 +1,6 @@ +set(POINT_CLOUDS_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/benchmarks/point_clouds) + +add_executable(ArborX_Benchmark_TriangulatedSurfaceDistance.exe triangulated_surface_distance.cpp) +target_include_directories(ArborX_Benchmark_TriangulatedSurfaceDistance.exe PRIVATE ${POINT_CLOUDS_INCLUDE_DIR}) +target_link_libraries(ArborX_Benchmark_TriangulatedSurfaceDistance.exe ArborX::ArborX Boost::program_options) +add_test(NAME ArborX_Benchmark_TriangulatedSurfaceDistance COMMAND ArborX_Benchmark_TriangulatedSurfaceDistance.exe) diff --git a/benchmarks/triangulated_surface_distance/triangulated_surface_distance.cpp b/benchmarks/triangulated_surface_distance/triangulated_surface_distance.cpp new file mode 100644 index 000000000..e03a74af8 --- /dev/null +++ b/benchmarks/triangulated_surface_distance/triangulated_surface_distance.cpp @@ -0,0 +1,312 @@ +/**************************************************************************** + * Copyright (c) 2023 by the ArborX authors * + * All rights reserved. * + * * + * This file is part of the ArborX library. ArborX is * + * distributed under a BSD 3-clause license. For the licensing terms see * + * the LICENSE file in the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +using Point = ArborX::ExperimentalHyperGeometry::Point<3>; +using Triangle = ArborX::ExperimentalHyperGeometry::Triangle<3>; + +auto icosahedron() +{ + auto a = Kokkos::numbers::phi_v; + auto b = 1.f; + + std::vector vertices; + + vertices.push_back(Point{0, b, -a}); + vertices.push_back(Point{b, a, 0}); + vertices.push_back(Point{-b, a, 0}); + vertices.push_back(Point{0, b, a}); + vertices.push_back(Point{0, -b, a}); + vertices.push_back(Point{-a, 0, b}); + vertices.push_back(Point{0, -b, -a}); + vertices.push_back(Point{a, 0, -b}); + vertices.push_back(Point{a, 0, b}); + vertices.push_back(Point{-a, 0, -b}); + vertices.push_back(Point{b, -a, 0}); + vertices.push_back(Point{-b, -a, 0}); + + std::vector triangles; + + triangles.insert(triangles.end(), {2, 1, 0}); + triangles.insert(triangles.end(), {1, 2, 3}); + triangles.insert(triangles.end(), {5, 4, 3}); + triangles.insert(triangles.end(), {4, 8, 3}); + triangles.insert(triangles.end(), {7, 6, 0}); + triangles.insert(triangles.end(), {6, 9, 0}); + triangles.insert(triangles.end(), {11, 10, 4}); + triangles.insert(triangles.end(), {10, 11, 6}); + triangles.insert(triangles.end(), {9, 5, 2}); + triangles.insert(triangles.end(), {5, 9, 11}); + triangles.insert(triangles.end(), {8, 7, 1}); + triangles.insert(triangles.end(), {7, 8, 10}); + triangles.insert(triangles.end(), {2, 5, 3}); + triangles.insert(triangles.end(), {8, 1, 3}); + triangles.insert(triangles.end(), {9, 2, 0}); + triangles.insert(triangles.end(), {1, 7, 0}); + triangles.insert(triangles.end(), {11, 9, 6}); + triangles.insert(triangles.end(), {7, 10, 6}); + triangles.insert(triangles.end(), {5, 11, 4}); + triangles.insert(triangles.end(), {10, 8, 4}); + + return std::make_pair(vertices, triangles); +} + +/* Subdivide every triangle into four + /\ /\ + / \ / \ + / \ ---> /____\ + / \ /\ /\ + / \ / \ / \ + /__________\ /____\/____\ +*/ +void subdivide(std::vector &vertices, std::vector &triangles) +{ + std::map, int> hash; + + int const num_triangles = triangles.size() / 3; + + int vindex = vertices.size(); + std::vector new_triangles; + for (int i = 0; i < num_triangles; ++i) + { + int v[3] = {triangles[3 * i + 0], triangles[3 * i + 1], + triangles[3 * i + 2]}; + int vmid[3]; + + for (int j = 0; j < 3; ++j) + { + int vmin = std::min(v[j], v[(j + 1) % 3]); + int vmax = std::max(v[j], v[(j + 1) % 3]); + + auto it = hash.find(std::make_pair(vmin, vmax)); + if (it != hash.end()) + { + vmid[j] = it->second; + } + else + { + vmid[j] = vindex; + hash[std::make_pair(vmin, vmax)] = vindex; + ++vindex; + + vertices.push_back(Point{(vertices[vmin][0] + vertices[vmax][0]) / 2, + (vertices[vmin][1] + vertices[vmax][1]) / 2, + (vertices[vmin][2] + vertices[vmax][2]) / 2}); + } + } + + new_triangles.insert(new_triangles.end(), {v[0], vmid[0], vmid[2]}); + new_triangles.insert(new_triangles.end(), {v[1], vmid[0], vmid[1]}); + new_triangles.insert(new_triangles.end(), {v[2], vmid[2], vmid[1]}); + new_triangles.insert(new_triangles.end(), {vmid[0], vmid[1], vmid[2]}); + } + triangles = new_triangles; +} + +void projectVerticesToSphere(std::vector &vertices, float radius) +{ + for (auto &v : vertices) + { + auto norm = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + v[0] *= radius / norm; + v[1] *= radius / norm; + v[2] *= radius / norm; + } +} + +auto buildTriangles(float radius, int num_refinements) +{ + Kokkos::Profiling::ScopedRegion guard("Benchmark::build_triangles"); + + auto [vertices, triangles] = icosahedron(); + for (int i = 1; i <= num_refinements; ++i) + subdivide(vertices, triangles); + + projectVerticesToSphere(vertices, radius); + + return std::make_pair(vertices, triangles); +} + +template +auto vec2view(std::vector const &in, std::string const &label = "") +{ + Kokkos::View out( + Kokkos::view_alloc(label, Kokkos::WithoutInitializing), in.size()); + Kokkos::deep_copy(out, Kokkos::View>{ + in.data(), in.size()}); + return out; +} + +template +struct Triangles +{ + Kokkos::View _points; + Kokkos::View _triangle_vertices; +}; + +template +class ArborX::AccessTraits, ArborX::PrimitivesTag> +{ + using Self = Triangles; + +public: + using memory_space = MemorySpace; + + static KOKKOS_FUNCTION auto size(Self const &self) + { + return self._triangle_vertices.size() / 3; + } + static KOKKOS_FUNCTION auto get(Self const &self, int i) + { + auto const &vertices = self._triangle_vertices; + return Triangle{self._points(vertices(3 * i + 0)), + self._points(vertices(3 * i + 1)), + self._points(vertices(3 * i + 2))}; + } +}; + +struct DistanceCallback +{ + template + KOKKOS_FUNCTION void operator()(Predicate const &predicate, + Value const &value, + OutputFunctor const &out) const + { + using ArborX::Details::distance; + out(distance(ArborX::getGeometry(predicate), value)); + } +}; + +template +void writeVtk(std::string const &filename, Points const &vertices, + Triangles const &triangles) +{ + int const num_vertices = vertices.size(); + int const num_elements = triangles.size() / 3; + + constexpr int DIM = + ArborX::GeometryTraits::dimension_v; + + auto vertices_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, vertices); + auto triangles_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, triangles); + + std::ofstream out(filename); + + out << "# vtk DataFile Version 2.0\n"; + out << "Mesh example\n"; + out << "ASCII\n"; + out << "DATASET POLYDATA\n\n"; + out << "POINTS " << num_vertices << " float\n"; + for (int i = 0; i < num_vertices; ++i) + { + for (int d = 0; d < DIM; ++d) + out << " " << vertices(i)[d]; + out << '\n'; + } + + int const num_cell_vertices = 3; + out << "\nPOLYGONS " << num_elements << " " + << (num_elements * (1 + num_cell_vertices)) << '\n'; + for (int i = 0; i < num_elements; ++i) + { + out << num_cell_vertices; + for (int j = 0; j < num_cell_vertices; ++j) + out << " " << triangles(i * num_cell_vertices + j); + out << '\n'; + } +} + +int main(int argc, char *argv[]) +{ + Kokkos::ScopeGuard guard(argc, argv); + + using ExecutionSpace = Kokkos::DefaultExecutionSpace; + using MemorySpace = typename ExecutionSpace::memory_space; + + std::cout << "ArborX version : " << ArborX::version() << '\n'; + std::cout << "ArborX hash : " << ArborX::gitCommitHash() << '\n'; + std::cout << "Kokkos version : " << ArborX::Details::KokkosExt::version() + << '\n'; + + namespace bpo = boost::program_options; + + bpo::options_description desc("Allowed options"); + int n; + int num_refinements; + float radius; + std::string vtk_filename; + // clang-format off + desc.add_options() + ( "help", "help message" ) + ( "n", bpo::value(&n)->default_value(1000), "number of points" ) + ( "radius", bpo::value(&radius)->default_value(1.f), "sphere radius" ) + ( "refinements", bpo::value(&num_refinements)->default_value(5), "number of icosahedron refinements" ) + ( "vtk-filename", bpo::value(&vtk_filename), "filename to dump mesh to in VTK format" ) + ; + // clang-format on + bpo::variables_map vm; + bpo::store(bpo::command_line_parser(argc, argv).options(desc).run(), vm); + bpo::notify(vm); + + if (vm.count("help") > 0) + { + std::cout << desc << '\n'; + return 1; + } + + auto [vertices_v, triangles_v] = buildTriangles(radius, num_refinements); + + auto vertices = vec2view(vertices_v); + auto triangles = vec2view(triangles_v); + + if (!vtk_filename.empty()) + writeVtk(vtk_filename, vertices, triangles); + + ExecutionSpace space; + + Kokkos::Profiling::pushRegion("Benchmark::build_points"); + Kokkos::View random_points( + Kokkos::view_alloc(space, Kokkos::WithoutInitializing, + "Benchmark::points"), + n); + generatePointCloud(PointCloudType::filled_box, std::cbrt(n), random_points); + Kokkos::Profiling::popRegion(); + + std::cout << "#triangles : " << triangles.size() << '\n'; + std::cout << "#queries : " << random_points.size() << '\n'; + + ArborX::BVH index( + space, Triangles{vertices, triangles}); + + Kokkos::View offset("Benchmark::offsets", 0); + Kokkos::View distances("Benchmark::distances", 0); + index.query(space, ArborX::Experimental::make_nearest(random_points, 1), + DistanceCallback{}, distances, offset); + + return 0; +}