From a79d4d49ae5779d30b5dc5a321392fc767584a31 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 17:16:35 -0400 Subject: [PATCH 01/10] Add shared device utilities for the ipc::cuda implementations Introduce src/ipc/utils/cuda/device_utils.cuh, the header every ipc::cuda translation unit needs before it can launch anything: - IPC_TOOLKIT_CUDA_CHECK, which turns a cudaError_t into a std::runtime_error naming the file and line. - KERNEL_BLOCK_SIZE and kernel_grid_size(), the single definition of the launch geometry, so the block size is not repeated per call site. - global_dof_index(), mirroring the index math of local_gradient_to_global_gradient() for device-side gradient scatter. - A compile-time guard rejecting compute capability < 6.0, where atomicAdd(double*, double) does not exist. Include directly: global_dof_index() compares VERTEX_DERIVATIVE_LAYOUT against Eigen::RowMajor, and config.hpp deliberately defines its own Eigen-free layout constants rather than pulling in Eigen, so the header would otherwise only compile when an includer happened to have included Eigen first. The header is CUDA-only and included from .cu files exclusively; it is wired in under IPC_TOOLKIT_WITH_CUDA so a non-CUDA build never sees it. --- src/ipc/utils/CMakeLists.txt | 4 ++ src/ipc/utils/cuda/CMakeLists.txt | 5 +++ src/ipc/utils/cuda/device_utils.cuh | 59 +++++++++++++++++++++++++++++ 3 files changed, 68 insertions(+) create mode 100644 src/ipc/utils/cuda/CMakeLists.txt create mode 100644 src/ipc/utils/cuda/device_utils.cuh 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..7d8376ad2 --- /dev/null +++ b/src/ipc/utils/cuda/CMakeLists.txt @@ -0,0 +1,5 @@ +set(SOURCES + device_utils.cuh +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) diff --git a/src/ipc/utils/cuda/device_utils.cuh b/src/ipc/utils/cuda/device_utils.cuh new file mode 100644 index 000000000..65efa3e46 --- /dev/null +++ b/src/ipc/utils/cuda/device_utils.cuh @@ -0,0 +1,59 @@ +// 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 // Eigen::RowMajor, for VERTEX_DERIVATIVE_LAYOUT + +#include +#include + +// atomicAdd(double*, double) requires compute capability 6.0+. +#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600 +#error "ipc::cuda requires compute capability 6.0+ (atomicAdd on double)." +#endif + +/// @brief Throw a std::runtime_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) { \ + throw std::runtime_error( \ + std::string("CUDA error at " __FILE__ ":") \ + + std::to_string(__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); +} + +/// @brief Global DOF index of component @p d of vertex @p vertex_id. +/// Mirrors the index math of local_gradient_to_global_gradient() +/// (see src/ipc/utils/local_to_global.hpp) for dim=3. +__device__ inline index_t global_dof_index( + const index_t vertex_id, const int d, const index_t n_total_vertices) +{ + if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) { + return 3 * vertex_id + d; + } else { + return n_total_vertices * d + vertex_id; + } +} + +} // namespace ipc::cuda + +#endif // IPC_TOOLKIT_WITH_CUDA From b2b35ca404c1a6b9302ef3ab173233909e36a04a Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 24 Jul 2026 17:52:28 -0700 Subject: [PATCH 02/10] Add ipc::cuda::LBVH: GPU-native LBVH broad phase A first-class GPU counterpart to ipc::LBVH (not a CPU-upload adapter): builds vertex/edge/face AABBs and their BVHs entirely on the device (Morton codes + Apetrei 2014 single-pass bottom-up construction, reusing the 32-byte ipc::LBVH::Node layout for host validation/interop), then runs candidate detection with the BVH descent and mesh-connectivity (shared-vertex) exclusion both on the device. The user vertex filter is honored on the device for the common accept-all case (new CollisionFilter::accepts_all()); a non-trivial filter falls back to a host pass over the device-emitted, connectivity-filtered candidates. Either path matches the CPU ipc::LBVH's candidate set exactly. Adds DeviceCandidateView + detect_*_candidates_device() so candidates can stay device-resident for a future GPU-native pipeline (e.g. device Additive CCD) instead of always materializing to host vectors. Supporting changes: ipc::math::morton_2D/3D and expand_bits_1/2 are now IPC_TOOLKIT_HOST_DEVICE so the device Morton codes reuse the exact CPU implementation; the Morton-normalization reciprocal is now precomputed once per build and multiplied per box instead of divided (CPU and GPU changed identically so their Morton codes stay bit-matched to each other). Validation: build + detect + custom-filter-fallback GPU-run-validated on an RTX 3070 (artemis): 150517 assertions across 3 test cases, plus exact candidate-set parity against the CPU LBVH for all 6 candidate types. Benchmarked against the CPU LBVH (edge-edge detection): 1.1-1.7x faster on every real mesh tested except a trivial two-cube case. The Morton reciprocal-multiply optimization and code cleanup (Eigen::Array3d in place of a hand-rolled Vec3d, .min()/.max() in place of manual fminf/fmaxf loops) landed after artemis went offline and are Docker-compile-validated only; pending a GPU re-run. Not yet done: ipc::cuda::LBVH is not registered in BroadPhaseMethod / create_broad_phase (deferred until the device-resident candidate path is consumed by something), and the connectivity/user-filter split does not yet support device-side patch/connected-component filters (would need a label-data CollisionFilter descriptor). Co-Authored-By: Claude Sonnet 5 --- src/ipc/broad_phase/CMakeLists.txt | 4 + src/ipc/broad_phase/cuda/CMakeLists.txt | 7 + src/ipc/broad_phase/cuda/lbvh.cu | 1360 ++++++++++++++++++ src/ipc/broad_phase/cuda/lbvh.hpp | 175 +++ src/ipc/broad_phase/cuda/lbvh_impl.cuh | 94 ++ src/ipc/collision_filter.hpp | 16 +- tests/src/tests/broad_phase/CMakeLists.txt | 6 + tests/src/tests/broad_phase/test_gpu_lbvh.cu | 309 ++++ tests/src/tests/broad_phase/test_lbvh.cpp | 37 + 9 files changed, 2007 insertions(+), 1 deletion(-) create mode 100644 src/ipc/broad_phase/cuda/CMakeLists.txt create mode 100644 src/ipc/broad_phase/cuda/lbvh.cu create mode 100644 src/ipc/broad_phase/cuda/lbvh.hpp create mode 100644 src/ipc/broad_phase/cuda/lbvh_impl.cuh create mode 100644 tests/src/tests/broad_phase/test_gpu_lbvh.cu diff --git a/src/ipc/broad_phase/CMakeLists.txt b/src/ipc/broad_phase/CMakeLists.txt index e0613ebcf..2ae7503e5 100644 --- a/src/ipc/broad_phase/CMakeLists.txt +++ b/src/ipc/broad_phase/CMakeLists.txt @@ -23,3 +23,7 @@ set(SOURCES ) target_sources(ipc_toolkit PRIVATE ${SOURCES}) + +if(IPC_TOOLKIT_WITH_CUDA) + add_subdirectory(cuda) +endif() 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..d2ae174da --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -0,0 +1,1360 @@ +#include "lbvh.hpp" + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace ipc::cuda { + +namespace { + + // Eigen::Array3d is passed to kernels by value, so it must be exactly three + // packed doubles (no vectorization padding) to have a stable layout. + static_assert( + sizeof(Eigen::Array3d) == 24, + "Eigen::Array3d must be 24 bytes (3 packed doubles)"); + + /// @brief Per-internal-node scratch used by the bottom-up build. The device + /// analog of ipc::LBVH::ConstructionInfo, kept separate on purpose: that + /// struct's visitation_count is a std::atomic, which cannot be used + /// here (atomicAdd needs an int*, and std::atomic is non-copyable so it + /// cannot be a thrust::device_vector element). A plain int suffices because + /// atomicAdd provides the atomicity the CPU gets from std::atomic. + struct DeviceConstructionInfo { + int left_range; + int right_range; + int visitation_count; + }; + + /// @brief Min/max domain accumulator for the Morton-normalization reduction. + struct Domain { + double mn[3]; + double mx[3]; + }; + + 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.mn[k] = fmin(a.mn[k], b.mn[k]); + r.mx[k] = fmax(a.mx[k], b.mx[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.mn[k] = box_min[3 * i + k]; + d.mx[k] = box_max[3 * i + k]; + } + return d; + } + }; + + // -- Box building ------------------------------------------------------- + // Matches ipc::build_*_boxes + AABB::conservative_inflation exactly: the + // double bounds are nudged outward with nextafter so the box is + // conservative. (The leaf nodes later apply a second float-nextafter in + // build_hierarchy_kernel, matching assign_inflated_aabb.) + + __global__ void build_vertex_boxes_static_kernel( + const double* __restrict__ vertices, // 3 * n, row-major + const int n, + 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) { + const double v = vertices[3 * i + k]; + box_min[3 * i + k] = nextafter(v - inflation_radius, -INFINITY); + box_max[3 * i + k] = nextafter(v + inflation_radius, INFINITY); + } + } + + __global__ void build_vertex_boxes_dynamic_kernel( + const double* __restrict__ vertices_t0, // 3 * n, row-major + const double* __restrict__ vertices_t1, // 3 * n, row-major + const int n, + 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) { + const double a = vertices_t0[3 * i + k]; + const double b = vertices_t1[3 * i + k]; + // union of the two inflated point boxes; nextafter is monotonic so + // min(nextafter(a),nextafter(b)) == nextafter(min(a,b)). + box_min[3 * i + k] = + nextafter(fmin(a, b) - inflation_radius, -INFINITY); + box_max[3 * i + k] = + nextafter(fmax(a, b) + inflation_radius, INFINITY); + } + } + + __global__ void build_edge_boxes_kernel( + const double* __restrict__ vbox_min, + const double* __restrict__ vbox_max, + const index_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 index_t e0 = edges[2 * i + 0]; + const index_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]); + } + } + + __global__ void build_face_boxes_kernel( + const double* __restrict__ vbox_min, + const double* __restrict__ vbox_max, + const index_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 index_t f0 = faces[3 * i + 0]; + const index_t f1 = faces[3 * i + 1]; + const index_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 Number of common leading bits between Morton codes at sorted + /// positions i and j (device port of the CPU delta()). Duplicate codes fall + /// back to the CLZ of the index XOR (offset by 32 so it sorts after any + /// code-level difference). + /// @param sorted_codes The Morton codes in sorted order. + /// @param n The number of codes. + /// @param i The first sorted position. + /// @param code_i The code at position i (passed to avoid a redundant look-up). + /// @param j The second sorted position. + /// @return The common-prefix length, or -1 when j is out of bounds. + __device__ inline int delta_device( + const uint64_t* __restrict__ sorted_codes, + const int n, + const int i, + const uint64_t code_i, + const int j) + { + if (j < 0 || j >= n) { + return -1; + } + const uint64_t code_j = sorted_codes[j]; + if (code_i == code_j) { + return 32 + __clz(i ^ j); + } + return __clzll(static_cast(code_i ^ code_j)); + } + + /// @brief Compute one Morton code per box from its (normalized) center. + /// Mirrors the compute_morton_codes block of ipc::LBVH::init_bvh. + __global__ void compute_morton_codes_kernel( + const double* __restrict__ box_min, // 3 * n, row-major + const double* __restrict__ box_max, // 3 * n, row-major + const int n, + const Eigen::Array3d mesh_min, + const Eigen::Array3d mesh_width_inv, + const int dim, + uint64_t* __restrict__ codes, + index_t* __restrict__ box_ids) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } + + const double cx = 0.5 * (box_min[3 * i + 0] + box_max[3 * i + 0]); + const double cy = 0.5 * (box_min[3 * i + 1] + box_max[3 * i + 1]); + const double cz = 0.5 * (box_min[3 * i + 2] + box_max[3 * i + 2]); + + // (center - mesh_min) * mesh_width_inv -- the reciprocal is + // precomputed once per build (see compute_domain) and multiplied here + // instead of dividing per box, matching the CPU (ipc::LBVH::init_bvh) + // bit-for-bit. + const double mx = (cx - mesh_min.x()) * mesh_width_inv.x(); + const double my = (cy - mesh_min.y()) * mesh_width_inv.y(); + const double mz = (cz - mesh_min.z()) * mesh_width_inv.z(); + + codes[i] = (dim == 2) ? morton_2D(mx, my) : morton_3D(mx, my, mz); + box_ids[i] = i; + } + + /// @brief Single-pass bottom-up hierarchy + AABB build (Apetrei 2014). + /// One thread per leaf. Direct port of the build_hierarchy_and_boxes block + /// of ipc::LBVH::init_bvh, with atomicAdd + __threadfence replacing the + /// std::atomic arrival gate. + __global__ void build_hierarchy_kernel( + const double* __restrict__ box_min, + const double* __restrict__ box_max, + const uint64_t* __restrict__ sorted_codes, + const index_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 int LEAF_OFFSET = N_LEAVES - 1; + + // --- Initialize leaf node --- + { + const index_t bid = sorted_box_ids[i]; + ipc::LBVH::Node leaf; +#pragma unroll + for (int k = 0; k < 3; ++k) { + // Round the float AABB out (matches assign_inflated_aabb). + leaf.aabb_min[k] = nextafterf( + static_cast(box_min[3 * bid + k]), -INFINITY); + leaf.aabb_max[k] = nextafterf( + static_cast(box_max[3 * bid + k]), INFINITY); + } + leaf.primitive_id = static_cast(bid); + leaf.is_inner_marker = 0; + nodes[LEAF_OFFSET + i] = leaf; + // A leaf's rightmost leaf is itself. + rightmost[LEAF_OFFSET + i] = i; + } + + // Single-node tree: the leaf is the root; no internal nodes to build. + if (N_LEAVES == 1) { + if (i == 0) { + *root_idx = 0; + } + return; + } + + // --- Bottom-up walk (Apetrei 2014, Fig. 2) --- + int left_key = i; + int right_key = i; + int current_node = LEAF_OFFSET + i; + + while (true) { + // Choose parent (see the CPU comment in ipc::LBVH::init_bvh). + const bool is_child_a = (left_key == 0) + || (right_key != N_LEAVES - 1 + && delta_device( + sorted_codes, N_LEAVES, right_key, + sorted_codes[right_key], right_key + 1) + > delta_device( + sorted_codes, N_LEAVES, left_key - 1, + sorted_codes[left_key - 1], left_key)); + const int parent = is_child_a ? right_key : left_key - 1; + + // Write the child pointer + range onto the parent. + 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; + } + + // Publish this child's node data and range to all threads before + // signaling arrival, so the second thread reads consistent state. + __threadfence(); + + // Atomic arrival gate: first thread stops; second proceeds knowing + // both children are complete. + if (atomicAdd(&infos[parent].visitation_count, 1) == 0) { + break; // first thread to arrive -> finished + } + + // Second thread: compute the parent AABB union and rightmost leaf. + const ipc::LBVH::Node& child_a = nodes[nodes[parent].left]; + const ipc::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); + rightmost[parent] = ::max( + rightmost[nodes[parent].left], rightmost[nodes[parent].right]); + + // 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) { + // Only one thread reaches the root. + *root_idx = current_node; + break; + } + } + } + + /// @brief Swap the node and rightmost-leaf entries at indices 0 and root + /// (runs on a single thread). + /// @param nodes The BVH nodes. + /// @param rightmost The per-node rightmost-leaf indices. + /// @param root The index to swap with index 0. + __global__ void swap_root_kernel( + ipc::LBVH::Node* __restrict__ nodes, + int32_t* __restrict__ rightmost, + const int root) + { + if (blockIdx.x == 0 && threadIdx.x == 0) { + const ipc::LBVH::Node tmp = nodes[0]; + nodes[0] = nodes[root]; + nodes[root] = tmp; + const int32_t t = rightmost[0]; + rightmost[0] = rightmost[root]; + rightmost[root] = t; + } + } + + /// @brief After the root swap, rewrite left pointers that referenced the + /// old node 0 to its new location. See the CPU swap_root_to_zero comment: + /// the old node 0 was only ever a left child, so only .left needs patching. + /// is_inner_marker aliases .right and is nonzero iff internal. + /// @param nodes The BVH nodes. + /// @param num_nodes The number of nodes. + /// @param root The new location of the old node 0. + __global__ void patch_left_kernel( + ipc::LBVH::Node* __restrict__ nodes, + const int num_nodes, + const int root) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= num_nodes) { + return; + } + if (nodes[i].is_inner_marker != 0 && nodes[i].left == 0) { + nodes[i].left = root; + } + } + + /// @brief Build one BVH on the device from device-resident box corners. + /// Mirrors ipc::LBVH::init_bvh; the output BVH is resized and filled in + /// place. + /// @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 mesh_min The Morton-normalization domain minimum. + /// @param mesh_width_inv The reciprocal of the Morton-normalization domain + /// extent (precomputed once per build; see compute_domain). + /// @param dim The simulation dimension (2 or 3). + /// @param bvh The BVH to build (output). + void build_tree( + const double* d_box_min, + const double* d_box_max, + const int n, + const Eigen::Array3d& mesh_min, + const Eigen::Array3d& mesh_width_inv, + const int dim, + LBVH::Impl::DeviceBVH& bvh) + { + bvh.n_leaves = n; + if (n == 0) { + bvh.nodes.clear(); + bvh.rightmost_leaves.clear(); + return; + } + + const size_t num_nodes = size_t(2) * n - 1; + bvh.nodes.resize(num_nodes); + bvh.rightmost_leaves.resize(num_nodes); + + thrust::device_vector morton_codes(n); + thrust::device_vector box_ids(n); + // Value-initialized to zero => visitation_count starts at 0. + thrust::device_vector infos(num_nodes); + thrust::device_vector d_root(1, -1); + + compute_morton_codes_kernel<<>>( + d_box_min, d_box_max, n, mesh_min, mesh_width_inv, dim, + thrust::raw_pointer_cast(morton_codes.data()), + thrust::raw_pointer_cast(box_ids.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + thrust::sort_by_key( + morton_codes.begin(), morton_codes.end(), box_ids.begin()); + + build_hierarchy_kernel<<>>( + d_box_min, d_box_max, thrust::raw_pointer_cast(morton_codes.data()), + thrust::raw_pointer_cast(box_ids.data()), n, + thrust::raw_pointer_cast(bvh.nodes.data()), + thrust::raw_pointer_cast(bvh.rightmost_leaves.data()), + thrust::raw_pointer_cast(infos.data()), + thrust::raw_pointer_cast(d_root.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + const int root = d_root[0]; // device->host read + if (root > 0) { + swap_root_kernel<<<1, 1>>>( + thrust::raw_pointer_cast(bvh.nodes.data()), + thrust::raw_pointer_cast(bvh.rightmost_leaves.data()), root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + patch_left_kernel<<< + kernel_grid_size(num_nodes), KERNEL_BLOCK_SIZE>>>( + thrust::raw_pointer_cast(bvh.nodes.data()), + static_cast(num_nodes), root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + } + + /// @brief Compute the Morton-normalization domain (min of mins, max of + /// maxs) over device-resident vertex box corners, and its reciprocal + /// extent. The reciprocal is computed once here (per build) and multiplied + /// per box in compute_morton_codes_kernel instead of dividing per box, + /// matching the CPU (ipc::LBVH::init_bvh) bit-for-bit. + void compute_domain( + const thrust::device_vector& vbox_min, + const thrust::device_vector& vbox_max, + const int n_vertices, + Eigen::Array3d& mesh_min, + Eigen::Array3d& mesh_width_inv) + { + Domain init; + for (int k = 0; k < 3; ++k) { + init.mn[k] = std::numeric_limits::max(); + init.mx[k] = std::numeric_limits::lowest(); + } + const Domain dom = thrust::transform_reduce( + thrust::counting_iterator(0), + thrust::counting_iterator(n_vertices), + MakeDomain { thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data()) }, + init, DomainReduce {}); + + mesh_min = Eigen::Array3d(dom.mn[0], dom.mn[1], dom.mn[2]); + const Eigen::Array3d mesh_width( + dom.mx[0] - dom.mn[0], dom.mx[1] - dom.mn[1], + dom.mx[2] - dom.mn[2]); + mesh_width_inv = 1.0 / mesh_width; + } + + // Upload an integer connectivity matrix (rowwise) as a flat row-major + // index_t device array. + template + thrust::device_vector + upload_connectivity(Eigen::ConstRef M) + { + const size_t n = M.rows(); + std::vector h(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)); + } + } + return thrust::device_vector(h); + } + + 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()); + thrust::copy(bvh.nodes.begin(), bvh.nodes.end(), nodes.begin()); + thrust::copy( + bvh.rightmost_leaves.begin(), bvh.rightmost_leaves.end(), + rightmost_leaves.begin()); + } + + /// @brief Given device-resident vertex boxes, 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 vbox_min The vertex box min corners (3 * n_vertices, device). + /// @param vbox_max The vertex box max corners (3 * n_vertices, device). + /// @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 thrust::device_vector& vbox_min, + const thrust::device_vector& vbox_max, + 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 to the device, and keep a host copy for the + // host-side can_*_collide filters. + impl.edges = upload_connectivity<2>(edges); + impl.faces = upload_connectivity<3>(faces); + + impl.h_edge_vertex_ids.resize(n_edges); + for (int i = 0; i < n_edges; ++i) { + impl.h_edge_vertex_ids[i] = { { static_cast(edges(i, 0)), + static_cast( + edges(i, 1)) } }; + } + impl.h_face_vertex_ids.resize(n_faces); + for (int i = 0; i < n_faces; ++i) { + impl.h_face_vertex_ids[i] = { { static_cast(faces(i, 0)), + static_cast(faces(i, 1)), + static_cast( + faces(i, 2)) } }; + } + + // Build edge/face boxes on the device from the vertex boxes. + thrust::device_vector ebox_min(3 * size_t(n_edges)); + thrust::device_vector ebox_max(3 * size_t(n_edges)); + if (n_edges > 0) { + build_edge_boxes_kernel<<< + kernel_grid_size(n_edges), KERNEL_BLOCK_SIZE>>>( + thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data()), + thrust::raw_pointer_cast(impl.edges.data()), n_edges, + thrust::raw_pointer_cast(ebox_min.data()), + thrust::raw_pointer_cast(ebox_max.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + + thrust::device_vector fbox_min(3 * size_t(n_faces)); + thrust::device_vector fbox_max(3 * size_t(n_faces)); + if (n_faces > 0) { + build_face_boxes_kernel<<< + kernel_grid_size(n_faces), KERNEL_BLOCK_SIZE>>>( + thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data()), + thrust::raw_pointer_cast(impl.faces.data()), n_faces, + thrust::raw_pointer_cast(fbox_min.data()), + thrust::raw_pointer_cast(fbox_max.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } + + // The CPU normalizes all three BVHs by the vertex box domain. + Eigen::Array3d mesh_min, mesh_width_inv; + compute_domain( + vbox_min, vbox_max, n_vertices, mesh_min, mesh_width_inv); + + build_tree( + thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data()), n_vertices, mesh_min, + mesh_width_inv, dim, impl.vertex_bvh); + build_tree( + thrust::raw_pointer_cast(ebox_min.data()), + thrust::raw_pointer_cast(ebox_max.data()), n_edges, mesh_min, + mesh_width_inv, dim, impl.edge_bvh); + build_tree( + thrust::raw_pointer_cast(fbox_min.data()), + thrust::raw_pointer_cast(fbox_max.data()), n_faces, mesh_min, + mesh_width_inv, dim, impl.face_bvh); + + IPC_TOOLKIT_CUDA_CHECK(cudaDeviceSynchronize()); + } + + // -- Traversal ---------------------------------------------------------- + + __device__ inline bool + aabb_intersects(const ipc::LBVH::Node& a, const ipc::LBVH::Node& b) + { + return a.aabb_min[0] <= b.aabb_max[0] && b.aabb_min[0] <= a.aabb_max[0] + && a.aabb_min[1] <= b.aabb_max[1] && b.aabb_min[1] <= a.aabb_max[1] + && a.aabb_min[2] <= b.aabb_max[2] && b.aabb_min[2] <= a.aabb_max[2]; + } + + /// @brief Whether two primitives share a vertex id (the device connectivity + /// filter). A vertex primitive's id set is {itself}; an edge's is its 2 + /// endpoints; a face's is its 3 vertices. This is exactly the + /// shared-endpoint exclusion in ipc::LBVH's can_*_collide (for + /// vertex-vertex it reduces to p_a == p_b). + /// @param p_a The first primitive id. + /// @param conn_a The first primitive's connectivity, or null for a vertex. + /// @param count_a The number of vertex ids per first primitive (1, 2, or 3). + /// @param p_b The second primitive id. + /// @param conn_b The second primitive's connectivity, or null for a vertex. + /// @param count_b The number of vertex ids per second primitive. + /// @return Whether the two primitives share any vertex id. + __device__ inline bool prim_shares_vertex( + const int p_a, + const index_t* __restrict__ conn_a, + const int count_a, + const int p_b, + const index_t* __restrict__ conn_b, + const int count_b) + { + index_t ids_a[3]; + index_t ids_b[3]; + if (conn_a == nullptr) { + ids_a[0] = p_a; + } else { + for (int k = 0; k < count_a; ++k) { + ids_a[k] = conn_a[count_a * p_a + k]; + } + } + if (conn_b == nullptr) { + ids_b[0] = p_b; + } else { + for (int k = 0; k < count_b; ++k) { + ids_b[k] = conn_b[count_b * p_b + k]; + } + } + for (int i = 0; i < count_a; ++i) { + for (int j = 0; j < count_b; ++j) { + if (ids_a[i] == ids_b[j]) { + return true; + } + } + } + return false; + } + + /// @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. + template + __device__ inline void emit_pair( + const int query_prim, + const int node_prim, + int32_t* __restrict__ out_a, + int32_t* __restrict__ out_b, + int* __restrict__ counter, + const int capacity) + { + int a = query_prim, b = node_prim; + if constexpr (swap_order) { + const int t = a; + a = b; + b = t; + } + const int slot = atomicAdd(counter, 1); + 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. Descent is a direct port of traverse_lbvh() + /// in lbvh.cpp (scalar path); the connectivity (shared-vertex) exclusion is + /// applied here on the device. The remaining user vertex filter (if any) is + /// applied on the host, so the final set matches the CPU ipc::LBVH. + /// @tparam triangular Self-collision: skip subtrees fully left of the query. + /// @tparam swap_order Emit (target_prim, source_prim) instead. + template + __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 index_t* __restrict__ source_conn, // null for vertex primitives + const int source_count, // ids per source primitive + const index_t* __restrict__ target_conn, // null for vertex primitives + const int target_count, // ids per target primitive + int32_t* __restrict__ out_a, + int32_t* __restrict__ out_b, + int* __restrict__ counter, + const int 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]; + const int query_leaf_idx = s; + + constexpr int MAX_STACK_SIZE = 64; + int stack[MAX_STACK_SIZE]; + int stack_ptr = 0; + stack[stack_ptr++] = ipc::LBVH::Node::INVALID_POINTER; // 0 + + int node_idx = 0; // root + do { + const ipc::LBVH::Node& node = target[node_idx]; + + if (target_size == 1) { // single node (only root, which is a leaf) + if constexpr (triangular) { + break; // no self-collision with a single primitive + } + if (aabb_intersects(node, query) + && !prim_shares_vertex( + query.primitive_id, source_conn, source_count, + node.primitive_id, target_conn, target_count)) { + emit_pair( + query.primitive_id, node.primitive_id, out_a, out_b, + counter, capacity); + } + break; + } + + const ipc::LBVH::Node& child_l = target[node.left]; + const ipc::LBVH::Node& child_r = target[node.right]; + bool intersects_l = aabb_intersects(child_l, query); + bool intersects_r = aabb_intersects(child_r, query); + + // Skip subtrees fully on the query's left (triangular only). + if constexpr (triangular) { + 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; + } + } + + // is_inner_marker aliases .right; it is 0 iff the node is a leaf. + const bool l_leaf = (child_l.is_inner_marker == 0); + const bool r_leaf = (child_r.is_inner_marker == 0); + + if (intersects_l && l_leaf + && !prim_shares_vertex( + query.primitive_id, source_conn, source_count, + child_l.primitive_id, target_conn, target_count)) { + emit_pair( + query.primitive_id, child_l.primitive_id, out_a, out_b, + counter, capacity); + } + if (intersects_r && r_leaf + && !prim_shares_vertex( + query.primitive_id, source_conn, source_count, + child_r.primitive_id, target_conn, target_count)) { + emit_pair( + query.primitive_id, child_r.primitive_id, out_a, out_b, + counter, capacity); + } + + const bool traverse_l = intersects_l && !l_leaf; + const bool traverse_r = intersects_r && !r_leaf; + + if (!traverse_l && !traverse_r) { + node_idx = stack[--stack_ptr]; + } else { + node_idx = traverse_l ? node.left : node.right; + if (traverse_l && traverse_r) { + stack[stack_ptr++] = node.right; + } + } + } while (node_idx != ipc::LBVH::Node::INVALID_POINTER); + } + + /// @brief Run the device traversal of the target BVH by the source leaves, + /// leaving the connectivity-filtered candidate pairs device-resident in the + /// output buffers (resized to the exact count). Grows the buffer and + /// re-runs once if the first pass overflows. + /// @tparam triangular Self-collision: skip subtrees fully left of the query. + /// @tparam swap_order Emit (target_prim, source_prim) instead. + /// @param source The BVH whose leaves are the queries. + /// @param target The BVH to descend. + /// @param source_conn The source primitives' connectivity (null for vertices). + /// @param source_count The vertex ids per source primitive (1, 2, or 3). + /// @param target_conn The target primitives' connectivity (null for vertices). + /// @param target_count The vertex ids per target primitive (1, 2, or 3). + /// @param d_a The first ids of each emitted pair (output, device). + /// @param d_b The second ids of each emitted pair (output, device). + /// @return The number of candidate pairs emitted. + template + size_t run_traversal( + const LBVH::Impl::DeviceBVH& source, + const LBVH::Impl::DeviceBVH& target, + const index_t* source_conn, + const int source_count, + const index_t* target_conn, + const int target_count, + thrust::device_vector& d_a, + thrust::device_vector& d_b) + { + const int n_source_leaves = source.n_leaves; + const int target_size = static_cast(target.nodes.size()); + if (n_source_leaves == 0 || target_size == 0) { + d_a.clear(); + d_b.clear(); + return 0; + } + const int source_leaf_offset = n_source_leaves - 1; + + int capacity = std::max(1024, 8 * n_source_leaves); + thrust::device_vector d_counter(1); + + int count = 0; + while (true) { + d_a.resize(capacity); + d_b.resize(capacity); + d_counter[0] = 0; + + traverse_kernel + <<>>( + thrust::raw_pointer_cast(source.nodes.data()), + n_source_leaves, source_leaf_offset, + thrust::raw_pointer_cast(target.nodes.data()), target_size, + thrust::raw_pointer_cast(target.rightmost_leaves.data()), + source_conn, source_count, target_conn, target_count, + thrust::raw_pointer_cast(d_a.data()), + thrust::raw_pointer_cast(d_b.data()), + thrust::raw_pointer_cast(d_counter.data()), capacity); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + count = d_counter[0]; // device->host read (also synchronizes) + if (count <= capacity) { + break; // everything fit + } + capacity = count; // exact size now known; the re-run will fit + } + + d_a.resize(count); // shrink to the exact candidate count (keeps data) + d_b.resize(count); + return static_cast(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 d_a The first ids of each candidate pair (device). + /// @param d_b The second ids of each candidate pair (device). + /// @param count The number of candidate pairs. + /// @param accepts_all Whether the user vertex filter accepts every pair. + /// @param can_collide The predicate applied when accepts_all is false. + /// @param out The materialized candidates (appended to). + template + void materialize( + const thrust::device_vector& d_a, + const thrust::device_vector& d_b, + const size_t count, + const bool accepts_all, + const std::function& can_collide, + std::vector& out) + { + if (count == 0) { + return; + } + std::vector h_a(count), h_b(count); + thrust::copy(d_a.begin(), d_a.begin() + count, h_a.begin()); + thrust::copy(d_b.begin(), d_b.begin() + count, h_b.begin()); + + out.reserve(out.size() + count); + if (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]); + } + } + } + } + +} // namespace + +LBVH::LBVH() : ipc::BroadPhase(), m_impl(std::make_unique()) { } + +LBVH::~LBVH() = default; + +LBVH::LBVH(LBVH&&) noexcept = default; +LBVH& LBVH::operator=(LBVH&&) noexcept = default; + +const LBVH::Impl& LBVH::impl() const { return *m_impl; } + +void LBVH::build( + Eigen::ConstRef vertices, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius) +{ + clear(); + + if (vertices.cols() != 3) { + log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); + } + dim = 3; + + const int n_vertices = static_cast(vertices.rows()); + if (n_vertices == 0) { + return; + } + + // Upload vertices as a flat row-major array. + std::vector h_verts(3 * size_t(n_vertices)); + for (int i = 0; i < n_vertices; ++i) { + for (int k = 0; k < 3; ++k) { + h_verts[3 * size_t(i) + k] = vertices(i, k); + } + } + const thrust::device_vector d_verts(h_verts); + + // Build vertex boxes on the device. + thrust::device_vector vbox_min(3 * size_t(n_vertices)); + thrust::device_vector vbox_max(3 * size_t(n_vertices)); + build_vertex_boxes_static_kernel<<< + kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( + thrust::raw_pointer_cast(d_verts.data()), n_vertices, inflation_radius, + thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + build_from_vertex_boxes( + *m_impl, dim, vbox_min, vbox_max, n_vertices, edges, faces); +} + +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()); + + clear(); + + if (vertices_t0.cols() != 3) { + log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); + } + dim = 3; + + const int n_vertices = static_cast(vertices_t0.rows()); + if (n_vertices == 0) { + return; + } + + std::vector h_v0(3 * size_t(n_vertices)); + std::vector h_v1(3 * size_t(n_vertices)); + for (int i = 0; i < n_vertices; ++i) { + for (int k = 0; k < 3; ++k) { + h_v0[3 * size_t(i) + k] = vertices_t0(i, k); + h_v1[3 * size_t(i) + k] = vertices_t1(i, k); + } + } + const thrust::device_vector d_v0(h_v0); + const thrust::device_vector d_v1(h_v1); + + thrust::device_vector vbox_min(3 * size_t(n_vertices)); + thrust::device_vector vbox_max(3 * size_t(n_vertices)); + build_vertex_boxes_dynamic_kernel<<< + kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( + thrust::raw_pointer_cast(d_v0.data()), + thrust::raw_pointer_cast(d_v1.data()), n_vertices, inflation_radius, + thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(vbox_max.data())); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + build_from_vertex_boxes( + *m_impl, dim, vbox_min, vbox_max, n_vertices, edges, faces); +} + +void LBVH::build( + const AABBs& vertex_boxes, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const uint8_t _dim) +{ + clear(); + + if (_dim != 3) { + log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); + } + dim = 3; + + 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]; + } + } + const thrust::device_vector vbox_min(h_min); + const thrust::device_vector vbox_max(h_max); + + build_from_vertex_boxes( + *m_impl, dim, vbox_min, vbox_max, n_vertices, edges, faces); +} + +void LBVH::clear() +{ + ipc::BroadPhase::clear(); + if (m_impl) { + 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. + +namespace { + // Raw device pointer to a connectivity array, or nullptr if empty (a + // vertex primitive has no connectivity array). + const index_t* conn_ptr(const thrust::device_vector& v) + { + return v.empty() ? nullptr : thrust::raw_pointer_cast(v.data()); + } + + // Fill buf with the device connectivity-filtered candidate pairs, then + // materialize them (host) into out, trimming with can_collide when the user + // filter is not accept-all. + template + void detect_host( + const LBVH::Impl::DeviceBVH& source, + const LBVH::Impl::DeviceBVH& target, + const index_t* source_conn, + const int source_count, + const index_t* target_conn, + const int target_count, + LBVH::Impl::DeviceCandidates& buf, + const bool accepts_all, + const std::function& can_collide, + std::vector& out) + { + const size_t count = run_traversal( + source, target, source_conn, source_count, target_conn, + target_count, buf.a, buf.b); + materialize( + buf.a, buf.b, count, accepts_all, can_collide, out); + } + + // Fill buf on the device and return a view of it. + template + LBVH::DeviceCandidateView detect_device( + const LBVH::Impl::DeviceBVH& source, + const LBVH::Impl::DeviceBVH& target, + const index_t* source_conn, + const int source_count, + const index_t* target_conn, + const int target_count, + LBVH::Impl::DeviceCandidates& buf) + { + const size_t count = run_traversal( + source, target, source_conn, source_count, target_conn, + target_count, buf.a, buf.b); + return LBVH::DeviceCandidateView { + count ? thrust::raw_pointer_cast(buf.a.data()) : nullptr, + count ? thrust::raw_pointer_cast(buf.b.data()) : nullptr, count + }; + } +} // namespace + +void LBVH::detect_vertex_vertex_candidates( + std::vector& candidates) const +{ + if (m_impl->vertex_bvh.n_leaves <= 1) { + return; // need at least 2 vertices for a collision + } + detect_host( + m_impl->vertex_bvh, m_impl->vertex_bvh, nullptr, 1, nullptr, 1, + m_impl->vv_candidates, can_vertices_collide.accepts_all(), + [this](size_t a, size_t b) { return can_vertices_collide(a, b); }, + candidates); +} + +void LBVH::detect_edge_vertex_candidates( + std::vector& candidates) const +{ + if (m_impl->edge_bvh.n_leaves == 0 || m_impl->vertex_bvh.n_leaves == 0) { + return; + } + detect_host( + m_impl->edge_bvh, m_impl->vertex_bvh, conn_ptr(m_impl->edges), 2, + nullptr, 1, m_impl->ev_candidates, can_vertices_collide.accepts_all(), + [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 +{ + if (m_impl->edge_bvh.n_leaves <= 1) { + return; // need at least 2 edges for a collision + } + detect_host( + m_impl->edge_bvh, m_impl->edge_bvh, conn_ptr(m_impl->edges), 2, + conn_ptr(m_impl->edges), 2, m_impl->ee_candidates, + can_vertices_collide.accepts_all(), + [this](size_t a, size_t b) { return can_edges_collide(a, b); }, + candidates); +} + +void LBVH::detect_face_vertex_candidates( + std::vector& candidates) const +{ + if (m_impl->face_bvh.n_leaves == 0 || m_impl->vertex_bvh.n_leaves == 0) { + return; + } + // Iterate over the vertices (source) and query the face BVH (target), + // swapping so the emitted pair is (face, vertex). Mirrors ipc::LBVH. + detect_host( + m_impl->vertex_bvh, m_impl->face_bvh, nullptr, 1, + conn_ptr(m_impl->faces), 3, m_impl->fv_candidates, + can_vertices_collide.accepts_all(), + [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 +{ + if (m_impl->edge_bvh.n_leaves == 0 || m_impl->face_bvh.n_leaves == 0) { + return; + } + // Iterate over the faces (source) and query the edge BVH (target), + // swapping so the emitted pair is (edge, face). Mirrors ipc::LBVH. + detect_host( + m_impl->face_bvh, m_impl->edge_bvh, conn_ptr(m_impl->faces), 3, + conn_ptr(m_impl->edges), 2, m_impl->ef_candidates, + can_vertices_collide.accepts_all(), + [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 +{ + if (m_impl->face_bvh.n_leaves <= 1) { + return; // need at least 2 faces for a collision + } + detect_host( + m_impl->face_bvh, m_impl->face_bvh, conn_ptr(m_impl->faces), 3, + conn_ptr(m_impl->faces), 3, m_impl->ff_candidates, + can_vertices_collide.accepts_all(), + [this](size_t a, size_t b) { return can_faces_collide(a, b); }, + candidates); +} + +// --------------------------------------------------------------------------- +// Device-resident candidate accessors. Run the traversal and return a view of +// the connectivity-filtered pairs left on the device (valid until the next +// call on the same type or clear()). For the accept-all filter this is the +// exact candidate set; otherwise it is a superset the caller must trim with +// the user vertex filter. + +LBVH::DeviceCandidateView LBVH::detect_vertex_vertex_candidates_device() const +{ + if (m_impl->vertex_bvh.n_leaves <= 1) { + m_impl->vv_candidates.clear(); + return {}; + } + return detect_device( + m_impl->vertex_bvh, m_impl->vertex_bvh, nullptr, 1, nullptr, 1, + m_impl->vv_candidates); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_vertex_candidates_device() const +{ + if (m_impl->edge_bvh.n_leaves == 0 || m_impl->vertex_bvh.n_leaves == 0) { + m_impl->ev_candidates.clear(); + return {}; + } + return detect_device( + m_impl->edge_bvh, m_impl->vertex_bvh, conn_ptr(m_impl->edges), 2, + nullptr, 1, m_impl->ev_candidates); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_edge_candidates_device() const +{ + if (m_impl->edge_bvh.n_leaves <= 1) { + m_impl->ee_candidates.clear(); + return {}; + } + return detect_device( + m_impl->edge_bvh, m_impl->edge_bvh, conn_ptr(m_impl->edges), 2, + conn_ptr(m_impl->edges), 2, m_impl->ee_candidates); +} + +LBVH::DeviceCandidateView LBVH::detect_face_vertex_candidates_device() const +{ + if (m_impl->face_bvh.n_leaves == 0 || m_impl->vertex_bvh.n_leaves == 0) { + m_impl->fv_candidates.clear(); + return {}; + } + return detect_device( + m_impl->vertex_bvh, m_impl->face_bvh, nullptr, 1, + conn_ptr(m_impl->faces), 3, m_impl->fv_candidates); +} + +LBVH::DeviceCandidateView LBVH::detect_edge_face_candidates_device() const +{ + if (m_impl->edge_bvh.n_leaves == 0 || m_impl->face_bvh.n_leaves == 0) { + m_impl->ef_candidates.clear(); + return {}; + } + return detect_device( + m_impl->face_bvh, m_impl->edge_bvh, conn_ptr(m_impl->faces), 3, + conn_ptr(m_impl->edges), 2, m_impl->ef_candidates); +} + +LBVH::DeviceCandidateView LBVH::detect_face_face_candidates_device() const +{ + if (m_impl->face_bvh.n_leaves <= 1) { + m_impl->ff_candidates.clear(); + return {}; + } + return detect_device( + m_impl->face_bvh, m_impl->face_bvh, conn_ptr(m_impl->faces), 3, + conn_ptr(m_impl->faces), 3, m_impl->ff_candidates); +} + +// --------------------------------------------------------------------------- +// Host-side can_*_collide filters (mesh connectivity + user vertex filter). +// Mirror ipc::LBVH's overrides, backed by the host connectivity copies. + +bool LBVH::can_edge_vertex_collide(size_t ei, size_t vi) const +{ + const auto& [e0i, e1i] = m_impl->h_edge_vertex_ids[ei]; + return vi != e0i && vi != e1i + && (can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i)); +} + +bool LBVH::can_edges_collide(size_t eai, size_t ebi) const +{ + const auto& [ea0i, ea1i] = m_impl->h_edge_vertex_ids[eai]; + const auto& [eb0i, eb1i] = m_impl->h_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)); +} + +bool LBVH::can_face_vertex_collide(size_t fi, size_t vi) const +{ + const auto& [f0i, f1i, f2i] = m_impl->h_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)); +} + +bool LBVH::can_edge_face_collide(size_t ei, size_t fi) const +{ + const auto& [e0i, e1i] = m_impl->h_edge_vertex_ids[ei]; + const auto& [f0i, f1i, f2i] = m_impl->h_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)); +} + +bool LBVH::can_faces_collide(size_t fai, size_t fbi) const +{ + const auto& [fa0i, fa1i, fa2i] = m_impl->h_face_vertex_ids[fai]; + const auto& [fb0i, fb1i, fb2i] = m_impl->h_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)); +} + +size_t LBVH::num_vertex_nodes() const +{ + return m_impl->vertex_bvh.nodes.size(); +} + +size_t LBVH::num_edge_nodes() const { return m_impl->edge_bvh.nodes.size(); } + +size_t LBVH::num_face_nodes() const { return m_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(m_impl->vertex_bvh, nodes, rightmost_leaves); +} + +void LBVH::edge_nodes_to_host( + ipc::LBVH::Nodes& nodes, ipc::LBVH::RightmostLeaves& rightmost_leaves) const +{ + to_host(m_impl->edge_bvh, nodes, rightmost_leaves); +} + +void LBVH::face_nodes_to_host( + ipc::LBVH::Nodes& nodes, ipc::LBVH::RightmostLeaves& rightmost_leaves) const +{ + to_host(m_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..bda7b444e --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh.hpp @@ -0,0 +1,175 @@ +#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. +/// +/// 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. +class LBVH : public ipc::BroadPhase { +public: + LBVH(); + ~LBVH(); + + LBVH(LBVH&&) noexcept; + LBVH& operator=(LBVH&&) 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 are valid until + /// the next detect_*_device() call on the same type or clear(). + 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. + /// @param vertices Vertex positions (rowwise, |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. + /// @param vertices_t0 Starting vertex positions (rowwise, |V| × 3). + /// @param vertices_t1 Ending vertex positions (rowwise, |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_t0, + Eigen::ConstRef vertices_t1, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const double inflation_radius = 0) override; + + /// @brief Build the broad phase from precomputed host vertex AABBs. + /// The vertex boxes are uploaded; edge/face boxes and all BVHs are built on + /// the device. + /// @param vertex_boxes Precomputed vertex AABBs. + /// @param edges Collision mesh edges. + /// @param faces Collision mesh faces. + /// @param dim Dimension of the simulation (must be 3). + void build( + const AABBs& vertex_boxes, + Eigen::ConstRef edges, + Eigen::ConstRef faces, + const uint8_t dim) override; + + /// @brief Clear any built data. + 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. + + 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. + + 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; + + // Pimpl pattern to keep CUDA types out of this header. The Impl is + // defined in lbvh_impl.cuh for use by the ipc::cuda implementation files + // (.cu) only. + struct Impl; + const Impl& impl() const; + +protected: + // 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. + 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: + std::unique_ptr m_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..b896540a5 --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh_impl.cuh @@ -0,0 +1,94 @@ +// 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 + +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 { + thrust::device_vector nodes; + thrust::device_vector rightmost_leaves; + int n_leaves = 0; + + void clear() + { + nodes.clear(); + rightmost_leaves.clear(); + n_leaves = 0; + } + }; + + 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. + struct DeviceCandidates { + thrust::device_vector a; + thrust::device_vector b; + + void clear() + { + a.clear(); + b.clear(); + } + }; + + DeviceCandidates vv_candidates; + DeviceCandidates ev_candidates; + DeviceCandidates ee_candidates; + DeviceCandidates fv_candidates; + DeviceCandidates ef_candidates; + DeviceCandidates ff_candidates; + + // Mesh connectivity, uploaded once and used by the device traversal's + // shared-vertex (connectivity) filter. Flat row-major: + // edges = 2 * n_edges, faces = 3 * n_faces. + thrust::device_vector edges; + thrust::device_vector faces; + + // Host copies of the connectivity, used by the host-side can_*_collide + // filters applied to the device-emitted candidate pairs. + std::vector> h_edge_vertex_ids; + std::vector> h_face_vertex_ids; + + void clear() + { + vertex_bvh.clear(); + edge_bvh.clear(); + face_bvh.clear(); + edges.clear(); + faces.clear(); + h_edge_vertex_ids.clear(); + h_face_vertex_ids.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/collision_filter.hpp b/src/ipc/collision_filter.hpp index 16fe1cd13..40c46a745 100644 --- a/src/ipc/collision_filter.hpp +++ b/src/ipc/collision_filter.hpp @@ -32,7 +32,11 @@ class CollisionFilter { // ── Construction ───────────────────────────────────────────────────────── /// @brief Default filter: accept all pairs. - CollisionFilter() : m_fn([](size_t, size_t) { return true; }) { } + CollisionFilter() + : m_fn([](size_t, size_t) { return true; }) + , m_accepts_all(true) + { + } /// @brief Construct from any callable bool(size_t, size_t). /// @note Disabled when Fn is CollisionFilter itself to avoid shadowing @@ -59,6 +63,13 @@ class CollisionFilter { /// @brief Implicit conversion to std::function. operator std::function() const { return m_fn; } + /// @brief Whether this filter trivially accepts every pair. + /// @return true only for the default-constructed (accept-all) filter; + /// conservatively false for any user-supplied or composed filter. + /// @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_accepts_all; } + // ── Composition ────────────────────────────────────────────────────────── /// @brief Union: accept if EITHER filter passes. @@ -100,6 +111,9 @@ class CollisionFilter { private: std::function m_fn; + /// @brief True only for the default (accept-all) filter. Any callable- or + /// composition-constructed filter leaves this false (conservative). + bool m_accepts_all = false; }; // ───────────────────────────────────────────────────────────────────────────── diff --git a/tests/src/tests/broad_phase/CMakeLists.txt b/tests/src/tests/broad_phase/CMakeLists.txt index 1806f2a49..9cb26926b 100644 --- a/tests/src/tests/broad_phase/CMakeLists.txt +++ b/tests/src/tests/broad_phase/CMakeLists.txt @@ -15,6 +15,12 @@ set(SOURCES brute_force_comparison.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/test_gpu_lbvh.cu b/tests/src/tests/broad_phase/test_gpu_lbvh.cu new file mode 100644 index 000000000..028d1afd0 --- /dev/null +++ b/tests/src/tests/broad_phase/test_gpu_lbvh.cu @@ -0,0 +1,309 @@ +// 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). + +#include + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include + +#include +#include + +#include +#include + +#include + +#include +#include + +using namespace ipc; + +namespace { + +bool has_cuda_device() +{ + int n = 0; + return cudaGetDeviceCount(&n) == cudaSuccess && n > 0; +} + +bool is_aabb_union( + const LBVH::Node& parent, + const LBVH::Node& child_a, + const LBVH::Node& child_b) +{ + const Eigen::Array3d cmin = + child_a.aabb_min.min(child_b.aabb_min).cast(); + const Eigen::Array3d cmax = + child_a.aabb_max.max(child_b.aabb_max).cast(); + constexpr float EPS = 1e-4f; + return (abs(parent.aabb_max.cast() - cmax) < EPS).all() + && (abs(parent.aabb_min.cast() - cmin) < EPS).all(); +} + +// Recursively verify reachability (each node visited exactly once) and that +// every internal node's AABB is the union of its children's. Collects the leaf +// primitive ids that are reached. +void traverse_and_check( + 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()); + CHECK(!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); + CHECK(is_aabb_union(node, child_a, child_b)); + } + traverse_and_check(nodes, node.left, visited, reached_leaves); + traverse_and_check(nodes, node.right, visited, reached_leaves); +} + +// Validate one device-built tree (copied back to the host) against the +// corresponding CPU-built node array. +void check_tree(const LBVH::Nodes& nodes, const LBVH::Nodes& cpu_nodes) +{ + if (nodes.size() <= 1) { + return; // single-node trees are not exercised here + } + REQUIRE(nodes.size() == cpu_nodes.size()); + REQUIRE(nodes.size() % 2 == 1); // 2n - 1 + const size_t n_leaves = (nodes.size() + 1) / 2; + + // -- Structural validity: reachable-once + AABB unions. -- + std::vector visited(nodes.size(), false); + std::vector reached_leaves; + traverse_and_check(nodes, 0, visited, reached_leaves); + CHECK( + std::all_of(visited.begin(), visited.end(), [](bool v) { return v; })); + + // -- Leaf set must be exactly {0, ..., n_leaves - 1}. -- + 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)); + } + + // -- Root AABB must equal the CPU root AABB (an order-independent union of + // identically-inflated boxes). -- + constexpr float EPS = 1e-4f; + CHECK((abs(nodes[0].aabb_min.cast() + - cpu_nodes[0].aabb_min.cast()) + < EPS) + .all()); + CHECK((abs(nodes[0].aabb_max.cast() + - cpu_nodes[0].aabb_max.cast()) + < EPS) + .all()); +} + +} // namespace + +TEST_CASE("GPU LBVH build", "[broad_phase][lbvh][cuda][gpu]") +{ + if (!has_cuda_device()) { + SKIP("No CUDA device available; kernels compiled but not executed."); + } + + 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); + check_tree(nodes, cpu_lbvh.vertex_nodes()); + } + SECTION("edges") + { + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + check_tree(nodes, cpu_lbvh.edge_nodes()); + } + SECTION("faces") + { + gpu_lbvh.face_nodes_to_host(nodes, rightmost); + check_tree(nodes, cpu_lbvh.face_nodes()); + } + + // clear() empties the device trees. + gpu_lbvh.clear(); + gpu_lbvh.vertex_nodes_to_host(nodes, rightmost); + CHECK(nodes.empty()); +} + +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 detect candidates", "[broad_phase][lbvh][cuda][gpu]") +{ + if (!has_cuda_device()) { + SKIP("No CUDA device available; kernels compiled but not executed."); + } + + 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). + 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); + } + { + 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]") +{ + if (!has_cuda_device()) { + SKIP("No CUDA device available; kernels compiled but not executed."); + } + + 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); + } + { + 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); + } +} + +#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..aa802bf12 100644 --- a/tests/src/tests/broad_phase/test_lbvh.cpp +++ b/tests/src/tests/broad_phase/test_lbvh.cpp @@ -7,6 +7,10 @@ #include #include +#ifdef IPC_TOOLKIT_WITH_CUDA +#include +#endif + #include #include @@ -344,6 +348,24 @@ TEST_CASE( lbvh->detect_edge_edge_candidates(ee_candidates); return ee_candidates.size(); }; + +#ifdef IPC_TOOLKIT_WITH_CUDA + 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(); + }; +#endif } TEST_CASE("Benchmark LBVH::build", "[!benchmark][broad_phase][lbvh]") @@ -388,5 +410,20 @@ TEST_CASE("Benchmark LBVH::build", "[!benchmark][broad_phase][lbvh]") vertices_t0, vertices_t1, edges, faces, inflation_radius); return lbvh->edge_nodes().size(); }; + +#ifdef IPC_TOOLKIT_WITH_CUDA + 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(fmt::format("cuda::LBVH::build [{}]", scene)) + { + gpu_lbvh.build( + vertices_t0, vertices_t1, edges, faces, inflation_radius); + return gpu_lbvh.num_edge_nodes(); + }; +#endif } } \ No newline at end of file From dcc32a701d13330d5e7c7463fff5b3e3ecd6cf8e Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 24 Jul 2026 18:33:25 -0700 Subject: [PATCH 03/10] Cache candidate-buffer capacity across calls in ipc::cuda::LBVH detect_*_candidates sized its output buffer from a fresh max(1024, 8 * n_source_leaves) guess on every single call, with nothing remembering what a prior call actually needed. Checking the real candidate counts on the CPU LBVH (proven exactly equal to the GPU's) showed 5 of 7 benchmarked meshes overflow that guess by up to 17x, so nearly every real mesh silently paid for two full kernel dispatches on every call: one that discovers the buffer is too small, then a full re-traversal at the corrected size. Add a predicted_capacity field to LBVH::Impl::DeviceCandidates (one per candidate type) that persists the largest count ever observed and seeds the next call's guess. It is deliberately not reset by clear(), since build() calls clear() every timestep and the hint must survive that or it never helps; it only ever grows for the object's lifetime, mirroring the predicted_*_candidates_size pattern already used by the (Slang) vulkan branch's LBVH. Also add a logger().warn() on overflow, matching that same branch, so a retry is no longer silent. Validated on artemis (RTX 3070): [lbvh][cuda] unchanged at 150517 assertions. Re-benchmarked detect_edge_edge_candidates against the CPU LBVH: the 2 meshes that never overflowed are byte-identical before/after as expected; the 5 that did are 13-29% faster (e.g. Rod-Twist 15.8ms -> 11.2ms, Puffer-Ball 1.097s -> 0.912s), widening the GPU's margin over the CPU across the board (e.g. Rod-Twist 1.23x -> 1.73x, Puffer-Ball 1.47x -> 1.76x faster than CPU). Co-Authored-By: Claude Sonnet 5 --- src/ipc/broad_phase/cuda/lbvh.cu | 53 ++++++++++++++++---------- src/ipc/broad_phase/cuda/lbvh_impl.cuh | 11 ++++++ 2 files changed, 43 insertions(+), 21 deletions(-) diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index d2ae174da..194ea7f13 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -794,9 +794,12 @@ namespace { } /// @brief Run the device traversal of the target BVH by the source leaves, - /// leaving the connectivity-filtered candidate pairs device-resident in the - /// output buffers (resized to the exact count). Grows the buffer and - /// re-runs once if the first pass overflows. + /// leaving the connectivity-filtered candidate pairs device-resident in + /// buf.a/buf.b (resized to the exact count). The initial buffer size is + /// seeded from buf.predicted_capacity (the largest count ever observed for + /// this type on this object), so only the first call -- or a call whose + /// count exceeds every prior call -- pays the overflow-and-retry cost; + /// every other call fits on the first pass. /// @tparam triangular Self-collision: skip subtrees fully left of the query. /// @tparam swap_order Emit (target_prim, source_prim) instead. /// @param source The BVH whose leaves are the queries. @@ -805,8 +808,7 @@ namespace { /// @param source_count The vertex ids per source primitive (1, 2, or 3). /// @param target_conn The target primitives' connectivity (null for vertices). /// @param target_count The vertex ids per target primitive (1, 2, or 3). - /// @param d_a The first ids of each emitted pair (output, device). - /// @param d_b The second ids of each emitted pair (output, device). + /// @param buf The output candidate buffer and capacity hint (in/out). /// @return The number of candidate pairs emitted. template size_t run_traversal( @@ -816,25 +818,26 @@ namespace { const int source_count, const index_t* target_conn, const int target_count, - thrust::device_vector& d_a, - thrust::device_vector& d_b) + LBVH::Impl::DeviceCandidates& buf) { const int n_source_leaves = source.n_leaves; const int target_size = static_cast(target.nodes.size()); if (n_source_leaves == 0 || target_size == 0) { - d_a.clear(); - d_b.clear(); + buf.a.clear(); + buf.b.clear(); return 0; } const int source_leaf_offset = n_source_leaves - 1; - int capacity = std::max(1024, 8 * n_source_leaves); + size_t capacity = std::max( + buf.predicted_capacity, + static_cast(std::max(1024, 8 * n_source_leaves))); thrust::device_vector d_counter(1); int count = 0; while (true) { - d_a.resize(capacity); - d_b.resize(capacity); + buf.a.resize(capacity); + buf.b.resize(capacity); d_counter[0] = 0; traverse_kernel @@ -844,20 +847,28 @@ namespace { thrust::raw_pointer_cast(target.nodes.data()), target_size, thrust::raw_pointer_cast(target.rightmost_leaves.data()), source_conn, source_count, target_conn, target_count, - thrust::raw_pointer_cast(d_a.data()), - thrust::raw_pointer_cast(d_b.data()), - thrust::raw_pointer_cast(d_counter.data()), capacity); + thrust::raw_pointer_cast(buf.a.data()), + thrust::raw_pointer_cast(buf.b.data()), + thrust::raw_pointer_cast(d_counter.data()), + static_cast(capacity)); IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); count = d_counter[0]; // device->host read (also synchronizes) - if (count <= capacity) { + if (static_cast(count) <= capacity) { break; // everything fit } - capacity = count; // exact size now known; the re-run will 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 } - d_a.resize(count); // shrink to the exact candidate count (keeps data) - d_b.resize(count); + buf.predicted_capacity = std::max(buf.predicted_capacity, capacity); + buf.a.resize(count); // shrink to the exact candidate count (keeps data) + buf.b.resize(count); return static_cast(count); } @@ -1072,7 +1083,7 @@ namespace { { const size_t count = run_traversal( source, target, source_conn, source_count, target_conn, - target_count, buf.a, buf.b); + target_count, buf); materialize( buf.a, buf.b, count, accepts_all, can_collide, out); } @@ -1090,7 +1101,7 @@ namespace { { const size_t count = run_traversal( source, target, source_conn, source_count, target_conn, - target_count, buf.a, buf.b); + target_count, buf); return LBVH::DeviceCandidateView { count ? thrust::raw_pointer_cast(buf.a.data()) : nullptr, count ? thrust::raw_pointer_cast(buf.b.data()) : nullptr, count diff --git a/src/ipc/broad_phase/cuda/lbvh_impl.cuh b/src/ipc/broad_phase/cuda/lbvh_impl.cuh index b896540a5..30dd20152 100644 --- a/src/ipc/broad_phase/cuda/lbvh_impl.cuh +++ b/src/ipc/broad_phase/cuda/lbvh_impl.cuh @@ -46,10 +46,21 @@ struct LBVH::Impl { thrust::device_vector a; thrust::device_vector b; + /// @brief Largest candidate count ever observed for this type on this + /// object, used to size the next traversal's output buffer so repeated + /// calls (e.g. one per Newton iteration, or one per build() at a new + /// timestep) don't pay the overflow-and-retry cost every time -- only + /// the first time, or when the count grows past every prior call. + /// Deliberately NOT reset by clear() (see below): build() calls + /// clear() every timestep, and this hint must survive that so the + /// learned size doesn't need re-discovering each time. + size_t predicted_capacity = 0; + void clear() { a.clear(); b.clear(); + // predicted_capacity is intentionally left untouched. } }; From abd5150f44eca3c988f3af09e11713f8130a3752 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 24 Jul 2026 18:39:58 -0700 Subject: [PATCH 04/10] Register ipc::cuda::LBVH in create_broad_phase Add BroadPhaseMethod::LBVH_CUDA, appended after SWEEP_AND_TINIEST_QUEUE to keep the existing enum values stable (the "Create broad phase" test casts consecutive integers to BroadPhaseMethod). The factory case mirrors the SWEEP_AND_TINIEST_QUEUE case exactly: returns ipc::cuda::LBVH under IPC_TOOLKIT_WITH_CUDA, otherwise throws with a message naming the CMake option to enable. Not added to tests/src/tests/utils.cpp's broad_phases() / BroadPhaseGenerator (used by most generic cross-broad-phase comparison tests): several of those exercise 2D meshes, and ipc::cuda::LBVH::build() currently throws on non-3D input (v1 scope), unlike SweepAndTiniestQueue which silently upgrades 2D to 3D via to_X3d() before building. Adding it there would break those tests immediately; left for a follow-up if 2D parity is wanted. Validated: host (non-CUDA) build passes "Create broad phase" (5 assertions, count unchanged). Artemis (CUDA, RTX 3070): same test passes with the bumped count (7 assertions); [lbvh][cuda] suite unaffected (150517 assertions, no regression). Co-Authored-By: Claude Sonnet 5 --- src/ipc/broad_phase/create_broad_phase.cpp | 8 ++++++++ src/ipc/broad_phase/create_broad_phase.hpp | 3 ++- tests/src/tests/broad_phase/test_broad_phase.cpp | 2 +- 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/ipc/broad_phase/create_broad_phase.cpp b/src/ipc/broad_phase/create_broad_phase.cpp index c812a1170..0a12b9e4e 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 @@ -29,6 +30,13 @@ create_broad_phase(const BroadPhaseMethod& broad_phase_method) #else 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 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..3ae624caa 100644 --- a/src/ipc/broad_phase/create_broad_phase.hpp +++ b/src/ipc/broad_phase/create_broad_phase.hpp @@ -11,7 +11,8 @@ enum class BroadPhaseMethod : uint8_t { SPATIAL_HASH, LBVH, SWEEP_AND_PRUNE, - SWEEP_AND_TINIEST_QUEUE + SWEEP_AND_TINIEST_QUEUE, + LBVH_CUDA }; std::shared_ptr diff --git a/tests/src/tests/broad_phase/test_broad_phase.cpp b/tests/src/tests/broad_phase/test_broad_phase.cpp index 37ab71545..451b2779f 100644 --- a/tests/src/tests/broad_phase/test_broad_phase.cpp +++ b/tests/src/tests/broad_phase/test_broad_phase.cpp @@ -296,7 +296,7 @@ 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; + uint8_t n_broad_phase_methods = 7; #else uint8_t n_broad_phase_methods = 5; #endif From d6b7df0566a5189b4b5e9eb04ce61a705912bf56 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Sun, 26 Jul 2026 13:19:47 -0700 Subject: [PATCH 05/10] Add 2D support to ipc::cuda::LBVH, matching the CPU LBVH exactly build() previously hard-coded dim = 3 and threw for non-3-column input. The vertex upload also unconditionally read 3 components per vertex, which would read out of bounds on a 2-column matrix -- dimension support was blocked below the validation check, not just at it. Mirror the CPU ipc::LBVH's actual semantics instead of the simpler upgrade-to-3D-via-to_X3d approach SweepAndTiniestQueue uses. The key subtlety: ipc::AABB's constructor zero-initializes its 3-wide array and only assigns the first `dim` components from the already-inflated input, so a 2D box's z bound is an exact, uninflated 0.0 -- not nextafter(0 +/- inflation_radius, ...). build_vertex_boxes_{static, dynamic}_kernel now take dim and, for components past it, write a hard 0.0 instead of running the inflation formula, matching that exactly. Vertex upload now sizes to dim * n instead of a fixed 3 * n. All three build() overloads relax to assert(dim == 2 || dim == 3) (matching the CPU's debug-only assert, not a throw) and set dim from the real input. The Morton-code kernel's dim == 2 branch already existed (copied from the CPU when first written) and needed no change; the edge/face box union kernels and the Apetrei hierarchy build are dim-agnostic and untouched. Add "GPU LBVH 2D build and detect" using the same mesh-2D CSV data as the CPU's own 2D test: checks vertex/edge BVH structural and root-AABB parity, plus exact detect_edge_vertex_candidates parity against the CPU LBVH (the only candidate type meaningful in 2D). Validated on artemis (RTX 3070): [lbvh][cuda] now 152785 assertions across 4 test cases (was 150517/3) -- the existing 3D paths are unregressed and the new 2D path matches the CPU exactly. Co-Authored-By: Claude Sonnet 5 --- src/ipc/broad_phase/cuda/lbvh.cu | 99 +++++++++++--------- src/ipc/broad_phase/cuda/lbvh.hpp | 2 +- tests/src/tests/broad_phase/test_gpu_lbvh.cu | 45 +++++++++ 3 files changed, 103 insertions(+), 43 deletions(-) diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index 194ea7f13..892d03b1b 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -80,12 +80,20 @@ namespace { // double bounds are nudged outward with nextafter so the box is // conservative. (The leaf nodes later apply a second float-nextafter in // build_hierarchy_kernel, matching assign_inflated_aabb.) + // + // For dim == 2 input, ipc::AABB always stores a 3-wide array whose z + // component is zero-initialized and never touched by conservative_inflation + // (only the first `dim` components of the constructor argument are + // assigned) -- so the z bound is an exact, uninflated 0.0, not + // nextafter(0 +/- inflation_radius, ...). Replicate that exactly: for + // k >= dim, write a hard 0.0 instead of inflating. __global__ void build_vertex_boxes_static_kernel( - const double* __restrict__ vertices, // 3 * n, row-major + const double* __restrict__ vertices, // dim * n, row-major const int n, + const int dim, const double inflation_radius, - double* __restrict__ box_min, + double* __restrict__ box_min, // always 3 * n, row-major double* __restrict__ box_max) { const int i = blockIdx.x * blockDim.x + threadIdx.x; @@ -94,18 +102,24 @@ namespace { } #pragma unroll for (int k = 0; k < 3; ++k) { - const double v = vertices[3 * i + k]; - box_min[3 * i + k] = nextafter(v - inflation_radius, -INFINITY); - box_max[3 * i + k] = nextafter(v + inflation_radius, INFINITY); + if (k < dim) { + const double v = vertices[dim * i + k]; + box_min[3 * i + k] = nextafter(v - inflation_radius, -INFINITY); + box_max[3 * i + k] = nextafter(v + inflation_radius, INFINITY); + } else { + box_min[3 * i + k] = 0.0; + box_max[3 * i + k] = 0.0; + } } } __global__ void build_vertex_boxes_dynamic_kernel( - const double* __restrict__ vertices_t0, // 3 * n, row-major - const double* __restrict__ vertices_t1, // 3 * n, row-major + const double* __restrict__ vertices_t0, // dim * n, row-major + const double* __restrict__ vertices_t1, // dim * n, row-major const int n, + const int dim, const double inflation_radius, - double* __restrict__ box_min, + double* __restrict__ box_min, // always 3 * n, row-major double* __restrict__ box_max) { const int i = blockIdx.x * blockDim.x + threadIdx.x; @@ -114,14 +128,20 @@ namespace { } #pragma unroll for (int k = 0; k < 3; ++k) { - const double a = vertices_t0[3 * i + k]; - const double b = vertices_t1[3 * i + k]; - // union of the two inflated point boxes; nextafter is monotonic so - // min(nextafter(a),nextafter(b)) == nextafter(min(a,b)). - box_min[3 * i + k] = - nextafter(fmin(a, b) - inflation_radius, -INFINITY); - box_max[3 * i + k] = - nextafter(fmax(a, b) + inflation_radius, INFINITY); + if (k < dim) { + const double a = vertices_t0[dim * i + k]; + const double b = vertices_t1[dim * i + k]; + // union of the two inflated point boxes; nextafter is + // monotonic so min(nextafter(a),nextafter(b)) == + // nextafter(min(a,b)). + box_min[3 * i + k] = + nextafter(fmin(a, b) - inflation_radius, -INFINITY); + box_max[3 * i + k] = + nextafter(fmax(a, b) + inflation_radius, INFINITY); + } else { + box_min[3 * i + k] = 0.0; + box_max[3 * i + k] = 0.0; + } } } @@ -931,32 +951,31 @@ void LBVH::build( { clear(); - if (vertices.cols() != 3) { - log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); - } - dim = 3; + assert(vertices.cols() == 2 || vertices.cols() == 3); + dim = static_cast(vertices.cols()); const int n_vertices = static_cast(vertices.rows()); if (n_vertices == 0) { return; } - // Upload vertices as a flat row-major array. - std::vector h_verts(3 * size_t(n_vertices)); + // Upload vertices as a flat row-major array (dim components per vertex; + // no padding -- the box kernel below fills the unused z for 2D input). + std::vector h_verts(size_t(dim) * size_t(n_vertices)); for (int i = 0; i < n_vertices; ++i) { - for (int k = 0; k < 3; ++k) { - h_verts[3 * size_t(i) + k] = vertices(i, k); + for (int k = 0; k < dim; ++k) { + h_verts[size_t(dim) * size_t(i) + k] = vertices(i, k); } } const thrust::device_vector d_verts(h_verts); - // Build vertex boxes on the device. + // Build vertex boxes on the device (always 3-wide storage). thrust::device_vector vbox_min(3 * size_t(n_vertices)); thrust::device_vector vbox_max(3 * size_t(n_vertices)); build_vertex_boxes_static_kernel<<< kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( - thrust::raw_pointer_cast(d_verts.data()), n_vertices, inflation_radius, - thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(d_verts.data()), n_vertices, dim, + inflation_radius, thrust::raw_pointer_cast(vbox_min.data()), thrust::raw_pointer_cast(vbox_max.data())); IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); @@ -976,22 +995,20 @@ void LBVH::build( clear(); - if (vertices_t0.cols() != 3) { - log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); - } - dim = 3; + 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; } - std::vector h_v0(3 * size_t(n_vertices)); - std::vector h_v1(3 * size_t(n_vertices)); + std::vector h_v0(size_t(dim) * size_t(n_vertices)); + std::vector h_v1(size_t(dim) * size_t(n_vertices)); for (int i = 0; i < n_vertices; ++i) { - for (int k = 0; k < 3; ++k) { - h_v0[3 * size_t(i) + k] = vertices_t0(i, k); - h_v1[3 * size_t(i) + k] = vertices_t1(i, k); + for (int k = 0; k < dim; ++k) { + h_v0[size_t(dim) * size_t(i) + k] = vertices_t0(i, k); + h_v1[size_t(dim) * size_t(i) + k] = vertices_t1(i, k); } } const thrust::device_vector d_v0(h_v0); @@ -1002,8 +1019,8 @@ void LBVH::build( build_vertex_boxes_dynamic_kernel<<< kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( thrust::raw_pointer_cast(d_v0.data()), - thrust::raw_pointer_cast(d_v1.data()), n_vertices, inflation_radius, - thrust::raw_pointer_cast(vbox_min.data()), + thrust::raw_pointer_cast(d_v1.data()), n_vertices, dim, + inflation_radius, thrust::raw_pointer_cast(vbox_min.data()), thrust::raw_pointer_cast(vbox_max.data())); IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); @@ -1019,10 +1036,8 @@ void LBVH::build( { clear(); - if (_dim != 3) { - log_and_throw_error("ipc::cuda::LBVH currently supports 3D only!"); - } - dim = 3; + assert(_dim == 2 || _dim == 3); + dim = _dim; const int n_vertices = static_cast(vertex_boxes.size()); if (n_vertices == 0) { diff --git a/src/ipc/broad_phase/cuda/lbvh.hpp b/src/ipc/broad_phase/cuda/lbvh.hpp index bda7b444e..3a19d5c31 100644 --- a/src/ipc/broad_phase/cuda/lbvh.hpp +++ b/src/ipc/broad_phase/cuda/lbvh.hpp @@ -80,7 +80,7 @@ class LBVH : public ipc::BroadPhase { /// @param vertex_boxes Precomputed vertex AABBs. /// @param edges Collision mesh edges. /// @param faces Collision mesh faces. - /// @param dim Dimension of the simulation (must be 3). + /// @param dim Dimension of the simulation (2 or 3). void build( const AABBs& vertex_boxes, Eigen::ConstRef edges, diff --git a/tests/src/tests/broad_phase/test_gpu_lbvh.cu b/tests/src/tests/broad_phase/test_gpu_lbvh.cu index 028d1afd0..a05a49f57 100644 --- a/tests/src/tests/broad_phase/test_gpu_lbvh.cu +++ b/tests/src/tests/broad_phase/test_gpu_lbvh.cu @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -306,4 +307,48 @@ TEST_CASE( } } +// 2D input has no faces; ipc::AABB zero-pads the unused z component without +// inflating it (see build_vertex_boxes_{static,dynamic}_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]") +{ + if (!has_cuda_device()) { + SKIP("No CUDA device available; kernels compiled but not executed."); + } + + 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); + check_tree(nodes, cpu_lbvh.vertex_nodes()); + gpu_lbvh.edge_nodes_to_host(nodes, rightmost); + check_tree(nodes, cpu_lbvh.edge_nodes()); + + // -- 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()); +} + #endif // IPC_TOOLKIT_WITH_CUDA From 2ba97811ff4a690298d48e4ad506703d1f66744f Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 17:16:45 -0400 Subject: [PATCH 06/10] Fix CUDA LBVH errors on MSVC nextafter(double, float) does not exist on the device, so use nextafter(double, double) with a constexpr for positive and negative infinity. INFINITY is a float macro, so nextafter(double, INFINITY) resolves to the host-only std::nextafter promotion template instead of CUDA's __device__ nextafter(double, double). --- src/ipc/broad_phase/cuda/lbvh.cu | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index 892d03b1b..0a7b71ee9 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -87,6 +87,13 @@ namespace { // assigned) -- so the z bound is an exact, uninflated 0.0, not // nextafter(0 +/- inflation_radius, ...). Replicate that exactly: for // k >= dim, write a hard 0.0 instead of inflating. + // + // The direction arguments must be doubles: INFINITY is a float macro, so + // nextafter(double, INFINITY) resolves to the host-only + // std::nextafter promotion template instead of CUDA's + // __device__ nextafter(double, double). + constexpr double POS_INF = std::numeric_limits::infinity(); + constexpr double NEG_INF = -POS_INF; __global__ void build_vertex_boxes_static_kernel( const double* __restrict__ vertices, // dim * n, row-major @@ -104,8 +111,8 @@ namespace { for (int k = 0; k < 3; ++k) { if (k < dim) { const double v = vertices[dim * i + k]; - box_min[3 * i + k] = nextafter(v - inflation_radius, -INFINITY); - box_max[3 * i + k] = nextafter(v + inflation_radius, INFINITY); + box_min[3 * i + k] = nextafter(v - inflation_radius, NEG_INF); + box_max[3 * i + k] = nextafter(v + inflation_radius, POS_INF); } else { box_min[3 * i + k] = 0.0; box_max[3 * i + k] = 0.0; @@ -135,9 +142,9 @@ namespace { // monotonic so min(nextafter(a),nextafter(b)) == // nextafter(min(a,b)). box_min[3 * i + k] = - nextafter(fmin(a, b) - inflation_radius, -INFINITY); + nextafter(fmin(a, b) - inflation_radius, NEG_INF); box_max[3 * i + k] = - nextafter(fmax(a, b) + inflation_radius, INFINITY); + nextafter(fmax(a, b) + inflation_radius, POS_INF); } else { box_min[3 * i + k] = 0.0; box_max[3 * i + k] = 0.0; From 419a82cf05bfbd1437792475c95d3afe7b5a02a4 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 14:43:30 -0400 Subject: [PATCH 07/10] Fix the GPU LBVH traversal stack corruption on sm_120 prim_shares_vertex used two runtime-indexed index_t[3] locals. Runtime indexing forces them into local memory, and on sm_120 ptxas sized traverse_kernel's frame at 0x110 bytes while basing those arrays at frame+0x100 -- 16 bytes of room for 24 bytes of object, on top of the 0x100-byte traversal stack based at frame+0. Writes landed on stack[0..1] and destroyed the INVALID_POINTER sentinel the descent loop terminates on, so the traversal popped past the bottom of the stack and read stack[-1]. The result was cudaErrorIllegalAddress, which surfaces as an apparent hang: the driver spins in the candidate-counter readback, and the poisoned context makes every later GPU test look stuck too. Hold the vertex ids in scalars instead, filling unused slots from slot 0 so every comparison stays well defined. The frame drops to 0x100 (exactly the stack) and local traffic to the 3 stack accesses. Scope of the miscompile: sm_120 only. sm_75/86/89 allocate 0x120 as expected, identically with -rdc=true and -rdc=false, and the driver's own JIT (CUDA 13.3) reproduces the 272 vs 288 split, so it is neither an -rdc nor a 12.8 artifact. Building the unfixed source as compute_89 PTX and JIT-ing onto the sm_120 device passes clean. A provably bounded index does not help, so this is not licensed by the latent UB. Tests: [gpu] ~[!benchmark] passes (158705 assertions, 28 cases) and compute-sanitizer memcheck reports 0 errors on [lbvh][gpu]; both faulted before. Co-Authored-By: Claude Opus 5 (1M context) --- src/ipc/broad_phase/cuda/lbvh.cu | 41 ++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index 0a7b71ee9..3c6ad8fbb 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -661,30 +661,35 @@ namespace { const index_t* __restrict__ conn_b, const int count_b) { - index_t ids_a[3]; - index_t ids_b[3]; + // Use scalars, not arrays. Runtime-indexed local arrays force local + // memory allocation, causing stack corruption on ptxas (sm_120) when + // the frame overflows into the traversal stack sentinel. Keeping values + // in registers limits the frame size to 0x100 and prevents invalid + // memory writes. Unused slots are filled from slot 0 for well-defined + // comparisons. + index_t a0, a1, a2; if (conn_a == nullptr) { - ids_a[0] = p_a; + a0 = a1 = a2 = p_a; } else { - for (int k = 0; k < count_a; ++k) { - ids_a[k] = conn_a[count_a * p_a + k]; - } + const index_t* row = conn_a + count_a * p_a; + a0 = row[0]; + a1 = count_a > 1 ? row[1] : a0; + a2 = count_a > 2 ? row[2] : a0; } + + index_t b0, b1, b2; if (conn_b == nullptr) { - ids_b[0] = p_b; + b0 = b1 = b2 = p_b; } else { - for (int k = 0; k < count_b; ++k) { - ids_b[k] = conn_b[count_b * p_b + k]; - } + const index_t* row = conn_b + count_b * p_b; + b0 = row[0]; + b1 = count_b > 1 ? row[1] : b0; + b2 = count_b > 2 ? row[2] : b0; } - for (int i = 0; i < count_a; ++i) { - for (int j = 0; j < count_b; ++j) { - if (ids_a[i] == ids_b[j]) { - return true; - } - } - } - return false; + + return a0 == b0 || a0 == b1 || a0 == b2 // + || a1 == b0 || a1 == b1 || a1 == b2 // + || a2 == b0 || a2 == b1 || a2 == b2; } /// @brief Append a (source_prim, target_prim) pair (post-swap) via an From a16ade8cc41a54204f99293dd05a17b42db9a293 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 15:26:53 -0400 Subject: [PATCH 08/10] Share the duplicated LBVH code between the CPU and CUDA broad phases ipc::LBVH and ipc::cuda::LBVH held line-for-line ports of the same algorithm, with agreement asserted only in comments and checked only by tests. Hoist the pieces that need no parallel abstraction into shared host/device code so the agreement is structural. New shared code: - ipc::morton_code() computes a box's Morton code from its center, normalizing by a domain whose width is passed as a reciprocal so the host and device multiply rather than divide. - ipc::count_leading_zeros() and ipc::morton_common_prefix() replace the per-platform CLZ dispatch and the duplicate-code fallback rule (Apetrei 2014's delta). - ipc::details::can_*_collide() hold the five mesh-connectivity filters. These were duplicated three times, not two: ipc::BroadPhase carries the same logic over AABB::vertex_ids. LBVH::Node's is_inner/is_leaf/is_valid/intersects are now IPC_TOOLKIT_HOST_DEVICE, so the traversal kernel calls the same predicates as the CPU instead of open-coding is_inner_marker == 0 and reimplementing the AABB overlap test. 167 duplicated lines collapse into 111 shared ones. Node::intersects() also generates better SASS than the hand-expanded aabb_intersects it replaces: traverse_kernel drops 320 -> 304 instructions, 33 -> 29 global loads, 16 -> 12 float compares and 5 -> 3 reconvergence pairs, with the register count (35-36) and the 0x100 local frame unchanged. Holding that frame is a hard requirement here -- the sm_120 miscompile fixed in c03e546d was frame-size sensitive. Morton codes are unchanged bit-for-bit. check_tree only compares root AABBs within 1e-4, so the suite cannot establish this; a standalone harness comparing the shared function against both prior forms over 200,000 random 2D and 3D cases found zero differences, and the compute_morton_codes_kernel opcode histogram is unchanged. Tested: full suite (4,348,911 assertions in 353 cases), compute-sanitizer memcheck on [lbvh][gpu] with 0 errors, clang-format and clang-tidy clean. Co-Authored-By: Claude Opus 5 (1M context) --- src/ipc/broad_phase/CMakeLists.txt | 1 + src/ipc/broad_phase/broad_phase.cpp | 44 ++---- src/ipc/broad_phase/cuda/lbvh.cu | 95 ++++-------- src/ipc/broad_phase/details/CMakeLists.txt | 6 + .../details/connectivity_filters.hpp | 138 ++++++++++++++++++ src/ipc/broad_phase/lbvh.cpp | 86 +++-------- src/ipc/broad_phase/lbvh.hpp | 15 +- src/ipc/math/morton.hpp | 94 ++++++++++++ 8 files changed, 308 insertions(+), 171 deletions(-) create mode 100644 src/ipc/broad_phase/details/CMakeLists.txt create mode 100644 src/ipc/broad_phase/details/connectivity_filters.hpp diff --git a/src/ipc/broad_phase/CMakeLists.txt b/src/ipc/broad_phase/CMakeLists.txt index 2ae7503e5..1c2346764 100644 --- a/src/ipc/broad_phase/CMakeLists.txt +++ b/src/ipc/broad_phase/CMakeLists.txt @@ -24,6 +24,7 @@ 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/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/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index 3c6ad8fbb..7d3537e8f 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -3,6 +3,7 @@ #ifdef IPC_TOOLKIT_WITH_CUDA #include +#include #include #include #include @@ -223,11 +224,7 @@ namespace { if (j < 0 || j >= n) { return -1; } - const uint64_t code_j = sorted_codes[j]; - if (code_i == code_j) { - return 32 + __clz(i ^ j); - } - return __clzll(static_cast(code_i ^ code_j)); + return ipc::morton_common_prefix(code_i, i, sorted_codes[j], j); } /// @brief Compute one Morton code per box from its (normalized) center. @@ -247,19 +244,16 @@ namespace { return; } - const double cx = 0.5 * (box_min[3 * i + 0] + box_max[3 * i + 0]); - const double cy = 0.5 * (box_min[3 * i + 1] + box_max[3 * i + 1]); - const double cz = 0.5 * (box_min[3 * i + 2] + box_max[3 * i + 2]); - - // (center - mesh_min) * mesh_width_inv -- the reciprocal is - // precomputed once per build (see compute_domain) and multiplied here - // instead of dividing per box, matching the CPU (ipc::LBVH::init_bvh) + // mesh_width_inv is the reciprocal of the domain width, computed + // once per build (see compute_domain), so ipc::morton_code() multiplies + // rather than divides and matches the CPU (ipc::LBVH::init_bvh) // bit-for-bit. - const double mx = (cx - mesh_min.x()) * mesh_width_inv.x(); - const double my = (cy - mesh_min.y()) * mesh_width_inv.y(); - const double mz = (cz - mesh_min.z()) * mesh_width_inv.z(); + 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] = (dim == 2) ? morton_2D(mx, my) : morton_3D(mx, my, mz); + codes[i] = ipc::morton_code(center, mesh_min, mesh_width_inv, dim); box_ids[i] = i; } @@ -392,7 +386,6 @@ namespace { /// @brief After the root swap, rewrite left pointers that referenced the /// old node 0 to its new location. See the CPU swap_root_to_zero comment: /// the old node 0 was only ever a left child, so only .left needs patching. - /// is_inner_marker aliases .right and is nonzero iff internal. /// @param nodes The BVH nodes. /// @param num_nodes The number of nodes. /// @param root The new location of the old node 0. @@ -405,7 +398,7 @@ namespace { if (i >= num_nodes) { return; } - if (nodes[i].is_inner_marker != 0 && nodes[i].left == 0) { + if (nodes[i].is_inner() && nodes[i].left == 0) { nodes[i].left = root; } } @@ -633,18 +626,10 @@ namespace { // -- Traversal ---------------------------------------------------------- - __device__ inline bool - aabb_intersects(const ipc::LBVH::Node& a, const ipc::LBVH::Node& b) - { - return a.aabb_min[0] <= b.aabb_max[0] && b.aabb_min[0] <= a.aabb_max[0] - && a.aabb_min[1] <= b.aabb_max[1] && b.aabb_min[1] <= a.aabb_max[1] - && a.aabb_min[2] <= b.aabb_max[2] && b.aabb_min[2] <= a.aabb_max[2]; - } - /// @brief Whether two primitives share a vertex id (the device connectivity /// filter). A vertex primitive's id set is {itself}; an edge's is its 2 /// endpoints; a face's is its 3 vertices. This is exactly the - /// shared-endpoint exclusion in ipc::LBVH's can_*_collide (for + /// shared-endpoint exclusion in ipc::details::can_*_collide (for /// vertex-vertex it reduces to p_a == p_b). /// @param p_a The first primitive id. /// @param conn_a The first primitive's connectivity, or null for a vertex. @@ -762,7 +747,7 @@ namespace { if constexpr (triangular) { break; // no self-collision with a single primitive } - if (aabb_intersects(node, query) + if (node.intersects(query) && !prim_shares_vertex( query.primitive_id, source_conn, source_count, node.primitive_id, target_conn, target_count)) { @@ -775,8 +760,8 @@ namespace { const ipc::LBVH::Node& child_l = target[node.left]; const ipc::LBVH::Node& child_r = target[node.right]; - bool intersects_l = aabb_intersects(child_l, query); - bool intersects_r = aabb_intersects(child_r, query); + bool intersects_l = child_l.intersects(query); + bool intersects_r = child_r.intersects(query); // Skip subtrees fully on the query's left (triangular only). if constexpr (triangular) { @@ -790,9 +775,8 @@ namespace { } } - // is_inner_marker aliases .right; it is 0 iff the node is a leaf. - const bool l_leaf = (child_l.is_inner_marker == 0); - const bool r_leaf = (child_r.is_inner_marker == 0); + const bool l_leaf = child_l.is_leaf(); + const bool r_leaf = child_r.is_leaf(); if (intersects_l && l_leaf && !prim_shares_vertex( @@ -1302,8 +1286,9 @@ LBVH::DeviceCandidateView LBVH::detect_face_face_candidates_device() const bool LBVH::can_edge_vertex_collide(size_t ei, size_t vi) const { const auto& [e0i, e1i] = m_impl->h_edge_vertex_ids[ei]; - return vi != e0i && vi != e1i - && (can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i)); + + return ipc::details::can_edge_vertex_collide( + e0i, e1i, vi, can_vertices_collide); } bool LBVH::can_edges_collide(size_t eai, size_t ebi) const @@ -1311,21 +1296,16 @@ bool LBVH::can_edges_collide(size_t eai, size_t ebi) const const auto& [ea0i, ea1i] = m_impl->h_edge_vertex_ids[eai]; const auto& [eb0i, eb1i] = m_impl->h_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 ipc::details::can_edges_collide( + ea0i, ea1i, eb0i, eb1i, can_vertices_collide); } bool LBVH::can_face_vertex_collide(size_t fi, size_t vi) const { const auto& [f0i, f1i, f2i] = m_impl->h_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 ipc::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 @@ -1333,14 +1313,8 @@ bool LBVH::can_edge_face_collide(size_t ei, size_t fi) const const auto& [e0i, e1i] = m_impl->h_edge_vertex_ids[ei]; const auto& [f0i, f1i, f2i] = m_impl->h_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 ipc::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 @@ -1348,19 +1322,8 @@ bool LBVH::can_faces_collide(size_t fai, size_t fbi) const const auto& [fa0i, fa1i, fa2i] = m_impl->h_face_vertex_ids[fai]; const auto& [fb0i, fb1i, fb2i] = m_impl->h_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 ipc::details::can_faces_collide( + fa0i, fa1i, fa2i, fb0i, fb1i, fb2i, can_vertices_collide); } size_t LBVH::num_vertex_nodes() const diff --git a/src/ipc/broad_phase/details/CMakeLists.txt b/src/ipc/broad_phase/details/CMakeLists.txt new file mode 100644 index 000000000..1fbd0bfac --- /dev/null +++ b/src/ipc/broad_phase/details/CMakeLists.txt @@ -0,0 +1,6 @@ +set(SOURCES + connectivity_filters.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..d0c3eeb8e --- /dev/null +++ b/src/ipc/broad_phase/details/connectivity_filters.hpp @@ -0,0 +1,138 @@ +#pragma once + +#include +#include + +#include + +namespace ipc::details { + +// Mesh-connectivity collision filters shared by every broad phase. +// +// ipc::BroadPhase, ipc::LBVH, 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. + +/// @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) +{ + return vi != e0i && vi != e1i + && (can_vertices_collide(vi, e0i) || can_vertices_collide(vi, e1i)); +} + +/// @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 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)); +} + +/// @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) +{ + return vi != f0i && vi != f1i && vi != f2i + && (can_vertices_collide(vi, f0i) || can_vertices_collide(vi, f1i) + || can_vertices_collide(vi, f2i)); +} + +/// @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 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)); +} + +/// @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 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)); +} + +} // namespace ipc::details diff --git a/src/ipc/broad_phase/lbvh.cpp b/src/ipc/broad_phase/lbvh.cpp index 97ba8a9c6..35ff8c5e7 100644 --- a/src/ipc/broad_phase/lbvh.cpp +++ b/src/ipc/broad_phase/lbvh.cpp @@ -1,5 +1,6 @@ #include "lbvh.hpp" +#include #include #include #include @@ -96,11 +97,11 @@ void LBVH::build( } 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). + /// 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. code_i is the + /// Morton code at position i, passed explicitly to avoid a redundant + /// lookup. The prefix itself is computed by ipc::morton_common_prefix(), + /// shared with the device build in ipc::cuda::LBVH. int delta( const LBVH::MortonCodeElements& sorted_morton_codes, int i, @@ -110,24 +111,8 @@ namespace { 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 + return morton_common_prefix( + code_i, i, sorted_morton_codes[j].morton_code, j); } } // namespace @@ -154,17 +139,8 @@ void LBVH::init_bvh( 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; }); } @@ -808,8 +784,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 +794,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 +803,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 +814,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 +825,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..f166d0156 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(); diff --git a/src/ipc/math/morton.hpp b/src/ipc/math/morton.hpp index d43c06bac..acd3b76db 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(WIN32) +#include // for __lzcnt / __lzcnt64 +#endif + namespace ipc { /// @brief Expands a 32-bit integer into 64 bits by inserting 1 zero after each bit. @@ -64,4 +71,91 @@ IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_3D(double x, double y, double z) return (xx << 2) | (yy << 1) | zz; } +/// @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 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(WIN32) + return static_cast(__lzcnt(v)); +#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(WIN32) + return static_cast(__lzcnt64(v)); +#else +#error "count_leading_zeros: no leading-zero-count intrinsic for this compiler" +#endif +} + +/// @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 32 so that any code-level difference always compares as the shorter +/// prefix. +/// +/// @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. +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 32 + count_leading_zeros(static_cast(i ^ j)); + } + return count_leading_zeros(code_i ^ code_j); +} + } // namespace ipc \ No newline at end of file From 6ce6abc0217ffc81555114f6833b33ab30c04c6b Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 15:55:43 -0400 Subject: [PATCH 09/10] Share the LBVH build and traversal between the CPU and CUDA broad phases The Apetrei 2014 bottom-up build and the BVH descent were line-for-line ports between ipc::LBVH and ipc::cuda::LBVH. Hoist both into shared host/device code, leaving each platform only what it genuinely owns: the parallel launch, the sort, and a small policy per difference. ipc::details::build_hierarchy_from_leaf() takes the sorted-code accessor (the host stores an array of structs, the device a flat array) and the atomic arrival gate. ipc::details::traverse_lbvh() takes what to do on an overlap, which is the whole of the host/device difference there: the host filters and appends to a std::vector, the device filters against the mesh connectivity and appends through an atomic counter. Also shared: set_inflated_aabb(), init_leaf_node(), delta(), is_left_child(), swap_root_to_zero() and patch_left_pointer(). LBVH::ConstructionInfo is now a template over its counter type, so the host uses std::atomic and the device a plain int, from one layout. 434 lines leave the two implementations for 197 lines of shared code. Fixes a latent race in the device build. The arrival gate had a __threadfence() on the release side but none on the acquire side, then read the sibling's child pointer, range endpoint and rightmost leaf with ordinary loads, which may be served from a stale L1 on another SM. The shared gate's contract requires both halves, and the device policy now fences after an increment that returns nonzero. The kernel's SASS gains exactly one MEMBAR.SC.GPU, giving MEMBAR.SC.GPU / ATOMG.E.ADD.STRONG.GPU / MEMBAR.SC.GPU with the paired CCTL.IVALL that invalidates L1. This would have corrupted internal-node AABBs and rightmost[] without breaking the tree structure, so check_tree's structural checks could not have caught it. The single-leaf build case is now explicit on both sides. The host previously relied on writing nodes[0].left = 0 over the lone leaf's primitive_id, which was only correct because a one-box sort always yields box_id 0. traverse_kernel's SASS is bit-identical after the change -- the templated descent and its lambda inline away completely -- and every kernel's register count is unchanged from before this series. The 0x100 frame that the sm_120 miscompile in c03e546d turned on is preserved; that was the acceptance gate for touching this kernel at all. Adds coverage for single-primitive BVHs, which no existing mesh reaches. Only face-vertex and edge-face put a one-node BVH in the traversal target position, so the test builds one face and one disjoint edge and checks both against BruteForce, and the device against the host. Verified non-vacuous by mutation: suppressing the emit in the single-node branch fails 4 of its assertions. Tested: full suite (4,348,947 assertions in 354 cases), compute-sanitizer memcheck on [lbvh][gpu] with 0 errors, clang-format and clang-tidy clean. Co-Authored-By: Claude Opus 5 (1M context) --- src/ipc/broad_phase/cuda/lbvh.cu | 288 +++++------------- src/ipc/broad_phase/details/CMakeLists.txt | 2 + src/ipc/broad_phase/details/lbvh_build.hpp | 259 ++++++++++++++++ src/ipc/broad_phase/details/lbvh_traverse.hpp | 124 ++++++++ src/ipc/broad_phase/lbvh.cpp | 264 +++------------- src/ipc/broad_phase/lbvh.hpp | 22 +- tests/src/tests/broad_phase/test_lbvh.cpp | 116 +++++++ 7 files changed, 634 insertions(+), 441 deletions(-) create mode 100644 src/ipc/broad_phase/details/lbvh_build.hpp create mode 100644 src/ipc/broad_phase/details/lbvh_traverse.hpp diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index 7d3537e8f..ea8be0b75 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -4,6 +4,8 @@ #include #include +#include +#include #include #include #include @@ -29,17 +31,12 @@ namespace { sizeof(Eigen::Array3d) == 24, "Eigen::Array3d must be 24 bytes (3 packed doubles)"); - /// @brief Per-internal-node scratch used by the bottom-up build. The device - /// analog of ipc::LBVH::ConstructionInfo, kept separate on purpose: that - /// struct's visitation_count is a std::atomic, which cannot be used - /// here (atomicAdd needs an int*, and std::atomic is non-copyable so it - /// cannot be a thrust::device_vector element). A plain int suffices because - /// atomicAdd provides the atomicity the CPU gets from std::atomic. - struct DeviceConstructionInfo { - int left_range; - int right_range; - int visitation_count; - }; + /// @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*, std::atomic is non-copyable and so cannot be + /// a thrust::device_vector element, and atomicAdd() supplies the same + /// atomicity. The layout is otherwise identical to the host's. + using DeviceConstructionInfo = ipc::LBVH::ConstructionInfo; /// @brief Min/max domain accumulator for the Morton-normalization reduction. struct Domain { @@ -204,29 +201,6 @@ namespace { // -- Tree building ------------------------------------------------------ - /// @brief Number of common leading bits between Morton codes at sorted - /// positions i and j (device port of the CPU delta()). Duplicate codes fall - /// back to the CLZ of the index XOR (offset by 32 so it sorts after any - /// code-level difference). - /// @param sorted_codes The Morton codes in sorted order. - /// @param n The number of codes. - /// @param i The first sorted position. - /// @param code_i The code at position i (passed to avoid a redundant look-up). - /// @param j The second sorted position. - /// @return The common-prefix length, or -1 when j is out of bounds. - __device__ inline int delta_device( - const uint64_t* __restrict__ sorted_codes, - const int n, - const int i, - const uint64_t code_i, - const int j) - { - if (j < 0 || j >= n) { - return -1; - } - return ipc::morton_common_prefix(code_i, i, sorted_codes[j], j); - } - /// @brief Compute one Morton code per box from its (normalized) center. /// Mirrors the compute_morton_codes block of ipc::LBVH::init_bvh. __global__ void compute_morton_codes_kernel( @@ -258,9 +232,18 @@ namespace { } /// @brief Single-pass bottom-up hierarchy + AABB build (Apetrei 2014). - /// One thread per leaf. Direct port of the build_hierarchy_and_boxes block - /// of ipc::LBVH::init_bvh, with atomicAdd + __threadfence replacing the - /// std::atomic arrival gate. + /// 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, @@ -277,94 +260,43 @@ namespace { return; } - const int LEAF_OFFSET = N_LEAVES - 1; - - // --- Initialize leaf node --- - { - const index_t bid = sorted_box_ids[i]; - ipc::LBVH::Node leaf; -#pragma unroll - for (int k = 0; k < 3; ++k) { - // Round the float AABB out (matches assign_inflated_aabb). - leaf.aabb_min[k] = nextafterf( - static_cast(box_min[3 * bid + k]), -INFINITY); - leaf.aabb_max[k] = nextafterf( - static_cast(box_max[3 * bid + k]), INFINITY); - } - leaf.primitive_id = static_cast(bid); - leaf.is_inner_marker = 0; - nodes[LEAF_OFFSET + i] = leaf; - // A leaf's rightmost leaf is itself. - rightmost[LEAF_OFFSET + i] = i; - } - - // Single-node tree: the leaf is the root; no internal nodes to build. - if (N_LEAVES == 1) { - if (i == 0) { - *root_idx = 0; - } - return; - } - - // --- Bottom-up walk (Apetrei 2014, Fig. 2) --- - int left_key = i; - int right_key = i; - int current_node = LEAF_OFFSET + i; - - while (true) { - // Choose parent (see the CPU comment in ipc::LBVH::init_bvh). - const bool is_child_a = (left_key == 0) - || (right_key != N_LEAVES - 1 - && delta_device( - sorted_codes, N_LEAVES, right_key, - sorted_codes[right_key], right_key + 1) - > delta_device( - sorted_codes, N_LEAVES, left_key - 1, - sorted_codes[left_key - 1], left_key)); - const int parent = is_child_a ? right_key : left_key - 1; - - // Write the child pointer + range onto the parent. - 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; - } - - // Publish this child's node data and range to all threads before - // signaling arrival, so the second thread reads consistent state. - __threadfence(); - - // Atomic arrival gate: first thread stops; second proceeds knowing - // both children are complete. - if (atomicAdd(&infos[parent].visitation_count, 1) == 0) { - break; // first thread to arrive -> finished - } + const index_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; + }); - // Second thread: compute the parent AABB union and rightmost leaf. - const ipc::LBVH::Node& child_a = nodes[nodes[parent].left]; - const ipc::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); - rightmost[parent] = ::max( - rightmost[nodes[parent].left], rightmost[nodes[parent].right]); - - // 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) { - // Only one thread reaches the root. - *root_idx = current_node; - break; - } + if (root >= 0) { + *root_idx = root; // only one thread reaches the root } } - /// @brief Swap the node and rightmost-leaf entries at indices 0 and root - /// (runs on a single thread). + /// @brief Swap the node and rightmost-leaf entries at index 0 and the root + /// (single thread). See ipc::details::swap_root_to_zero(). /// @param nodes The BVH nodes. /// @param rightmost The per-node rightmost-leaf indices. /// @param root The index to swap with index 0. @@ -374,18 +306,12 @@ namespace { const int root) { if (blockIdx.x == 0 && threadIdx.x == 0) { - const ipc::LBVH::Node tmp = nodes[0]; - nodes[0] = nodes[root]; - nodes[root] = tmp; - const int32_t t = rightmost[0]; - rightmost[0] = rightmost[root]; - rightmost[root] = t; + 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 the CPU swap_root_to_zero comment: - /// the old node 0 was only ever a left child, so only .left needs patching. + /// 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. @@ -398,9 +324,7 @@ namespace { if (i >= num_nodes) { return; } - if (nodes[i].is_inner() && nodes[i].left == 0) { - nodes[i].left = root; - } + ipc::details::patch_left_pointer(nodes[i], root); } /// @brief Build one BVH on the device from device-resident box corners. @@ -704,12 +628,27 @@ namespace { /// @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. Descent is a direct port of traverse_lbvh() - /// in lbvh.cpp (scalar path); the connectivity (shared-vertex) exclusion is - /// applied here on the device. The remaining user vertex filter (if any) is - /// applied on the host, so the final set matches the CPU ipc::LBVH. + /// pair to the output arrays. The descent is ipc::details::traverse_lbvh(), + /// shared with the CPU ipc::LBVH; the connectivity (shared-vertex) + /// exclusion is applied here on the device. The remaining user vertex + /// filter (if any) is applied on the host, so the final set matches the CPU + /// ipc::LBVH. /// @tparam triangular Self-collision: skip subtrees fully left of the query. /// @tparam swap_order Emit (target_prim, source_prim) instead. + /// @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 source_count The vertex ids per source primitive (1, 2, or 3). + /// @param target_conn The target connectivity (null for vertices). + /// @param target_count The vertex ids per target primitive (1, 2, or 3). + /// @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 __global__ void traverse_kernel( const ipc::LBVH::Node* __restrict__ source, @@ -732,81 +671,18 @@ namespace { return; } const ipc::LBVH::Node query = source[source_leaf_offset + s]; - const int query_leaf_idx = s; - - constexpr int MAX_STACK_SIZE = 64; - int stack[MAX_STACK_SIZE]; - int stack_ptr = 0; - stack[stack_ptr++] = ipc::LBVH::Node::INVALID_POINTER; // 0 - - int node_idx = 0; // root - do { - const ipc::LBVH::Node& node = target[node_idx]; - if (target_size == 1) { // single node (only root, which is a leaf) - if constexpr (triangular) { - break; // no self-collision with a single primitive - } - if (node.intersects(query) - && !prim_shares_vertex( + ipc::details::traverse_lbvh( + query, s, target, target_size, target_rightmost, + [&](const ipc::LBVH::Node& leaf) { + if (!prim_shares_vertex( query.primitive_id, source_conn, source_count, - node.primitive_id, target_conn, target_count)) { + leaf.primitive_id, target_conn, target_count)) { emit_pair( - query.primitive_id, node.primitive_id, out_a, out_b, + query.primitive_id, leaf.primitive_id, out_a, out_b, counter, capacity); } - break; - } - - const ipc::LBVH::Node& child_l = target[node.left]; - const ipc::LBVH::Node& child_r = target[node.right]; - bool intersects_l = child_l.intersects(query); - bool intersects_r = child_r.intersects(query); - - // Skip subtrees fully on the query's left (triangular only). - if constexpr (triangular) { - 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; - } - } - - const bool l_leaf = child_l.is_leaf(); - const bool r_leaf = child_r.is_leaf(); - - if (intersects_l && l_leaf - && !prim_shares_vertex( - query.primitive_id, source_conn, source_count, - child_l.primitive_id, target_conn, target_count)) { - emit_pair( - query.primitive_id, child_l.primitive_id, out_a, out_b, - counter, capacity); - } - if (intersects_r && r_leaf - && !prim_shares_vertex( - query.primitive_id, source_conn, source_count, - child_r.primitive_id, target_conn, target_count)) { - emit_pair( - query.primitive_id, child_r.primitive_id, out_a, out_b, - counter, capacity); - } - - const bool traverse_l = intersects_l && !l_leaf; - const bool traverse_r = intersects_r && !r_leaf; - - if (!traverse_l && !traverse_r) { - node_idx = stack[--stack_ptr]; - } else { - node_idx = traverse_l ? node.left : node.right; - if (traverse_l && traverse_r) { - stack[stack_ptr++] = node.right; - } - } - } while (node_idx != ipc::LBVH::Node::INVALID_POINTER); + }); } /// @brief Run the device traversal of the target BVH by the source leaves, diff --git a/src/ipc/broad_phase/details/CMakeLists.txt b/src/ipc/broad_phase/details/CMakeLists.txt index 1fbd0bfac..942162c1f 100644 --- a/src/ipc/broad_phase/details/CMakeLists.txt +++ b/src/ipc/broad_phase/details/CMakeLists.txt @@ -1,5 +1,7 @@ set(SOURCES connectivity_filters.hpp + lbvh_build.hpp + lbvh_traverse.hpp spatial_hash_impl.hpp ) 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..502a52d1c --- /dev/null +++ b/src/ipc/broad_phase/details/lbvh_build.hpp @@ -0,0 +1,259 @@ +#pragma once + +#include +#include +#include + +#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. +/// +/// @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..784570346 --- /dev/null +++ b/src/ipc/broad_phase/details/lbvh_traverse.hpp @@ -0,0 +1,124 @@ +#pragma once + +#include +#include + +#include +#include + +namespace ipc::details { + +/// @brief Descends a target BVH for one query leaf, reporting every target leaf +/// whose AABB overlaps the query. +/// +/// 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 and ipc::cuda::LBVH. What differs between them +/// is only what happens on an overlap, which is why that is a policy: 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 Emit Callable (const LBVH::Node& leaf) -> void, invoked for every +/// overlapping target leaf. It owns both the collision filtering and the +/// recording of the pair. +/// +/// @param query The querying leaf node. +/// @param query_leaf_idx The query's position in its own Morton-sorted leaf +/// order. 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 emit The per-overlap callback. +template +IPC_TOOLKIT_HOST_DEVICE void traverse_lbvh( + const LBVH::Node& query, + const int query_leaf_idx, + const LBVH::Node* target, + const int target_size, + const int32_t* target_rightmost, + Emit&& emit) +{ + // A fixed-size stack keeps the descent free of dynamic allocation. + 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 = 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 + } + if (node.intersects(query)) { + emit(node); + } + 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]; + bool intersects_l = child_l.intersects(query); + bool intersects_r = child_r.intersects(query); + + // Ignore a subtree lying entirely to the query's left; that pair is + // reported when the other primitive is the query instead. + if constexpr (triangular) { + 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; + } + } + + const bool l_leaf = child_l.is_leaf(); + const bool r_leaf = child_r.is_leaf(); + + // An overlapped leaf is a candidate. + if (intersects_l && l_leaf) { + emit(child_l); + } + if (intersects_r && r_leaf) { + emit(child_r); + } + + // An overlapped inner node is descended into. + const bool traverse_l = intersects_l && !l_leaf; + const bool traverse_r = intersects_r && !r_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) { + assert(stack_ptr < MAX_STACK_SIZE); + 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/lbvh.cpp b/src/ipc/broad_phase/lbvh.cpp index 35ff8c5e7..283543aaa 100644 --- a/src/ipc/broad_phase/lbvh.cpp +++ b/src/ipc/broad_phase/lbvh.cpp @@ -1,6 +1,8 @@ #include "lbvh.hpp" #include +#include +#include #include #include #include @@ -25,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( @@ -96,26 +81,6 @@ void LBVH::build( face_boxes.clear(); } -namespace { - /// 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. code_i is the - /// Morton code at position i, passed explicitly to avoid a redundant - /// lookup. The prefix itself is computed by ipc::morton_common_prefix(), - /// shared with the device build in ipc::cuda::LBVH. - 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; - } - return morton_common_prefix( - code_i, i, sorted_morton_codes[j].morton_code, j); - } -} // namespace - void LBVH::init_bvh( const AABBs& boxes, Nodes& lbvh, RightmostLeaves& rightmost_leaves) const { @@ -156,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()); @@ -172,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); }); } } @@ -342,6 +216,9 @@ 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 ipc::cuda::LBVH. template void traverse_lbvh( const LBVH::Node& query, @@ -351,81 +228,12 @@ 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( + query, int(query_leaf_idx), lbvh.data(), int(lbvh.size()), + rightmost_leaves.data(), [&](const LBVH::Node& leaf) { 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 diff --git a/src/ipc/broad_phase/lbvh.hpp b/src/ipc/broad_phase/lbvh.hpp index f166d0156..4274e61d8 100644 --- a/src/ipc/broad_phase/lbvh.hpp +++ b/src/ipc/broad_phase/lbvh.hpp @@ -91,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/tests/src/tests/broad_phase/test_lbvh.cpp b/tests/src/tests/broad_phase/test_lbvh.cpp index aa802bf12..d68b254ae 100644 --- a/tests/src/tests/broad_phase/test_lbvh.cpp +++ b/tests/src/tests/broad_phase/test_lbvh.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -289,6 +290,121 @@ 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)); + } + +#ifdef IPC_TOOLKIT_WITH_CUDA + // The device build has its own single-leaf branch, so check it agrees with + // the host on exactly these trees. + cuda::LBVH gpu_lbvh; + gpu_lbvh.build(vertices, edges, faces, inflation_radius); + + REQUIRE(gpu_lbvh.num_face_nodes() == 1); + REQUIRE(gpu_lbvh.num_edge_nodes() == 1); + + { + std::vector gpu_candidates, cpu_candidates; + gpu_lbvh.detect_face_vertex_candidates(gpu_candidates); + lbvh.detect_face_vertex_candidates(cpu_candidates); + REQUIRE(!cpu_candidates.empty()); + CHECK(gpu_candidates.size() == cpu_candidates.size()); + CHECK(contains_all_candidates(gpu_candidates, cpu_candidates)); + } + + { + std::vector gpu_candidates, cpu_candidates; + gpu_lbvh.detect_edge_face_candidates(gpu_candidates); + lbvh.detect_edge_face_candidates(cpu_candidates); + REQUIRE(!cpu_candidates.empty()); + CHECK(gpu_candidates.size() == cpu_candidates.size()); + CHECK(contains_all_candidates(gpu_candidates, cpu_candidates)); + } + + { + std::vector gpu_candidates, cpu_candidates; + gpu_lbvh.detect_edge_vertex_candidates(gpu_candidates); + lbvh.detect_edge_vertex_candidates(cpu_candidates); + REQUIRE(!cpu_candidates.empty()); + CHECK(gpu_candidates.size() == cpu_candidates.size()); + CHECK(contains_all_candidates(gpu_candidates, cpu_candidates)); + } +#endif +} + TEST_CASE( "Benchmark LBVH::detect_edge_edge_candidates", "[!benchmark][broad_phase][lbvh]") From abfa013052c3d0070995b22c8cf803d3b9fa57d0 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Wed, 9 Sep 2026 15:32:32 -0400 Subject: [PATCH 10/10] Remove warning from spdlog.cmake - Warning falsly triggers on our own dependencies because spdlog.cmake takes precedence over downstream spdlog.cmake scripts. --- cmake/recipes/spdlog.cmake | 12 ------------ 1 file changed, 12 deletions(-) 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()