diff --git a/cmake/recipes/spdlog.cmake b/cmake/recipes/spdlog.cmake index 5e083f888..d3be8c70c 100644 --- a/cmake/recipes/spdlog.cmake +++ b/cmake/recipes/spdlog.cmake @@ -1,18 +1,6 @@ # spdlog (https://github.com/gabime/spdlog) # License: MIT if(TARGET spdlog::spdlog) - # Someone else created the target, so the fmt patch below never runs. That - # is fine without CUDA, but the patch is what makes the bundled fmt - # compile under nvcc at all, so warn rather than fail at the first .cu. - if(IPC_TOOLKIT_WITH_CUDA) - message(WARNING - "spdlog::spdlog was provided by an enclosing project, so IPC " - "Toolkit's cmake/patches/fmt-nvcc-compat.patch was not applied. " - "The bundled fmt does not compile under nvcc unpatched: its " - "literal-encoding probe misfires and a char32_t table uses hex " - "escapes with the high bit set. Apply the same patch to your " - "spdlog, or let IPC Toolkit fetch its own.") - endif() return() endif() diff --git a/src/ipc/broad_phase/CMakeLists.txt b/src/ipc/broad_phase/CMakeLists.txt index e0613ebcf..1c2346764 100644 --- a/src/ipc/broad_phase/CMakeLists.txt +++ b/src/ipc/broad_phase/CMakeLists.txt @@ -23,3 +23,8 @@ set(SOURCES ) target_sources(ipc_toolkit PRIVATE ${SOURCES}) + +add_subdirectory(details) +if(IPC_TOOLKIT_WITH_CUDA) + add_subdirectory(cuda) +endif() diff --git a/src/ipc/broad_phase/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/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/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..ea8be0b75 --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -0,0 +1,1237 @@ +#include "lbvh.hpp" + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include +#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 + /// 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 { + 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.) + // + // 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. + // + // 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 + const int n, + const int dim, + const double inflation_radius, + double* __restrict__ box_min, // always 3 * n, row-major + double* __restrict__ box_max) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } +#pragma unroll + for (int k = 0; k < 3; ++k) { + if (k < dim) { + const double v = vertices[dim * i + k]; + 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; + } + } + } + + __global__ void build_vertex_boxes_dynamic_kernel( + 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, // always 3 * n, row-major + double* __restrict__ box_max) + { + const int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i >= n) { + return; + } +#pragma unroll + for (int k = 0; k < 3; ++k) { + if (k < dim) { + const double a = vertices_t0[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, NEG_INF); + box_max[3 * i + k] = + nextafter(fmax(a, b) + inflation_radius, POS_INF); + } else { + box_min[3 * i + k] = 0.0; + box_max[3 * i + k] = 0.0; + } + } + } + + __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 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; + } + + // 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 Eigen::Array3d center( + 0.5 * (box_min[3 * i + 0] + box_max[3 * i + 0]), + 0.5 * (box_min[3 * i + 1] + box_max[3 * i + 1]), + 0.5 * (box_min[3 * i + 2] + box_max[3 * i + 2])); + + codes[i] = ipc::morton_code(center, mesh_min, mesh_width_inv, dim); + box_ids[i] = i; + } + + /// @brief Single-pass bottom-up hierarchy + AABB build (Apetrei 2014). + /// One thread per leaf, driving ipc::details::build_hierarchy_from_leaf() + /// -- the same walk the CPU build runs -- with atomicAdd() and the fences + /// around it standing in for the host's std::atomic arrival. + /// @param box_min The box min corners (3 * n, row-major). + /// @param box_max The box max corners (3 * n, row-major). + /// @param sorted_codes The Morton codes in sorted order. + /// @param sorted_box_ids The box ids in Morton-sorted order. + /// @param N_LEAVES The number of leaves. + /// @param[out] nodes The BVH nodes. + /// @param[out] rightmost The per-node rightmost-leaf indices. + /// @param[in,out] infos The per-node construction scratch (zeroed). + /// @param[out] root_idx The root's index. + __global__ void build_hierarchy_kernel( + const double* __restrict__ box_min, + const double* __restrict__ box_max, + const uint64_t* __restrict__ sorted_codes, + const 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 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; + }); + + if (root >= 0) { + *root_idx = root; // only one thread reaches the root + } + } + + /// @brief Swap the node and rightmost-leaf entries at index 0 and the root + /// (single thread). See ipc::details::swap_root_to_zero(). + /// @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) { + ipc::details::swap_root_to_zero(nodes, rightmost, root); + } + } + + /// @brief After the root swap, rewrite left pointers that referenced the + /// old node 0 to its new location. See ipc::details::patch_left_pointer(). + /// @param nodes The BVH nodes. + /// @param num_nodes The number of nodes. + /// @param root The new location of the old node 0. + __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; + } + ipc::details::patch_left_pointer(nodes[i], 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 ---------------------------------------------------------- + + /// @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::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. + /// @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) + { + // 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) { + a0 = a1 = a2 = p_a; + } else { + 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) { + b0 = b1 = b2 = p_b; + } else { + 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; + } + + 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 + /// 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. 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, + 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]; + + 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, + leaf.primitive_id, target_conn, target_count)) { + emit_pair( + query.primitive_id, leaf.primitive_id, out_a, out_b, + counter, capacity); + } + }); + } + + /// @brief Run the device traversal of the target BVH by the source leaves, + /// 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. + /// @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 buf The output candidate buffer and capacity hint (in/out). + /// @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, + 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) { + buf.a.clear(); + buf.b.clear(); + return 0; + } + const int source_leaf_offset = n_source_leaves - 1; + + 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) { + buf.a.resize(capacity); + buf.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(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 (static_cast(count) <= capacity) { + break; // everything fit + } + + logger().warn( + "[ipc::cuda::LBVH] candidate count {} exceeded preallocated " + "capacity {}; re-running with the exact size (this cost is " + "amortized: later calls reuse the learned capacity)", + count, capacity); + capacity = static_cast(count); // exact size now known + } + + 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); + } + + /// @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(); + + 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 (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 < 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 (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, dim, + 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(); + + 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(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 < 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); + 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, dim, + 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(); + + assert(_dim == 2 || _dim == 3); + dim = _dim; + + 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); + 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); + 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 ipc::details::can_edge_vertex_collide( + e0i, e1i, vi, can_vertices_collide); +} + +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]; + + 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 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 +{ + const auto& [e0i, e1i] = m_impl->h_edge_vertex_ids[ei]; + const auto& [f0i, f1i, f2i] = m_impl->h_face_vertex_ids[fi]; + + 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 +{ + const auto& [fa0i, fa1i, fa2i] = m_impl->h_face_vertex_ids[fai]; + const auto& [fb0i, fb1i, fb2i] = m_impl->h_face_vertex_ids[fbi]; + + return ipc::details::can_faces_collide( + fa0i, fa1i, fa2i, fb0i, fb1i, fb2i, can_vertices_collide); +} + +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..3a19d5c31 --- /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 (2 or 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..30dd20152 --- /dev/null +++ b/src/ipc/broad_phase/cuda/lbvh_impl.cuh @@ -0,0 +1,105 @@ +// 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; + + /// @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. + } + }; + + 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/broad_phase/details/CMakeLists.txt b/src/ipc/broad_phase/details/CMakeLists.txt new file mode 100644 index 000000000..942162c1f --- /dev/null +++ b/src/ipc/broad_phase/details/CMakeLists.txt @@ -0,0 +1,8 @@ +set(SOURCES + connectivity_filters.hpp + lbvh_build.hpp + lbvh_traverse.hpp + spatial_hash_impl.hpp +) + +target_sources(ipc_toolkit PRIVATE ${SOURCES}) diff --git a/src/ipc/broad_phase/details/connectivity_filters.hpp b/src/ipc/broad_phase/details/connectivity_filters.hpp new file mode 100644 index 000000000..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/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 97ba8a9c6..283543aaa 100644 --- a/src/ipc/broad_phase/lbvh.cpp +++ b/src/ipc/broad_phase/lbvh.cpp @@ -1,5 +1,8 @@ #include "lbvh.hpp" +#include +#include +#include #include #include #include @@ -24,23 +27,6 @@ using namespace std::placeholders; namespace ipc { -namespace { - // Helper to safely convert double AABB to float AABB - inline void assign_inflated_aabb(const AABB& box, LBVH::Node& node) - { - // Round Min down - node.aabb_min = box.min.unaryExpr([](double val) { - return std::nextafter( - float(val), -std::numeric_limits::infinity()); - }); - // Round Max up - node.aabb_max = box.max.unaryExpr([](double val) { - return std::nextafter( - float(val), std::numeric_limits::infinity()); - }); - } -} // namespace - LBVH::LBVH() : BroadPhase() { static_assert( @@ -95,42 +81,6 @@ void LBVH::build( face_boxes.clear(); } -namespace { - /// Returns the number of common leading bits (CLZ of XOR) between sorted - /// Morton codes at positions i and j. code_i is the Morton code at position - /// i, passed explicitly to avoid a redundant lookup. Returns -1 when j is - /// out of bounds. Duplicate codes fall back to CLZ of the index XOR - /// (offset by 32 so it sorts after any code-level difference). - int delta( - const LBVH::MortonCodeElements& sorted_morton_codes, - int i, - uint64_t code_i, - int j) - { - if (j < 0 || j >= sorted_morton_codes.size()) { - return -1; - } - uint64_t code_j = sorted_morton_codes[j].morton_code; - if (code_i == code_j) { - // handle duplicate morton codes - int element_idx_i = i; - int element_idx_j = j; - - // add 32 for common prefix of code_i ^ code_j -#if defined(__GNUC__) || defined(__clang__) - return 32 + __builtin_clz(element_idx_i ^ element_idx_j); -#elif defined(WIN32) - return 32 + __lzcnt(element_idx_i ^ element_idx_j); -#endif - } -#if defined(__GNUC__) || defined(__clang__) - return __builtin_clzll(code_i ^ code_j); -#elif defined(WIN32) - return __lzcnt64(code_i ^ code_j); -#endif - } -} // namespace - void LBVH::init_bvh( const AABBs& boxes, Nodes& lbvh, RightmostLeaves& rightmost_leaves) const { @@ -154,17 +104,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; }); } @@ -180,7 +121,6 @@ void LBVH::init_bvh( assert(boxes.size() <= std::numeric_limits::max()); const int N_LEAVES = int(boxes.size()); - const int LEAF_OFFSET = N_LEAVES - 1; if (rightmost_leaves.size() != lbvh.size()) { rightmost_leaves.resize(lbvh.size()); @@ -196,136 +136,46 @@ void LBVH::init_bvh( } // Apetrei 2014: single bottom-up pass that simultaneously builds the - // hierarchy and computes bounding boxes. Each leaf thread walks toward the - // root, choosing its parent in O(1) by comparing the CLZ-delta values at - // the two ends of its current key range. - // - // In this layout internal node j always splits between sorted keys j and - // j+1. The root is NOT necessarily at index 0, so after construction we - // swap the root into position 0 to match the traversal code's expectation. + // hierarchy and computes bounding boxes. See + // ipc::details::build_hierarchy_from_leaf(), shared with the device build + // in ipc::cuda::LBVH. std::atomic root_idx(-1); { IPC_TOOLKIT_PROFILE_BLOCK("build_hierarchy_and_boxes"); tbb::parallel_for(0, N_LEAVES, [&](int i) { - // --- Initialize leaf node --- - { - const auto& box = boxes[morton_codes[i].box_id]; - - Node leaf_node; // Create leaf node - assign_inflated_aabb(box, leaf_node); - leaf_node.primitive_id = morton_codes[i].box_id; - leaf_node.is_inner_marker = 0; - lbvh[LEAF_OFFSET + i] = leaf_node; // Store leaf - // A leaf's rightmost leaf is itself - rightmost_leaves[LEAF_OFFSET + i] = i; - } - - // --- Bottom-up walk (Apetrei 2014, Fig. 2) --- - // Invariant: the current subtree covers the sorted-key range - // [left_key, right_key]. - int left_key = i; - int right_key = i; - int current_node = LEAF_OFFSET + i; - - while (true) { - // Choose parent. Candidates are internal node right_key - // (current becomes its left / childA) or internal node - // left_key-1 (current becomes its right / childB). Our delta() - // returns CLZ (higher = more-similar = finer split), so the - // CLOSER ancestor has the LARGER delta — hence ">". - // - // Boundary rules: - // left_key == 0 → must be childA (no node -1) - // right_key == n-1 → must be childB (no node n-1) - const bool is_child_a = (left_key == 0) - || (right_key != N_LEAVES - 1 - && delta( - morton_codes, right_key, - morton_codes[right_key].morton_code, - right_key + 1) - > delta( - morton_codes, left_key - 1, - morton_codes[left_key - 1].morton_code, - left_key)); - const int parent = is_child_a ? right_key : left_key - 1; - - auto& info = construction_infos[parent]; - - // Write the child pointer on the parent node. - // childA writes .left; childB writes .right. - if (is_child_a) { - lbvh[parent].left = current_node; - info.left_range = left_key; - } else { - lbvh[parent].right = current_node; - info.right_range = right_key; - } - - // Atomic arrival gate: the first thread to reach this parent - // stops; the second thread proceeds (it now knows both children - // are complete). - - if (info.visitation_count++ == 0) { - // this is the first thread that arrived at this - // node -> finished - break; - } - // this is the second thread that arrived at this node, - // both children are computed -> compute aabb union and - // continue - assert(lbvh[parent].is_inner()); - const Node& child_a = lbvh[lbvh[parent].left]; - const Node& child_b = lbvh[lbvh[parent].right]; - lbvh[parent].aabb_min = child_a.aabb_min.min(child_b.aabb_min); - lbvh[parent].aabb_max = child_a.aabb_max.max(child_b.aabb_max); - - // Compute rightmost leaf: max of children's rightmost - rightmost_leaves[parent] = std::max( - rightmost_leaves[lbvh[parent].left], - rightmost_leaves[lbvh[parent].right]); - - // Reconstruct the full key range for the parent. - left_key = construction_infos[parent].left_range; - right_key = construction_infos[parent].right_range; - current_node = parent; - - if (left_key == 0 && right_key == N_LEAVES - 1) { - // only one thread should reach the root - int expected = -1; - [[maybe_unused]] bool set = - root_idx.compare_exchange_strong( - expected, current_node); - assert(set); - break; // root AABB is complete - } + const size_t box_id = morton_codes[i].box_id; + details::init_leaf_node( + i, N_LEAVES, box_id, boxes[box_id].min, boxes[box_id].max, + lbvh.data(), rightmost_leaves.data()); + + const int root = details::build_hierarchy_from_leaf( + i, N_LEAVES, [&](int k) { return morton_codes[k].morton_code; }, + lbvh.data(), rightmost_leaves.data(), construction_infos.data(), + // std::atomic's post-increment is sequentially consistent, so + // it already orders this thread's writes before the arrival + // and the arrival before its later reads. + [](std::atomic& count) { return count++; }); + + if (root >= 0) { + // Only one thread should ever reach the root. + int expected = -1; + [[maybe_unused]] const bool set = + root_idx.compare_exchange_strong(expected, root); + assert(set); } }); } // --- Move the root to index 0 so traversal can start there. --- // In the Apetrei layout the root's index equals the global split position, - // which is generally != 0. We swap the root node into position 0 and patch - // up the single affected child pointer. - // - // Key invariant (Apetrei): node 0's subtree always has left_key=0, so it is - // only ever written as a LEFT child — meaning no internal node ever has - // right==0. Therefore swapping node 0 cannot create a spurious - // is_inner_marker==0 (which would look like a leaf). + // which is generally != 0. const int root = root_idx.load(); if (root > 0) { IPC_TOOLKIT_PROFILE_BLOCK("swap_root_to_zero"); - std::swap(lbvh[0], lbvh[root]); - std::swap(rightmost_leaves[0], rightmost_leaves[root]); - - // The root (now at 0) is never any node's child, so no pointer - // references R that needs rewriting to 0. The only pointers that - // referenced 0 (the old node-0) must be rewritten to R. And since old - // node-0 was only ever a LEFT child (see invariant above), we only need - // to patch .left pointers. + details::swap_root_to_zero(lbvh.data(), rightmost_leaves.data(), root); + tbb::parallel_for(size_t(0), lbvh.size(), [&](size_t i) { - if (lbvh[i].is_inner() && lbvh[i].left == 0) { - lbvh[i].left = root; - } + details::patch_left_pointer(lbvh[i], root); }); } } @@ -366,6 +216,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, @@ -375,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()) { - attempt_add_candidate( - query, child_l, can_collide, candidates); - } - if (intersects_r && child_r.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_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 @@ -808,8 +592,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 +602,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 +611,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 +622,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 +633,8 @@ bool LBVH::can_faces_collide(size_t fai, size_t fbi) const assert(fbi < face_vertex_ids.size()); const auto& [fb0i, fb1i, fb2i] = face_vertex_ids[fbi]; - const bool share_endpoint = fa0i == fb0i || fa0i == fb1i || fa0i == fb2i - || fa1i == fb0i || fa1i == fb1i || fa1i == fb2i || fa2i == fb0i - || fa2i == fb1i || fa2i == fb2i; - - return !share_endpoint - && (can_vertices_collide(fa0i, fb0i) || can_vertices_collide(fa0i, fb1i) - || can_vertices_collide(fa0i, fb2i) - || can_vertices_collide(fa1i, fb0i) - || can_vertices_collide(fa1i, fb1i) - || can_vertices_collide(fa1i, fb2i) - || can_vertices_collide(fa2i, fb0i) - || can_vertices_collide(fa2i, fb1i) - || can_vertices_collide(fa2i, fb2i)); + return details::can_faces_collide( + fa0i, fa1i, fa2i, fb0i, fb1i, fb2i, can_vertices_collide); } } // namespace ipc \ No newline at end of file diff --git a/src/ipc/broad_phase/lbvh.hpp b/src/ipc/broad_phase/lbvh.hpp index 64f37de77..4274e61d8 100644 --- a/src/ipc/broad_phase/lbvh.hpp +++ b/src/ipc/broad_phase/lbvh.hpp @@ -1,5 +1,6 @@ #pragma once +#include // for IPC_TOOLKIT_HOST_DEVICE #include #include @@ -48,14 +49,20 @@ class LBVH : public BroadPhase { #pragma GCC diagnostic pop + // These are host/device so the CUDA broad phase (ipc::cuda::LBVH) + // traverses the same node with the same predicates as the CPU one. + /// @brief Check if this node is an inner node. - bool is_inner() const { return is_inner_marker; } + IPC_TOOLKIT_HOST_DEVICE bool is_inner() const + { + return is_inner_marker; + } /// @brief Check if this node is a leaf node. - bool is_leaf() const { return !is_inner(); } + IPC_TOOLKIT_HOST_DEVICE bool is_leaf() const { return !is_inner(); } /// @brief Check if this node is valid. - bool is_valid() const + IPC_TOOLKIT_HOST_DEVICE bool is_valid() const { return is_inner() ? (left != INVALID_POINTER && right != INVALID_POINTER) @@ -63,7 +70,7 @@ class LBVH : public BroadPhase { } /// @brief Check if this node's AABB intersects with another node's AABB. - bool intersects(const Node& other) const + IPC_TOOLKIT_HOST_DEVICE bool intersects(const Node& other) const { return (aabb_min <= other.aabb_max).all() && (other.aabb_min <= aabb_max).all(); @@ -84,18 +91,26 @@ class LBVH : public BroadPhase { /// Used to skip subtrees during triangular (self-collision) traversal. using RightmostLeaves = std::vector>; -private: - struct ConstructionInfo { + /// @brief Per-internal-node scratch for the bottom-up build. + /// @tparam Counter The visitation counter's type: std::atomic here, + /// and a plain int for the device build in ipc::cuda::LBVH, where + /// atomicAdd() supplies the atomicity. + /// @see ipc::details::build_hierarchy_from_leaf + template struct ConstructionInfo { /// @brief Left range endpoint passed up by the left child. - int32_t left_range; + int left_range; /// @brief Right range endpoint passed up by the right child. - int32_t right_range; + int right_range; /// @brief Number of threads that arrived at this node. - std::atomic visitation_count; + Counter visitation_count; }; - using ConstructionInfos = - std::vector>; +private: + using HostConstructionInfo = ConstructionInfo>; + + using ConstructionInfos = std::vector< + HostConstructionInfo, + DefaultInitAllocator>; public: LBVH(); diff --git a/src/ipc/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/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 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 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_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 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..a05a49f57 --- /dev/null +++ b/tests/src/tests/broad_phase/test_gpu_lbvh.cu @@ -0,0 +1,354 @@ +// 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 +#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); + } +} + +// 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 diff --git a/tests/src/tests/broad_phase/test_lbvh.cpp b/tests/src/tests/broad_phase/test_lbvh.cpp index 08d4b39eb..d68b254ae 100644 --- a/tests/src/tests/broad_phase/test_lbvh.cpp +++ b/tests/src/tests/broad_phase/test_lbvh.cpp @@ -3,10 +3,15 @@ #include #include +#include #include #include #include +#ifdef IPC_TOOLKIT_WITH_CUDA +#include +#endif + #include #include @@ -285,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]") @@ -344,6 +464,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 +526,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