diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2450f2973..a9c5442a1 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -11,4 +11,4 @@ repos: - -i - --style=file - --verbose - files: '\.(c|cc|cpp|h|hpp|tpp|cxx|hh|inl|ipp)$' + files: '\.(c|cc|cpp|cu|cuh|h|hpp|tpp|cxx|hh|inl|ipp)$' diff --git a/CMakeLists.txt b/CMakeLists.txt index 68fd32051..66b98594c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -344,14 +344,13 @@ if(IPC_TOOLKIT_WITH_SIMD) # compiled twice and the nvcc pass owns the float and double # instantiations (see ipc_toolkit_shared_device_sources.cmake), so without # this they would be the only scalar code in the library built without - # SIMD support. Hand the flags to the host compiler nvcc drives rather - # than to nvcc itself. - set(SIMD_CUDA_HOST_FLAGS "") - foreach(simd_flag IN LISTS SIMD_CXX_FLAGS) - list(APPEND SIMD_CUDA_HOST_FLAGS "-Xcompiler=${simd_flag}") - endforeach() - target_compile_options(ipc_toolkit PUBLIC - "$<$,$>:${SIMD_CUDA_HOST_FLAGS}>") + # SIMD support. ipc_toolkit_filter_nvcc_flags() hands the flags to the + # host compiler nvcc drives (via -Xcompiler) rather than to nvcc itself, + # keeping only those that pass compiles with. + include(ipc_toolkit_filter_flags) + set(SIMD_CUDA_HOST_FLAGS ${SIMD_CXX_FLAGS}) + ipc_toolkit_filter_nvcc_flags(SIMD_CUDA_HOST_FLAGS) + target_compile_options(ipc_toolkit PUBLIC ${SIMD_CUDA_HOST_FLAGS}) endif() # Link against cross-platform xsimd library diff --git a/cmake/ipc_toolkit/ipc_toolkit_filter_flags.cmake b/cmake/ipc_toolkit/ipc_toolkit_filter_flags.cmake index b6dbcd0ff..36583d36e 100644 --- a/cmake/ipc_toolkit/ipc_toolkit_filter_flags.cmake +++ b/cmake/ipc_toolkit/ipc_toolkit_filter_flags.cmake @@ -22,4 +22,27 @@ function(ipc_toolkit_filter_flags flags) endif() endforeach() set(${flags} ${output_flags} PARENT_SCOPE) -endfunction() \ No newline at end of file +endfunction() + +# The nvcc counterpart of ipc_toolkit_filter_flags(): keep the flags of `flags` +# that nvcc's host compiler accepts, checked by actually compiling with nvcc and +# `-Xcompiler=`, and wrap each so it applies only to CUDA sources compiled +# by nvcc. `-Xcompiler` is required: nvcc parses some host flags itself with a +# different meaning (`-Werror` takes nvcc's own diagnostic names, `-march=...` +# is read as an input file), so a bare host flag on the nvcc command line is +# unsafe. Requires the CUDA language to be enabled. +function(ipc_toolkit_filter_nvcc_flags flags) + include(CheckCompilerFlag) + set(output_flags) + foreach(FLAG IN ITEMS ${${flags}}) + string(REPLACE "=" "-" FLAG_VAR "${FLAG}") + if(NOT DEFINED IS_SUPPORTED_NVCC_HOST_${FLAG_VAR}) + check_compiler_flag(CUDA "-Xcompiler=${FLAG}" IS_SUPPORTED_NVCC_HOST_${FLAG_VAR}) + endif() + if(IS_SUPPORTED_NVCC_HOST_${FLAG_VAR}) + list(APPEND output_flags + "$<$,$>:-Xcompiler=${FLAG}>") + endif() + endforeach() + set(${flags} ${output_flags} PARENT_SCOPE) +endfunction() diff --git a/cmake/ipc_toolkit/ipc_toolkit_warnings.cmake b/cmake/ipc_toolkit/ipc_toolkit_warnings.cmake index 22928c570..ff799263b 100644 --- a/cmake/ipc_toolkit/ipc_toolkit_warnings.cmake +++ b/cmake/ipc_toolkit/ipc_toolkit_warnings.cmake @@ -195,3 +195,54 @@ add_library(ipc::toolkit::warnings ALIAS ipc_toolkit_warnings) include(ipc_toolkit_filter_flags) ipc_toolkit_filter_flags(IPC_TOOLKIT_WARNING_FLAGS) target_compile_options(ipc_toolkit_warnings INTERFACE ${IPC_TOOLKIT_WARNING_FLAGS}) + +# nvcc forwards none of the flags above to the host compiler it drives, so a +# .cu is otherwise compiled with no warnings at all. Hand the host pass the +# flags that catch bugs in OUR code (the device pass has no use for them). +# +# The full set is deliberately not forwarded: nvcc's generated host code and the +# CUDA headers are not clean under -Wpedantic, -Wold-style-cast or -Wsign-promo +# (thousands of hits in its stubs and in crt/device_functions.hpp), and its +# rewriting of aggregate initializers trips -Werror=missing-braces. +# +# ipc_toolkit_filter_nvcc_flags() checks each flag against nvcc and scopes it +# to CUDA sources compiled by nvcc (clang as the CUDA compiler has no -Xcompiler +# and takes the C++ flags above directly). +if(IPC_TOOLKIT_WITH_CUDA AND NOT MSVC) + set(IPC_TOOLKIT_NVCC_HOST_WARNING_FLAGS + -Wall + -Wextra + -Wshadow + -Woverloaded-virtual + -Wuninitialized + -Wcast-qual + -Wpointer-arith + -Werror=return-type + -Werror=non-virtual-dtor + -Werror=delete-non-virtual-dtor + -Wno-unused-parameter + -Wno-sign-compare + -Wno-unknown-pragmas # nvcc's own pragmas (unroll) reach the host pass + ) + ipc_toolkit_filter_nvcc_flags(IPC_TOOLKIT_NVCC_HOST_WARNING_FLAGS) + target_compile_options(ipc_toolkit_warnings INTERFACE ${IPC_TOOLKIT_NVCC_HOST_WARNING_FLAGS}) +endif() + +# The device pass has no -Wall; nvcc's own diagnostics are few, and these are +# the ones with a bug-finding record here. They go to nvcc directly (not via +# -Xcompiler): nvcc's -Werror takes nvcc's own diagnostic names. +# * all-warnings: the front end's own warnings (e.g. "variable used before its +# value is set") become errors. +# * cross-execution-space-call: calling a host-only function from device code +# is an error, the typical mistake in code shared via IPC_TOOLKIT_HOST_DEVICE. +# Not used: -Xptxas -warn-spills. Under separable compilation the shared device +# sources are relocatable __device__ functions that follow the ABI, and ptxas +# reports their callee-saved register traffic as spills (64 hits in the barrier +# functions alone), so it is noise here; check kernels with --resource-usage on +# a non-rdc compile instead. -warn-lmem-usage likewise: the traversal stack is +# intentional local memory. +if(IPC_TOOLKIT_WITH_CUDA) + target_compile_options(ipc_toolkit_warnings INTERFACE + "$<$,$>:SHELL:-Werror all-warnings>" + "$<$,$>:SHELL:-Werror cross-execution-space-call>") +endif() diff --git a/cmake/recipes/spdlog.cmake b/cmake/recipes/spdlog.cmake index 5e083f888..d3be8c70c 100644 --- a/cmake/recipes/spdlog.cmake +++ b/cmake/recipes/spdlog.cmake @@ -1,18 +1,6 @@ # spdlog (https://github.com/gabime/spdlog) # License: MIT if(TARGET spdlog::spdlog) - # Someone else created the target, so the fmt patch below never runs. That - # is fine without CUDA, but the patch is what makes the bundled fmt - # compile under nvcc at all, so warn rather than fail at the first .cu. - if(IPC_TOOLKIT_WITH_CUDA) - message(WARNING - "spdlog::spdlog was provided by an enclosing project, so IPC " - "Toolkit's cmake/patches/fmt-nvcc-compat.patch was not applied. " - "The bundled fmt does not compile under nvcc unpatched: its " - "literal-encoding probe misfires and a char32_t table uses hex " - "escapes with the high bit set. Apply the same patch to your " - "spdlog, or let IPC Toolkit fetch its own.") - endif() return() endif() diff --git a/docs/source/Doxyfile b/docs/source/Doxyfile index 6a7d7d1d2..4c881b0da 100644 --- a/docs/source/Doxyfile +++ b/docs/source/Doxyfile @@ -2315,7 +2315,7 @@ INCLUDE_FILE_PATTERNS = # recursively expanded use the := operator instead of the = operator. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -PREDEFINED = IPC_TOOLKIT_WITH_INEXACT_CCD IPC_TOOLKIT_WITH_ROBIN_MAP IPC_TOOLKIT_WITH_ABSEIL IPC_TOOLKIT_WITH_FILIB IPC_TOOLKIT_WITH_MESHFEM_SPARSE +PREDEFINED = IPC_TOOLKIT_WITH_CUDA IPC_TOOLKIT_WITH_INEXACT_CCD IPC_TOOLKIT_WITH_ROBIN_MAP IPC_TOOLKIT_WITH_ABSEIL IPC_TOOLKIT_WITH_FILIB IPC_TOOLKIT_WITH_MESHFEM_SPARSE # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this # tag can be used to specify a list of macro names that should be expanded. The diff --git a/docs/source/about/release_notes.rst b/docs/source/about/release_notes.rst index 08793fa8b..ee1acda7c 100644 --- a/docs/source/about/release_notes.rst +++ b/docs/source/about/release_notes.rst @@ -25,10 +25,17 @@ New Features |:rocket:| - Expose the intersection coordinates of an edge–triangle intersection through a new :cpp:func:`ipc::edge_triangle_intersection` overload, which reports the barycentric coordinates :math:`(u, v)` on the triangle and the parameter :math:`t` along the edge (`#245 `_). - Add :cpp:func:`ipc::CollisionMesh::face_normals`, computing the unit normal of each face for a given set of vertex positions (3D only) (`#245 `_). +- Add :cpp:class:`ipc::cuda::LBVH`, a GPU broad phase that builds the vertex/edge/face AABBs and BVHs and runs the traversal and shared-vertex filtering on the device, producing the same candidates as :cpp:class:`ipc::LBVH` for any vertex filter (`#260 `_). Requires ``IPC_TOOLKIT_WITH_CUDA``; selectable through ``BroadPhaseMethod::LBVH_CUDA``. + + - The Apetrei :cite:p:`Apetrei2014FastAS` bottom-up build, the BVH descent, and the shared-vertex exclusion are one ``ipc::details`` implementation shared by the CPU and CUDA broad phases; each platform supplies only its parallel launch, its sort, and its atomics. + - ``detect_*_candidates_device()`` return a view of the candidate pairs left on the device for a GPU-native pipeline. API Changes |:wrench:| ~~~~~~~~~~~~~~~~~~~~~~ +- ``BroadPhase::detect_*_candidates()`` now uniformly **clear** their output vector first on every broad phase, so it holds exactly that detection's result (`#260 `_). Previously half the implementations overwrote and half appended; all in-library callers pass an empty vector, so their results are unchanged. +- Add :cpp:func:`ipc::CollisionFilter::accepts_all`, true for a filter that holds no predicate (the default), so a broad phase can skip per-pair filtering entirely (`#260 `_). Composing with an accept-all filter now short-circuits: ``f | accept_all`` is accept-all and ``f & accept_all`` is ``f``. +- Add ``BroadPhaseMethod::NUM_BROAD_PHASE_METHODS`` as a sentinel for the number of methods (`#260 `_). - Update Tight Inclusion from ``1.0.6`` to ``1.1.0`` (`#248 `_). - Adds a ``BUCKET_DEPTH_FIRST_SEARCH`` root-finding method, which upstream makes the default for ``edgeEdgeCCD`` and ``vertexFaceCCD``. @@ -186,6 +193,7 @@ Python |:snake:| - Validate preconditions in the bindings instead of relying on the C++ ``assert``\ s, which are compiled out under ``NDEBUG`` and would let a release build silently accept a bad value (`#247 `_). ``BarrierPotential`` now raises ``ValueError`` for a non-positive or NaN ``dhat``/``stiffness`` and for a null barrier. - Bind ``edge_triangle_intersection()``, returning an ``(intersects, u, v, t)`` tuple since Python has no out-parameters, and ``CollisionMesh.face_normals()``, returning an (#F × 3) array to match the other per-element accessors (`#245 `_). ``face_normals()`` raises ``ValueError`` on a 2D mesh rather than invoking undefined behavior. +- Add the ``ipctk.cuda`` submodule, mirroring the C++ ``ipc::cuda`` namespace, with ``ipctk.cuda.LBVH`` (``ipc::cuda::LBVH``) in CUDA builds (`#260 `_). Documentation ~~~~~~~~~~~~~ diff --git a/docs/source/cpp-api/broad_phase.rst b/docs/source/cpp-api/broad_phase.rst index 70c4f5278..ec0e28ad7 100644 --- a/docs/source/cpp-api/broad_phase.rst +++ b/docs/source/cpp-api/broad_phase.rst @@ -40,7 +40,18 @@ Sweep and Prune Sweep and Tiniest Queue ----------------------- -.. .. doxygenclass:: ipc::SweepAndTiniestQueueGPU +Requires CUDA (``IPC_TOOLKIT_WITH_CUDA``). + +.. doxygenclass:: ipc::SweepAndTiniestQueue + :allow-dot-graphs: + +LBVH (CUDA) +----------- + +Requires CUDA (``IPC_TOOLKIT_WITH_CUDA``). + +.. doxygenclass:: ipc::cuda::LBVH + :allow-dot-graphs: AABB ---- diff --git a/docs/source/python-api/broad_phase.rst b/docs/source/python-api/broad_phase.rst index 022a01668..a55fe3efe 100644 --- a/docs/source/python-api/broad_phase.rst +++ b/docs/source/python-api/broad_phase.rst @@ -46,8 +46,24 @@ Sweep and Prune Sweep and Tiniest Queue ----------------------- -.. .. autoclass:: ipctk.SweepAndTiniestQueueGPU -.. :members: +``ipctk.SweepAndTiniestQueue`` is available only when ``ipctk`` is built with +CUDA (``IPC_TOOLKIT_WITH_CUDA``), which the documentation build is not. + +.. .. autoclass:: ipctk.SweepAndTiniestQueue +.. +.. .. autoclasstoc:: + +LBVH (CUDA) +----------- + +``ipctk.cuda.LBVH`` is the GPU counterpart of ``ipctk.LBVH`` (C++: +``ipc::cuda::LBVH``). The ``ipctk.cuda`` submodule mirrors the C++ ``ipc::cuda`` +namespace; its classes exist only when ``ipctk`` is built with CUDA +(``IPC_TOOLKIT_WITH_CUDA``), which the documentation build is not. + +.. .. autoclass:: ipctk.cuda.LBVH +.. +.. .. autoclasstoc:: AABB ---- diff --git a/docs/source/tutorials/getting_started.rst b/docs/source/tutorials/getting_started.rst index efc83a52a..5487ef413 100644 --- a/docs/source/tutorials/getting_started.rst +++ b/docs/source/tutorials/getting_started.rst @@ -539,7 +539,7 @@ The ``Candidates`` class represents the culled set of candidate pairs and is bui collision_mesh, vertices_t0, vertices_t1, broad_phase=ipctk.LBVH()) -Possible values for ``broad_phase`` are: ``BruteForce`` (parallel brute force culling), ``HashGrid``, ``SpatialHash`` (implementation from the original IPC codebase), ``LBVH`` (CPU implementation of :cite:t:`Karras2012HPG` using TBB), ``SweepAndPrune`` (a.k.a. Sort-and-Sweep from :cite:t:`Baraff1992PhD`), or ``SweepAndTiniestQueue`` (method of :cite:t:`Belgrod2023Time`; requires CUDA). The default is ``LBVH``. +Possible values for ``broad_phase`` are: ``BruteForce`` (parallel brute force culling), ``HashGrid``, ``SpatialHash`` (implementation from the original IPC codebase), ``LBVH`` (CPU implementation of :cite:t:`Karras2012HPG` using TBB), ``SweepAndPrune`` (a.k.a. Sort-and-Sweep from :cite:t:`Baraff1992PhD`), ``SweepAndTiniestQueue`` (method of :cite:t:`Belgrod2023Time`; requires CUDA), or ``cuda.LBVH`` (``ipc::cuda::LBVH`` in C++; the GPU counterpart of ``LBVH``, building and traversing the same hierarchy on the device; requires CUDA). The default is ``LBVH``. Narrow-Phase ^^^^^^^^^^^^ diff --git a/python/src/bindings.cpp b/python/src/bindings.cpp index 8ab40e974..06206814d 100644 --- a/python/src/bindings.cpp +++ b/python/src/bindings.cpp @@ -32,6 +32,12 @@ PYBIND11_MODULE(ipctk, m) define_sweep_and_tiniest_queue(m); define_voxel_size_heuristic(m); + // GPU implementations, mirroring the C++ ipc::cuda namespace. The submodule + // always exists; its classes are defined only in CUDA builds. + py::module_ cuda = m.def_submodule( + "cuda", "GPU (CUDA) implementations; populated only in CUDA builds."); + define_cuda_lbvh(cuda); + // candidates define_candidates(m); define_collision_stencil(m); diff --git a/python/src/broad_phase/CMakeLists.txt b/python/src/broad_phase/CMakeLists.txt index dd0db17f2..72b6700d4 100644 --- a/python/src/broad_phase/CMakeLists.txt +++ b/python/src/broad_phase/CMakeLists.txt @@ -2,6 +2,7 @@ set(SOURCES aabb.cpp broad_phase.cpp brute_force.cpp + cuda_lbvh.cpp hash_grid.cpp lbvh.cpp spatial_hash.cpp diff --git a/python/src/broad_phase/bindings.hpp b/python/src/broad_phase/bindings.hpp index 13edeeeb5..89cbd2dc2 100644 --- a/python/src/broad_phase/bindings.hpp +++ b/python/src/broad_phase/bindings.hpp @@ -5,6 +5,7 @@ void define_aabb(py::module_& m); void define_broad_phase(py::module_& m); void define_brute_force(py::module_& m); +void define_cuda_lbvh(py::module_& cuda); // into the ipctk.cuda submodule void define_hash_grid(py::module_& m); void define_lbvh(py::module_& m); void define_spatial_hash(py::module_& m); diff --git a/python/src/broad_phase/cuda_lbvh.cpp b/python/src/broad_phase/cuda_lbvh.cpp new file mode 100644 index 000000000..948e034f5 --- /dev/null +++ b/python/src/broad_phase/cuda_lbvh.cpp @@ -0,0 +1,32 @@ +#include + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA +using namespace ipc; // not defined if IPC_TOOLKIT_WITH_CUDA is not defined +#endif + +void define_cuda_lbvh(py::module_& m) // m is the ipctk.cuda submodule +{ +#ifdef IPC_TOOLKIT_WITH_CUDA + py::class_>( + m, "LBVH", + R"ipc_Qu8mg5v7( + GPU Linear Bounding Volume Hierarchy (LBVH) broad phase (ipc::cuda::LBVH). + + Builds the vertex/edge/face AABBs and BVHs, and runs the traversal and + mesh-connectivity filtering, on the device. Produces the same candidates + as ipctk.LBVH for any vertex filter. Available only in CUDA builds. + )ipc_Qu8mg5v7") + .def(py::init()) + .def_property_readonly( + "num_vertex_nodes", &cuda::LBVH::num_vertex_nodes, + "Number of nodes in the vertex BVH (2 * n_leaves - 1, or 0).") + .def_property_readonly( + "num_edge_nodes", &cuda::LBVH::num_edge_nodes, + "Number of nodes in the edge BVH (2 * n_leaves - 1, or 0).") + .def_property_readonly( + "num_face_nodes", &cuda::LBVH::num_face_nodes, + "Number of nodes in the face BVH (2 * n_leaves - 1, or 0)."); +#endif +} diff --git a/python/tests/utils.py b/python/tests/utils.py index 10d89fb8d..0dfec8065 100644 --- a/python/tests/utils.py +++ b/python/tests/utils.py @@ -45,6 +45,8 @@ def broad_phases(): yield ipctk.SpatialHash() yield ipctk.LBVH() yield ipctk.SweepAndPrune() + if hasattr(ipctk.cuda, "LBVH"): # only in CUDA builds + yield ipctk.cuda.LBVH() def finite_jacobian(x, f, h=1e-8): diff --git a/src/ipc/broad_phase/CMakeLists.txt b/src/ipc/broad_phase/CMakeLists.txt index e0613ebcf..1c2346764 100644 --- a/src/ipc/broad_phase/CMakeLists.txt +++ b/src/ipc/broad_phase/CMakeLists.txt @@ -23,3 +23,8 @@ set(SOURCES ) target_sources(ipc_toolkit PRIVATE ${SOURCES}) + +add_subdirectory(details) +if(IPC_TOOLKIT_WITH_CUDA) + add_subdirectory(cuda) +endif() diff --git a/src/ipc/broad_phase/aabb.cpp b/src/ipc/broad_phase/aabb.cpp index 4716988a4..2a5edefcd 100644 --- a/src/ipc/broad_phase/aabb.cpp +++ b/src/ipc/broad_phase/aabb.cpp @@ -44,13 +44,11 @@ void AABB::conservative_inflation( // Nudge the bounds outward to ensure conservativity. min = min.unaryExpr([inflation_radius](double v) { - return std::nextafter( - v - inflation_radius, -std::numeric_limits::infinity()); + return conservative_lower_bound(v, inflation_radius); }); max = max.unaryExpr([inflation_radius](double v) { - return std::nextafter( - v + inflation_radius, std::numeric_limits::infinity()); + return conservative_upper_bound(v, inflation_radius); }); } diff --git a/src/ipc/broad_phase/aabb.hpp b/src/ipc/broad_phase/aabb.hpp index ab03b6492..bc38bfcd4 100644 --- a/src/ipc/broad_phase/aabb.hpp +++ b/src/ipc/broad_phase/aabb.hpp @@ -3,8 +3,10 @@ #include #include #include +#include // for infinity() #include +#include // for nextafter namespace ipc { @@ -60,6 +62,33 @@ class alignas(64) AABB { Eigen::Ref max, const double inflation_radius); + /// @brief Conservatively inflate one lower bound. + /// + /// The single-coordinate policy behind conservative_inflation(): the bound + /// is moved out by the radius, then nudged to the next representable double + /// away from the box so rounding can never shrink it. Host/device so the + /// GPU broad phases build bit-identical boxes. + /// + /// @param v The coordinate to bound from below. + /// @param inflation_radius The radius to inflate by. + /// @return The conservative lower bound. + IPC_TOOLKIT_HOST_DEVICE static double + conservative_lower_bound(const double v, const double inflation_radius) + { + return nextafter(v - inflation_radius, -infinity()); + } + + /// @brief Conservatively inflate one upper bound. + /// @see conservative_lower_bound + /// @param v The coordinate to bound from above. + /// @param inflation_radius The radius to inflate by. + /// @return The conservative upper bound. + IPC_TOOLKIT_HOST_DEVICE static double + conservative_upper_bound(const double v, const double inflation_radius) + { + return nextafter(v + inflation_radius, infinity()); + } + public: /// @brief Minimum corner of the AABB. Eigen::Array3d min; diff --git a/src/ipc/broad_phase/broad_phase.cpp b/src/ipc/broad_phase/broad_phase.cpp index 114dfbc8f..5f00b3dda 100644 --- a/src/ipc/broad_phase/broad_phase.cpp +++ b/src/ipc/broad_phase/broad_phase.cpp @@ -1,6 +1,7 @@ #include "broad_phase.hpp" #include +#include #include #include @@ -131,8 +132,7 @@ bool BroadPhase::can_edge_vertex_collide(size_t ei, size_t vi) const assert(ei < edge_boxes.size()); const auto& [e0i, e1i, _] = edge_boxes[ei].vertex_ids; - return vi != e0i && vi != e1i - && (can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i)); + return details::can_edge_vertex_collide(e0i, e1i, vi, can_vertices_collide); } bool BroadPhase::can_edges_collide(size_t eai, size_t ebi) const @@ -142,13 +142,8 @@ bool BroadPhase::can_edges_collide(size_t eai, size_t ebi) const assert(ebi < edge_boxes.size()); const auto& [eb0i, eb1i, __] = edge_boxes[ebi].vertex_ids; - const bool share_endpoint = - ea0i == eb0i || ea0i == eb1i || ea1i == eb0i || ea1i == eb1i; - - return !share_endpoint - && (can_vertices_collide(ea0i, eb0i) || can_vertices_collide(ea0i, eb1i) - || can_vertices_collide(ea1i, eb0i) - || can_vertices_collide(ea1i, eb1i)); + return details::can_edges_collide( + ea0i, ea1i, eb0i, eb1i, can_vertices_collide); } bool BroadPhase::can_face_vertex_collide(size_t fi, size_t vi) const @@ -156,9 +151,8 @@ bool BroadPhase::can_face_vertex_collide(size_t fi, size_t vi) const assert(fi < face_boxes.size()); const auto& [f0i, f1i, f2i] = face_boxes[fi].vertex_ids; - return vi != f0i && vi != f1i && vi != f2i - && (can_vertices_collide(vi, f0i) || can_vertices_collide(vi, f1i) - || can_vertices_collide(vi, f2i)); + return details::can_face_vertex_collide( + f0i, f1i, f2i, vi, can_vertices_collide); } bool BroadPhase::can_edge_face_collide(size_t ei, size_t fi) const @@ -168,14 +162,8 @@ bool BroadPhase::can_edge_face_collide(size_t ei, size_t fi) const assert(fi < face_boxes.size()); const auto& [f0i, f1i, f2i] = face_boxes[fi].vertex_ids; - const bool share_endpoint = e0i == f0i || e0i == f1i || e0i == f2i - || e1i == f0i || e1i == f1i || e1i == f2i; - - return !share_endpoint - && (can_vertices_collide(e0i, f0i) || can_vertices_collide(e0i, f1i) - || can_vertices_collide(e0i, f2i) || can_vertices_collide(e1i, f0i) - || can_vertices_collide(e1i, f1i) - || can_vertices_collide(e1i, f2i)); + return details::can_edge_face_collide( + e0i, e1i, f0i, f1i, f2i, can_vertices_collide); } bool BroadPhase::can_faces_collide(size_t fai, size_t fbi) const @@ -185,20 +173,8 @@ bool BroadPhase::can_faces_collide(size_t fai, size_t fbi) const assert(fbi < face_boxes.size()); const auto& [fb0i, fb1i, fb2i] = face_boxes[fbi].vertex_ids; - const bool share_endpoint = fa0i == fb0i || fa0i == fb1i || fa0i == fb2i - || fa1i == fb0i || fa1i == fb1i || fa1i == fb2i || fa2i == fb0i - || fa2i == fb1i || fa2i == fb2i; - - return !share_endpoint - && (can_vertices_collide(fa0i, fb0i) // - || can_vertices_collide(fa0i, fb1i) - || can_vertices_collide(fa0i, fb2i) - || can_vertices_collide(fa1i, fb0i) - || can_vertices_collide(fa1i, fb1i) - || can_vertices_collide(fa1i, fb2i) - || can_vertices_collide(fa2i, fb0i) - || can_vertices_collide(fa2i, fb1i) - || can_vertices_collide(fa2i, fb2i)); + return details::can_faces_collide( + fa0i, fa1i, fa2i, fb0i, fb1i, fb2i, can_vertices_collide); } } // namespace ipc diff --git a/src/ipc/broad_phase/broad_phase.hpp b/src/ipc/broad_phase/broad_phase.hpp index 1b7ea96b5..5f333f6c6 100644 --- a/src/ipc/broad_phase/broad_phase.hpp +++ b/src/ipc/broad_phase/broad_phase.hpp @@ -63,36 +63,36 @@ class BroadPhase { virtual void clear(); /// @brief Detect all collision candidates needed for a given dimensional simulation. - /// @param candidates The detected collision candidates. + /// @param candidates The detected collision candidates (cleared first). void detect_collision_candidates(Candidates& candidates) const; /// @brief Find the candidate vertex-vertex collisions. - /// @param[out] candidates The candidate vertex-vertex collisions. + /// @param[out] candidates The candidate vertex-vertex collisions (cleared first). virtual void detect_vertex_vertex_candidates( std::vector& candidates) const = 0; /// @brief Find the candidate edge-vertex collisions. - /// @param[out] candidates The candidate edge-vertex collisions. + /// @param[out] candidates The candidate edge-vertex collisions (cleared first). virtual void detect_edge_vertex_candidates( std::vector& candidates) const = 0; /// @brief Find the candidate edge-edge collisions. - /// @param[out] candidates The candidate edge-edge collisions. + /// @param[out] candidates The candidate edge-edge collisions (cleared first). virtual void detect_edge_edge_candidates( std::vector& candidates) const = 0; /// @brief Find the candidate face-vertex collisions. - /// @param[out] candidates The candidate face-vertex collisions. + /// @param[out] candidates The candidate face-vertex collisions (cleared first). virtual void detect_face_vertex_candidates( std::vector& candidates) const = 0; /// @brief Find the candidate edge-face intersections. - /// @param[out] candidates The candidate edge-face intersections. + /// @param[out] candidates The candidate edge-face intersections (cleared first). virtual void detect_edge_face_candidates( std::vector& candidates) const = 0; /// @brief Find the candidate face-face collisions. - /// @param[out] candidates The candidate face-face collisions. + /// @param[out] candidates The candidate face-face collisions (cleared first). virtual void detect_face_face_candidates( std::vector& candidates) const = 0; diff --git a/src/ipc/broad_phase/brute_force.cpp b/src/ipc/broad_phase/brute_force.cpp index 166380c55..9e9425254 100644 --- a/src/ipc/broad_phase/brute_force.cpp +++ b/src/ipc/broad_phase/brute_force.cpp @@ -63,6 +63,7 @@ void BruteForce::detect_candidates( void BruteForce::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( vertex_boxes, vertex_boxes, can_vertices_collide, candidates); } @@ -70,6 +71,7 @@ void BruteForce::detect_vertex_vertex_candidates( void BruteForce::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_boxes, vertex_boxes, std::bind(&BruteForce::can_edge_vertex_collide, this, _1, _2), @@ -79,6 +81,7 @@ void BruteForce::detect_edge_vertex_candidates( void BruteForce::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_boxes, edge_boxes, std::bind(&BruteForce::can_edges_collide, this, _1, _2), candidates); @@ -87,6 +90,7 @@ void BruteForce::detect_edge_edge_candidates( void BruteForce::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( face_boxes, vertex_boxes, std::bind(&BruteForce::can_face_vertex_collide, this, _1, _2), @@ -96,6 +100,7 @@ void BruteForce::detect_face_vertex_candidates( void BruteForce::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_boxes, face_boxes, std::bind(&BruteForce::can_edge_face_collide, this, _1, _2), @@ -105,6 +110,7 @@ void BruteForce::detect_edge_face_candidates( void BruteForce::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( face_boxes, face_boxes, std::bind(&BruteForce::can_faces_collide, this, _1, _2), candidates); diff --git a/src/ipc/broad_phase/create_broad_phase.cpp b/src/ipc/broad_phase/create_broad_phase.cpp index c812a1170..2cc58e622 100644 --- a/src/ipc/broad_phase/create_broad_phase.cpp +++ b/src/ipc/broad_phase/create_broad_phase.cpp @@ -1,6 +1,7 @@ #include "create_broad_phase.hpp" #include +#include #include #include #include @@ -30,6 +31,14 @@ create_broad_phase(const BroadPhaseMethod& broad_phase_method) log_and_throw_error( "Sweep and Tiniest Queue broad phase requires CUDA! Enable it through CMake option IPC_TOOLKIT_WITH_CUDA."); #endif + case BroadPhaseMethod::LBVH_CUDA: +#ifdef IPC_TOOLKIT_WITH_CUDA + return std::make_shared(); +#else + log_and_throw_error( + "CUDA LBVH broad phase requires CUDA! Enable it through CMake option IPC_TOOLKIT_WITH_CUDA."); +#endif + case BroadPhaseMethod::NUM_BROAD_PHASE_METHODS: default: log_and_throw_error("Unknown broad phase type!"); } diff --git a/src/ipc/broad_phase/create_broad_phase.hpp b/src/ipc/broad_phase/create_broad_phase.hpp index ef6b2ee89..e80a92b65 100644 --- a/src/ipc/broad_phase/create_broad_phase.hpp +++ b/src/ipc/broad_phase/create_broad_phase.hpp @@ -5,15 +5,24 @@ namespace ipc { +/// @brief The broad phase methods create_broad_phase() can construct. enum class BroadPhaseMethod : uint8_t { BRUTE_FORCE, HASH_GRID, SPATIAL_HASH, LBVH, SWEEP_AND_PRUNE, - SWEEP_AND_TINIEST_QUEUE + SWEEP_AND_TINIEST_QUEUE, ///< Requires CUDA (IPC_TOOLKIT_WITH_CUDA). + LBVH_CUDA, ///< Requires CUDA (IPC_TOOLKIT_WITH_CUDA). + /// @brief The number of methods; not a method itself. + NUM_BROAD_PHASE_METHODS }; +/// @brief Construct a broad phase of the given method. +/// @param broad_phase_method The method to construct. +/// @return The broad phase. +/// @throws std::runtime_error if the method requires CUDA and the library was +/// built without it, or if the method is not a valid BroadPhaseMethod. std::shared_ptr create_broad_phase(const BroadPhaseMethod& broad_phase_method); diff --git a/src/ipc/broad_phase/cuda/CMakeLists.txt b/src/ipc/broad_phase/cuda/CMakeLists.txt new file mode 100644 index 000000000..2a375b892 --- /dev/null +++ b/src/ipc/broad_phase/cuda/CMakeLists.txt @@ -0,0 +1,7 @@ +set(SOURCES + lbvh.cu + lbvh.hpp + lbvh_impl.cuh +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu new file mode 100644 index 000000000..b5b62f340 --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -0,0 +1,1342 @@ +#include "lbvh.hpp" + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace ipc::cuda { + +namespace { + + // ipc::LBVH::Node is shared through device memory and copied back to the + // host bytewise, so both sides must agree on its layout. ipc::LBVH asserts + // the size in a host-only constructor no device path instantiates; this + // one is evaluated by nvcc's front end, which is what lays the type out for + // the device. + static_assert( + sizeof(ipc::LBVH::Node) == 32, + "ipc::LBVH::Node must be 32 bytes to share it with the host"); + static_assert( + alignof(ipc::LBVH::Node) == 32, + "ipc::LBVH::Node must be 32-byte aligned to share it with the host"); + + /// @brief Per-internal-node scratch used by the bottom-up build. The + /// counter is a plain int rather than the host build's std::atomic: + /// atomicAdd() needs an int*, and supplies the same atomicity. The layout + /// is otherwise identical to the host's. + using DeviceConstructionInfo = ipc::LBVH::ConstructionInfo; + + using Domain = LBVH::Impl::Domain; + static_assert( + std::is_trivially_copyable_v, + "Domain is copied into kernel parameter space and reduced by CUB"); + + struct DomainReduce { + __host__ __device__ Domain + operator()(const Domain& a, const Domain& b) const + { + Domain r; +#pragma unroll + for (int k = 0; k < 3; ++k) { + r.min[k] = fmin(a.min[k], b.min[k]); + r.max[k] = fmax(a.max[k], b.max[k]); + } + return r; + } + }; + + struct MakeDomain { + const double* box_min; + const double* box_max; + __host__ __device__ Domain operator()(const int i) const + { + Domain d; +#pragma unroll + for (int k = 0; k < 3; ++k) { + d.min[k] = box_min[3 * i + k]; + d.max[k] = box_max[3 * i + k]; + } + return d; + } + }; + + // -- Box building ------------------------------------------------------- + + /// @brief One vertex box from the vertex's positions at t0 and t1; pass + /// the same array twice for a static box. + /// + /// Mirrors AABB::from_point(p_t0, p_t1, r), i.e. the union of the two + /// inflated points, bit-for-bit: the per-coordinate inflation is + /// AABB::conservative_{lower,upper}_bound() -- the very function the host + /// calls -- and taking the union before inflating equals inflating before + /// the union because nextafter is monotone: + /// min(nextafter(a - r), nextafter(b - r)) == nextafter(min(a, b) - r). + /// + /// For dim == 2 input, ipc::AABB always stores a 3-wide array whose z + /// component is zero-initialized and never touched by the inflation (only + /// the first `dim` components of the constructor argument are assigned) -- + /// so the z bound is an exact, uninflated 0.0. Replicate that exactly. + /// + /// @param vertices_t0 Positions at t0, column-major (dim * n). + /// @param vertices_t1 Positions at t1, column-major (dim * n). + /// @param n The number of vertices. + /// @param dim The simulation dimension (2 or 3). + /// @param inflation_radius The inflation radius. + /// @param[out] box_min The box min corners (always 3 * n, row-major). + /// @param[out] box_max The box max corners (always 3 * n, row-major). + __global__ void build_vertex_boxes_kernel( + const double* vertices_t0, // not __restrict__: may alias vertices_t1 + const double* vertices_t1, // (the static build passes one array twice) + const int n, + const int dim, + const double inflation_radius, + double* __restrict__ box_min, + double* __restrict__ box_max) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } +#pragma unroll + for (int k = 0; k < 3; ++k) { + if (k < dim) { + const double a = vertices_t0[k * n + i]; + const double b = vertices_t1[k * n + i]; + box_min[3 * i + k] = AABB::conservative_lower_bound( + fmin(a, b), inflation_radius); + box_max[3 * i + k] = AABB::conservative_upper_bound( + fmax(a, b), inflation_radius); + } else { + box_min[3 * i + k] = 0.0; + box_max[3 * i + k] = 0.0; + } + } + } + + /// @brief One edge box as the union of its two vertex boxes, i.e. the + /// AABB(aabb1, aabb2) constructor ipc::build_edge_boxes uses. + __global__ void build_edge_boxes_kernel( + const double* __restrict__ vbox_min, + const double* __restrict__ vbox_max, + const int32_t* __restrict__ edges, // 2 * n, row-major + const int n, + double* __restrict__ box_min, + double* __restrict__ box_max) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } + const int32_t e0 = edges[2 * i + 0]; + const int32_t e1 = edges[2 * i + 1]; +#pragma unroll + for (int k = 0; k < 3; ++k) { + box_min[3 * i + k] = + fmin(vbox_min[3 * e0 + k], vbox_min[3 * e1 + k]); + box_max[3 * i + k] = + fmax(vbox_max[3 * e0 + k], vbox_max[3 * e1 + k]); + } + } + + /// @brief One face box as the union of its three vertex boxes, i.e. the + /// AABB(aabb1, aabb2, aabb3) constructor ipc::build_face_boxes uses. + __global__ void build_face_boxes_kernel( + const double* __restrict__ vbox_min, + const double* __restrict__ vbox_max, + const int32_t* __restrict__ faces, // 3 * n, row-major + const int n, + double* __restrict__ box_min, + double* __restrict__ box_max) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } + const int32_t f0 = faces[3 * i + 0]; + const int32_t f1 = faces[3 * i + 1]; + const int32_t f2 = faces[3 * i + 2]; +#pragma unroll + for (int k = 0; k < 3; ++k) { + box_min[3 * i + k] = fmin( + vbox_min[3 * f0 + k], + fmin(vbox_min[3 * f1 + k], vbox_min[3 * f2 + k])); + box_max[3 * i + k] = fmax( + vbox_max[3 * f0 + k], + fmax(vbox_max[3 * f1 + k], vbox_max[3 * f2 + k])); + } + } + + // -- Tree building ------------------------------------------------------ + + /// @brief Compute one Morton code per box from its (normalized) center. + /// Mirrors the compute_morton_codes block of ipc::LBVH::init_bvh. + /// @param box_min The box min corners (3 * n, row-major). + /// @param box_max The box max corners (3 * n, row-major). + /// @param n The number of boxes. + /// @param domain The Morton-normalization domain (device resident). + /// @param dim The simulation dimension (2 or 3). + /// @param[out] codes The Morton codes. + /// @param[out] box_ids The box ids (the identity, to be sorted with codes). + __global__ void compute_morton_codes_kernel( + const double* __restrict__ box_min, + const double* __restrict__ box_max, + const int n, + const Domain* __restrict__ domain, + const int dim, + uint64_t* __restrict__ codes, + int32_t* __restrict__ box_ids) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } + + const Eigen::Array3d mesh_min( + domain->min[0], domain->min[1], domain->min[2]); + const Eigen::Array3d mesh_max( + domain->max[0], domain->max[1], domain->max[2]); + // The same derivation the CPU (ipc::LBVH::init_bvh) evaluates once per + // build; IEEE division is deterministic, so evaluating it per thread + // gives the identical reciprocal and so identical codes. + const Eigen::Array3d mesh_width_inv = + ipc::morton_domain_width_inv(mesh_min, mesh_max); + + const Eigen::Array3d center( + 0.5 * (box_min[3 * i + 0] + box_max[3 * i + 0]), + 0.5 * (box_min[3 * i + 1] + box_max[3 * i + 1]), + 0.5 * (box_min[3 * i + 2] + box_max[3 * i + 2])); + + codes[i] = ipc::morton_code(center, mesh_min, mesh_width_inv, dim); + box_ids[i] = i; + } + + /// @brief Single-pass bottom-up hierarchy + AABB build (Apetrei 2014). + /// One thread per leaf, driving ipc::details::build_hierarchy_from_leaf() + /// -- the same walk the CPU build runs -- with atomicAdd() and the fences + /// around it standing in for the host's std::atomic arrival. + /// @param box_min The box min corners (3 * n, row-major). + /// @param box_max The box max corners (3 * n, row-major). + /// @param sorted_codes The Morton codes in sorted order. + /// @param sorted_box_ids The box ids in Morton-sorted order. + /// @param n_leaves The number of leaves. + /// @param[out] nodes The BVH nodes. + /// @param[out] rightmost The per-node rightmost-leaf indices. + /// @param[in,out] infos The per-node construction scratch (zeroed). + /// @param[out] root_idx The root's index. + __global__ void build_hierarchy_kernel( + const double* __restrict__ box_min, + const double* __restrict__ box_max, + const uint64_t* __restrict__ sorted_codes, + const int32_t* __restrict__ sorted_box_ids, + const int n_leaves, + ipc::LBVH::Node* __restrict__ nodes, + int32_t* __restrict__ rightmost, + DeviceConstructionInfo* __restrict__ infos, + int* __restrict__ root_idx) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n_leaves) { + return; + } + + const int32_t bid = sorted_box_ids[i]; + ipc::details::init_leaf_node( + i, n_leaves, bid, + Eigen::Array3d( + box_min[3 * bid + 0], box_min[3 * bid + 1], + box_min[3 * bid + 2]), + Eigen::Array3d( + box_max[3 * bid + 0], box_max[3 * bid + 1], + box_max[3 * bid + 2]), + nodes, rightmost); + + const int root = ipc::details::build_hierarchy_from_leaf( + i, n_leaves, [sorted_codes](int k) { return sorted_codes[k]; }, + nodes, rightmost, infos, + [](int& count) { + // Release: publish this thread's child pointer, range endpoint + // and leaf/subtree AABB before announcing arrival, so whoever + // continues sees a complete child. + __threadfence(); + const int previous = atomicAdd(&count, 1); + if (previous != 0) { + // Acquire: this thread continues and reads the sibling's + // node, range endpoint and rightmost leaf. Those are + // ordinary loads, so without this fence they may be served + // from a stale L1 on another SM. + __threadfence(); + } + return previous; + }); + + if (root >= 0) { + *root_idx = root; // only one thread reaches the root + } + } + + /// @brief Swap the node and rightmost-leaf entries at index 0 and the root + /// (single thread). See ipc::details::swap_root_to_zero(). Reads the root + /// from device memory so the build never round-trips through the host. + /// @param nodes The BVH nodes. + /// @param rightmost The per-node rightmost-leaf indices. + /// @param root The root's index; a no-op if it is already 0 (or unset). + __global__ void swap_root_kernel( + ipc::LBVH::Node* __restrict__ nodes, + int32_t* __restrict__ rightmost, + const int* __restrict__ root) + { + if (blockIdx.x == 0 && threadIdx.x == 0 && *root > 0) { + ipc::details::swap_root_to_zero(nodes, rightmost, *root); + } + } + + /// @brief After the root swap, rewrite left pointers that referenced the + /// old node 0 to its new location. See ipc::details::patch_left_pointer(). + /// @param nodes The BVH nodes. + /// @param num_nodes The number of nodes. + /// @param root The new location of the old node 0; a no-op if the root was + /// already 0 (or unset). + __global__ void patch_left_kernel( + ipc::LBVH::Node* __restrict__ nodes, + const int num_nodes, + const int* __restrict__ root) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + const int r = *root; + if (i >= num_nodes || r <= 0) { + return; + } + ipc::details::patch_left_pointer(nodes[i], r); + } + + /// @brief Build one BVH on the device from device-resident box corners. + /// Mirrors ipc::LBVH::init_bvh. Fully asynchronous: nothing here waits on + /// the device; the caller synchronizes once after all three trees. + /// @param impl The pimpl whose scratch buffers to use. + /// @param d_box_min The box min corners (3 * n, row-major, device). + /// @param d_box_max The box max corners (3 * n, row-major, device). + /// @param n The number of boxes (leaves). + /// @param dim The simulation dimension (2 or 3). + /// @param[out] bvh The BVH to build. + /// @param[out] d_root Where to write the root's index (device); -1 if the + /// build never reaches a root. + void build_tree( + LBVH::Impl& impl, + const double* d_box_min, + const double* d_box_max, + const int n, + const int dim, + LBVH::Impl::DeviceBVH& bvh, + int* d_root) + { + if (n == 0) { + bvh.clear(); + return; + } + + const size_t num_nodes = size_t(2) * n - 1; + bvh.nodes.resize(num_nodes); + bvh.rightmost_leaves.resize(num_nodes); + + for (auto& codes : impl.morton_codes) { + codes.resize(n); + } + for (auto& ids : impl.box_ids) { + ids.resize(n); + } + // Only the visitation counts need zeroing; a memset is the cheapest + // way to do it. + impl.construction_infos.resize(num_nodes); + impl.construction_infos.zero(); + IPC_TOOLKIT_CUDA_CHECK(cudaMemsetAsync(d_root, 0xFF, sizeof(int))); + + compute_morton_codes_kernel<<>>( + d_box_min, d_box_max, n, impl.domain.data(), dim, + impl.morton_codes[0].data(), impl.box_ids[0].data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + // Radix sort the (code, id) pairs by code. + // + // A radix sort is out of place: each pass reads one array and writes + // another, so CUB is given two buffers per sequence (a DoubleBuffer) + // and ping-pongs between them, pass after pass. The unsorted input is + // in buffer [0], where the codes kernel wrote it; after the sort, + // Current() says which of the two holds the result -- it depends on + // the number of passes, so it must be read back rather than assumed. + cub::DoubleBuffer keys( + impl.morton_codes[0].data(), impl.morton_codes[1].data()); + cub::DoubleBuffer values( + impl.box_ids[0].data(), impl.box_ids[1].data()); + + // CUB's two-phase convention: with a null storage pointer the call + // only reports the scratch bytes it needs; the second call sorts. + // Keeping that scratch in a persistent buffer is what makes this + // allocation-free per build (thrust::sort_by_key would cudaMalloc and + // cudaFree it every call). + size_t temp_bytes = 0; + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceRadixSort::SortPairs( + nullptr, temp_bytes, keys, values, n)); + impl.sort_temp.resize(temp_bytes); + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceRadixSort::SortPairs( + impl.sort_temp.data(), temp_bytes, keys, values, n)); + + // Current() is the sorted half of each ping-pong pair. + build_hierarchy_kernel<<>>( + d_box_min, d_box_max, keys.Current(), values.Current(), n, + bvh.nodes.data(), bvh.rightmost_leaves.data(), + impl.construction_infos.data(), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + swap_root_kernel<<<1, 1>>>( + bvh.nodes.data(), bvh.rightmost_leaves.data(), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + patch_left_kernel<<>>( + bvh.nodes.data(), static_cast(num_nodes), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + + /// @brief Reduce the Morton-normalization domain (min of mins, max of + /// maxs) over the device-resident vertex box corners into impl.domain. + /// The same quantity BroadPhase::compute_mesh_aabb() computes on the host + /// from the same boxes, with the same seeds; min/max are order-independent, + /// so the two agree bit-for-bit. Stays on the device: the codes kernel + /// reads it directly, with no host round-trip. + void compute_domain(LBVH::Impl& impl, const int n_vertices) + { + Domain init; + for (int k = 0; k < 3; ++k) { + init.min[k] = std::numeric_limits::max(); + init.max[k] = std::numeric_limits::lowest(); + } + + const auto domains = thrust::make_transform_iterator( + thrust::counting_iterator(0), + MakeDomain { impl.vbox_min.data(), impl.vbox_max.data() }); + + impl.domain.resize(1); + size_t temp_bytes = 0; + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceReduce::Reduce( + nullptr, temp_bytes, domains, impl.domain.data(), n_vertices, + DomainReduce {}, init)); + impl.reduce_temp.resize(temp_bytes); + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceReduce::Reduce( + impl.reduce_temp.data(), temp_bytes, domains, + impl.domain.data(), n_vertices, DomainReduce {}, init)); + } + + /// @brief Upload vertex positions column-major, straight from the matrix + /// when it is contiguous (no host transpose). + void upload_vertices( + Eigen::ConstRef vertices, DeviceBuffer& d) + { + const size_t n = static_cast(vertices.size()); + if (vertices.innerStride() == 1 + && vertices.outerStride() == vertices.rows()) { + d.upload(vertices.data(), n); + } else { + const Eigen::MatrixXd contiguous = vertices; + d.upload(contiguous.data(), n); + } + } + + /// @brief Flatten an integer connectivity matrix (rowwise) to row-major + /// 32-bit ids on the host, and upload the same array to the device. + template + void upload_connectivity( + Eigen::ConstRef M, + std::vector& h, + DeviceBuffer& d) + { + const size_t n = M.rows(); + h.resize(Cols * n); + for (size_t i = 0; i < n; ++i) { + for (int k = 0; k < Cols; ++k) { + h[Cols * i + k] = static_cast(M(i, k)); + } + } + d.upload(h.data(), h.size()); + } + + /// @brief Copy a device BVH to the host. Bytewise: ipc::LBVH::Node holds + /// only floats and ints, so a memcpy is its copy. + void to_host( + const LBVH::Impl::DeviceBVH& bvh, + ipc::LBVH::Nodes& nodes, + ipc::LBVH::RightmostLeaves& rightmost_leaves) + { + nodes.resize(bvh.nodes.size()); + rightmost_leaves.resize(bvh.rightmost_leaves.size()); + bvh.nodes.download(nodes.data()); + bvh.rightmost_leaves.download(rightmost_leaves.data()); + } + + /// @brief Given device-resident vertex boxes (impl.vbox_min/max), build the + /// edge/face boxes and all three BVHs. Shared by every build() overload. + /// @param impl The pimpl to fill (output). + /// @param dim The simulation dimension (2 or 3). + /// @param n_vertices The number of vertices. + /// @param edges The mesh edges. + /// @param faces The mesh faces. + void build_from_vertex_boxes( + LBVH::Impl& impl, + const int dim, + const int n_vertices, + Eigen::ConstRef edges, + Eigen::ConstRef faces) + { + assert(edges.size() == 0 || edges.cols() == 2); + assert(faces.size() == 0 || faces.cols() == 3); + + const int n_edges = static_cast(edges.rows()); + const int n_faces = static_cast(faces.rows()); + + upload_connectivity<2>(edges, impl.h_edges, impl.edges); + upload_connectivity<3>(faces, impl.h_faces, impl.faces); + + // Build edge/face boxes on the device from the vertex boxes. + impl.ebox_min.resize(3 * size_t(n_edges)); + impl.ebox_max.resize(3 * size_t(n_edges)); + if (n_edges > 0) { + build_edge_boxes_kernel<<< + kernel_grid_size(n_edges), KERNEL_BLOCK_SIZE>>>( + impl.vbox_min.data(), impl.vbox_max.data(), impl.edges.data(), + n_edges, impl.ebox_min.data(), impl.ebox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + + impl.fbox_min.resize(3 * size_t(n_faces)); + impl.fbox_max.resize(3 * size_t(n_faces)); + if (n_faces > 0) { + build_face_boxes_kernel<<< + kernel_grid_size(n_faces), KERNEL_BLOCK_SIZE>>>( + impl.vbox_min.data(), impl.vbox_max.data(), impl.faces.data(), + n_faces, impl.fbox_min.data(), impl.fbox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + + // The CPU normalizes all three BVHs by the vertex box domain. + compute_domain(impl, n_vertices); + + impl.roots.resize(3); + build_tree( + impl, impl.vbox_min.data(), impl.vbox_max.data(), n_vertices, dim, + impl.vertex_bvh, impl.roots.data() + 0); + build_tree( + impl, impl.ebox_min.data(), impl.ebox_max.data(), n_edges, dim, + impl.edge_bvh, impl.roots.data() + 1); + build_tree( + impl, impl.fbox_min.data(), impl.fbox_max.data(), n_faces, dim, + impl.face_bvh, impl.roots.data() + 2); + + // The one synchronization of the build: surfaces any kernel fault and + // lets the roots be read back -- all three at once. + IPC_TOOLKIT_CUDA_CHECK(cudaDeviceSynchronize()); + + // A hierarchy build that never reaches its root leaves -1 behind. The + // traversal always starts at node 0, which would then be an arbitrary + // interior node, silently dropping every candidate outside its + // subtree. Refuse to hand out such a tree. + int roots[3]; + impl.roots.download(roots); + const int n_leaves[3] = { n_vertices, n_edges, n_faces }; + for (int k = 0; k < 3; ++k) { + if (n_leaves[k] > 0 && roots[k] < 0) { + log_and_throw_error( + "ipc::cuda::LBVH: the device hierarchy build did not " + "reach a root (tree {} of 3, {} leaves); the BVH is " + "malformed", + k, n_leaves[k]); + } + } + } + + // -- Traversal ---------------------------------------------------------- + + /// @brief Load a primitive's vertex ids into registers: a vertex's id is + /// itself; an edge's or a face's come from its connectivity row. Fully + /// unrolled, so the array is only ever indexed by constants and stays in + /// registers (a runtime-indexed local array would go to local memory). + /// @tparam N The number of vertex ids per primitive (1, 2, or 3). + /// @param prim The primitive id. + /// @param conn The connectivity (N ids per primitive, row-major); unused + /// and may be null for N == 1. + /// @param[out] ids The primitive's vertex ids. + template + __device__ inline void load_vertex_ids( + const int32_t prim, const int32_t* __restrict__ conn, int32_t (&ids)[N]) + { + if constexpr (N == 1) { + ids[0] = prim; + } else { + assert(conn != nullptr); +#pragma unroll + for (int k = 0; k < N; ++k) { + ids[k] = conn[N * prim + k]; + } + } + } + + /// @brief Append a (source_prim, target_prim) pair (post-swap) via an + /// atomic counter. Writes only if the slot is within capacity; the counter + /// still advances on overflow so the caller learns the required size. The + /// counter and capacity are 64-bit so neither can wrap or truncate. + template + __device__ inline void emit_pair( + const int32_t query_prim, + const int32_t node_prim, + int32_t* __restrict__ out_a, + int32_t* __restrict__ out_b, + unsigned long long* __restrict__ counter, + const unsigned long long capacity) + { + int32_t a = query_prim, b = node_prim; + if constexpr (swap_order) { + const int32_t t = a; + a = b; + b = t; + } + const unsigned long long slot = atomicAdd(counter, 1ULL); + if (slot < capacity) { + out_a[slot] = a; + out_b[slot] = b; + } + } + + /// @brief One thread per source leaf: descend the target BVH and append + /// every AABB-overlapping, connectivity-passing (source_prim, target_prim) + /// pair to the output arrays. The descent is ipc::details::traverse_lbvh() + /// and the shared-vertex exclusion ipc::details::share_vertex(), both + /// shared with the CPU ipc::LBVH. The remaining user vertex filter (if + /// any) is applied on the host, so the final set matches the CPU. + /// @tparam triangular Self-collision: skip subtrees left of the query. + /// @tparam swap_order Emit (target_prim, source_prim) instead. + /// @tparam SourceCount The vertex ids per source primitive (1, 2, or 3). + /// @tparam TargetCount The vertex ids per target primitive (1, 2, or 3). + /// @param source The BVH whose leaves are the queries. + /// @param n_source_leaves The number of source leaves. + /// @param source_leaf_offset The index of the source BVH's first leaf. + /// @param target The BVH to descend. + /// @param target_size The number of nodes in the target BVH. + /// @param target_rightmost The target's per-node rightmost-leaf indices. + /// @param source_conn The source connectivity (null for vertices). + /// @param target_conn The target connectivity (null for vertices). + /// @param[out] out_a The first ids of the emitted pairs. + /// @param[out] out_b The second ids of the emitted pairs. + /// @param[in,out] counter The emitted-pair counter. + /// @param capacity The output arrays' capacity. + template < + bool triangular, + bool swap_order, + int SourceCount, + int TargetCount> + __global__ void traverse_kernel( + const ipc::LBVH::Node* __restrict__ source, + const int n_source_leaves, + const int source_leaf_offset, + const ipc::LBVH::Node* __restrict__ target, + const int target_size, + const int32_t* __restrict__ target_rightmost, + const int32_t* __restrict__ source_conn, + const int32_t* __restrict__ target_conn, + int32_t* __restrict__ out_a, + int32_t* __restrict__ out_b, + unsigned long long* __restrict__ counter, + const unsigned long long capacity) + { + const int s = blockIdx.x * blockDim.x + threadIdx.x; + if (s >= n_source_leaves) { + return; + } + const ipc::LBVH::Node query = source[source_leaf_offset + s]; + + int32_t query_ids[SourceCount]; + load_vertex_ids( + query.primitive_id, source_conn, query_ids); + + ipc::details::traverse_lbvh( + s, target, target_size, target_rightmost, + [&](const ipc::LBVH::Node& node) { return node.intersects(query); }, + [&](const ipc::LBVH::Node& leaf, const int /*leaf_idx*/, + const bool /*intersects*/) { + int32_t leaf_ids[TargetCount]; + load_vertex_ids( + leaf.primitive_id, target_conn, leaf_ids); + if (!ipc::details::share_vertex(query_ids, leaf_ids)) { + emit_pair( + query.primitive_id, leaf.primitive_id, out_a, out_b, + counter, capacity); + } + }); + } + + /// @brief How one candidate type is traversed: which BVH is the source + /// (its leaves are the queries), which is the target (descended), how many + /// vertex ids each primitive has, whether the pair is triangular (a BVH + /// against itself) and whether it is emitted swapped. Each tuple is stated + /// exactly once here and drives both the host-materializing and the + /// device-view detect paths, so the two cannot disagree. + template struct Traversal; + + template <> struct Traversal { + static constexpr bool triangular = true; + static constexpr bool swap_order = false; + static constexpr int source_count = 1; + static constexpr int target_count = 1; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.vertex_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.vertex_bvh; + } + static const int32_t* source_conn(const LBVH::Impl&) { return nullptr; } + static const int32_t* target_conn(const LBVH::Impl&) { return nullptr; } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.vv_candidates; + } + }; + + // In 2D and for codimensional edge-vertex collisions there are more + // vertices than edges, so iterate over the edges. Mirrors ipc::LBVH. + template <> struct Traversal { + static constexpr bool triangular = false; + static constexpr bool swap_order = false; + static constexpr int source_count = 2; + static constexpr int target_count = 1; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.edge_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.vertex_bvh; + } + static const int32_t* source_conn(const LBVH::Impl& impl) + { + return impl.edges.data(); + } + static const int32_t* target_conn(const LBVH::Impl&) { return nullptr; } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.ev_candidates; + } + }; + + template <> struct Traversal { + static constexpr bool triangular = true; + static constexpr bool swap_order = false; + static constexpr int source_count = 2; + static constexpr int target_count = 2; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.edge_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.edge_bvh; + } + static const int32_t* source_conn(const LBVH::Impl& impl) + { + return impl.edges.data(); + } + static const int32_t* target_conn(const LBVH::Impl& impl) + { + return impl.edges.data(); + } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.ee_candidates; + } + }; + + // The ratio vertices:faces is 1:2, so iterate over the vertices and query + // the face BVH, swapping so the emitted pair is (face, vertex). Mirrors + // ipc::LBVH. + template <> struct Traversal { + static constexpr bool triangular = false; + static constexpr bool swap_order = true; + static constexpr int source_count = 1; + static constexpr int target_count = 3; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.vertex_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.face_bvh; + } + static const int32_t* source_conn(const LBVH::Impl&) { return nullptr; } + static const int32_t* target_conn(const LBVH::Impl& impl) + { + return impl.faces.data(); + } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.fv_candidates; + } + }; + + // The ratio edges:faces is 3:2, so iterate over the faces and query the + // edge BVH, swapping so the emitted pair is (edge, face). Mirrors + // ipc::LBVH. + template <> struct Traversal { + static constexpr bool triangular = false; + static constexpr bool swap_order = true; + static constexpr int source_count = 3; + static constexpr int target_count = 2; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.face_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.edge_bvh; + } + static const int32_t* source_conn(const LBVH::Impl& impl) + { + return impl.faces.data(); + } + static const int32_t* target_conn(const LBVH::Impl& impl) + { + return impl.edges.data(); + } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.ef_candidates; + } + }; + + template <> struct Traversal { + static constexpr bool triangular = true; + static constexpr bool swap_order = false; + static constexpr int source_count = 3; + static constexpr int target_count = 3; + static const LBVH::Impl::DeviceBVH& source(const LBVH::Impl& impl) + { + return impl.face_bvh; + } + static const LBVH::Impl::DeviceBVH& target(const LBVH::Impl& impl) + { + return impl.face_bvh; + } + static const int32_t* source_conn(const LBVH::Impl& impl) + { + return impl.faces.data(); + } + static const int32_t* target_conn(const LBVH::Impl& impl) + { + return impl.faces.data(); + } + static LBVH::Impl::DeviceCandidates& buffer(LBVH::Impl& impl) + { + return impl.ff_candidates; + } + }; + + /// @brief Run the device traversal for one candidate type, leaving the + /// connectivity-filtered pairs device-resident in its buffer. + /// + /// The first pass is sized from the buffer's high-water mark (the largest + /// count any earlier call on this object needed), or a guess of 8 pairs + /// per source leaf, whichever is larger. The kernel counts every pair it + /// finds but writes only those that fit, so if the first pass overflows + /// the exact count is known and a second pass always fits: at most two + /// passes, and only the first call -- or a call whose count exceeds every + /// prior one -- pays for the second. + /// + /// @tparam Candidate The candidate type; see Traversal. + /// @param impl The pimpl; the caller must hold impl.mutex. + /// @return The number of candidate pairs emitted. + template size_t run_traversal(LBVH::Impl& impl) + { + using T = Traversal; + const LBVH::Impl::DeviceBVH& source = T::source(impl); + const LBVH::Impl::DeviceBVH& target = T::target(impl); + LBVH::Impl::DeviceCandidates& buf = T::buffer(impl); + + buf.clear(); + + const int n_source_leaves = source.n_leaves(); + const int target_size = static_cast(target.nodes.size()); + // A triangular traversal is a BVH against itself, and a lone primitive + // cannot collide with itself. + if (n_source_leaves == 0 || target_size == 0 + || (T::triangular && n_source_leaves < 2)) { + return 0; + } + const int source_leaf_offset = n_source_leaves - 1; + + size_t capacity = std::max( + { buf.a.capacity(), size_t(1024), + size_t(8) * static_cast(n_source_leaves) }); + impl.counter.resize(1); + + unsigned long long count = 0; + for (int pass = 0; pass < 2; ++pass) { + buf.a.reserve(capacity); + buf.b.reserve(capacity); + impl.counter.zero(); + + traverse_kernel< + T::triangular, T::swap_order, T::source_count, T::target_count> + <<>>( + source.nodes.data(), n_source_leaves, source_leaf_offset, + target.nodes.data(), target_size, + target.rightmost_leaves.data(), T::source_conn(impl), + T::target_conn(impl), buf.a.data(), buf.b.data(), + impl.counter.data(), + static_cast(capacity)); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + impl.counter.download(&count); // synchronizes + if (count <= capacity) { + break; // everything fit + } + + logger().warn( + "[ipc::cuda::LBVH] candidate count {} exceeded preallocated " + "capacity {}; re-running with the exact size (this cost is " + "amortized: later calls reuse the learned capacity)", + count, capacity); + capacity = static_cast(count); // exact size now known + } + if (count > capacity) { + // The count is a deterministic function of the two trees, so the + // second pass must fit; anything else is a device fault. + log_and_throw_error( + "ipc::cuda::LBVH: the candidate count changed between two " + "identical traversals ({} > capacity {})", + count, capacity); + } + + buf.count = static_cast(count); + buf.a.resize(buf.count); + buf.b.resize(buf.count); + return buf.count; + } + + /// @brief Copy the device-resident candidate pairs to host Candidate + /// objects. For the accept-all filter every pair is kept (the device set + /// is already exact); otherwise the user vertex filter trims the + /// connectivity-filtered superset. + /// @param buf The device pairs. + /// @param filter The user vertex filter. + /// @param can_collide The full predicate applied when the filter is not accept-all. + /// @param[out] out The materialized candidates (cleared first). + template + void materialize( + const LBVH::Impl::DeviceCandidates& buf, + const CollisionFilter& filter, + const CanCollide& can_collide, + std::vector& out) + { + out.clear(); + const size_t count = buf.count; + if (count == 0) { + return; + } + std::vector h_a(count), h_b(count); + buf.a.download(h_a.data()); + buf.b.download(h_b.data()); + + out.reserve(count); + if (filter.accepts_all()) { + for (size_t k = 0; k < count; ++k) { + out.emplace_back(h_a[k], h_b[k]); + } + } else { + for (size_t k = 0; k < count; ++k) { + if (can_collide(h_a[k], h_b[k])) { + out.emplace_back(h_a[k], h_b[k]); + } + } + } + } + + /// @brief Traverse on the device, then materialize on the host. + template + void detect_host( + LBVH::Impl& impl, + const CollisionFilter& filter, + const CanCollide& can_collide, + std::vector& out) + { + // The detect_*() methods are const on the BroadPhase interface, and + // ipc::LBVH's really are read-only, so a caller may legitimately run + // two of them concurrently on one shared object. Here every one of + // them writes shared device state -- the pair counter, and the + // candidate buffer of its type -- so without this lock two detections + // would race on the counter and on buffer reallocation. Held for the + // materialize too, so another call cannot overwrite the buffer while + // it is being copied to the host. + const std::lock_guard lock(impl.mutex); + run_traversal(impl); + materialize( + Traversal::buffer(impl), filter, can_collide, out); + } + + /// @brief Traverse on the device and return a view of the result. + template + LBVH::DeviceCandidateView detect_device(LBVH::Impl& impl) + { + // Same reason as detect_host(): const on the interface, but writes the + // shared counter and this type's candidate buffer. The lock ends with + // the call, so the returned view is only as safe as the caller's own + // ordering of later detect_*() calls (see DeviceCandidateView). + const std::lock_guard lock(impl.mutex); + const size_t count = run_traversal(impl); + const LBVH::Impl::DeviceCandidates& buf = + Traversal::buffer(impl); + return LBVH::DeviceCandidateView { count ? buf.a.data() : nullptr, + count ? buf.b.data() : nullptr, + count }; + } + +} // namespace + +// --------------------------------------------------------------------------- + +LBVH::LBVH() : ipc::BroadPhase(), m_impl(std::make_unique()) { } + +LBVH::~LBVH() = default; + +// ipc::BroadPhase declares a destructor and so has no move operations: moving +// it as a whole would invoke its copy, which copies a std::function and may +// throw. Its members are moved individually instead, which cannot. The +// moved-from object is left in the cleared state (no Impl, dim 0, default +// filter); impl() re-seeds it on its next use. + +LBVH::LBVH(LBVH&& other) noexcept + : ipc::BroadPhase() + , m_impl(std::move(other.m_impl)) +{ + can_vertices_collide = std::move(other.can_vertices_collide); + other.can_vertices_collide = CollisionFilter(); + dim = other.dim; + other.dim = 0; +} + +LBVH& LBVH::operator=(LBVH&& other) noexcept +{ + if (this != &other) { + m_impl = std::move(other.m_impl); + can_vertices_collide = std::move(other.can_vertices_collide); + other.can_vertices_collide = CollisionFilter(); + dim = other.dim; + other.dim = 0; + } + return *this; +} + +LBVH::Impl& LBVH::impl() const +{ + if (!m_impl) { + m_impl = std::make_unique(); + } + return *m_impl; +} + +void LBVH::build( + Eigen::ConstRef vertices, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius) +{ + // A static box is a temporal one whose endpoints coincide. The dynamic + // build notices the two refs alias and uploads the vertices only once. + build(vertices, vertices, edges, faces, inflation_radius); +} + +void LBVH::build( + Eigen::ConstRef vertices_t0, + Eigen::ConstRef vertices_t1, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius) +{ + assert(vertices_t0.rows() == vertices_t1.rows()); + assert(vertices_t0.cols() == vertices_t1.cols()); + assert(vertices_t0.rows() <= std::numeric_limits::max()); + + clear(); + + assert(vertices_t0.cols() == 2 || vertices_t0.cols() == 3); + dim = static_cast(vertices_t0.cols()); + + const int n_vertices = static_cast(vertices_t0.rows()); + if (n_vertices == 0) { + return; + } + + Impl& device = impl(); + upload_vertices(vertices_t0, device.vertices_t0); + + // The static build passes the same matrix twice; upload it once and point + // the kernel's t1 at the t0 copy rather than paying a second transfer. + const bool same_vertices = vertices_t0.data() == vertices_t1.data() + && vertices_t0.outerStride() == vertices_t1.outerStride(); + if (!same_vertices) { + upload_vertices(vertices_t1, device.vertices_t1); + } + const double* d_vertices_t1 = + same_vertices ? device.vertices_t0.data() : device.vertices_t1.data(); + + // Build vertex boxes on the device (always 3-wide storage). + device.vbox_min.resize(3 * size_t(n_vertices)); + device.vbox_max.resize(3 * size_t(n_vertices)); + build_vertex_boxes_kernel<<< + kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( + device.vertices_t0.data(), d_vertices_t1, n_vertices, dim, + inflation_radius, device.vbox_min.data(), device.vbox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + build_from_vertex_boxes(device, dim, n_vertices, edges, faces); +} + +void LBVH::build( + Eigen::ConstRef edges, + Eigen::ConstRef faces) +{ + // BroadPhase::build(const AABBs&, edges, faces, dim) has cleared us and + // filled vertex_boxes and dim. + assert(dim == 2 || dim == 3); + assert(vertex_boxes.size() <= std::numeric_limits::max()); + + const int n_vertices = static_cast(vertex_boxes.size()); + if (n_vertices == 0) { + return; + } + + // Upload the precomputed vertex boxes. + std::vector h_min(3 * size_t(n_vertices)); + std::vector h_max(3 * size_t(n_vertices)); + for (int i = 0; i < n_vertices; ++i) { + for (int k = 0; k < 3; ++k) { + h_min[3 * size_t(i) + k] = vertex_boxes[i].min[k]; + h_max[3 * size_t(i) + k] = vertex_boxes[i].max[k]; + } + } + Impl& device = impl(); + device.vbox_min.upload(h_min.data(), h_min.size()); + device.vbox_max.upload(h_max.data(), h_max.size()); + + build_from_vertex_boxes(device, dim, n_vertices, edges, faces); + + // As in ipc::LBVH: the host boxes are redundant once the trees exist. + vertex_boxes.clear(); +} + +void LBVH::clear() +{ + ipc::BroadPhase::clear(); + if (m_impl) { // a moved-from object has nothing to clear + m_impl->clear(); + } +} + +// --------------------------------------------------------------------------- +// BroadPhase interface. Device BVH descent + device connectivity filter; the +// user vertex filter is applied on the host only when it is not accept-all. + +void LBVH::detect_vertex_vertex_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_vertices_collide(a, b); }, + candidates); +} + +void LBVH::detect_edge_vertex_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_edge_vertex_collide(a, b); }, + candidates); +} + +void LBVH::detect_edge_edge_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_edges_collide(a, b); }, + candidates); +} + +void LBVH::detect_face_vertex_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_face_vertex_collide(a, b); }, + candidates); +} + +void LBVH::detect_edge_face_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_edge_face_collide(a, b); }, + candidates); +} + +void LBVH::detect_face_face_candidates( + std::vector& candidates) const +{ + detect_host( + impl(), can_vertices_collide, + [this](size_t a, size_t b) { return can_faces_collide(a, b); }, + candidates); +} + +// --------------------------------------------------------------------------- +// Device-resident candidate accessors. + +LBVH::DeviceCandidateView LBVH::detect_vertex_vertex_candidates_device() const +{ + return detect_device(impl()); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_vertex_candidates_device() const +{ + return detect_device(impl()); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_edge_candidates_device() const +{ + return detect_device(impl()); +} + +LBVH::DeviceCandidateView LBVH::detect_face_vertex_candidates_device() const +{ + return detect_device(impl()); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_face_candidates_device() const +{ + return detect_device(impl()); +} + +LBVH::DeviceCandidateView LBVH::detect_face_face_candidates_device() const +{ + return detect_device(impl()); +} + +// --------------------------------------------------------------------------- +// Host-side can_*_collide filters (mesh connectivity + user vertex filter). +// Mirror ipc::LBVH's overrides, backed by the host connectivity copies. The +// ids arrive from device memory, so the bounds asserts are the tripwire that +// localizes a host/device desync to the traversal. + +bool LBVH::can_edge_vertex_collide(size_t ei, size_t vi) const +{ + const std::vector& edges = impl().h_edges; + assert(2 * ei + 1 < edges.size()); + + return ipc::details::can_edge_vertex_collide( + edges[2 * ei], edges[2 * ei + 1], vi, can_vertices_collide); +} + +bool LBVH::can_edges_collide(size_t eai, size_t ebi) const +{ + const std::vector& edges = impl().h_edges; + assert(2 * eai + 1 < edges.size()); + assert(2 * ebi + 1 < edges.size()); + + return ipc::details::can_edges_collide( + edges[2 * eai], edges[2 * eai + 1], edges[2 * ebi], edges[2 * ebi + 1], + can_vertices_collide); +} + +bool LBVH::can_face_vertex_collide(size_t fi, size_t vi) const +{ + const std::vector& faces = impl().h_faces; + assert(3 * fi + 2 < faces.size()); + + return ipc::details::can_face_vertex_collide( + faces[3 * fi], faces[3 * fi + 1], faces[3 * fi + 2], vi, + can_vertices_collide); +} + +bool LBVH::can_edge_face_collide(size_t ei, size_t fi) const +{ + const std::vector& edges = impl().h_edges; + const std::vector& faces = impl().h_faces; + assert(2 * ei + 1 < edges.size()); + assert(3 * fi + 2 < faces.size()); + + return ipc::details::can_edge_face_collide( + edges[2 * ei], edges[2 * ei + 1], faces[3 * fi], faces[3 * fi + 1], + faces[3 * fi + 2], can_vertices_collide); +} + +bool LBVH::can_faces_collide(size_t fai, size_t fbi) const +{ + const std::vector& faces = impl().h_faces; + assert(3 * fai + 2 < faces.size()); + assert(3 * fbi + 2 < faces.size()); + + return ipc::details::can_faces_collide( + faces[3 * fai], faces[3 * fai + 1], faces[3 * fai + 2], faces[3 * fbi], + faces[3 * fbi + 1], faces[3 * fbi + 2], can_vertices_collide); +} + +size_t LBVH::num_vertex_nodes() const { return impl().vertex_bvh.nodes.size(); } + +size_t LBVH::num_edge_nodes() const { return impl().edge_bvh.nodes.size(); } + +size_t LBVH::num_face_nodes() const { return impl().face_bvh.nodes.size(); } + +// --------------------------------------------------------------------------- +// Debug / validation. + +void LBVH::vertex_nodes_to_host( + ipc::LBVH::Nodes& nodes, ipc::LBVH::RightmostLeaves& rightmost_leaves) const +{ + to_host(impl().vertex_bvh, nodes, rightmost_leaves); +} + +void LBVH::edge_nodes_to_host( + ipc::LBVH::Nodes& nodes, ipc::LBVH::RightmostLeaves& rightmost_leaves) const +{ + to_host(impl().edge_bvh, nodes, rightmost_leaves); +} + +void LBVH::face_nodes_to_host( + ipc::LBVH::Nodes& nodes, ipc::LBVH::RightmostLeaves& rightmost_leaves) const +{ + to_host(impl().face_bvh, nodes, rightmost_leaves); +} + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/broad_phase/cuda/lbvh.hpp b/src/ipc/broad_phase/cuda/lbvh.hpp new file mode 100644 index 000000000..af974d721 --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh.hpp @@ -0,0 +1,219 @@ +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include // ipc::LBVH::Node / Nodes / RightmostLeaves + +#include + +namespace ipc::cuda { + +/// @brief GPU Linear Bounding Volume Hierarchy (LBVH) broad phase. +/// +/// A first-class GPU counterpart to ipc::LBVH: it builds the vertex/edge/face +/// AABBs and their BVHs, and runs the traversal and mesh-connectivity +/// filtering, entirely on the device. Construction uses Morton codes + the +/// Apetrei [2014] single-pass bottom-up build and reuses the 32-byte +/// ipc::LBVH::Node layout, so the device tree can be copied back to the host +/// and validated against -- or traversed by -- the CPU code. The build, the +/// scalar descent, and the shared-vertex exclusion are the same +/// ipc::details functions ipc::LBVH runs; only the parallel launches, the +/// sort, and the atomics are CUDA's. +/// +/// Detection runs the BVH descent (AABB overlap + triangular dedup) and the +/// connectivity (shared-vertex) exclusion on the device. The user vertex filter +/// (can_vertices_collide) is honored on the device only when it is the default +/// accept-all filter; a non-trivial filter is applied on the host while +/// materializing the device-emitted (connectivity-filtered) candidates. Either +/// way the output matches the CPU ipc::LBVH exactly for any filter. +/// +/// All device-side primitive and vertex ids are 32-bit, independent of +/// ipc::index_t: ipc::LBVH::Node stores its primitive id as an int32_t, which +/// already caps every LBVH at 2^31 - 1 primitives, and the connectivity and +/// candidate buffers use the same width so a downstream kernel reads one +/// consistent id type (see DeviceCandidateView). +/// +/// The detect_*() methods are const, as the BroadPhase interface requires, but +/// share device buffers; concurrent calls on one object are serialized by an +/// internal mutex. Device memory is retained across clear() and build() (the +/// candidate buffers keep their high-water-mark capacity so a per-frame call +/// allocates nothing) and released by the destructor. +class LBVH : public ipc::BroadPhase { +public: + LBVH(); + ~LBVH(); + + /// @brief Move; the moved-from object is left cleared and usable. + LBVH(LBVH&& other) noexcept; + /// @brief Move-assign; the moved-from object is left cleared and usable. + LBVH& operator=(LBVH&& other) noexcept; + LBVH(const LBVH&) = delete; + LBVH& operator=(const LBVH&) = delete; + + /// @brief Non-owning view of device-resident candidate pairs (SoA). + /// + /// The pointers address device memory owned by this LBVH and stay valid + /// until the next detect_*() or detect_*_device() call of the SAME + /// candidate type (both variants share one buffer per type), or until + /// clear(), build(), or destruction. The ids are int32_t, not index_t: + /// see the class comment. + struct DeviceCandidateView { + const int32_t* a = nullptr; ///< Device pointer to the first ids. + const int32_t* b = nullptr; ///< Device pointer to the second ids. + size_t size = 0; ///< Number of candidate pairs. + }; + + /// @brief Get the name of the broad phase method. + std::string name() const override { return "LBVH (CUDA)"; } + + using ipc::BroadPhase::build; + + /// @brief Build the broad phase for static collision detection. + /// The vertex boxes and everything after them are built on the device. + /// @param vertices Vertex positions (rowwise, |V| × 2 or |V| × 3). + /// @param edges Collision mesh edges. + /// @param faces Collision mesh faces. + /// @param inflation_radius Radius of inflation around all elements. + void build( + Eigen::ConstRef vertices, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius = 0) override; + + /// @brief Build the broad phase for continuous collision detection. + /// The vertex boxes and everything after them are built on the device. + /// @param vertices_t0 Starting vertex positions (rowwise). + /// @param vertices_t1 Ending vertex positions (rowwise). + /// @param edges Collision mesh edges. + /// @param faces Collision mesh faces. + /// @param inflation_radius Radius of inflation around all elements. + void build( + Eigen::ConstRef vertices_t0, + Eigen::ConstRef vertices_t1, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius = 0) override; + + // BroadPhase::build(const AABBs&, edges, faces, dim) is inherited: it + // copies the boxes into BroadPhase::vertex_boxes and calls the protected + // build(edges, faces) below, which uploads them. + + /// @brief Clear any built data. Device memory is retained for reuse. + void clear() override; + + // ------------------------------------------------------------------ + // BroadPhase interface (host-materializing). The BVH descent (AABB overlap + // + triangular dedup) and the mesh-connectivity (shared-vertex) exclusion + // both run on the device. The user vertex filter is applied on the host + // only when it is not accept-all (see can_*_collide); the output matches + // the CPU ipc::LBVH exactly for any filter. Like every BroadPhase, these + // clear the output vector first. Each invalidates any DeviceCandidateView + // of the same type. + + void detect_vertex_vertex_candidates( + std::vector& candidates) const override; + void detect_edge_vertex_candidates( + std::vector& candidates) const override; + void detect_edge_edge_candidates( + std::vector& candidates) const override; + void detect_face_vertex_candidates( + std::vector& candidates) const override; + void detect_edge_face_candidates( + std::vector& candidates) const override; + void detect_face_face_candidates( + std::vector& candidates) const override; + + // ------------------------------------------------------------------ + // Device-resident candidate accessors (GPU-native pipeline). Each runs the + // filtered traversal and returns a view of the connectivity-filtered pairs + // left on the device. For the default (accept-all) vertex filter the view + // is the exact candidate set; otherwise it is a connectivity-filtered + // superset the caller must trim with the user vertex filter. See + // DeviceCandidateView for the view's lifetime. + + DeviceCandidateView detect_vertex_vertex_candidates_device() const; + DeviceCandidateView detect_edge_vertex_candidates_device() const; + DeviceCandidateView detect_edge_edge_candidates_device() const; + DeviceCandidateView detect_face_vertex_candidates_device() const; + DeviceCandidateView detect_edge_face_candidates_device() const; + DeviceCandidateView detect_face_face_candidates_device() const; + + // ------------------------------------------------------------------ + // Sizes (cheap; no device->host node copy). + + /// @brief Number of nodes in the vertex BVH (2 * n_leaves - 1, or 0). + size_t num_vertex_nodes() const; + /// @brief Number of nodes in the edge BVH (2 * n_leaves - 1, or 0). + size_t num_edge_nodes() const; + /// @brief Number of nodes in the face BVH (2 * n_leaves - 1, or 0). + size_t num_face_nodes() const; + + // ------------------------------------------------------------------ + // Debug / validation: copy the device trees back to the host. + + /// @brief Copy the vertex BVH back to the host. + void vertex_nodes_to_host( + ipc::LBVH::Nodes& nodes, + ipc::LBVH::RightmostLeaves& rightmost_leaves) const; + /// @brief Copy the edge BVH back to the host. + void edge_nodes_to_host( + ipc::LBVH::Nodes& nodes, + ipc::LBVH::RightmostLeaves& rightmost_leaves) const; + /// @brief Copy the face BVH back to the host. + void face_nodes_to_host( + ipc::LBVH::Nodes& nodes, + ipc::LBVH::RightmostLeaves& rightmost_leaves) const; + + // ------------------------------------------------------------------ + + /// @brief Opaque implementation type (pimpl), keeping CUDA types out of + /// this header. Defined in lbvh_impl.cuh for the ipc::cuda implementation + /// files (.cu) only; it is not reachable through this interface and is + /// not part of it. Declared here (rather than privately) only so that the + /// implementation's file-local helpers can name it. + struct Impl; + +protected: + /// @brief Build the device trees from BroadPhase::vertex_boxes. + /// + /// Reached through the inherited BroadPhase::build(const AABBs&, edges, + /// faces, dim), which has copied the caller's boxes into vertex_boxes. The + /// boxes are uploaded and then, as in ipc::LBVH, cleared from the host: + /// the device trees make them redundant. + /// + /// @param edges Collision mesh edges. + /// @param faces Collision mesh faces. + void build( + Eigen::ConstRef edges, + Eigen::ConstRef faces) override; + + // Host-side collision filters, used to trim the device-emitted candidates + // only when the user vertex filter is not accept-all (the device already + // excludes shared-vertex pairs). Mirror ipc::LBVH, backed by host copies + // of the connectivity. + bool can_edge_vertex_collide(size_t ei, size_t vi) const override; + bool can_edges_collide(size_t eai, size_t ebi) const override; + bool can_face_vertex_collide(size_t fi, size_t vi) const override; + bool can_edge_face_collide(size_t ei, size_t fi) const override; + bool can_faces_collide(size_t fai, size_t fbi) const override; + +private: + /// @brief The implementation, created on first use. + /// + /// A moved-from object hands its Impl over and is left with none; rather + /// than guard every method against that, this re-seeds it on demand, so a + /// moved-from object behaves exactly like a cleared one. Allocating here + /// instead of in the move is what lets the move be noexcept. Mutable + /// because the const detect_*() methods, which the BroadPhase interface + /// requires, do use (and mutate) the device buffers; see the class comment. + Impl& impl() const; + + mutable std::unique_ptr m_impl; ///< See impl(). +}; + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/broad_phase/cuda/lbvh_impl.cuh b/src/ipc/broad_phase/cuda/lbvh_impl.cuh new file mode 100644 index 000000000..4beda3c1b --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh_impl.cuh @@ -0,0 +1,147 @@ +// Definition of the pimpl struct of ipc::cuda::LBVH. This header is CUDA-only +// and must be included from the ipc::cuda implementation files (.cu) +// exclusively. + +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include + +#include +#include +#include + +namespace ipc::cuda { + +struct LBVH::Impl { + /// @brief A single BVH: the node array (root at index 0, same 32-byte + /// ipc::LBVH::Node layout as the CPU path) plus the per-node Morton-sorted + /// rightmost-leaf index used to skip subtrees in triangular traversal. + struct DeviceBVH { + DeviceBuffer nodes; + DeviceBuffer rightmost_leaves; + + /// @brief The number of leaves, derived from the node count (a BVH + /// over n primitives has 2n - 1 nodes) so the two cannot disagree. + int n_leaves() const + { + return nodes.empty() ? 0 : static_cast((nodes.size() + 1) / 2); + } + + void clear() + { + nodes.clear(); + rightmost_leaves.clear(); + } + }; + + DeviceBVH vertex_bvh; + DeviceBVH edge_bvh; + DeviceBVH face_bvh; + + /// @brief Device-resident candidate pairs (SoA) for one collision type, + /// connectivity-filtered on the device. For the default (accept-all) vertex + /// filter this is already the exact candidate set; otherwise it is a + /// superset the host trims with the user filter. + /// + /// The buffers are grow-only pools. Their capacity is the largest count any + /// traversal of this type on this object has needed, and it survives + /// clear() -- which build() calls every timestep -- so a per-frame call + /// sizes its first pass from the previous frame's count and pays the + /// overflow-and-retry only when the count grows past every prior call. + /// Nothing is initialized: the kernel writes exactly the slots it fills. + /// The memory is released only by the destructor. + struct DeviceCandidates { + DeviceBuffer a; + DeviceBuffer b; + size_t count = 0; ///< The number of valid pairs in a/b. + + void clear() + { + count = 0; + a.clear(); + b.clear(); + } + }; + + DeviceCandidates vv_candidates; + DeviceCandidates ev_candidates; + DeviceCandidates ee_candidates; + DeviceCandidates fv_candidates; + DeviceCandidates ef_candidates; + DeviceCandidates ff_candidates; + + /// @brief Mesh connectivity, flat row-major (2 ids per edge, 3 per face). + /// On the device for the traversal's shared-vertex filter; mirrored on the + /// host for the can_*_collide filters, which run only when the user vertex + /// filter is not accept-all. + DeviceBuffer edges; + DeviceBuffer faces; + std::vector h_edges; + std::vector h_faces; + + /// @brief The Morton-normalization domain: the union of the vertex boxes. + /// A plain aggregate so it is trivially copyable and CUB can reduce it. + struct Domain { + double min[3]; + double max[3]; + }; + + // -- Build scratch ------------------------------------------------------ + // Persistent and grow-only so a per-frame rebuild allocates nothing (every + // cudaFree synchronizes the whole device). Nothing here is initialized + // except where a kernel needs it: the construction infos are zeroed and + // the roots set to -1 on every build. + + DeviceBuffer vertices_t0; ///< Uploaded positions (column-major). + DeviceBuffer vertices_t1; ///< Uploaded positions (column-major). + DeviceBuffer vbox_min; ///< Vertex box min corners (3 per box). + DeviceBuffer vbox_max; ///< Vertex box max corners (3 per box). + DeviceBuffer ebox_min; ///< Edge box min corners (3 per box). + DeviceBuffer ebox_max; ///< Edge box max corners (3 per box). + DeviceBuffer fbox_min; ///< Face box min corners (3 per box). + DeviceBuffer fbox_max; ///< Face box max corners (3 per box). + DeviceBuffer domain; ///< The normalization domain (1). + /// @brief The Morton codes (keys) and box ids (values) to sort, as the two + /// ping-pong buffers each an out-of-place radix sort needs: the codes + /// kernel fills [0], and CUB alternates between [0] and [1] per pass (see + /// build_tree()). + DeviceBuffer morton_codes[2]; + DeviceBuffer box_ids[2]; + DeviceBuffer> construction_infos; + DeviceBuffer roots; ///< One root index per BVH (3). + DeviceBuffer sort_temp; ///< CUB radix-sort storage. + DeviceBuffer reduce_temp; ///< CUB reduce storage. + DeviceBuffer counter; ///< Emitted-pair counter (1). + + /// @brief Serializes the detect_*() calls, which are const on the + /// BroadPhase interface but share the candidate buffers and the counter. + std::mutex mutex; + + /// @brief Forget the built trees, connectivity, and candidates. Every + /// allocation is kept for the next build (see DeviceCandidates). + void clear() + { + vertex_bvh.clear(); + edge_bvh.clear(); + face_bvh.clear(); + edges.clear(); + faces.clear(); + h_edges.clear(); + h_faces.clear(); + vv_candidates.clear(); + ev_candidates.clear(); + ee_candidates.clear(); + fv_candidates.clear(); + ef_candidates.clear(); + ff_candidates.clear(); + } +}; + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/broad_phase/details/CMakeLists.txt b/src/ipc/broad_phase/details/CMakeLists.txt new file mode 100644 index 000000000..942162c1f --- /dev/null +++ b/src/ipc/broad_phase/details/CMakeLists.txt @@ -0,0 +1,8 @@ +set(SOURCES + connectivity_filters.hpp + lbvh_build.hpp + lbvh_traverse.hpp + spatial_hash_impl.hpp +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) diff --git a/src/ipc/broad_phase/details/connectivity_filters.hpp b/src/ipc/broad_phase/details/connectivity_filters.hpp new file mode 100644 index 000000000..6fd4f08a4 --- /dev/null +++ b/src/ipc/broad_phase/details/connectivity_filters.hpp @@ -0,0 +1,193 @@ +#pragma once + +#include +#include + +#include + +namespace ipc::details { + +// Mesh-connectivity collision filters shared by every broad phase. +// +// ipc::BroadPhase, ipc::LBVH, ipc::SweepAndTiniestQueue, and ipc::cuda::LBVH +// each store the connectivity differently -- in the AABBs' vertex_ids, in a +// dedicated host copy, or in a host mirror of the device arrays -- but they all +// apply the same rule: exclude primitive pairs that share a vertex, then accept +// the pair only if the user vertex filter passes for at least one of the +// remaining vertex pairs. These take the vertex ids directly so each broad +// phase can supply them from whatever storage it has. +// +// The rule is split in two so its halves can live where they are needed: +// share_vertex() has no dependency on the user filter and is host/device, so +// the device traversal of ipc::cuda::LBVH applies exactly the same exclusion +// the host does; any_vertex_pair_can_collide() is the user-filter half, which +// only ever runs on the host. + +/// @brief Whether two primitives share a vertex id. +/// +/// A vertex primitive has one id, an edge two, a face three. Fully unrolled at +/// compile time, so every id stays in a register on the device. +/// +/// @tparam NA The number of vertex ids of the first primitive (1, 2, or 3). +/// @tparam NB The number of vertex ids of the second primitive (1, 2, or 3). +/// @tparam Id The vertex id type. +/// @param a The first primitive's vertex ids. +/// @param b The second primitive's vertex ids. +/// @return Whether any id of @p a equals any id of @p b. +template +IPC_TOOLKIT_HOST_DEVICE inline bool +share_vertex(const Id (&a)[NA], const Id (&b)[NB]) +{ + static_assert(NA >= 1 && NA <= 3 && NB >= 1 && NB <= 3); + bool shared = false; + // The unroll pragma is nvcc's; the host compilers unroll these constant + // trip counts on their own and would only warn about the unknown pragma. +#ifdef __CUDA_ARCH__ +#pragma unroll +#endif + for (int i = 0; i < NA; ++i) { +#ifdef __CUDA_ARCH__ +#pragma unroll +#endif + for (int j = 0; j < NB; ++j) { + shared |= a[i] == b[j]; + } + } + return shared; +} + +/// @brief Whether the user vertex filter passes for at least one pair of +/// vertices drawn from two primitives. +/// @tparam NA The number of vertex ids of the first primitive (1, 2, or 3). +/// @tparam NB The number of vertex ids of the second primitive (1, 2, or 3). +/// @tparam Id The vertex id type. +/// @param a The first primitive's vertex ids. +/// @param b The second primitive's vertex ids. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether any vertex pair passes the filter. Always true for the +/// accept-all filter, without evaluating anything. +template +inline bool any_vertex_pair_can_collide( + const Id (&a)[NA], + const Id (&b)[NB], + const CollisionFilter& can_vertices_collide) +{ + static_assert(NA >= 1 && NA <= 3 && NB >= 1 && NB <= 3); + if (can_vertices_collide.accepts_all()) { + return true; + } + for (int i = 0; i < NA; ++i) { + for (int j = 0; j < NB; ++j) { + if (can_vertices_collide(a[i], b[j])) { + return true; + } + } + } + return false; +} + +/// @brief Whether an edge and a vertex can collide. +/// @param e0i The first vertex of the edge. +/// @param e1i The second vertex of the edge. +/// @param vi The vertex. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether the pair should be considered for collision. +inline bool can_edge_vertex_collide( + const index_t e0i, + const index_t e1i, + const size_t vi, + const CollisionFilter& can_vertices_collide) +{ + const index_t e[2] = { e0i, e1i }; + const index_t v[1] = { static_cast(vi) }; + return !share_vertex(e, v) + && any_vertex_pair_can_collide(v, e, can_vertices_collide); +} + +/// @brief Whether two edges can collide. +/// @param ea0i The first vertex of the first edge. +/// @param ea1i The second vertex of the first edge. +/// @param eb0i The first vertex of the second edge. +/// @param eb1i The second vertex of the second edge. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether the pair should be considered for collision. +inline bool can_edges_collide( + const index_t ea0i, + const index_t ea1i, + const index_t eb0i, + const index_t eb1i, + const CollisionFilter& can_vertices_collide) +{ + const index_t ea[2] = { ea0i, ea1i }; + const index_t eb[2] = { eb0i, eb1i }; + return !share_vertex(ea, eb) + && any_vertex_pair_can_collide(ea, eb, can_vertices_collide); +} + +/// @brief Whether a face and a vertex can collide. +/// @param f0i The first vertex of the face. +/// @param f1i The second vertex of the face. +/// @param f2i The third vertex of the face. +/// @param vi The vertex. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether the pair should be considered for collision. +inline bool can_face_vertex_collide( + const index_t f0i, + const index_t f1i, + const index_t f2i, + const size_t vi, + const CollisionFilter& can_vertices_collide) +{ + const index_t f[3] = { f0i, f1i, f2i }; + const index_t v[1] = { static_cast(vi) }; + return !share_vertex(f, v) + && any_vertex_pair_can_collide(v, f, can_vertices_collide); +} + +/// @brief Whether an edge and a face can intersect. +/// @param e0i The first vertex of the edge. +/// @param e1i The second vertex of the edge. +/// @param f0i The first vertex of the face. +/// @param f1i The second vertex of the face. +/// @param f2i The third vertex of the face. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether the pair should be considered for collision. +inline bool can_edge_face_collide( + const index_t e0i, + const index_t e1i, + const index_t f0i, + const index_t f1i, + const index_t f2i, + const CollisionFilter& can_vertices_collide) +{ + const index_t e[2] = { e0i, e1i }; + const index_t f[3] = { f0i, f1i, f2i }; + return !share_vertex(e, f) + && any_vertex_pair_can_collide(e, f, can_vertices_collide); +} + +/// @brief Whether two faces can collide. +/// @param fa0i The first vertex of the first face. +/// @param fa1i The second vertex of the first face. +/// @param fa2i The third vertex of the first face. +/// @param fb0i The first vertex of the second face. +/// @param fb1i The second vertex of the second face. +/// @param fb2i The third vertex of the second face. +/// @param can_vertices_collide The user vertex filter. +/// @return Whether the pair should be considered for collision. +inline bool can_faces_collide( + const index_t fa0i, + const index_t fa1i, + const index_t fa2i, + const index_t fb0i, + const index_t fb1i, + const index_t fb2i, + const CollisionFilter& can_vertices_collide) +{ + const index_t fa[3] = { fa0i, fa1i, fa2i }; + const index_t fb[3] = { fb0i, fb1i, fb2i }; + return !share_vertex(fa, fb) + && any_vertex_pair_can_collide(fa, fb, can_vertices_collide); +} + +} // namespace ipc::details diff --git a/src/ipc/broad_phase/details/lbvh_build.hpp b/src/ipc/broad_phase/details/lbvh_build.hpp new file mode 100644 index 000000000..ebe393338 --- /dev/null +++ b/src/ipc/broad_phase/details/lbvh_build.hpp @@ -0,0 +1,264 @@ +#pragma once + +#include +#include +#include +#include // for infinity() + +#include + +#include +#include +#include + +namespace ipc::details { + +// The LBVH construction of Apetrei [2014], shared by ipc::LBVH and +// ipc::cuda::LBVH. +// +// Everything here is host/device and addresses the tree through raw pointers, +// so the CPU can drive it from a tbb::parallel_for over std::vectors and the +// GPU from a kernel over thrust::device_vectors. What each platform still owns +// is only the parallel launch, the sort, and the two policies these take: how +// to read a sorted Morton code, and how to perform the atomic arrival. + +/// @brief Rounds a double AABB outward to the smallest enclosing float AABB. +/// +/// Each corner is nudged to the next representable float away from the box, so +/// the float AABB always encloses the double one and never clips a primitive. +/// +/// @note The directions are typed floats (not the INFINITY macro, whose type +/// is unspecified) so this resolves to the same nextafterf on host and device. +/// +/// @param box_min The minimum corner of the double AABB. +/// @param box_max The maximum corner of the double AABB. +/// @param[out] node The node whose AABB is set. +IPC_TOOLKIT_HOST_DEVICE inline void set_inflated_aabb( + const Eigen::Array3d& box_min, + const Eigen::Array3d& box_max, + LBVH::Node& node) +{ + for (int k = 0; k < 3; ++k) { + node.aabb_min[k] = + nextafterf(static_cast(box_min[k]), -infinity()); + node.aabb_max[k] = + nextafterf(static_cast(box_max[k]), infinity()); + } +} + +/// @brief Initializes the leaf node for sorted position i. +/// +/// Leaves occupy the upper half of the node array, at n_leaves - 1 + i. +/// +/// @param i The leaf's position in the Morton-sorted order. +/// @param n_leaves The number of leaves. +/// @param box_id The id of the primitive this leaf holds. +/// @param box_min The minimum corner of the primitive's AABB. +/// @param box_max The maximum corner of the primitive's AABB. +/// @param[out] nodes The BVH nodes. +/// @param[out] rightmost_leaves The per-node rightmost-leaf indices. +IPC_TOOLKIT_HOST_DEVICE inline void init_leaf_node( + const int i, + const int n_leaves, + const index_t box_id, + const Eigen::Array3d& box_min, + const Eigen::Array3d& box_max, + LBVH::Node* nodes, + int32_t* rightmost_leaves) +{ + LBVH::Node leaf; + set_inflated_aabb(box_min, box_max, leaf); + leaf.primitive_id = static_cast(box_id); + leaf.is_inner_marker = 0; + + const int leaf_idx = n_leaves - 1 + i; + nodes[leaf_idx] = leaf; + rightmost_leaves[leaf_idx] = i; // a leaf's rightmost leaf is itself +} + +/// @brief Returns the length of the common leading-bit prefix of the sorted +/// Morton codes at positions i and j, or -1 when j is out of bounds. +/// @tparam CodeAt Callable (int) -> uint64_t returning a sorted Morton code. +/// @param code_at The sorted Morton code accessor. +/// @param n_leaves The number of codes. +/// @param i The first sorted position. +/// @param j The second sorted position. +/// @return The common-prefix length, or -1 when j is out of bounds. +template +IPC_TOOLKIT_HOST_DEVICE inline int +delta(CodeAt&& code_at, const int n_leaves, const int i, const int j) +{ + if (j < 0 || j >= n_leaves) { + return -1; + } + return morton_common_prefix(code_at(i), i, code_at(j), j); +} + +/// @brief Whether the subtree spanning [left_key, right_key] is its parent's +/// left child. +/// +/// The two candidate parents are internal node right_key (which would make +/// this the left child) and internal node left_key - 1 (the right child). +/// delta() grows with the similarity of the codes, so the nearer ancestor is +/// the one with the LARGER delta -- hence ">". +/// +/// At the boundaries only one candidate exists: a range starting at 0 has no +/// node -1 to its left, and a range ending at n_leaves - 1 has no node +/// n_leaves - 1 to its right. +/// +/// @tparam CodeAt Callable (int) -> uint64_t returning a sorted Morton code. +/// @param code_at The sorted Morton code accessor. +/// @param n_leaves The number of leaves. +/// @param left_key The left endpoint of the subtree's sorted-key range. +/// @param right_key The right endpoint of the subtree's sorted-key range. +/// @return Whether this subtree is its parent's left child. +template +IPC_TOOLKIT_HOST_DEVICE inline bool is_left_child( + CodeAt&& code_at, + const int n_leaves, + const int left_key, + const int right_key) +{ + return left_key == 0 + || (right_key != n_leaves - 1 + && delta(code_at, n_leaves, right_key, right_key + 1) + > delta(code_at, n_leaves, left_key - 1, left_key)); +} + +/// @brief Walks one leaf up to the root, building the hierarchy, the internal +/// AABBs and the rightmost-leaf indices (Apetrei [2014], Fig. 2). +/// +/// Each leaf's thread climbs toward the root, choosing its parent in O(1) from +/// the delta values at the two ends of its current key range. At every parent +/// the first of the two arriving threads stops and the second continues, so +/// whoever continues knows both children are complete. +/// +/// In this layout internal node j always splits between sorted keys j and +/// j + 1, so the root is generally NOT at index 0; swap_root_to_zero() moves +/// it there afterwards, which is what the traversal expects. +/// +/// @tparam Counter The visitation counter's type (see +/// ipc::LBVH::ConstructionInfo). +/// @tparam CodeAt Callable (int) -> uint64_t returning a sorted Morton code. +/// @tparam Arrive Callable (Counter&) -> int that atomically increments the +/// counter and returns its previous value. It must order this thread's earlier +/// writes before the increment, and -- when it returns nonzero, so this thread +/// continues -- order the increment before this thread's later reads. Without +/// both halves the continuing thread can read a stale sibling. +/// +/// @param i The leaf's position in the Morton-sorted order. +/// @param n_leaves The number of leaves. +/// @param code_at The sorted Morton code accessor. +/// @param[in,out] nodes The BVH nodes; the leaves must already be initialized. +/// @param[in,out] rightmost_leaves The per-node rightmost-leaf indices. +/// @param[in,out] infos The per-node construction scratch, zero-initialized. +/// @param arrive The atomic arrival gate. +/// @return The root's index if this leaf's walk reached the root, else -1. +template +IPC_TOOLKIT_HOST_DEVICE int build_hierarchy_from_leaf( + const int i, + const int n_leaves, + CodeAt&& code_at, + LBVH::Node* nodes, + int32_t* rightmost_leaves, + LBVH::ConstructionInfo* infos, + Arrive&& arrive) +{ + // A single-leaf tree is its own root and has no internal nodes to build. + if (n_leaves == 1) { + return i == 0 ? 0 : -1; + } + + // Invariant: the current subtree covers the sorted-key range + // [left_key, right_key]. + int left_key = i; + int right_key = i; + int current_node = n_leaves - 1 + i; + + while (true) { + const bool is_child_a = + is_left_child(code_at, n_leaves, left_key, right_key); + const int parent = is_child_a ? right_key : left_key - 1; + + // Write the child pointer and the range endpoint onto the parent. The + // left child writes .left and the left endpoint, the right child + // writes .right and the right endpoint. + if (is_child_a) { + nodes[parent].left = current_node; + infos[parent].left_range = left_key; + } else { + nodes[parent].right = current_node; + infos[parent].right_range = right_key; + } + + if (arrive(infos[parent].visitation_count) == 0) { + return -1; // first thread to arrive here -> done + } + + // Second thread to arrive: both children are complete, so their AABBs + // and rightmost leaves can be combined into the parent's. + assert(nodes[parent].is_inner()); + const LBVH::Node& child_a = nodes[nodes[parent].left]; + const LBVH::Node& child_b = nodes[nodes[parent].right]; + nodes[parent].aabb_min = child_a.aabb_min.min(child_b.aabb_min); + nodes[parent].aabb_max = child_a.aabb_max.max(child_b.aabb_max); + + const int32_t rightmost_a = rightmost_leaves[nodes[parent].left]; + const int32_t rightmost_b = rightmost_leaves[nodes[parent].right]; + rightmost_leaves[parent] = + rightmost_a > rightmost_b ? rightmost_a : rightmost_b; + + // Reconstruct the parent's full key range and continue upward. + left_key = infos[parent].left_range; + right_key = infos[parent].right_range; + current_node = parent; + + if (left_key == 0 && right_key == n_leaves - 1) { + return current_node; // the root's AABB is complete + } + } +} + +/// @brief Swaps the node and rightmost-leaf entries at index 0 and the root, so +/// that traversal can start at index 0. +/// +/// The root is never any node's child, so no pointer needs rewriting to reach +/// its new home at 0. Pointers that referenced the old node 0 do, which is what +/// patch_left_pointer() handles. +/// +/// @param[in,out] nodes The BVH nodes. +/// @param[in,out] rightmost_leaves The per-node rightmost-leaf indices. +/// @param root The root's index, which must be greater than 0. +IPC_TOOLKIT_HOST_DEVICE inline void +swap_root_to_zero(LBVH::Node* nodes, int32_t* rightmost_leaves, const int root) +{ + assert(root > 0); + + const LBVH::Node node = nodes[0]; + nodes[0] = nodes[root]; + nodes[root] = node; + + const int32_t rightmost = rightmost_leaves[0]; + rightmost_leaves[0] = rightmost_leaves[root]; + rightmost_leaves[root] = rightmost; +} + +/// @brief Rewrites a left pointer that referenced the old node 0 to the root's +/// new location, after swap_root_to_zero(). +/// +/// Apetrei's layout guarantees node 0's subtree has left_key == 0, so node 0 is +/// only ever written as a LEFT child. Two things follow: only .left pointers +/// need patching, and swapping node 0 away cannot leave a node whose .right -- +/// which aliases is_inner_marker -- is 0 and so reads as a leaf. +/// +/// @param[in,out] node The node to patch. +/// @param root The new location of the old node 0. +IPC_TOOLKIT_HOST_DEVICE inline void +patch_left_pointer(LBVH::Node& node, const int root) +{ + if (node.is_inner() && node.left == 0) { + node.left = root; + } +} + +} // namespace ipc::details diff --git a/src/ipc/broad_phase/details/lbvh_traverse.hpp b/src/ipc/broad_phase/details/lbvh_traverse.hpp new file mode 100644 index 000000000..f4ba48c92 --- /dev/null +++ b/src/ipc/broad_phase/details/lbvh_traverse.hpp @@ -0,0 +1,178 @@ +#pragma once + +#include +#include +#include // for MORTON_KEY_BITS +#include // for log_and_throw_error (host only) +#include // for any(bool) + +#include +#include +#include + +namespace ipc::details { + +/// @brief Descends a target BVH for one query -- or one batch of queries -- +/// reporting every target leaf whose AABB overlaps. +/// +/// A stackless-style descent with an explicit stack: at each inner node the +/// overlapping children are handled immediately if they are leaves, descended +/// into if only one is inner, and the right one postponed on the stack if both +/// are. The root lives at index 0, and LBVH::Node::INVALID_POINTER (which is +/// also 0) doubles as the stack's bottom sentinel -- popping it ends the walk, +/// because no node other than the root ever lives at index 0. +/// +/// This is shared by ipc::LBVH -- both its scalar and its SIMD traversal -- and +/// ipc::cuda::LBVH. What differs between them is how a node is tested against +/// the query and what happens on an overlap, which is why both are policies: +/// the scalar host and the device test one query and get a bool; the SIMD host +/// tests a batch of queries at once and gets a lane mask. On an overlap the +/// host filters and appends to a std::vector, while the device filters against +/// the mesh connectivity and appends through an atomic counter. +/// +/// @tparam triangular Self-collision: skip any subtree lying entirely to the +/// left of the query, so each unordered pair is reported exactly once. +/// @tparam Intersects Callable (const LBVH::Node&) -> Mask, where Mask is bool +/// for a single query or a lane mask for a batch of queries. Either must +/// support any(mask) -- ipc::any for bool, xsimd::any by ADL for a batch -- +/// and Mask(false). +/// @tparam Emit Callable (const LBVH::Node& leaf, int leaf_idx, const Mask&) +/// -> void, invoked for every target leaf overlapping at least one query. It +/// owns both the collision filtering and the recording of the pair. For a +/// batch, it must re-check the triangular skip per lane against +/// target_rightmost[leaf_idx], because the skip below is conservative for the +/// batch as a whole. +/// +/// @param query_leaf_idx The query's position in its own Morton-sorted leaf +/// order -- for a batch, the smallest position in it. Used only by the +/// triangular skip. +/// @param target The target BVH's nodes, root at index 0. +/// @param target_size The number of nodes in the target BVH. +/// @param target_rightmost The target's per-node rightmost-leaf indices. Used +/// only by the triangular skip. +/// @param intersects The per-node overlap test. +/// @param emit The per-overlap callback. +template +IPC_TOOLKIT_HOST_DEVICE void traverse_lbvh( + const int query_leaf_idx, + const LBVH::Node* target, + const int target_size, + const int32_t* target_rightmost, + Intersects&& intersects, + Emit&& emit) +{ + using Mask = std::decay_t; + + // A fixed-size stack keeps the descent free of dynamic allocation. Every + // internal node on a root-to- leaf path splits at a distinct, strictly + // increasing prefix length of the MORTON_KEY_BITS-bit key, so a path holds + // at most MORTON_KEY_BITS internal nodes, and at most one right child per + // internal node on the current path is ever pending -- plus the sentinel. + // The overflow check below can therefore never fire on a well-formed tree; + // it turns a malformed one into a hard failure instead of an out-of-bounds + // write. + constexpr int MAX_STACK_SIZE = MORTON_KEY_BITS + 1; + int stack[MAX_STACK_SIZE]; + int stack_ptr = 0; + stack[stack_ptr++] = LBVH::Node::INVALID_POINTER; + + int node_idx = 0; // root + do { + const LBVH::Node& node = target[node_idx]; + + if (target_size == 1) { // only the root, which is therefore a leaf + assert(node.is_leaf()); + if constexpr (triangular) { + break; // a lone primitive cannot collide with itself + } + const Mask mask = intersects(node); + if (any(mask)) { + emit(node, node_idx, mask); + } + break; + } + + assert(node.is_inner()); // so .left and .right are valid pointers + +#if !defined(__CUDA_ARCH__) && (defined(__GNUC__) || defined(__clang__)) + // Prefetch the children to reduce cache misses. The device needs no + // equivalent; it hides the latency with its other resident warps. + __builtin_prefetch(&target[node.left], 0, 1); + __builtin_prefetch(&target[node.right], 0, 1); +#endif + + const LBVH::Node& child_l = target[node.left]; + const LBVH::Node& child_r = target[node.right]; + Mask intersects_l = intersects(child_l); + Mask intersects_r = intersects(child_r); + + // Ignore a subtree lying entirely to the query's left; that pair is + // reported when the other primitive is the query instead. For a batch + // this uses its smallest query position, so it is conservative: a + // subtree left of every query is skipped here, and emit() re-checks + // the rest per lane. + if constexpr (triangular) { + if constexpr (std::is_same_v) { + // A scalar overlap result is a free branch, so test it first + // and skip the rightmost-leaf load when nothing overlaps. + if (intersects_l + && target_rightmost[node.left] <= query_leaf_idx) { + intersects_l = false; + } + if (intersects_r + && target_rightmost[node.right] <= query_leaf_idx) { + intersects_r = false; + } + } else { + // For a batch, any(mask) is a movemask plus a second + // data-dependent branch per child, which costs more than the + // one load it would save -- so test the position alone. + if (target_rightmost[node.left] <= query_leaf_idx) { + intersects_l = Mask(false); + } + if (target_rightmost[node.right] <= query_leaf_idx) { + intersects_r = Mask(false); + } + } + } + + const bool any_l = any(intersects_l); + const bool any_r = any(intersects_r); + + // An overlapped leaf is a candidate. + if (any_l && child_l.is_leaf()) { + emit(child_l, node.left, intersects_l); + } + if (any_r && child_r.is_leaf()) { + emit(child_r, node.right, intersects_r); + } + + // An overlapped inner node is descended into. + const bool traverse_l = any_l && !child_l.is_leaf(); + const bool traverse_r = any_r && !child_r.is_leaf(); + + if (!traverse_l && !traverse_r) { + assert(stack_ptr > 0); + node_idx = stack[--stack_ptr]; + } else { + node_idx = traverse_l ? node.left : node.right; + if (traverse_l && traverse_r) { + if (stack_ptr >= MAX_STACK_SIZE) { + // Unreachable on a well-formed tree (see MAX_STACK_SIZE); + // fail hard rather than write past the stack. +#ifdef __CUDA_ARCH__ + __trap(); +#else + log_and_throw_error( + "ipc::details::traverse_lbvh: traversal stack " + "overflow; the BVH is deeper than the Morton key " + "width allows and so is malformed"); +#endif + } + stack[stack_ptr++] = node.right; // postpone the right child + } + } + } while (node_idx != LBVH::Node::INVALID_POINTER); +} + +} // namespace ipc::details diff --git a/src/ipc/broad_phase/hash_grid.cpp b/src/ipc/broad_phase/hash_grid.cpp index 89eeb86c2..58e1c41dd 100644 --- a/src/ipc/broad_phase/hash_grid.cpp +++ b/src/ipc/broad_phase/hash_grid.cpp @@ -286,6 +286,7 @@ void HashGrid::detect_candidates( void HashGrid::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( vertex_items, vertex_boxes, can_vertices_collide, candidates); } @@ -293,6 +294,7 @@ void HashGrid::detect_vertex_vertex_candidates( void HashGrid::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_items, vertex_items, edge_boxes, vertex_boxes, std::bind(&HashGrid::can_edge_vertex_collide, this, _1, _2), @@ -302,6 +304,7 @@ void HashGrid::detect_edge_vertex_candidates( void HashGrid::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_items, edge_boxes, std::bind(&HashGrid::can_edges_collide, this, _1, _2), candidates); @@ -310,6 +313,7 @@ void HashGrid::detect_edge_edge_candidates( void HashGrid::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( face_items, vertex_items, face_boxes, vertex_boxes, std::bind(&HashGrid::can_face_vertex_collide, this, _1, _2), @@ -319,6 +323,7 @@ void HashGrid::detect_face_vertex_candidates( void HashGrid::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( edge_items, face_items, edge_boxes, face_boxes, std::bind(&HashGrid::can_edge_face_collide, this, _1, _2), candidates); @@ -327,6 +332,7 @@ void HashGrid::detect_edge_face_candidates( void HashGrid::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); detect_candidates( face_items, face_boxes, std::bind(&HashGrid::can_faces_collide, this, _1, _2), candidates); diff --git a/src/ipc/broad_phase/lbvh.cpp b/src/ipc/broad_phase/lbvh.cpp index 97ba8a9c6..e39fee3a4 100644 --- a/src/ipc/broad_phase/lbvh.cpp +++ b/src/ipc/broad_phase/lbvh.cpp @@ -1,5 +1,8 @@ #include "lbvh.hpp" +#include +#include +#include #include #include #include @@ -24,23 +27,6 @@ using namespace std::placeholders; namespace ipc { -namespace { - // Helper to safely convert double AABB to float AABB - inline void assign_inflated_aabb(const AABB& box, LBVH::Node& node) - { - // Round Min down - node.aabb_min = box.min.unaryExpr([](double val) { - return std::nextafter( - float(val), -std::numeric_limits::infinity()); - }); - // Round Max up - node.aabb_max = box.max.unaryExpr([](double val) { - return std::nextafter( - float(val), std::numeric_limits::infinity()); - }); - } -} // namespace - LBVH::LBVH() : BroadPhase() { static_assert( @@ -95,42 +81,6 @@ void LBVH::build( face_boxes.clear(); } -namespace { - /// Returns the number of common leading bits (CLZ of XOR) between sorted - /// Morton codes at positions i and j. code_i is the Morton code at position - /// i, passed explicitly to avoid a redundant lookup. Returns -1 when j is - /// out of bounds. Duplicate codes fall back to CLZ of the index XOR - /// (offset by 32 so it sorts after any code-level difference). - int delta( - const LBVH::MortonCodeElements& sorted_morton_codes, - int i, - uint64_t code_i, - int j) - { - if (j < 0 || j >= sorted_morton_codes.size()) { - return -1; - } - uint64_t code_j = sorted_morton_codes[j].morton_code; - if (code_i == code_j) { - // handle duplicate morton codes - int element_idx_i = i; - int element_idx_j = j; - - // add 32 for common prefix of code_i ^ code_j -#if defined(__GNUC__) || defined(__clang__) - return 32 + __builtin_clz(element_idx_i ^ element_idx_j); -#elif defined(WIN32) - return 32 + __lzcnt(element_idx_i ^ element_idx_j); -#endif - } -#if defined(__GNUC__) || defined(__clang__) - return __builtin_clzll(code_i ^ code_j); -#elif defined(WIN32) - return __lzcnt64(code_i ^ code_j); -#endif - } -} // namespace - void LBVH::init_bvh( const AABBs& boxes, Nodes& lbvh, RightmostLeaves& rightmost_leaves) const { @@ -150,21 +100,12 @@ void LBVH::init_bvh( IPC_TOOLKIT_PROFILE_BLOCK("compute_morton_codes"); const Eigen::Array3d mesh_width_inv = - 1.0 / (mesh_aabb.max - mesh_aabb.min); + morton_domain_width_inv(mesh_aabb.min, mesh_aabb.max); tbb::parallel_for(size_t(0), boxes.size(), [&](size_t i) { const auto& box = boxes[i]; - const Eigen::Array3d center = 0.5 * (box.min + box.max); - const Eigen::Array3d mapped_center = - (center - mesh_aabb.min) * mesh_width_inv; - - if (dim == 2) { - morton_codes[i].morton_code = - morton_2D(mapped_center.x(), mapped_center.y()); - } else { - morton_codes[i].morton_code = morton_3D( - mapped_center.x(), mapped_center.y(), mapped_center.z()); - } + morton_codes[i].morton_code = morton_code( + 0.5 * (box.min + box.max), mesh_aabb.min, mesh_width_inv, dim); morton_codes[i].box_id = i; }); } @@ -180,7 +121,6 @@ void LBVH::init_bvh( assert(boxes.size() <= std::numeric_limits::max()); const int N_LEAVES = int(boxes.size()); - const int LEAF_OFFSET = N_LEAVES - 1; if (rightmost_leaves.size() != lbvh.size()) { rightmost_leaves.resize(lbvh.size()); @@ -196,136 +136,46 @@ void LBVH::init_bvh( } // Apetrei 2014: single bottom-up pass that simultaneously builds the - // hierarchy and computes bounding boxes. Each leaf thread walks toward the - // root, choosing its parent in O(1) by comparing the CLZ-delta values at - // the two ends of its current key range. - // - // In this layout internal node j always splits between sorted keys j and - // j+1. The root is NOT necessarily at index 0, so after construction we - // swap the root into position 0 to match the traversal code's expectation. + // hierarchy and computes bounding boxes. See + // ipc::details::build_hierarchy_from_leaf(), shared with the device build + // in ipc::cuda::LBVH. std::atomic root_idx(-1); { IPC_TOOLKIT_PROFILE_BLOCK("build_hierarchy_and_boxes"); tbb::parallel_for(0, N_LEAVES, [&](int i) { - // --- Initialize leaf node --- - { - const auto& box = boxes[morton_codes[i].box_id]; - - Node leaf_node; // Create leaf node - assign_inflated_aabb(box, leaf_node); - leaf_node.primitive_id = morton_codes[i].box_id; - leaf_node.is_inner_marker = 0; - lbvh[LEAF_OFFSET + i] = leaf_node; // Store leaf - // A leaf's rightmost leaf is itself - rightmost_leaves[LEAF_OFFSET + i] = i; - } - - // --- Bottom-up walk (Apetrei 2014, Fig. 2) --- - // Invariant: the current subtree covers the sorted-key range - // [left_key, right_key]. - int left_key = i; - int right_key = i; - int current_node = LEAF_OFFSET + i; - - while (true) { - // Choose parent. Candidates are internal node right_key - // (current becomes its left / childA) or internal node - // left_key-1 (current becomes its right / childB). Our delta() - // returns CLZ (higher = more-similar = finer split), so the - // CLOSER ancestor has the LARGER delta — hence ">". - // - // Boundary rules: - // left_key == 0 → must be childA (no node -1) - // right_key == n-1 → must be childB (no node n-1) - const bool is_child_a = (left_key == 0) - || (right_key != N_LEAVES - 1 - && delta( - morton_codes, right_key, - morton_codes[right_key].morton_code, - right_key + 1) - > delta( - morton_codes, left_key - 1, - morton_codes[left_key - 1].morton_code, - left_key)); - const int parent = is_child_a ? right_key : left_key - 1; - - auto& info = construction_infos[parent]; - - // Write the child pointer on the parent node. - // childA writes .left; childB writes .right. - if (is_child_a) { - lbvh[parent].left = current_node; - info.left_range = left_key; - } else { - lbvh[parent].right = current_node; - info.right_range = right_key; - } - - // Atomic arrival gate: the first thread to reach this parent - // stops; the second thread proceeds (it now knows both children - // are complete). - - if (info.visitation_count++ == 0) { - // this is the first thread that arrived at this - // node -> finished - break; - } - // this is the second thread that arrived at this node, - // both children are computed -> compute aabb union and - // continue - assert(lbvh[parent].is_inner()); - const Node& child_a = lbvh[lbvh[parent].left]; - const Node& child_b = lbvh[lbvh[parent].right]; - lbvh[parent].aabb_min = child_a.aabb_min.min(child_b.aabb_min); - lbvh[parent].aabb_max = child_a.aabb_max.max(child_b.aabb_max); - - // Compute rightmost leaf: max of children's rightmost - rightmost_leaves[parent] = std::max( - rightmost_leaves[lbvh[parent].left], - rightmost_leaves[lbvh[parent].right]); - - // Reconstruct the full key range for the parent. - left_key = construction_infos[parent].left_range; - right_key = construction_infos[parent].right_range; - current_node = parent; - - if (left_key == 0 && right_key == N_LEAVES - 1) { - // only one thread should reach the root - int expected = -1; - [[maybe_unused]] bool set = - root_idx.compare_exchange_strong( - expected, current_node); - assert(set); - break; // root AABB is complete - } + const size_t box_id = morton_codes[i].box_id; + details::init_leaf_node( + i, N_LEAVES, box_id, boxes[box_id].min, boxes[box_id].max, + lbvh.data(), rightmost_leaves.data()); + + const int root = details::build_hierarchy_from_leaf( + i, N_LEAVES, [&](int k) { return morton_codes[k].morton_code; }, + lbvh.data(), rightmost_leaves.data(), construction_infos.data(), + // std::atomic's post-increment is sequentially consistent, so + // it already orders this thread's writes before the arrival + // and the arrival before its later reads. + [](std::atomic& count) { return count++; }); + + if (root >= 0) { + // Only one thread should ever reach the root. + int expected = -1; + [[maybe_unused]] const bool set = + root_idx.compare_exchange_strong(expected, root); + assert(set); } }); } // --- Move the root to index 0 so traversal can start there. --- // In the Apetrei layout the root's index equals the global split position, - // which is generally != 0. We swap the root node into position 0 and patch - // up the single affected child pointer. - // - // Key invariant (Apetrei): node 0's subtree always has left_key=0, so it is - // only ever written as a LEFT child — meaning no internal node ever has - // right==0. Therefore swapping node 0 cannot create a spurious - // is_inner_marker==0 (which would look like a leaf). + // which is generally != 0. const int root = root_idx.load(); if (root > 0) { IPC_TOOLKIT_PROFILE_BLOCK("swap_root_to_zero"); - std::swap(lbvh[0], lbvh[root]); - std::swap(rightmost_leaves[0], rightmost_leaves[root]); - - // The root (now at 0) is never any node's child, so no pointer - // references R that needs rewriting to 0. The only pointers that - // referenced 0 (the old node-0) must be rewritten to R. And since old - // node-0 was only ever a LEFT child (see invariant above), we only need - // to patch .left pointers. + details::swap_root_to_zero(lbvh.data(), rightmost_leaves.data(), root); + tbb::parallel_for(size_t(0), lbvh.size(), [&](size_t i) { - if (lbvh[i].is_inner() && lbvh[i].left == 0) { - lbvh[i].left = root; - } + details::patch_left_pointer(lbvh[i], root); }); } } @@ -366,6 +216,10 @@ namespace { candidates.emplace_back(i, j); } + /// Scalar traversal: descend the target BVH for one query leaf and record + /// every overlapping, filter-passing pair. The descent itself is + /// ipc::details::traverse_lbvh(), shared with the SIMD traversal below and + /// with ipc::cuda::LBVH; only the overlap test and the emission are ours. template void traverse_lbvh( const LBVH::Node& query, @@ -375,86 +229,23 @@ namespace { const std::function& can_collide, std::vector& candidates) { - // Use a fixed-size array as a stack to avoid dynamic allocations - constexpr int MAX_STACK_SIZE = 64; - int stack[MAX_STACK_SIZE]; - int stack_ptr = 0; - stack[stack_ptr++] = LBVH::Node::INVALID_POINTER; - - int node_idx = 0; // root - do { - const LBVH::Node& node = lbvh[node_idx]; - - if (lbvh.size() == 1) { // Single node case (only root) - assert(node.is_leaf()); // Only one node, so it must be a leaf - if constexpr (triangular) { - break; // No self-collision if only one node - } - if (node.intersects(query)) { - attempt_add_candidate( - query, node, can_collide, candidates); - } - break; - } - - // Check left and right are valid pointers - assert(node.is_inner()); - -#if defined(__GNUC__) || defined(__clang__) - // Prefetch child nodes to reduce cache misses - __builtin_prefetch(&lbvh[node.left], 0, 1); - __builtin_prefetch(&lbvh[node.right], 0, 1); -#endif - - const LBVH::Node& child_l = lbvh[node.left]; - const LBVH::Node& child_r = lbvh[node.right]; - bool intersects_l = child_l.intersects(query); - bool intersects_r = child_r.intersects(query); - - // Ignore overlap if the subtree is fully on the - // left-hand side of the query (triangular traversal only). - if constexpr (triangular) { - if (intersects_l - && rightmost_leaves[node.left] <= query_leaf_idx) { - intersects_l = false; - } - if (intersects_r - && rightmost_leaves[node.right] <= query_leaf_idx) { - intersects_r = false; - } - } - - // Query overlaps a leaf node => report collision. - if (intersects_l && child_l.is_leaf()) { + details::traverse_lbvh( + int(query_leaf_idx), lbvh.data(), int(lbvh.size()), + rightmost_leaves.data(), + [&](const LBVH::Node& node) { return node.intersects(query); }, + [&](const LBVH::Node& leaf, const int /*leaf_idx*/, + const bool /*intersects*/) { attempt_add_candidate( - query, child_l, can_collide, candidates); - } - if (intersects_r && child_r.is_leaf()) { - attempt_add_candidate( - query, child_r, can_collide, candidates); - } - - // Query overlaps an internal node => traverse. - bool traverse_l = (intersects_l && !child_l.is_leaf()); - bool traverse_r = (intersects_r && !child_r.is_leaf()); - - if (!traverse_l && !traverse_r) { - assert(stack_ptr > 0); - node_idx = stack[--stack_ptr]; - } else { - node_idx = traverse_l ? node.left : node.right; - if (traverse_l && traverse_r) { - // Postpone traversal of the right child - assert(stack_ptr < MAX_STACK_SIZE); - stack[stack_ptr++] = node.right; - } - } - } while (node_idx != LBVH::Node::INVALID_POINTER); // Same as root + query, leaf, can_collide, candidates); + }); } #ifdef IPC_TOOLKIT_WITH_SIMD - // SIMD Traversal - // Traverses multiple queries simultaneously using SIMD. + /// SIMD traversal: descend the target BVH once for a batch of up to + /// xs::batch::size query leaves, testing every node against all of + /// them at once. The descent is the same ipc::details::traverse_lbvh() as + /// the scalar path; the overlap test returns a lane mask instead of a bool, + /// and the emission fans a leaf out to the lanes that overlap it. template void traverse_lbvh_simd( const LBVH::Node* queries, @@ -466,6 +257,7 @@ namespace { std::vector& candidates) { using batch_t = xs::batch; + using mask_t = xs::batch_bool; assert(n_queries >= 1 && n_queries <= batch_t::size); // Load queries into single registers @@ -502,136 +294,36 @@ namespace { const auto q_max_z = make_simd([&](int k) { return queries[k].aabb_max.z(); }); - // Use a fixed-size array as a stack to avoid dynamic allocations - constexpr int MAX_STACK_SIZE = 64; - int stack[MAX_STACK_SIZE]; - int stack_ptr = 0; - stack[stack_ptr++] = LBVH::Node::INVALID_POINTER; - - int node_idx = 0; // root - do { - const LBVH::Node& node = lbvh[node_idx]; - - if (lbvh.size() == 1) { // Single node case (only root) - assert(node.is_leaf()); // Only one node, so it must be a leaf - if constexpr (triangular) { - break; // No self-collision if only one node - } - // Check intersection with all queries simultaneously - const xs::batch_bool intersects = - (node.aabb_min.x() <= q_max_x) + details::traverse_lbvh( + int(first_query_leaf_idx), lbvh.data(), int(lbvh.size()), + rightmost_leaves.data(), + // Intersect all queries at once: + // (node.min <= query.max) && (query.min <= node.max) + [&](const LBVH::Node& node) -> mask_t { + return (node.aabb_min.x() <= q_max_x) & (node.aabb_min.y() <= q_max_y) & (node.aabb_min.z() <= q_max_z) & (q_min_x <= node.aabb_max.x()) & (q_min_y <= node.aabb_max.y()) & (q_min_z <= node.aabb_max.z()); - if (xs::any(intersects)) { - for (int k = 0; k < n_queries; ++k) { - if (intersects.get(k)) { - attempt_add_candidate( - queries[k], node, can_collide, candidates); - } - } - } - break; - } - - // Check left and right are valid pointers - assert(node.is_inner()); - -#if defined(__GNUC__) || defined(__clang__) - // Prefetch child nodes to reduce cache misses - __builtin_prefetch(&lbvh[node.left], 0, 1); - __builtin_prefetch(&lbvh[node.right], 0, 1); -#endif - - const LBVH::Node& child_l = lbvh[node.left]; - const LBVH::Node& child_r = lbvh[node.right]; - - // 1. Intersect multiple queries at once - // (child_l.min <= query.max) && (query.min <= child_l.max) - xs::batch_bool intersects_l = - (child_l.aabb_min.x() <= q_max_x) - & (child_l.aabb_min.y() <= q_max_y) - & (child_l.aabb_min.z() <= q_max_z) - & (q_min_x <= child_l.aabb_max.x()) - & (q_min_y <= child_l.aabb_max.y()) - & (q_min_z <= child_l.aabb_max.z()); - - // 2. Intersect multiple queries at once - // (child_r.min <= query.max) && (query.min <= child_r.max) - xs::batch_bool intersects_r = - (child_r.aabb_min.x() <= q_max_x) - & (child_r.aabb_min.y() <= q_max_y) - & (child_r.aabb_min.z() <= q_max_z) - & (q_min_x <= child_r.aabb_max.x()) - & (q_min_y <= child_r.aabb_max.y()) - & (q_min_z <= child_r.aabb_max.z()); - - // Ignore overlap if the subtree is fully on the left-hand side - // of all queries (triangular traversal only). - // We use first_query_leaf_idx (the smallest query leaf index - // in the SIMD batch) for a conservative check: if all leaves - // in the subtree are <= the smallest query, they are also <= - // every other query in the batch. - if constexpr (triangular) { - if (rightmost_leaves[node.left] <= first_query_leaf_idx) { - intersects_l = xs::batch_bool(false); - } - if (rightmost_leaves[node.right] <= first_query_leaf_idx) { - intersects_r = xs::batch_bool(false); - } - } - - const bool any_intersects_l = xs::any(intersects_l); - const bool any_intersects_r = xs::any(intersects_r); - - // Query overlaps a leaf node => report collision - if (any_intersects_l && child_l.is_leaf()) { - for (int k = 0; k < n_queries; ++k) { + }, + [&](const LBVH::Node& leaf, const int leaf_idx, + const mask_t& intersects) { + for (size_t k = 0; k < n_queries; ++k) { if constexpr (triangular) { - if (rightmost_leaves[node.left] + // The shared descent skipped subtrees left of the + // batch's FIRST query; finish the check per lane. + if (rightmost_leaves[leaf_idx] <= first_query_leaf_idx + k) { continue; } } - if (intersects_l.get(k)) { + if (intersects.get(k)) { attempt_add_candidate( - queries[k], child_l, can_collide, candidates); + queries[k], leaf, can_collide, candidates); } } - } - if (any_intersects_r && child_r.is_leaf()) { - for (int k = 0; k < n_queries; ++k) { - if constexpr (triangular) { - if (rightmost_leaves[node.right] - <= first_query_leaf_idx + k) { - continue; - } - } - if (intersects_r.get(k)) { - attempt_add_candidate( - queries[k], child_r, can_collide, candidates); - } - } - } - - // Query overlaps an internal node => traverse. - bool traverse_l = (any_intersects_l && !child_l.is_leaf()); - bool traverse_r = (any_intersects_r && !child_r.is_leaf()); - - if (!traverse_l && !traverse_r) { - assert(stack_ptr > 0); - node_idx = stack[--stack_ptr]; - } else { - node_idx = traverse_l ? node.left : node.right; - if (traverse_l && traverse_r) { - // Postpone traversal of the right child - assert(stack_ptr < MAX_STACK_SIZE); - stack[stack_ptr++] = node.right; - } - } - } while (node_idx != LBVH::Node::INVALID_POINTER); // Same as root + }); } #endif @@ -718,6 +410,7 @@ void LBVH::detect_candidates( void LBVH::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (vertex_bvh.size() <= 1) { // Need at least 2 vertices for a collision return; } @@ -731,6 +424,7 @@ void LBVH::detect_vertex_vertex_candidates( void LBVH::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (!has_edges() || !has_vertices()) { return; } @@ -747,6 +441,7 @@ void LBVH::detect_edge_vertex_candidates( void LBVH::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); if (edge_bvh.size() <= 1) { // Need at least 2 edges for a collision return; } @@ -761,6 +456,7 @@ void LBVH::detect_edge_edge_candidates( void LBVH::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (!has_faces() || !has_vertices()) { return; } @@ -776,6 +472,7 @@ void LBVH::detect_face_vertex_candidates( void LBVH::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); if (!has_edges() || !has_faces()) { return; } @@ -791,6 +488,7 @@ void LBVH::detect_edge_face_candidates( void LBVH::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); if (face_bvh.size() <= 1) { // Need at least 2 faces for a collision return; } @@ -808,8 +506,7 @@ bool LBVH::can_edge_vertex_collide(size_t ei, size_t vi) const assert(ei < edge_vertex_ids.size()); const auto& [e0i, e1i] = edge_vertex_ids[ei]; - return vi != e0i && vi != e1i - && (can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i)); + return details::can_edge_vertex_collide(e0i, e1i, vi, can_vertices_collide); } bool LBVH::can_edges_collide(size_t eai, size_t ebi) const @@ -819,13 +516,8 @@ bool LBVH::can_edges_collide(size_t eai, size_t ebi) const assert(ebi < edge_vertex_ids.size()); const auto& [eb0i, eb1i] = edge_vertex_ids[ebi]; - const bool share_endpoint = - ea0i == eb0i || ea0i == eb1i || ea1i == eb0i || ea1i == eb1i; - - return !share_endpoint - && (can_vertices_collide(ea0i, eb0i) || can_vertices_collide(ea0i, eb1i) - || can_vertices_collide(ea1i, eb0i) - || can_vertices_collide(ea1i, eb1i)); + return details::can_edges_collide( + ea0i, ea1i, eb0i, eb1i, can_vertices_collide); } bool LBVH::can_face_vertex_collide(size_t fi, size_t vi) const @@ -833,9 +525,8 @@ bool LBVH::can_face_vertex_collide(size_t fi, size_t vi) const assert(fi < face_vertex_ids.size()); const auto& [f0i, f1i, f2i] = face_vertex_ids[fi]; - return vi != f0i && vi != f1i && vi != f2i - && (can_vertices_collide(vi, f0i) || can_vertices_collide(vi, f1i) - || can_vertices_collide(vi, f2i)); + return details::can_face_vertex_collide( + f0i, f1i, f2i, vi, can_vertices_collide); } bool LBVH::can_edge_face_collide(size_t ei, size_t fi) const @@ -845,14 +536,8 @@ bool LBVH::can_edge_face_collide(size_t ei, size_t fi) const assert(fi < face_vertex_ids.size()); const auto& [f0i, f1i, f2i] = face_vertex_ids[fi]; - const bool share_endpoint = e0i == f0i || e0i == f1i || e0i == f2i - || e1i == f0i || e1i == f1i || e1i == f2i; - - return !share_endpoint - && (can_vertices_collide(e0i, f0i) || can_vertices_collide(e0i, f1i) - || can_vertices_collide(e0i, f2i) || can_vertices_collide(e1i, f0i) - || can_vertices_collide(e1i, f1i) - || can_vertices_collide(e1i, f2i)); + return details::can_edge_face_collide( + e0i, e1i, f0i, f1i, f2i, can_vertices_collide); } bool LBVH::can_faces_collide(size_t fai, size_t fbi) const @@ -862,19 +547,8 @@ bool LBVH::can_faces_collide(size_t fai, size_t fbi) const assert(fbi < face_vertex_ids.size()); const auto& [fb0i, fb1i, fb2i] = face_vertex_ids[fbi]; - const bool share_endpoint = fa0i == fb0i || fa0i == fb1i || fa0i == fb2i - || fa1i == fb0i || fa1i == fb1i || fa1i == fb2i || fa2i == fb0i - || fa2i == fb1i || fa2i == fb2i; - - return !share_endpoint - && (can_vertices_collide(fa0i, fb0i) || can_vertices_collide(fa0i, fb1i) - || can_vertices_collide(fa0i, fb2i) - || can_vertices_collide(fa1i, fb0i) - || can_vertices_collide(fa1i, fb1i) - || can_vertices_collide(fa1i, fb2i) - || can_vertices_collide(fa2i, fb0i) - || can_vertices_collide(fa2i, fb1i) - || can_vertices_collide(fa2i, fb2i)); + return details::can_faces_collide( + fa0i, fa1i, fa2i, fb0i, fb1i, fb2i, can_vertices_collide); } } // namespace ipc \ No newline at end of file diff --git a/src/ipc/broad_phase/lbvh.hpp b/src/ipc/broad_phase/lbvh.hpp index 64f37de77..4274e61d8 100644 --- a/src/ipc/broad_phase/lbvh.hpp +++ b/src/ipc/broad_phase/lbvh.hpp @@ -1,5 +1,6 @@ #pragma once +#include // for IPC_TOOLKIT_HOST_DEVICE #include #include @@ -48,14 +49,20 @@ class LBVH : public BroadPhase { #pragma GCC diagnostic pop + // These are host/device so the CUDA broad phase (ipc::cuda::LBVH) + // traverses the same node with the same predicates as the CPU one. + /// @brief Check if this node is an inner node. - bool is_inner() const { return is_inner_marker; } + IPC_TOOLKIT_HOST_DEVICE bool is_inner() const + { + return is_inner_marker; + } /// @brief Check if this node is a leaf node. - bool is_leaf() const { return !is_inner(); } + IPC_TOOLKIT_HOST_DEVICE bool is_leaf() const { return !is_inner(); } /// @brief Check if this node is valid. - bool is_valid() const + IPC_TOOLKIT_HOST_DEVICE bool is_valid() const { return is_inner() ? (left != INVALID_POINTER && right != INVALID_POINTER) @@ -63,7 +70,7 @@ class LBVH : public BroadPhase { } /// @brief Check if this node's AABB intersects with another node's AABB. - bool intersects(const Node& other) const + IPC_TOOLKIT_HOST_DEVICE bool intersects(const Node& other) const { return (aabb_min <= other.aabb_max).all() && (other.aabb_min <= aabb_max).all(); @@ -84,18 +91,26 @@ class LBVH : public BroadPhase { /// Used to skip subtrees during triangular (self-collision) traversal. using RightmostLeaves = std::vector>; -private: - struct ConstructionInfo { + /// @brief Per-internal-node scratch for the bottom-up build. + /// @tparam Counter The visitation counter's type: std::atomic here, + /// and a plain int for the device build in ipc::cuda::LBVH, where + /// atomicAdd() supplies the atomicity. + /// @see ipc::details::build_hierarchy_from_leaf + template struct ConstructionInfo { /// @brief Left range endpoint passed up by the left child. - int32_t left_range; + int left_range; /// @brief Right range endpoint passed up by the right child. - int32_t right_range; + int right_range; /// @brief Number of threads that arrived at this node. - std::atomic visitation_count; + Counter visitation_count; }; - using ConstructionInfos = - std::vector>; +private: + using HostConstructionInfo = ConstructionInfo>; + + using ConstructionInfos = std::vector< + HostConstructionInfo, + DefaultInitAllocator>; public: LBVH(); diff --git a/src/ipc/broad_phase/spatial_hash.cpp b/src/ipc/broad_phase/spatial_hash.cpp index 6bd45a00e..4d66c597a 100644 --- a/src/ipc/broad_phase/spatial_hash.cpp +++ b/src/ipc/broad_phase/spatial_hash.cpp @@ -236,6 +236,7 @@ namespace { void SpatialHash::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (vertex_boxes.empty()) { return; } @@ -251,6 +252,7 @@ void SpatialHash::detect_vertex_vertex_candidates( void SpatialHash::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (edge_boxes.empty() || vertex_boxes.empty()) { return; } @@ -267,6 +269,7 @@ void SpatialHash::detect_edge_vertex_candidates( void SpatialHash::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); if (edge_boxes.empty()) { return; } @@ -282,6 +285,7 @@ void SpatialHash::detect_edge_edge_candidates( void SpatialHash::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); if (face_boxes.empty() || vertex_boxes.empty()) { return; } @@ -299,6 +303,7 @@ void SpatialHash::detect_face_vertex_candidates( void SpatialHash::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); if (edge_boxes.empty() || face_boxes.empty()) { return; } @@ -315,6 +320,7 @@ void SpatialHash::detect_edge_face_candidates( void SpatialHash::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); if (face_boxes.empty()) { return; } diff --git a/src/ipc/broad_phase/sweep_and_prune.cpp b/src/ipc/broad_phase/sweep_and_prune.cpp index eeb78cfb0..1088c1055 100644 --- a/src/ipc/broad_phase/sweep_and_prune.cpp +++ b/src/ipc/broad_phase/sweep_and_prune.cpp @@ -102,6 +102,7 @@ void SweepAndPrune::clear() void SweepAndPrune::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep(boxes->vertices, vv_sort_axis, overlaps); @@ -115,6 +116,7 @@ void SweepAndPrune::detect_vertex_vertex_candidates( void SweepAndPrune::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep( boxes->edges, boxes->vertices, ev_sort_axis, overlaps); @@ -129,6 +131,7 @@ void SweepAndPrune::detect_edge_vertex_candidates( void SweepAndPrune::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep(boxes->edges, ee_sort_axis, overlaps); @@ -142,6 +145,7 @@ void SweepAndPrune::detect_edge_edge_candidates( void SweepAndPrune::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep( boxes->faces, boxes->vertices, fv_sort_axis, overlaps); @@ -156,6 +160,7 @@ void SweepAndPrune::detect_face_vertex_candidates( void SweepAndPrune::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep( boxes->edges, boxes->faces, ef_sort_axis, overlaps); @@ -170,6 +175,7 @@ void SweepAndPrune::detect_edge_face_candidates( void SweepAndPrune::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); std::vector> overlaps; scalable_ccd::sort_and_sweep(boxes->faces, ff_sort_axis, overlaps); diff --git a/src/ipc/broad_phase/sweep_and_tiniest_queue.cu b/src/ipc/broad_phase/sweep_and_tiniest_queue.cu index bf68f3d93..78b21d4b5 100644 --- a/src/ipc/broad_phase/sweep_and_tiniest_queue.cu +++ b/src/ipc/broad_phase/sweep_and_tiniest_queue.cu @@ -2,6 +2,8 @@ #ifdef IPC_TOOLKIT_WITH_CUDA +#include + #include namespace ipc { @@ -71,7 +73,7 @@ void SweepAndTiniestQueue::build( } void SweepAndTiniestQueue::build( - const AABBs& vertex_boxes, + const AABBs& _vertex_boxes, // not the (unused) inherited member Eigen::ConstRef edges, Eigen::ConstRef faces, const uint8_t _dim) @@ -84,20 +86,20 @@ void SweepAndTiniestQueue::build( dim = _dim; // Convert from ipc::AABB to scalable_ccd::cuda::AABB - boxes->vertices.resize(vertex_boxes.size()); - for (int i = 0; i < vertex_boxes.size(); ++i) { - boxes->vertices[i].min.x = vertex_boxes[i].min.x(); - boxes->vertices[i].min.y = vertex_boxes[i].min.y(); - boxes->vertices[i].min.z = vertex_boxes[i].min.z(); + boxes->vertices.resize(_vertex_boxes.size()); + for (int i = 0; i < _vertex_boxes.size(); ++i) { + boxes->vertices[i].min.x = _vertex_boxes[i].min.x(); + boxes->vertices[i].min.y = _vertex_boxes[i].min.y(); + boxes->vertices[i].min.z = _vertex_boxes[i].min.z(); - boxes->vertices[i].max.x = vertex_boxes[i].max.x(); - boxes->vertices[i].max.y = vertex_boxes[i].max.y(); - boxes->vertices[i].max.z = vertex_boxes[i].max.z(); + boxes->vertices[i].max.x = _vertex_boxes[i].max.x(); + boxes->vertices[i].max.y = _vertex_boxes[i].max.y(); + boxes->vertices[i].max.z = _vertex_boxes[i].max.z(); // If vertex id == -1 it means this slot is not used. // But Scalable CCD does not have this kind of special value so we map // it to unique negative id. - const auto [vi, vj, vk] = vertex_boxes[i].vertex_ids; + const auto [vi, vj, vk] = _vertex_boxes[i].vertex_ids; assert(vi >= 0); boxes->vertices[i].vertex_ids.x = vi; boxes->vertices[i].vertex_ids.y = vj >= 0 ? vj : (-vi - 1); @@ -121,6 +123,7 @@ void SweepAndTiniestQueue::clear() void SweepAndTiniestQueue::detect_vertex_vertex_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_vertex_boxes broad_phase.build( @@ -136,6 +139,7 @@ void SweepAndTiniestQueue::detect_vertex_vertex_candidates( void SweepAndTiniestQueue::detect_edge_vertex_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_vertex_boxes and d_edge_boxes broad_phase.build( @@ -152,6 +156,7 @@ void SweepAndTiniestQueue::detect_edge_vertex_candidates( void SweepAndTiniestQueue::detect_edge_edge_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_edge_boxes broad_phase.build( @@ -167,6 +172,7 @@ void SweepAndTiniestQueue::detect_edge_edge_candidates( void SweepAndTiniestQueue::detect_face_vertex_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_vertex_boxes and d_face_boxes broad_phase.build( @@ -183,6 +189,7 @@ void SweepAndTiniestQueue::detect_face_vertex_candidates( void SweepAndTiniestQueue::detect_edge_face_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_face_boxes and d_edge_boxes broad_phase.build( @@ -199,6 +206,7 @@ void SweepAndTiniestQueue::detect_edge_face_candidates( void SweepAndTiniestQueue::detect_face_face_candidates( std::vector& candidates) const { + candidates.clear(); scalable_ccd::cuda::BroadPhase broad_phase; // TODO: Precompute d_face_boxes broad_phase.build( @@ -213,70 +221,66 @@ void SweepAndTiniestQueue::detect_face_face_candidates( // ---------------------------------------------------------------------------- +// Scalable CCD already excludes pairs that share a vertex, so these apply only +// the user-filter half of ipc::details' connectivity rule (and assert the +// other half held). + bool SweepAndTiniestQueue::can_edge_vertex_collide(size_t ei, size_t vi) const { const auto& [e0i, e1i, _] = boxes->edges[ei].vertex_ids; + const index_t e[2] = { e0i, e1i }; + const index_t v[1] = { static_cast(vi) }; - // Checked by scalable_ccd - assert(vi != e0i && vi != e1i); + assert(!details::share_vertex(e, v)); // Checked by scalable_ccd - return can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i); + return details::any_vertex_pair_can_collide(v, e, can_vertices_collide); } bool SweepAndTiniestQueue::can_edges_collide(size_t eai, size_t ebi) const { const auto& [ea0i, ea1i, _] = boxes->edges[eai].vertex_ids; const auto& [eb0i, eb1i, __] = boxes->edges[ebi].vertex_ids; + const index_t ea[2] = { ea0i, ea1i }; + const index_t eb[2] = { eb0i, eb1i }; - // Checked by scalable_ccd - assert(ea0i != eb0i && ea0i != eb1i && ea1i != eb0i && ea1i != eb1i); + assert(!details::share_vertex(ea, eb)); // Checked by scalable_ccd - return can_vertices_collide(ea0i, eb0i) || can_vertices_collide(ea0i, eb1i) - || can_vertices_collide(ea1i, eb0i) || can_vertices_collide(ea1i, eb1i); + return details::any_vertex_pair_can_collide(ea, eb, can_vertices_collide); } bool SweepAndTiniestQueue::can_face_vertex_collide(size_t fi, size_t vi) const { const auto& [f0i, f1i, f2i] = boxes->faces[fi].vertex_ids; + const index_t f[3] = { f0i, f1i, f2i }; + const index_t v[1] = { static_cast(vi) }; - // Checked by scalable_ccd - assert(vi != f0i && vi != f1i && vi != f2i); + assert(!details::share_vertex(f, v)); // Checked by scalable_ccd - return can_vertices_collide(vi, f0i) || can_vertices_collide(vi, f1i) - || can_vertices_collide(vi, f2i); + return details::any_vertex_pair_can_collide(v, f, can_vertices_collide); } bool SweepAndTiniestQueue::can_edge_face_collide(size_t ei, size_t fi) const { const auto& [e0i, e1i, _] = boxes->edges[ei].vertex_ids; const auto& [f0i, f1i, f2i] = boxes->faces[fi].vertex_ids; + const index_t e[2] = { e0i, e1i }; + const index_t f[3] = { f0i, f1i, f2i }; - // Checked by scalable_ccd - assert( - e0i != f0i && e0i != f1i && e0i != f2i && e1i != f0i && e1i != f1i - && e1i != f2i); + assert(!details::share_vertex(e, f)); // Checked by scalable_ccd - return can_vertices_collide(e0i, f0i) || can_vertices_collide(e0i, f1i) - || can_vertices_collide(e0i, f2i) || can_vertices_collide(e1i, f0i) - || can_vertices_collide(e1i, f1i) || can_vertices_collide(e1i, f2i); + return details::any_vertex_pair_can_collide(e, f, can_vertices_collide); } bool SweepAndTiniestQueue::can_faces_collide(size_t fai, size_t fbi) const { const auto& [fa0i, fa1i, fa2i] = boxes->faces[fai].vertex_ids; const auto& [fb0i, fb1i, fb2i] = boxes->faces[fbi].vertex_ids; + const index_t fa[3] = { fa0i, fa1i, fa2i }; + const index_t fb[3] = { fb0i, fb1i, fb2i }; + + assert(!details::share_vertex(fa, fb)); // Checked by scalable_ccd - // Checked by scalable_ccd - assert( - fa0i != fb0i && fa0i != fb1i && fa0i != fb2i && fa1i != fb0i - && fa1i != fb1i && fa1i != fb2i && fa2i != fb0i && fa2i != fb1i - && fa2i != fb2i); - - return can_vertices_collide(fa0i, fb0i) || can_vertices_collide(fa0i, fb1i) - || can_vertices_collide(fa0i, fb2i) || can_vertices_collide(fa1i, fb0i) - || can_vertices_collide(fa1i, fb1i) || can_vertices_collide(fa1i, fb2i) - || can_vertices_collide(fa2i, fb0i) || can_vertices_collide(fa2i, fb1i) - || can_vertices_collide(fa2i, fb2i); + return details::any_vertex_pair_can_collide(fa, fb, can_vertices_collide); } } // namespace ipc diff --git a/src/ipc/collision_filter.hpp b/src/ipc/collision_filter.hpp index 16fe1cd13..a4d65d60f 100644 --- a/src/ipc/collision_filter.hpp +++ b/src/ipc/collision_filter.hpp @@ -25,6 +25,12 @@ namespace ipc { // operator| → union (true if EITHER filter passes) // operator& → intersection (true if BOTH filters pass) // operator! → negation +// +// The accept-all filter is represented by an EMPTY std::function rather than by +// a callable that returns true, so accepts_all() is a property of the filter's +// state (nothing can be desynchronized from it) and the common no-filter path +// costs a null check instead of an indirect call per pair. Composition +// short-circuits on it: `f | accept_all` is accept_all, `f & accept_all` is f. // ───────────────────────────────────────────────────────────────────────────── class CollisionFilter { @@ -32,11 +38,12 @@ class CollisionFilter { // ── Construction ───────────────────────────────────────────────────────── /// @brief Default filter: accept all pairs. - CollisionFilter() : m_fn([](size_t, size_t) { return true; }) { } + CollisionFilter() = default; /// @brief Construct from any callable bool(size_t, size_t). /// @note Disabled when Fn is CollisionFilter itself to avoid shadowing /// the copy constructor. + /// @note An empty std::function is the accept-all filter. template < typename Fn, typename = std::enable_if_t< @@ -52,18 +59,42 @@ class CollisionFilter { /// @param vi Index of the first vertex. /// @param vj Index of the second vertex. /// @return true if the pair should be considered for collision. - bool operator()(size_t vi, size_t vj) const { return m_fn(vi, vj); } + bool operator()(size_t vi, size_t vj) const + { + return !m_fn || m_fn(vi, vj); + } // ── Implicit conversion ────────────────────────────────────────────────── /// @brief Implicit conversion to std::function. - operator std::function() const { return m_fn; } + /// @note Always returns a callable function, even for the accept-all + /// filter (whose stored function is empty). + operator std::function() const + { + if (accepts_all()) { + return [](size_t, size_t) { return true; }; + } + return m_fn; + } + + /// @brief Whether this filter accepts every pair. + /// @return true for the default filter, for one constructed from an empty + /// std::function, and for any composition that reduces to one + /// (e.g. the union of two accept-all filters); false for any filter + /// holding a user-supplied callable, even one that happens to + /// return true for every pair. + /// @note Used by GPU broad phases to skip host-side filtering entirely when + /// the device-emitted (connectivity-filtered) set is already exact. + bool accepts_all() const { return !m_fn; } // ── Composition ────────────────────────────────────────────────────────── /// @brief Union: accept if EITHER filter passes. friend CollisionFilter operator|(CollisionFilter lhs, CollisionFilter rhs) { + if (lhs.accepts_all() || rhs.accepts_all()) { + return CollisionFilter(); // accept-all absorbs the union + } return CollisionFilter([l = std::move(lhs.m_fn), r = std::move(rhs.m_fn)](size_t vi, size_t vj) { return l(vi, vj) || r(vi, vj); @@ -73,6 +104,12 @@ class CollisionFilter { /// @brief Intersection: accept only if BOTH filters pass. friend CollisionFilter operator&(CollisionFilter lhs, CollisionFilter rhs) { + if (lhs.accepts_all()) { + return rhs; // accept-all is the identity of the intersection + } + if (rhs.accepts_all()) { + return lhs; + } return CollisionFilter([l = std::move(lhs.m_fn), r = std::move(rhs.m_fn)](size_t vi, size_t vj) { return l(vi, vj) && r(vi, vj); @@ -82,6 +119,9 @@ class CollisionFilter { /// @brief Negation: accept only if this filter rejects. CollisionFilter operator!() const { + if (accepts_all()) { + return CollisionFilter([](size_t, size_t) { return false; }); + } return CollisionFilter( [f = m_fn](size_t vi, size_t vj) { return !f(vi, vj); }); } @@ -99,6 +139,7 @@ class CollisionFilter { } private: + /// @brief The predicate; empty for the accept-all filter. std::function m_fn; }; diff --git a/src/ipc/math/morton.hpp b/src/ipc/math/morton.hpp index d43c06bac..12928d15c 100644 --- a/src/ipc/math/morton.hpp +++ b/src/ipc/math/morton.hpp @@ -3,8 +3,15 @@ #include // for IPC_TOOLKIT_HOST_DEVICE #include // for clamp +#include + #include // for uint64_t +#if !defined(__CUDA_ARCH__) && !defined(__GNUC__) && !defined(__clang__) \ + && defined(_MSC_VER) +#include // for _BitScanReverse / _BitScanReverse64 +#endif + namespace ipc { /// @brief Expands a 32-bit integer into 64 bits by inserting 1 zero after each bit. @@ -64,4 +71,151 @@ IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_3D(double x, double y, double z) return (xx << 2) | (yy << 1) | zz; } +/// @brief Computes the reciprocal width of a Morton normalization domain. +/// +/// The one place this derivation lives, so the host and device LBVH builds -- +/// which must produce bit-identical codes -- cannot drift apart on it. +/// +/// A degenerate axis -- every box spanning the same coordinate, as for a +/// planar mesh with no inflation, the unused z of a 2D mesh, or a single box +/// -- has zero width. Its reciprocal is 0 rather than infinity, so every +/// center normalizes to 0 along it and the codes simply carry no information +/// on that axis (as they should); 1/0 would give infinity and then 0 * inf = +/// NaN in morton_code(), whose conversion to an integer is undefined. +/// +/// @param domain_min The minimum corner of the normalization domain. +/// @param domain_max The maximum corner of the normalization domain. +/// @return The reciprocal of the domain's width along each axis, or 0 along a +/// zero-width axis. +IPC_TOOLKIT_HOST_DEVICE inline Eigen::Array3d morton_domain_width_inv( + const Eigen::Array3d& domain_min, const Eigen::Array3d& domain_max) +{ + Eigen::Array3d width_inv; + for (int k = 0; k < 3; ++k) { + const double width = domain_max[k] - domain_min[k]; + width_inv[k] = width > 0 ? 1.0 / width : 0.0; + } + return width_inv; +} + +/// @brief Calculates the Morton code of a box from its center. +/// +/// The center is normalized into the unit square/cube by the given domain +/// before being encoded. The domain's width is passed as a reciprocal (see +/// morton_domain_width_inv()) so this multiplies rather than divides, letting +/// the host and device LBVH builds agree bit-for-bit. +/// +/// @param center The center of the box. +/// @param domain_min The minimum corner of the normalization domain. +/// @param domain_width_inv The reciprocal of the normalization domain's width. +/// @param dim The dimension of the simulation (2 or 3). +/// @return The Morton code of the normalized center. +IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_code( + const Eigen::Array3d& center, + const Eigen::Array3d& domain_min, + const Eigen::Array3d& domain_width_inv, + const int dim) +{ + const double x = (center.x() - domain_min.x()) * domain_width_inv.x(); + const double y = (center.y() - domain_min.y()) * domain_width_inv.y(); + if (dim == 2) { + return morton_2D(x, y); + } + const double z = (center.z() - domain_min.z()) * domain_width_inv.z(); + return morton_3D(x, y, z); +} + +/// @brief Counts the leading zero bits of a 32-bit value. +/// @note Undefined for v == 0, matching the underlying intrinsics. +/// @param v The value to count the leading zeros of. +/// @return The number of leading zero bits. +IPC_TOOLKIT_HOST_DEVICE inline int count_leading_zeros(const uint32_t v) +{ +#ifdef __CUDA_ARCH__ + return __clz(static_cast(v)); +#elif defined(__GNUC__) || defined(__clang__) + return __builtin_clz(v); +#elif defined(_MSC_VER) + // Not __lzcnt: without LZCNT/ABM support that encoding decodes as bsr and + // returns the index of the highest set bit instead. _BitScanReverse is the + // bsr itself, so the count is 31 minus the index on every x86/ARM target. + unsigned long index; // NOLINT(google-runtime-int) + _BitScanReverse(&index, v); + return 31 - static_cast(index); +#else +#error "count_leading_zeros: no leading-zero-count intrinsic for this compiler" +#endif +} + +/// @brief Counts the leading zero bits of a 64-bit value. +/// @note Undefined for v == 0, matching the underlying intrinsics. +/// @param v The value to count the leading zeros of. +/// @return The number of leading zero bits. +IPC_TOOLKIT_HOST_DEVICE inline int count_leading_zeros(const uint64_t v) +{ +#ifdef __CUDA_ARCH__ + return __clzll(static_cast(v)); +#elif defined(__GNUC__) || defined(__clang__) + return __builtin_clzll(v); +#elif defined(_MSC_VER) +#if defined(_M_X64) || defined(_M_ARM64) + unsigned long index; // NOLINT(google-runtime-int) + _BitScanReverse64(&index, v); + return 63 - static_cast(index); +#else + // _BitScanReverse64 does not exist on 32-bit targets: scan the halves. + const uint32_t hi = static_cast(v >> 32); + return hi != 0 ? count_leading_zeros(hi) + : 32 + count_leading_zeros(static_cast(v)); +#endif +#else +#error "count_leading_zeros: no leading-zero-count intrinsic for this compiler" +#endif +} + +/// @brief The number of bits in the Morton code (the width of its type). +constexpr int MORTON_CODE_BITS = 8 * sizeof(uint64_t); + +/// @brief The number of bits in the position used to break ties between +/// duplicate Morton codes (see morton_common_prefix()): the width of the +/// sorted-position type. +constexpr int MORTON_TIE_BREAK_BITS = 8 * sizeof(int); + +/// @brief The width of the augmented key morton_common_prefix() compares: the +/// Morton code followed by the sorted position. +/// +/// Every internal node of an LBVH splits its range at a distinct prefix length +/// of this key, and prefix lengths strictly increase from the root down, so no +/// root-to-leaf path has more than this many internal nodes. This bounds the +/// depth of the tree and so sizes the traversal stack. +constexpr int MORTON_KEY_BITS = MORTON_CODE_BITS + MORTON_TIE_BREAK_BITS; + +/// @brief Computes the length of the common leading-bit prefix of two sorted +/// Morton codes. +/// +/// This is the delta of Apetrei [2014]: a larger value means the two positions +/// are separated by a finer split, and so have a nearer common ancestor. +/// Duplicate codes fall back to the leading zeros of the positions' XOR, offset +/// by the code width so that any code-level difference always compares as the +/// shorter prefix -- as if the position were appended to the code (Karras +/// [2012]). +/// +/// @note The two positions must differ (i != j). This holds for every delta the +/// LBVH build evaluates, as it only ever compares adjacent positions. +/// +/// @param code_i The Morton code at sorted position i. +/// @param i The first sorted position. +/// @param code_j The Morton code at sorted position j. +/// @param j The second sorted position. +/// @return The length of the common leading-bit prefix, in [0, MORTON_KEY_BITS). +IPC_TOOLKIT_HOST_DEVICE inline int morton_common_prefix( + const uint64_t code_i, const int i, const uint64_t code_j, const int j) +{ + if (code_i == code_j) { + return MORTON_CODE_BITS + + count_leading_zeros(static_cast(i ^ j)); + } + return count_leading_zeros(code_i ^ code_j); +} + } // namespace ipc \ No newline at end of file diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt index ed1a23faf..0660e56a5 100644 --- a/src/ipc/utils/CMakeLists.txt +++ b/src/ipc/utils/CMakeLists.txt @@ -27,3 +27,7 @@ set(SOURCES ) target_sources(ipc_toolkit PRIVATE ${SOURCES}) + +if(IPC_TOOLKIT_WITH_CUDA) + add_subdirectory(cuda) +endif() diff --git a/src/ipc/utils/cuda/CMakeLists.txt b/src/ipc/utils/cuda/CMakeLists.txt new file mode 100644 index 000000000..26eb8ff94 --- /dev/null +++ b/src/ipc/utils/cuda/CMakeLists.txt @@ -0,0 +1,6 @@ +set(SOURCES + device_buffer.cuh + device_utils.cuh +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) diff --git a/src/ipc/utils/cuda/device_buffer.cuh b/src/ipc/utils/cuda/device_buffer.cuh new file mode 100644 index 000000000..ee42fee70 --- /dev/null +++ b/src/ipc/utils/cuda/device_buffer.cuh @@ -0,0 +1,146 @@ +// Raw device memory for the ipc::cuda implementation files. +// This header is CUDA-only and must be included from .cu files exclusively. + +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include + +#include + +#include +#include + +namespace ipc::cuda { + +/// @brief A device array that is never initialized, only ever grows, and is +/// released without throwing. +/// +/// This replaces thrust::device_vector for per-frame scratch and output +/// buffers, for two reasons: +/// +/// 1. device_vector::resize() value-initializes, launching a fill over every +/// appended element. For a buffer the next kernel writes from scratch that +/// is a full pass of wasted memory traffic, and for a candidate buffer sized +/// to its high-water mark it is the dominant per-call cost. +/// 2. device_vector's destructor throws if cudaFree fails. After a sticky +/// kernel fault every CUDA call fails, so unwinding past a live +/// device_vector calls std::terminate and the caller never sees the error +/// IPC_TOOLKIT_CUDA_CHECK raised. This buffer ignores cudaFree's result. +/// +/// The capacity is a high-water mark: clear() and a smaller resize() keep the +/// allocation, so a per-frame rebuild of the same mesh allocates nothing. +/// Reallocating does NOT preserve the contents; every user here refills from +/// scratch. +template class DeviceBuffer { +public: + DeviceBuffer() = default; + ~DeviceBuffer() { release(); } + + DeviceBuffer(const DeviceBuffer&) = delete; + DeviceBuffer& operator=(const DeviceBuffer&) = delete; + + DeviceBuffer(DeviceBuffer&& other) noexcept + : m_data(std::exchange(other.m_data, nullptr)) + , m_size(std::exchange(other.m_size, size_t(0))) + , m_capacity(std::exchange(other.m_capacity, size_t(0))) + { + } + + DeviceBuffer& operator=(DeviceBuffer&& other) noexcept + { + if (this != &other) { + release(); + m_data = std::exchange(other.m_data, nullptr); + m_size = std::exchange(other.m_size, size_t(0)); + m_capacity = std::exchange(other.m_capacity, size_t(0)); + } + return *this; + } + + /// @brief Ensure room for @p n elements. + /// @note Reallocating discards the contents; nothing is initialized. + void reserve(const size_t n) + { + if (n <= m_capacity) { + return; + } + release(); + IPC_TOOLKIT_CUDA_CHECK( + cudaMalloc(reinterpret_cast(&m_data), n * sizeof(T))); + m_capacity = n; + } + + /// @brief Set the size to @p n, reserving as needed (see reserve()). + void resize(const size_t n) + { + reserve(n); + m_size = n; + } + + /// @brief Set the size to zero, keeping the allocation. + void clear() { m_size = 0; } + + /// @brief Free the allocation. + /// @note Never throws: after a sticky device error cudaFree fails too, and + /// a destructor cannot propagate that without terminating the program. + void release() noexcept + { + if (m_data != nullptr) { + static_cast(cudaFree(m_data)); + m_data = nullptr; + } + m_size = 0; + m_capacity = 0; + } + + /// @brief Set every byte of the first size() elements to @p byte + /// (asynchronous, on the default stream). + void fill_bytes(const int byte) + { + if (m_size > 0) { + IPC_TOOLKIT_CUDA_CHECK( + cudaMemsetAsync(m_data, byte, m_size * sizeof(T))); + } + } + + /// @brief Zero the first size() elements (asynchronous). + void zero() { fill_bytes(0); } + + /// @brief Resize to @p n and copy @p n elements from host memory. + void upload(const T* host, const size_t n) + { + resize(n); + if (n > 0) { + IPC_TOOLKIT_CUDA_CHECK(cudaMemcpy( + m_data, host, n * sizeof(T), cudaMemcpyHostToDevice)); + } + } + + /// @brief Copy the first size() elements to host memory (synchronous). + void download(T* host) const + { + if (m_size > 0) { + IPC_TOOLKIT_CUDA_CHECK(cudaMemcpy( + host, m_data, m_size * sizeof(T), cudaMemcpyDeviceToHost)); + } + } + + T* data() { return m_data; } + const T* data() const { return m_data; } + size_t size() const { return m_size; } + size_t capacity() const { return m_capacity; } + bool empty() const { return m_size == 0; } + +private: + T* m_data = nullptr; + size_t m_size = 0; + size_t m_capacity = 0; +}; + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/utils/cuda/device_utils.cuh b/src/ipc/utils/cuda/device_utils.cuh new file mode 100644 index 000000000..ea84b83f1 --- /dev/null +++ b/src/ipc/utils/cuda/device_utils.cuh @@ -0,0 +1,39 @@ +// Device-side utilities shared by the ipc::cuda implementation files. +// This header is CUDA-only and must be included from .cu files exclusively. + +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include + +#include + +/// @brief Log and throw (ipc::log_and_throw_error) if a CUDA runtime call +/// fails. +#define IPC_TOOLKIT_CUDA_CHECK(expr) \ + do { \ + const cudaError_t ipc_cuda_check_err = (expr); \ + if (ipc_cuda_check_err != cudaSuccess) { \ + ::ipc::log_and_throw_error( \ + "CUDA error at {}:{}: {}", __FILE__, __LINE__, \ + cudaGetErrorString(ipc_cuda_check_err)); \ + } \ + } while (false) + +namespace ipc::cuda { + +/// @brief Number of threads per block used by the ipc::cuda kernels. +constexpr int KERNEL_BLOCK_SIZE = 256; + +/// @brief Compute the launch grid size for @p n threads. +inline int kernel_grid_size(const size_t n) +{ + return static_cast((n + KERNEL_BLOCK_SIZE - 1) / KERNEL_BLOCK_SIZE); +} + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/src/ipc/utils/simd.hpp b/src/ipc/utils/simd.hpp index 72c099635..970f5b3e6 100644 --- a/src/ipc/utils/simd.hpp +++ b/src/ipc/utils/simd.hpp @@ -66,6 +66,12 @@ IPC_TOOLKIT_HOST_DEVICE inline T select(const bool mask, const T& a, const T& b) return mask ? a : b; } +/// @brief Whether any lane of a mask is set. +/// +/// The scalar counterpart of `xsimd::any`, which ADL finds for a batch `mask`, +/// so one `any(mask)` compiles for both (see ipc::details::traverse_lbvh). +IPC_TOOLKIT_HOST_DEVICE constexpr bool any(const bool mask) { return mask; } + /// @brief Clamp `v` to `[lo, hi]`. /// /// Not `std::clamp`, which cannot be called from device code: MSVC's debug STL diff --git a/tests/src/tests/broad_phase/CMakeLists.txt b/tests/src/tests/broad_phase/CMakeLists.txt index 1806f2a49..572cf94ac 100644 --- a/tests/src/tests/broad_phase/CMakeLists.txt +++ b/tests/src/tests/broad_phase/CMakeLists.txt @@ -13,8 +13,15 @@ set(SOURCES # Utilities brute_force_comparison.cpp brute_force_comparison.hpp + lbvh_validation.hpp ) +if(IPC_TOOLKIT_WITH_CUDA) + list(APPEND SOURCES + test_gpu_lbvh.cu + ) +endif() + target_sources(ipc_toolkit_tests PRIVATE ${SOURCES}) ################################################################################ diff --git a/tests/src/tests/broad_phase/lbvh_validation.hpp b/tests/src/tests/broad_phase/lbvh_validation.hpp new file mode 100644 index 000000000..28ddeb5f8 --- /dev/null +++ b/tests/src/tests/broad_phase/lbvh_validation.hpp @@ -0,0 +1,109 @@ +#pragma once + +// Structural validators for ipc::LBVH node arrays, shared by the CPU tests +// (test_lbvh.cpp) and the GPU parity tests (test_gpu_lbvh.cu) so both trees +// are held to the same predicate. Host-only: no CUDA here. + +#include + +#include + +#include +#include +#include + +namespace ipc::tests { + +/// @brief Whether a parent's AABB is exactly the union of its children's. +/// +/// Exact, not approximate: the build computes the parent as the component-wise +/// float min/max of its children on the host and the device alike, so anything +/// short of bit equality means a parent was combined from the wrong children or +/// was written before both had arrived. +inline bool is_aabb_union( + const LBVH::Node& parent, + const LBVH::Node& child_a, + const LBVH::Node& child_b) +{ + return (parent.aabb_min == child_a.aabb_min.min(child_b.aabb_min)).all() + && (parent.aabb_max == child_a.aabb_max.max(child_b.aabb_max)).all(); +} + +/// @brief Recursively verify that every node below @p index is reached exactly +/// once and that every internal node's AABB is the union of its children's, +/// collecting the primitive ids of the leaves reached. +/// +/// A revisit is fatal (REQUIRE), not just a failure: a shared child would +/// otherwise be re-traversed exponentially, and a cycle forever. +inline void traverse_lbvh_nodes( + const LBVH::Nodes& nodes, + const int32_t index, + std::vector& visited, + std::vector& reached_leaves) +{ + REQUIRE(index >= 0); + REQUIRE(index < int32_t(nodes.size())); + const LBVH::Node& node = nodes[index]; + CHECK(node.is_valid()); + REQUIRE(!visited[index]); + visited[index] = true; + + if (node.is_leaf()) { + reached_leaves.push_back(node.primitive_id); + return; + } + + const LBVH::Node& child_a = nodes[node.left]; + const LBVH::Node& child_b = nodes[node.right]; + { + CAPTURE( + index, node.left, node.right, node.aabb_min.transpose(), + child_a.aabb_min.transpose(), child_b.aabb_min.transpose(), + node.aabb_max.transpose(), child_a.aabb_max.transpose(), + child_b.aabb_max.transpose()); + CHECK(is_aabb_union(node, child_a, child_b)); + } + traverse_lbvh_nodes(nodes, node.left, visited, reached_leaves); + traverse_lbvh_nodes(nodes, node.right, visited, reached_leaves); +} + +/// @brief Validate a built LBVH: 2n - 1 nodes for n leaves, every node +/// reachable exactly once from the root at index 0, every internal AABB the +/// union of its children's, and the leaves holding exactly the primitive ids +/// {0, ..., n - 1}. A single-node tree (one primitive) is valid. +inline void check_valid_lbvh_nodes(const LBVH::Nodes& nodes) +{ + REQUIRE(!nodes.empty()); + REQUIRE(nodes.size() % 2 == 1); + const size_t n_leaves = (nodes.size() + 1) / 2; + + std::vector visited(nodes.size(), false); + std::vector reached_leaves; + traverse_lbvh_nodes(nodes, 0, visited, reached_leaves); + CHECK( + std::all_of(visited.begin(), visited.end(), [](bool v) { return v; })); + + REQUIRE(reached_leaves.size() == n_leaves); + std::sort(reached_leaves.begin(), reached_leaves.end()); + for (size_t i = 0; i < reached_leaves.size(); ++i) { + CHECK(reached_leaves[i] == int32_t(i)); + } +} + +/// @brief Validate a tree built by one implementation against the same tree +/// built by another from the same boxes: it must be a valid LBVH of the same +/// size, and its root AABB -- an order-independent union of identically +/// inflated boxes -- must be exactly equal. +/// +/// The size check comes first, so an implementation that returns no nodes +/// fails rather than passing vacuously. +inline void +check_lbvh_nodes_match(const LBVH::Nodes& nodes, const LBVH::Nodes& reference) +{ + REQUIRE(nodes.size() == reference.size()); + check_valid_lbvh_nodes(nodes); + CHECK((nodes[0].aabb_min == reference[0].aabb_min).all()); + CHECK((nodes[0].aabb_max == reference[0].aabb_max).all()); +} + +} // namespace ipc::tests diff --git a/tests/src/tests/broad_phase/test_broad_phase.cpp b/tests/src/tests/broad_phase/test_broad_phase.cpp index 37ab71545..e30ea47e3 100644 --- a/tests/src/tests/broad_phase/test_broad_phase.cpp +++ b/tests/src/tests/broad_phase/test_broad_phase.cpp @@ -295,12 +295,23 @@ TEST_CASE("Broad phase build from boxes", "[broad_phase]") TEST_CASE("Create broad phase", "[broad_phase]") { -#ifdef IPC_TOOLKIT_WITH_CUDA - uint8_t n_broad_phase_methods = 6; -#else - uint8_t n_broad_phase_methods = 5; + using ipc::BroadPhaseMethod; + + constexpr auto NUM_METHODS = + static_cast(BroadPhaseMethod::NUM_BROAD_PHASE_METHODS); + for (uint8_t i = 0; i < NUM_METHODS; i++) { + const auto method = static_cast(i); + CAPTURE(i); +#ifndef IPC_TOOLKIT_WITH_CUDA + if (method == BroadPhaseMethod::SWEEP_AND_TINIEST_QUEUE + || method == BroadPhaseMethod::LBVH_CUDA) { + CHECK_THROWS(create_broad_phase(method)); + continue; + } #endif - for (uint8_t i = 0; i < n_broad_phase_methods; i++) { - CHECK(create_broad_phase(static_cast(i))); + CHECK(create_broad_phase(method)); } + + // The sentinel is not a method. + CHECK_THROWS(create_broad_phase(BroadPhaseMethod::NUM_BROAD_PHASE_METHODS)); } \ No newline at end of file diff --git a/tests/src/tests/broad_phase/test_gpu_lbvh.cu b/tests/src/tests/broad_phase/test_gpu_lbvh.cu new file mode 100644 index 000000000..a3ba945fa --- /dev/null +++ b/tests/src/tests/broad_phase/test_gpu_lbvh.cu @@ -0,0 +1,560 @@ +// Validates the GPU-built LBVH (ipc::cuda::LBVH) against the CPU ipc::LBVH: +// ipc::cuda::LBVH builds the vertex/edge/face AABBs and BVHs entirely on the +// device. The copied-back trees must be structurally valid (every node +// reachable exactly once, every internal AABB the union of its children, leaf +// set = {0..n-1}) and must agree with the CPU build on node count and root +// AABB (an order-independent union of identically-inflated boxes). The +// detected candidate sets must be exactly equal. + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +using namespace ipc; + +namespace { + +// The GPU and CPU candidate sets are determined by the (bit-identical) box +// overlaps + the same can_*_collide predicate, independent of tree structure, +// so they must be exactly equal as sets. +template +void compare_candidates_exact( + std::vector gpu, std::vector cpu) +{ + std::sort(gpu.begin(), gpu.end()); + std::sort(cpu.begin(), cpu.end()); + CHECK(gpu.size() == cpu.size()); + CHECK(gpu == cpu); +} + +} // namespace + +TEST_CASE("GPU LBVH build", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + constexpr double inflation_radius = 1e-3; + + const std::string mesh = GENERATE("cube.ply", "bunny.ply"); + CAPTURE(mesh); + + Eigen::MatrixXd vertices; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh(mesh, vertices, edges, faces)); + + // GPU build (boxes + BVHs all on the device). + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices, edges, faces, inflation_radius); + + // CPU reference. + LBVH cpu_lbvh; + cpu_lbvh.build(vertices, edges, faces, inflation_radius); + + LBVH::Nodes nodes; + LBVH::RightmostLeaves rightmost; + + SECTION("vertices") + { + gpu_lbvh.vertex_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.vertex_nodes()); + } + SECTION("edges") + { + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.edge_nodes()); + } + SECTION("faces") + { + gpu_lbvh.face_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.face_nodes()); + } + + // clear() empties the device trees. + gpu_lbvh.clear(); + CHECK(gpu_lbvh.num_vertex_nodes() == 0); + gpu_lbvh.vertex_nodes_to_host(nodes, rightmost); + CHECK(nodes.empty()); +} + +TEST_CASE("GPU LBVH detect candidates", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + constexpr double inflation_radius = 0; + + std::string mesh_t0, mesh_t1; + SECTION("Two cubes") + { + mesh_t0 = "two-cubes-far.ply"; + mesh_t1 = "two-cubes-intersecting.ply"; + } + SECTION("Cloth-Ball") + { + mesh_t0 = "cloth_ball92.ply"; + mesh_t1 = "cloth_ball93.ply"; + } + + Eigen::MatrixXd vertices_t0, vertices_t1; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh(mesh_t0, vertices_t0, edges, faces)); + REQUIRE(tests::load_mesh(mesh_t1, vertices_t1, edges, faces)); + + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices_t0, vertices_t1, edges, faces, inflation_radius); + + LBVH cpu_lbvh; + cpu_lbvh.build(vertices_t0, vertices_t1, edges, faces, inflation_radius); + + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_vertex_vertex_candidates(gpu_c); + cpu_lbvh.detect_vertex_vertex_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_vertex_candidates(gpu_c); + cpu_lbvh.detect_edge_vertex_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_edge_candidates(gpu_c); + cpu_lbvh.detect_edge_edge_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + + // With the default (accept-all) filter the device-resident buffer is + // already the exact set (no host trimming needed), and it holds the + // same pairs the host variant materialized. + const cuda::LBVH::DeviceCandidateView view = + gpu_lbvh.detect_edge_edge_candidates_device(); + REQUIRE(view.size == cpu_c.size()); + std::vector a(view.size), b(view.size); + REQUIRE_CUDA(cudaMemcpy( + a.data(), view.a, view.size * sizeof(int32_t), + cudaMemcpyDeviceToHost)); + REQUIRE_CUDA(cudaMemcpy( + b.data(), view.b, view.size * sizeof(int32_t), + cudaMemcpyDeviceToHost)); + std::vector view_c; + for (size_t k = 0; k < view.size; ++k) { + view_c.emplace_back(a[k], b[k]); + } + compare_candidates_exact(view_c, cpu_c); + + // Like every BroadPhase, detection clears its output first: a second + // call replaces the vector rather than doubling it. + gpu_lbvh.detect_edge_edge_candidates(gpu_c); + CHECK(gpu_c.size() == cpu_c.size()); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_face_vertex_candidates(gpu_c); + cpu_lbvh.detect_face_vertex_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_face_candidates(gpu_c); + cpu_lbvh.detect_edge_face_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_face_face_candidates(gpu_c); + cpu_lbvh.detect_face_face_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } +} + +// Exercises the host-fallback path: a non-trivial user vertex filter is not +// device-representable yet, so the device emits the connectivity-filtered +// superset and the host trims it. The result must still match the CPU exactly. +TEST_CASE( + "GPU LBVH detect candidates (custom filter)", + "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + Eigen::MatrixXd vertices_t0, vertices_t1; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh("two-cubes-far.ply", vertices_t0, edges, faces)); + REQUIRE( + tests::load_mesh( + "two-cubes-intersecting.ply", vertices_t1, edges, faces)); + + // An arbitrary (not device-representable) filter -> host fallback. + const auto filter = [](size_t a, size_t b) { return ((a + b) % 2) == 0; }; + + cuda::LBVH gpu_lbvh; + gpu_lbvh.can_vertices_collide = filter; + gpu_lbvh.build(vertices_t0, vertices_t1, edges, faces, 0); + REQUIRE_FALSE(gpu_lbvh.can_vertices_collide.accepts_all()); + + LBVH cpu_lbvh; + cpu_lbvh.can_vertices_collide = filter; + cpu_lbvh.build(vertices_t0, vertices_t1, edges, faces, 0); + + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_edge_candidates(gpu_c); + cpu_lbvh.detect_edge_edge_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + + // The device view is the connectivity-filtered superset the host + // trimmed: never smaller than the exact set. + CHECK( + gpu_lbvh.detect_edge_edge_candidates_device().size >= cpu_c.size()); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_face_vertex_candidates(gpu_c); + cpu_lbvh.detect_face_vertex_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_face_candidates(gpu_c); + cpu_lbvh.detect_edge_face_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + } +} + +// 2D input has no faces; ipc::AABB zero-pads the unused z component without +// inflating it (see build_vertex_boxes_kernel in lbvh.cu), so this also +// exercises that padding path against the CPU's exact behavior. +TEST_CASE("GPU LBVH 2D build and detect", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + Eigen::MatrixXd tmp; + REQUIRE(igl::readCSV((tests::DATA_DIR / "mesh-2D/V_t0.csv").string(), tmp)); + const Eigen::MatrixXd V0 = tmp.leftCols(2); + REQUIRE(igl::readCSV((tests::DATA_DIR / "mesh-2D/V_t1.csv").string(), tmp)); + const Eigen::MatrixXd V1 = tmp.leftCols(2); + Eigen::MatrixXi E; + REQUIRE(igl::readCSV((tests::DATA_DIR / "mesh-2D/E.csv").string(), E)); + E.array() -= 1; // Convert from OBJ format to 0-indexed + const Eigen::MatrixXi F(0, 3); + + constexpr double inflation_radius = 1e-3; + + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(V0, V1, E, F, inflation_radius); + + LBVH cpu_lbvh; + cpu_lbvh.build(V0, V1, E, F, inflation_radius); + + // -- Build parity (structure + root AABB, same checks as the 3D case). -- + LBVH::Nodes nodes; + LBVH::RightmostLeaves rightmost; + gpu_lbvh.vertex_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.vertex_nodes()); + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.edge_nodes()); + CHECK(gpu_lbvh.num_face_nodes() == 0); + + // -- Detection parity (only edge-vertex is meaningful in 2D; mirrors + // BroadPhase::detect_collision_candidates's dim == 2 branch). -- + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_vertex_candidates(gpu_c); + cpu_lbvh.detect_edge_vertex_candidates(cpu_c); + compare_candidates_exact(gpu_c, cpu_c); + CHECK(!gpu_c.empty()); +} + +// A BVH over a single primitive is one node, both root and leaf, and the +// shared descent takes a dedicated branch for such a TARGET. The CPU half of +// this scenario is checked against brute force in test_lbvh.cpp; here the +// device build -- which reaches the same branch through the same shared code +// but from a kernel -- must agree with the host on exactly these trees. +TEST_CASE("GPU LBVH single-primitive trees", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + // One face and one edge, sharing no vertices so the connectivity filter + // keeps the pair, and inflated enough that the AABBs actually overlap. + Eigen::MatrixXd vertices(5, 3); + vertices << 0.00, 0.00, 0.00, // 0 | + 1.00, 0.00, 0.00, // 1 |- the face + 0.00, 1.00, 0.00, // 2 | + 0.05, 0.05, 0.05, // 3 |- the edge + 0.15, 0.05, 0.05; // 4 | + + Eigen::MatrixXi edges(1, 2); + edges << 3, 4; + + Eigen::MatrixXi faces(1, 3); + faces << 0, 1, 2; + + constexpr double inflation_radius = 0.1; + + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices, edges, faces, inflation_radius); + + LBVH cpu_lbvh; + cpu_lbvh.build(vertices, edges, faces, inflation_radius); + + // The branch under test is only reached if these really are single nodes. + REQUIRE(gpu_lbvh.num_face_nodes() == 1); + REQUIRE(gpu_lbvh.num_edge_nodes() == 1); + REQUIRE(cpu_lbvh.face_nodes().size() == 1); + REQUIRE(cpu_lbvh.edge_nodes().size() == 1); + + LBVH::Nodes nodes; + LBVH::RightmostLeaves rightmost; + gpu_lbvh.face_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.face_nodes()); + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.edge_nodes()); + + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_face_vertex_candidates(gpu_c); + cpu_lbvh.detect_face_vertex_candidates(cpu_c); + // Without this the checks would pass on an empty set, which is + // exactly what a broken single-node branch would produce. + REQUIRE(!cpu_c.empty()); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_face_candidates(gpu_c); + cpu_lbvh.detect_edge_face_candidates(cpu_c); + REQUIRE(!cpu_c.empty()); + compare_candidates_exact(gpu_c, cpu_c); + } + { + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_vertex_candidates(gpu_c); + cpu_lbvh.detect_edge_vertex_candidates(cpu_c); + REQUIRE(!cpu_c.empty()); + compare_candidates_exact(gpu_c, cpu_c); + } +} + +// A planar mesh with no inflation makes the Morton normalization domain +// zero-width along z (see morton_domain_width_inv()). The device build must +// take the same guarded path as the host and produce the same trees. +TEST_CASE("GPU LBVH degenerate domain", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + constexpr int N = 4; + Eigen::MatrixXd vertices(N * N, 3); + for (int i = 0; i < N; ++i) { + for (int j = 0; j < N; ++j) { + vertices.row(N * i + j) << i, j, 0.0; + } + } + Eigen::MatrixXi faces(2 * (N - 1) * (N - 1), 3); + for (int i = 0, f = 0; i < N - 1; ++i) { + for (int j = 0; j < N - 1; ++j) { + const int v00 = N * i + j, v10 = v00 + N; + faces.row(f++) << v00, v10, v00 + 1; + faces.row(f++) << v10, v10 + 1, v00 + 1; + } + } + Eigen::MatrixXi edges; + igl::edges(faces, edges); + + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices, edges, faces, /*inflation_radius=*/0); + + LBVH cpu_lbvh; + cpu_lbvh.build(vertices, edges, faces, 0); + + LBVH::Nodes nodes; + LBVH::RightmostLeaves rightmost; + gpu_lbvh.vertex_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.vertex_nodes()); + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.edge_nodes()); + gpu_lbvh.face_nodes_to_host(nodes, rightmost); + tests::check_lbvh_nodes_match(nodes, cpu_lbvh.face_nodes()); + + std::vector gpu_c, cpu_c; + gpu_lbvh.detect_edge_edge_candidates(gpu_c); + cpu_lbvh.detect_edge_edge_candidates(cpu_c); + REQUIRE(!cpu_c.empty()); + compare_candidates_exact(gpu_c, cpu_c); +} + +// The moved-from object must stay usable: it is cleared, not left with a null +// implementation. +TEST_CASE("GPU LBVH move", "[broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + Eigen::MatrixXd vertices; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh("cube.ply", vertices, edges, faces)); + + cuda::LBVH a; + a.build(vertices, edges, faces, 1e-3); + const size_t num_nodes = a.num_face_nodes(); + REQUIRE(num_nodes > 1); + + cuda::LBVH b(std::move(a)); + CHECK(b.num_face_nodes() == num_nodes); + CHECK(a.num_face_nodes() == 0); // NOLINT(bugprone-use-after-move) + + std::vector candidates; + a.detect_face_face_candidates(candidates); // cleared, not null + CHECK(candidates.empty()); + + a.build(vertices, edges, faces, 1e-3); // and rebuildable + CHECK(a.num_face_nodes() == num_nodes); + + cuda::LBVH c; + c = std::move(b); + CHECK(c.num_face_nodes() == num_nodes); + CHECK(b.num_face_nodes() == 0); // NOLINT(bugprone-use-after-move) +} + +// --------------------------------------------------------------------------- +// Benchmarks. Hidden ([!benchmark]) and GPU-gated like every other case here, +// so a CUDA build without a device skips rather than fails them. + +TEST_CASE( + "Benchmark cuda::LBVH::detect_edge_edge_candidates", + "[!benchmark][broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + constexpr double inflation_radius = 0; + + std::string mesh_t0, mesh_t1; + SECTION("Two cubes") + { + mesh_t0 = "two-cubes-far.ply"; + mesh_t1 = "two-cubes-intersecting.ply"; + } + SECTION("Cloth-Ball") + { + mesh_t0 = "cloth_ball92.ply"; + mesh_t1 = "cloth_ball93.ply"; + } +#ifdef NDEBUG + SECTION("Armadillo-Rollers") + { + mesh_t0 = "armadillo-rollers/326.ply"; + mesh_t1 = "armadillo-rollers/327.ply"; + } + SECTION("Cloth-Funnel") + { + mesh_t0 = "cloth-funnel/227.ply"; + mesh_t1 = "cloth-funnel/228.ply"; + } + SECTION("N-Body-Simulation") + { + mesh_t0 = "n-body-simulation/balls16_18.ply"; + mesh_t1 = "n-body-simulation/balls16_19.ply"; + } + SECTION("Rod-Twist") + { + mesh_t0 = "rod-twist/3036.ply"; + mesh_t1 = "rod-twist/3037.ply"; + } +#endif + SECTION("Puffer-Ball") + { + mesh_t0 = "puffer-ball/20.ply"; + mesh_t1 = "puffer-ball/21.ply"; + } + + Eigen::MatrixXd vertices_t0, vertices_t1; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh(mesh_t0, vertices_t0, edges, faces)); + REQUIRE(tests::load_mesh(mesh_t1, vertices_t1, edges, faces)); + + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices_t0, vertices_t1, edges, faces, inflation_radius); + // Warm up the CUDA context so the first sample is not skewed by lazy + // context/allocation initialization. + { + std::vector warmup; + gpu_lbvh.detect_edge_edge_candidates(warmup); + } + + BENCHMARK("cuda::LBVH::detect_edge_edge_candidates") + { + std::vector ee_candidates; + gpu_lbvh.detect_edge_edge_candidates(ee_candidates); + return ee_candidates.size(); + }; +} + +TEST_CASE( + "Benchmark cuda::LBVH::build", "[!benchmark][broad_phase][lbvh][cuda][gpu]") +{ + tests::skip_if_no_cuda_device(); + + constexpr double inflation_radius = 0; + + struct Scene { + std::string name, mesh_t0, mesh_t1; + }; + +#ifdef NDEBUG + constexpr int NUM_SCENES = 6; +#else + constexpr int NUM_SCENES = 1; +#endif + + const std::array scenes = { { + Scene { "Cloth-Ball", "cloth_ball92.ply", "cloth_ball93.ply" }, +#ifdef NDEBUG + Scene { "Cloth-Funnel", "cloth-funnel/227.ply", + "cloth-funnel/228.ply" }, + Scene { "Armadillo-Rollers", "armadillo-rollers/326.ply", + "armadillo-rollers/327.ply" }, + Scene { "Rod-Twist", "rod-twist/3036.ply", "rod-twist/3037.ply" }, + Scene { "N-Body-Simulation", "n-body-simulation/balls16_18.ply", + "n-body-simulation/balls16_19.ply" }, + Scene { "Puffer-Ball", "puffer-ball/20.ply", "puffer-ball/21.ply" }, +#endif + } }; + + for (const auto& [scene, mesh_t0, mesh_t1] : scenes) { + Eigen::MatrixXd vertices_t0, vertices_t1; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh(mesh_t0, vertices_t0, edges, faces)); + REQUIRE(tests::load_mesh(mesh_t1, vertices_t1, edges, faces)); + + cuda::LBVH gpu_lbvh; + // Warm up the CUDA context so the first sample is not skewed by lazy + // context/allocation initialization. + gpu_lbvh.build( + vertices_t0, vertices_t1, edges, faces, inflation_radius); + + BENCHMARK("cuda::LBVH::build [" + scene + "]") + { + gpu_lbvh.build( + vertices_t0, vertices_t1, edges, faces, inflation_radius); + return gpu_lbvh.num_edge_nodes(); + }; + } +} + +#endif // IPC_TOOLKIT_WITH_CUDA diff --git a/tests/src/tests/broad_phase/test_lbvh.cpp b/tests/src/tests/broad_phase/test_lbvh.cpp index 08d4b39eb..f1a630823 100644 --- a/tests/src/tests/broad_phase/test_lbvh.cpp +++ b/tests/src/tests/broad_phase/test_lbvh.cpp @@ -3,76 +3,22 @@ #include #include +#include + +#include #include #include #include +#include + #include #include #include using namespace ipc; - -namespace { - -bool is_aabb_union( - const LBVH::Node& parent, - const LBVH::Node& childA, - const LBVH::Node& childB) -{ - AABB children; - children.min = childA.aabb_min.min(childB.aabb_min).cast(); - children.max = childA.aabb_max.max(childB.aabb_max).cast(); - constexpr float EPS = 1e-4f; - return (abs(parent.aabb_max.cast() - children.max) < EPS).all() - && (abs(parent.aabb_min.cast() - children.min) < EPS).all(); -} - -void traverse_lbvh( - const LBVH::Nodes& lbvh_nodes, - const uint32_t index, - std::vector& visited) -{ - const LBVH::Node& node = lbvh_nodes[index]; - CHECK(node.is_valid()); - - if (node.is_leaf()) { - // leaf - CHECK(!visited[index]); - visited[index] = true; - } else { - // inner node - CHECK(!visited[index]); - visited[index] = true; - - // verify aabbs - LBVH::Node childA = lbvh_nodes[node.left]; - LBVH::Node childB = lbvh_nodes[node.right]; - - { - CAPTURE( - index, node.left, node.right, node.aabb_min.transpose(), - childA.aabb_min.transpose(), childB.aabb_min.transpose(), - node.aabb_max.transpose(), childA.aabb_max.transpose(), - childB.aabb_max.transpose()); - CHECK(is_aabb_union(node, childA, childB)); - } - - // continue traversal - traverse_lbvh(lbvh_nodes, node.left, visited); - traverse_lbvh(lbvh_nodes, node.right, visited); - } -} - -void check_valid_lbvh_nodes(const LBVH::Nodes& lbvh_nodes) -{ - std::vector visited(lbvh_nodes.size(), false); - traverse_lbvh(lbvh_nodes, 0, visited); - REQUIRE( - std::all_of(visited.begin(), visited.end(), [](bool v) { return v; })); -} -} // namespace +using ipc::tests::check_valid_lbvh_nodes; TEST_CASE("LBVH::build", "[broad_phase][lbvh]") { @@ -285,6 +231,140 @@ TEST_CASE("LBVH::detect_*_candidates", "[broad_phase][lbvh]") #endif } +TEST_CASE("LBVH single-primitive trees", "[broad_phase][lbvh]") +{ + // A BVH over a single primitive is one node, which is both the root and a + // leaf. When such a BVH is the traversal TARGET the descent takes a + // dedicated branch, because the root cannot be descended into. Only two + // detections put a BVH there that can have one node -- face-vertex (the + // face BVH) and edge-face (the edge BVH) -- and the meshes the other tests + // load never reduce either to a single primitive. + // + // One face and one edge, sharing no vertices so the connectivity filter + // keeps the pair, and inflated enough that the AABBs actually overlap. + Eigen::MatrixXd vertices(5, 3); + vertices << 0.00, 0.00, 0.00, // 0 | + 1.00, 0.00, 0.00, // 1 |- the face + 0.00, 1.00, 0.00, // 2 | + 0.05, 0.05, 0.05, // 3 |- the edge + 0.15, 0.05, 0.05; // 4 | + + Eigen::MatrixXi edges(1, 2); + edges << 3, 4; + + Eigen::MatrixXi faces(1, 3); + faces << 0, 1, 2; + + constexpr double inflation_radius = 0.1; + + LBVH lbvh; + lbvh.build(vertices, edges, faces, inflation_radius); + + BruteForce brute_force; + brute_force.build(vertices, edges, faces, inflation_radius); + + // The branch under test is only reached if these really are single nodes. + REQUIRE(lbvh.face_nodes().size() == 1); + REQUIRE(lbvh.edge_nodes().size() == 1); + + // The LBVH rounds its AABBs outward to floats, so it may report a superset + // of the exact (double-precision) brute-force set, never a subset. + { + std::vector fv_candidates, expected; + lbvh.detect_face_vertex_candidates(fv_candidates); + brute_force.detect_face_vertex_candidates(expected); + + // Without this the checks below would pass on an empty set, which is + // exactly what a broken single-node branch would produce. + REQUIRE(!expected.empty()); + CHECK(fv_candidates.size() >= expected.size()); + CHECK(contains_all_candidates(fv_candidates, expected)); + } + + { + std::vector ef_candidates, expected; + lbvh.detect_edge_face_candidates(ef_candidates); + brute_force.detect_edge_face_candidates(expected); + + REQUIRE(!expected.empty()); + CHECK(ef_candidates.size() >= expected.size()); + CHECK(contains_all_candidates(ef_candidates, expected)); + } + + // The remaining types traverse multi-node targets here, but are cheap to + // check on a mesh this small. + { + std::vector vv_candidates, expected; + lbvh.detect_vertex_vertex_candidates(vv_candidates); + brute_force.detect_vertex_vertex_candidates(expected); + CHECK(contains_all_candidates(vv_candidates, expected)); + } + + { + std::vector ev_candidates, expected; + lbvh.detect_edge_vertex_candidates(ev_candidates); + brute_force.detect_edge_vertex_candidates(expected); + REQUIRE(!expected.empty()); + CHECK(contains_all_candidates(ev_candidates, expected)); + } + + // The device build reaches the same branch through the same shared code; + // its parity with these trees is checked in test_gpu_lbvh.cu, where a + // missing GPU skips rather than fails. +} + +// A planar mesh with no inflation has zero-width vertex boxes along z, so the +// Morton normalization domain is degenerate on that axis. The reciprocal width +// must be 0 there, not infinity: 0 * inf is NaN, and converting NaN to an +// integer code is undefined. The tree must still be valid and complete. +TEST_CASE("LBVH degenerate domain", "[broad_phase][lbvh]") +{ + // A 4x4 grid of vertices in the z = 0 plane, triangulated. + constexpr int N = 4; + Eigen::MatrixXd vertices(N * N, 3); + for (int i = 0; i < N; ++i) { + for (int j = 0; j < N; ++j) { + vertices.row(N * i + j) << i, j, 0.0; + } + } + Eigen::MatrixXi faces(2 * (N - 1) * (N - 1), 3); + for (int i = 0, f = 0; i < N - 1; ++i) { + for (int j = 0; j < N - 1; ++j) { + const int v00 = N * i + j, v10 = v00 + N; + faces.row(f++) << v00, v10, v00 + 1; + faces.row(f++) << v10, v10 + 1, v00 + 1; + } + } + Eigen::MatrixXi edges; + igl::edges(faces, edges); + + // Inflation 0 keeps the z extent of every box exactly zero. + LBVH lbvh; + lbvh.build(vertices, edges, faces, /*inflation_radius=*/0); + + check_valid_lbvh_nodes(lbvh.vertex_nodes()); + check_valid_lbvh_nodes(lbvh.edge_nodes()); + check_valid_lbvh_nodes(lbvh.face_nodes()); + + // And the result still matches brute force (a superset, as usual). + BruteForce brute_force; + brute_force.build(vertices, edges, faces, 0); + { + std::vector candidates, expected; + lbvh.detect_edge_edge_candidates(candidates); + brute_force.detect_edge_edge_candidates(expected); + REQUIRE(!expected.empty()); // coplanar neighbors do overlap + CHECK(contains_all_candidates(candidates, expected)); + } + { + std::vector candidates, expected; + lbvh.detect_face_vertex_candidates(candidates); + brute_force.detect_face_vertex_candidates(expected); + REQUIRE(!expected.empty()); + CHECK(contains_all_candidates(candidates, expected)); + } +} + TEST_CASE( "Benchmark LBVH::detect_edge_edge_candidates", "[!benchmark][broad_phase][lbvh]") @@ -344,6 +424,8 @@ TEST_CASE( lbvh->detect_edge_edge_candidates(ee_candidates); return ee_candidates.size(); }; + // The cuda::LBVH counterpart lives in test_gpu_lbvh.cu, behind the GPU + // skip. } TEST_CASE("Benchmark LBVH::build", "[!benchmark][broad_phase][lbvh]") @@ -388,5 +470,7 @@ TEST_CASE("Benchmark LBVH::build", "[!benchmark][broad_phase][lbvh]") vertices_t0, vertices_t1, edges, faces, inflation_radius); return lbvh->edge_nodes().size(); }; + // The cuda::LBVH counterpart lives in test_gpu_lbvh.cu, behind the + // GPU skip. } } \ No newline at end of file diff --git a/tests/src/tests/test_collision_filter.cpp b/tests/src/tests/test_collision_filter.cpp index d2b2bbe26..03cb5fb9a 100644 --- a/tests/src/tests/test_collision_filter.cpp +++ b/tests/src/tests/test_collision_filter.cpp @@ -234,3 +234,82 @@ TEST_CASE("CollisionFilter composition chain", "[collision_filter]") CHECK_FALSE(active(0, 1)); CHECK_FALSE(active(3, 5)); } + +// ───────────────────────────────────────────────────────────────────────────── + +TEST_CASE("CollisionFilter accepts_all", "[collision_filter]") +{ + // accepts_all() is what lets a GPU broad phase skip host-side filtering, + // so a wrong `true` silently drops the user's filter. Pin both directions. + + SECTION("default filter") + { + const CollisionFilter f; + CHECK(f.accepts_all()); + CHECK(f(0, 1)); + CHECK(f(7, 7)); + + // The conversion must still yield a callable function. + const std::function fn = f; + REQUIRE(fn); + CHECK(fn(3, 4)); + } + + SECTION("an empty std::function is the accept-all filter") + { + const std::function empty; + const CollisionFilter f(empty); + CHECK(f.accepts_all()); + CHECK(f(0, 1)); + } + + SECTION("a user callable is never accept-all, even if it returns true") + { + const CollisionFilter f([](size_t, size_t) { return true; }); + CHECK_FALSE(f.accepts_all()); + CHECK(f(0, 1)); + } + + SECTION("compositions") + { + const CollisionFilter all; + const CollisionFilter odd( + [](size_t vi, size_t vj) { return (vi + vj) % 2 == 1; }); + + // accept-all absorbs a union ... + CHECK((all | all).accepts_all()); + CHECK((odd | all).accepts_all()); + CHECK((all | odd).accepts_all()); + CHECK((odd | all)(0, 2)); + + // ... and is the identity of an intersection. + CHECK((all & all).accepts_all()); + CHECK_FALSE((odd & all).accepts_all()); + CHECK_FALSE((all & odd).accepts_all()); + CHECK((all & odd)(0, 1)); + CHECK_FALSE((all & odd)(0, 2)); + + // Negating accept-all rejects everything and is not accept-all. + const CollisionFilter none = !all; + CHECK_FALSE(none.accepts_all()); + CHECK_FALSE(none(0, 1)); + CHECK_FALSE((!odd).accepts_all()); + + // Compound assignment follows the same rules. + CollisionFilter g; + g &= odd; + CHECK_FALSE(g.accepts_all()); + CHECK(g(0, 1)); + g |= all; + CHECK(g.accepts_all()); + } + + SECTION("factories are never accept-all") + { + Eigen::VectorXi patches(2); + patches << 0, 1; + CHECK_FALSE(make_vertex_patches_filter(patches).accepts_all()); + CHECK_FALSE(make_static_obstacle_filter(1).accepts_all()); + CHECK_FALSE(make_codim_cross_filter(1).accepts_all()); + } +} diff --git a/tests/src/tests/utils.cpp b/tests/src/tests/utils.cpp index e05fb47b5..5e3a9dd34 100644 --- a/tests/src/tests/utils.cpp +++ b/tests/src/tests/utils.cpp @@ -8,6 +8,7 @@ #include #include #ifdef IPC_TOOLKIT_WITH_CUDA +#include #include #endif @@ -30,6 +31,7 @@ std::vector> broad_phases() std::make_shared(), #ifdef IPC_TOOLKIT_WITH_CUDA std::make_shared(), + std::make_shared(), #endif } }; }