diff --git a/include/uipc/constitution/embedded_collision_mesh.h b/include/uipc/constitution/embedded_collision_mesh.h new file mode 100644 index 000000000..b0641f142 --- /dev/null +++ b/include/uipc/constitution/embedded_collision_mesh.h @@ -0,0 +1,30 @@ +#pragma once +#include +#include + +namespace uipc::constitution +{ +// Barycentric embedding of a dense surface mesh onto a coarse tet mesh. +// +// Usage (in IIPCSolver::build_scene): +// EmbeddedCollisionMesh ecm; +// ecm.apply_to(tet_sc, surface_sc); +// +// The surface_sc must be passive (no FEM constitution). +// The tet_sc must already have a FEM constitution applied. +// +// At runtime the CUDA backend reads attributes written here: +// ecm_tet_index (IndexT per surface vertex) — containing tet index +// ecm_bary (Vector4 per surface vertex) — barycentric coords +// ecm_driven (IndexT meta, value=1) — marks the mesh as driven +// +// Reference: SOFA BarycentricMapping — apply() / applyJT() +class UIPC_CONSTITUTION_API EmbeddedCollisionMesh +{ + public: + // Build barycentric embedding CPU-side (point-in-tet query). + // Returns false if any surface vertex falls outside all tets (clamped). + bool apply_to(geometry::SimplicialComplex& tet_sc, + geometry::SimplicialComplex& surface_sc) const; +}; +} diff --git a/src/backends/cuda/engine/advance_ipc.cu b/src/backends/cuda/engine/advance_ipc.cu index c48792977..3062e6140 100644 --- a/src/backends/cuda/engine/advance_ipc.cu +++ b/src/backends/cuda/engine/advance_ipc.cu @@ -309,6 +309,9 @@ void SimEngine::advance() step_animation(); m_time_integrator_manager->predict_dof(); + // Forward pass: update passive surface positions from tet DOFs + event_before_collision_detection(); + // 3. Adaptive Parameter Calculation detect_dcd_candidates(); compute_adaptive_kappa(); @@ -340,6 +343,8 @@ void SimEngine::advance() m_state = SimEngineState::ComputeDyTopoEffect; compute_dytopo_effect(); + // Backward pass: scatter surface contact forces to tet DOFs + event_after_contact_assembly(); // 4) Solve Global Linear System => dx = A^-1 * b m_state = SimEngineState::SolveGlobalLinearSystem; diff --git a/src/backends/cuda/engine/sim_engine.cu b/src/backends/cuda/engine/sim_engine.cu index 23af7fec3..44ae46f7f 100644 --- a/src/backends/cuda/engine/sim_engine.cu +++ b/src/backends/cuda/engine/sim_engine.cu @@ -98,6 +98,18 @@ void SimEngine::event_write_scene() action(); } +void SimEngine::event_before_collision_detection() +{ + for(auto& action : m_on_before_collision_detection.view()) + action(); +} + +void SimEngine::event_after_contact_assembly() +{ + for(auto& action : m_on_after_contact_assembly.view()) + action(); +} + void SimEngine::dump_global_surface() { BackendPathTool tool{workspace()}; diff --git a/src/backends/cuda/finite_element/constraints/embedded_collision_mesh.cu b/src/backends/cuda/finite_element/constraints/embedded_collision_mesh.cu new file mode 100644 index 000000000..ce09b8928 --- /dev/null +++ b/src/backends/cuda/finite_element/constraints/embedded_collision_mesh.cu @@ -0,0 +1,442 @@ +// EmbeddedCollisionMesh — GPU backend +// +// Barycentric coupling: passive dense surface mesh driven by coarse FEM tet mesh. +// +// Forward (event_before_collision_detection): x_surface[i] = J * x_tet +// Backward (event_after_contact_assembly): grad_tet += Jᵀ * grad_surface +// +// Components: +// EmbeddedCollisionMeshVertexReporter — registers surface verts in GlobalVertexManager (UID=3) +// EmbeddedCollisionMeshSystem — holds GPU buffers, drives forward/backward passes +// EmbeddedCollisionMeshForceReporter — injects tet_extra_grad into FEM linear system +// +// Reference: SOFA CudaBarycentricMapping — apply() / applyJT() + +#include +#include +#include +#include +#include +#include +#include + +namespace uipc::backend::cuda +{ + +// ================================================================ +// Shared GPU buffers +// ================================================================ +struct EmbeddedCollisionMeshData +{ + muda::DeviceBuffer tet_index; // [n_surface] — which tet each surface vert belongs to + muda::DeviceBuffer bary; // [n_surface] — barycentric coords + muda::DeviceBuffer tet_cells; // [n_tets] — tet connectivity + muda::DeviceBuffer surface_pos_init; // [n_surface] — initial surface positions (for reporter) + muda::DeviceBuffer surface_grad_dense; // [n_surface] — scattered contact grads + muda::DeviceBuffer tet_extra_grad; // [n_tet_verts] — accumulated grad for FEM + + IndexT tet_vertex_offset = 0; + IndexT surface_vertex_offset = 0; + SizeT n_tet_verts = 0; + SizeT n_surface_verts = 0; + SizeT n_tets = 0; +}; + +// ================================================================ +// CUDA kernels — file-scope to satisfy nvcc __device__ lambda rules +// ================================================================ + +__global__ void kernel_ecm_forward( + const Vector3* tet_pos, + const Vector4i* tet_cells, + const IndexT* tet_index, + const Vector4* bary, + Vector3* surface_pos, + int n_surface) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if(i >= n_surface) + return; + + Vector4i cell = tet_cells[tet_index[i]]; + surface_pos[i] = bary[i](0) * tet_pos[cell[0]] + + bary[i](1) * tet_pos[cell[1]] + + bary[i](2) * tet_pos[cell[2]] + + bary[i](3) * tet_pos[cell[3]]; +} + +// atomicAdd on Float* flat — Vector3 layout not guaranteed contiguous for CUDA atomics +__global__ void kernel_ecm_backward( + const Vector3* surface_grad, + const Vector4i* tet_cells, + const IndexT* tet_index, + const Vector4* bary, + Float* tet_extra_grad_flat, + int n_surface) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if(i >= n_surface) + return; + + Vector3 f = surface_grad[i]; + Vector4i cell = tet_cells[tet_index[i]]; + for(int j = 0; j < 4; ++j) + { + Float w = bary[i](j); + int v = cell[j]; + atomicAdd(&tet_extra_grad_flat[v * 3 + 0], w * f(0)); + atomicAdd(&tet_extra_grad_flat[v * 3 + 1], w * f(1)); + atomicAdd(&tet_extra_grad_flat[v * 3 + 2], w * f(2)); + } +} + +// Scatter CDoubletVectorView COO entries → dense surface_grad_dense +__global__ void kernel_coo_to_dense( + const int* coo_indices, + const Float* coo_values_flat, // stride 3 (reinterpret of Segment[]) + int coo_count, + IndexT surface_off, + int n_surface, + Float* dense_flat) // stride 3 (reinterpret of Vector3[]) +{ + int I = blockIdx.x * blockDim.x + threadIdx.x; + if(I >= coo_count) + return; + + int l_i = coo_indices[I] - (int)surface_off; + if(l_i < 0 || l_i >= n_surface) + return; + + atomicAdd(&dense_flat[l_i * 3 + 0], coo_values_flat[I * 3 + 0]); + atomicAdd(&dense_flat[l_i * 3 + 1], coo_values_flat[I * 3 + 1]); + atomicAdd(&dense_flat[l_i * 3 + 2], coo_values_flat[I * 3 + 2]); +} + +// ================================================================ +// Forward declarations +// ================================================================ +class EmbeddedCollisionMeshSystem; +class EmbeddedCollisionMeshForceReporter; + +// ================================================================ +// EmbeddedCollisionMeshVertexReporter +// +// Registers the passive surface mesh in GlobalVertexManager. +// UID=3 — placed after ABD(0), FEM(1), HalfPlane(2) in the +// sorted reporter list, so surface verts are appended at the end. +// +// report_displacements() → zero: the surface is driven by the +// forward pass kernel writing directly into GlobalVertexManager +// positions, not by the dx system. CCD still sees the surface +// because AABB is built from [pos, pos+dx*alpha]; with dx=0 the +// AABB degenerates to a point but is still d_hat-expanded and +// inserted in the BVH — no silent skip. +// ================================================================ +class EmbeddedCollisionMeshVertexReporter final : public VertexReporter +{ + public: + using VertexReporter::VertexReporter; + + EmbeddedCollisionMeshSystem* m_ecm_system = nullptr; + + void do_build(BuildInfo&) override; + + void do_report_count(VertexCountInfo& info) override; + + void do_report_attributes(VertexAttributeInfo& info) override; + + void do_report_displacements(VertexDisplacementInfo& info) override + { + // Surface verts are repositioned each Newton iter by kernel_ecm_forward + // writing directly into the GlobalVertexManager positions buffer. + // The dx system is bypassed — displacement is always zero here. + info.displacements().fill(Vector3::Zero()); + } + + U64 get_uid() const noexcept override { return 3; } +}; + +// ================================================================ +// EmbeddedCollisionMeshSystem +// +// Owns GPU buffers and drives the two passes. +// All methods that access SimSystem API (world(), require()) +// are instance methods of this class, which inherits SimSystem. +// ================================================================ +class EmbeddedCollisionMeshSystem final : public SimSystem +{ + public: + using SimSystem::SimSystem; + + EmbeddedCollisionMeshData data; + + // CPU-side initial positions — needed by the vertex reporter at frame 0 + std::vector h_surface_pos_init; + + // Called by EmbeddedCollisionMeshForceReporter::do_assemble() + void assemble(FEMLinearSubsystemReporter::AssembleInfo& info) + { + auto& d = data; + if(d.n_tet_verts == 0) + return; + + muda::ParallelFor() + .file_line(__FILE__, __LINE__) + .apply( + d.n_tet_verts, + [grads = info.gradients().viewer().name("ecm_grad"), + extra = d.tet_extra_grad.cviewer().name("ecm_extra"), + t_off = (int)d.tet_vertex_offset] __device__(int i) mutable + { + Vector3 g = extra(i); + if(g.squaredNorm() > Float(0)) + grads(i).write(t_off + i, g); + }); + } + + protected: + void do_build() override + { + on_init_scene([this]() { _init(*this); }); + on_before_collision_detection([this]() { _forward(*this); }); + on_after_contact_assembly([this]() { _scatter_surface_grad(*this); }); + } + + private: + static void _init(EmbeddedCollisionMeshSystem& self) + { + auto& d = self.data; + auto geo_slots = self.world().scene().geometries(); + + const geometry::SimplicialComplex* surface_sc = nullptr; + const geometry::SimplicialComplex* tet_sc = nullptr; + + for(SizeT gi = 0; gi < geo_slots.size(); ++gi) + { + auto* sc = + geo_slots[gi]->geometry().as(); + if(!sc || !sc->meta().find("ecm_driven")) + continue; + + surface_sc = sc; + auto tet_id_attr = sc->meta().find("ecm_tet_geo_id"); + if(!tet_id_attr) + break; + + IndexT tet_gi = tet_id_attr->view()[0]; + if(tet_gi >= 0 && tet_gi < (IndexT)geo_slots.size()) + tet_sc = geo_slots[tet_gi] + ->geometry() + .as(); + break; + } + + if(!surface_sc || !tet_sc) + return; + + _init_buffers(self, *surface_sc, *tet_sc); + } + + static void _init_buffers(EmbeddedCollisionMeshSystem& self, + const geometry::SimplicialComplex& surface_sc, + const geometry::SimplicialComplex& tet_sc) + { + auto& d = self.data; + + auto idx_attr = surface_sc.vertices().find("ecm_tet_index"); + auto bary_attr = surface_sc.vertices().find("ecm_bary"); + if(!idx_attr || !bary_attr) + return; + + auto idx_view = idx_attr->view(); + auto bary_view = bary_attr->view(); + d.n_surface_verts = idx_view.size(); + + d.tet_index.resize(d.n_surface_verts); + d.tet_index.view().copy_from(idx_view.data()); + + d.bary.resize(d.n_surface_verts); + d.bary.view().copy_from(bary_view.data()); + + // global_vertex_offset is written on ALL geometries by their vertex reporters + // (both FEM and passive). It is available here because _init runs in + // on_init_scene, which fires after all reporters have initialized. + auto tet_off_slot = + tet_sc.meta().find(builtin::global_vertex_offset); + if(!tet_off_slot) + return; + d.tet_vertex_offset = tet_off_slot->view()[0]; + d.n_tet_verts = tet_sc.positions().view().size(); + + auto surf_off_slot = + surface_sc.meta().find(builtin::global_vertex_offset); + if(!surf_off_slot) + return; + d.surface_vertex_offset = surf_off_slot->view()[0]; + + auto& fem = self.require(); + auto tets_v = fem.tets(); + d.n_tets = tets_v.size(); + d.tet_cells.resize(d.n_tets); + d.tet_cells.view().copy_from(tets_v); + + d.surface_grad_dense.resize(d.n_surface_verts, Vector3::Zero()); + d.tet_extra_grad.resize(d.n_tet_verts, Vector3::Zero()); + } + + static void _forward(EmbeddedCollisionMeshSystem& self) + { + auto& d = self.data; + if(d.n_surface_verts == 0) + return; + + auto* gvm = &self.require(); + auto all_pos = gvm->positions(); + + // Write directly into the GlobalVertexManager positions buffer at the + // surface mesh offset — this is what the CCD and contact systems see. + Vector3* surf_ptr = const_cast(all_pos.data()) + d.surface_vertex_offset; + const Vector3* tet_ptr = all_pos.data() + d.tet_vertex_offset; + + int block = 256; + int grid = ((int)d.n_surface_verts + block - 1) / block; + kernel_ecm_forward<<>>( + tet_ptr, + d.tet_cells.data(), + d.tet_index.data(), + d.bary.data(), + surf_ptr, + (int)d.n_surface_verts); + } + + static void _scatter_surface_grad(EmbeddedCollisionMeshSystem& self) + { + auto& d = self.data; + if(d.n_surface_verts == 0) + return; + + auto* dytopo = &self.require(); + auto coo = dytopo->gradients(); + + d.surface_grad_dense.view().fill(Vector3::Zero()); + d.tet_extra_grad.view().fill(Vector3::Zero()); + + int coo_count = (int)coo.doublet_count(); + if(coo_count == 0) + return; + + { + int block = 256; + int grid = (coo_count + block - 1) / block; + kernel_coo_to_dense<<>>( + coo.indices().data(), + reinterpret_cast(coo.values().data()), + coo_count, + (IndexT)d.surface_vertex_offset, + (int)d.n_surface_verts, + reinterpret_cast(d.surface_grad_dense.data())); + } + + { + int block = 256; + int grid = ((int)d.n_surface_verts + block - 1) / block; + kernel_ecm_backward<<>>( + d.surface_grad_dense.data(), + d.tet_cells.data(), + d.tet_index.data(), + d.bary.data(), + reinterpret_cast(d.tet_extra_grad.data()), + (int)d.n_surface_verts); + } + } +}; + +REGISTER_SIM_SYSTEM(EmbeddedCollisionMeshSystem); + +// ================================================================ +// EmbeddedCollisionMeshVertexReporter — method bodies +// (defined after EmbeddedCollisionMeshSystem is complete) +// ================================================================ + +void EmbeddedCollisionMeshVertexReporter::do_build(BuildInfo&) +{ + m_ecm_system = &require(); + // VertexReporter::do_build() (the non-virtual base) calls add_reporter(this) + // automatically — no need to call it here. +} + +void EmbeddedCollisionMeshVertexReporter::do_report_count(VertexCountInfo& info) +{ + info.count(m_ecm_system->data.n_surface_verts); +} + +void EmbeddedCollisionMeshVertexReporter::do_report_attributes(VertexAttributeInfo& info) +{ + auto& d = m_ecm_system->data; + + if(info.frame() == 0) + { + // Write initial positions from CPU buffer + auto& h_pos = m_ecm_system->h_surface_pos_init; + if(!h_pos.empty()) + info.positions().copy_from(h_pos.data()); + + // The surface mesh has no contact element — use default (0) + // Thickness from d_hat default; no special body_id needed + + // Write global_vertex_offset on the surface geometry meta attribute + // so EmbeddedCollisionMeshSystem::_init_buffers() can read it. + auto global_offset = info.coindices().offset(); + auto geo_slots = world().scene().geometries(); + for(SizeT gi = 0; gi < geo_slots.size(); ++gi) + { + auto* sc = geo_slots[gi]->geometry().as(); + if(!sc || !sc->meta().find("ecm_driven")) + continue; + + auto gvo = sc->meta().find(builtin::global_vertex_offset); + if(!gvo) + gvo = sc->meta().create(builtin::global_vertex_offset); + view(*gvo)[0] = (IndexT)global_offset; + break; + } + } + // Frames > 0: positions are updated by kernel_ecm_forward directly into the + // GlobalVertexManager buffer — no copy needed here. +} + +REGISTER_SIM_SYSTEM(EmbeddedCollisionMeshVertexReporter); + +// ================================================================ +// EmbeddedCollisionMeshForceReporter +// ================================================================ +class EmbeddedCollisionMeshForceReporter final + : public FEMLinearSubsystemReporter +{ + public: + using FEMLinearSubsystemReporter::FEMLinearSubsystemReporter; + + SimSystemSlot m_ecm_system; + + void do_build(BuildInfo& info) override + { + m_ecm_system = require(); + } + + void do_init(InitInfo& info) override {} + + void do_report_extent(ReportExtentInfo& info) override + { + info.gradient_count(m_ecm_system->data.n_tet_verts); + info.hessian_count(0); + } + + void do_assemble(AssembleInfo& info) override + { + m_ecm_system->assemble(info); + } +}; + +REGISTER_SIM_SYSTEM(EmbeddedCollisionMeshForceReporter); + +} // namespace uipc::backend::cuda diff --git a/src/backends/cuda/sim_engine.h b/src/backends/cuda/sim_engine.h index ad832b70b..840c0736c 100644 --- a/src/backends/cuda/sim_engine.h +++ b/src/backends/cuda/sim_engine.h @@ -80,6 +80,10 @@ class SimEngine final : public backend::SimEngine void event_rebuild_scene(); SimActionCollection m_on_write_scene; void event_write_scene(); + SimActionCollection m_on_before_collision_detection; + void event_before_collision_detection(); + SimActionCollection m_on_after_contact_assembly; + void event_after_contact_assembly(); private: diff --git a/src/backends/cuda/sim_system.cpp b/src/backends/cuda/sim_system.cpp index 667cc3669..55c460b1e 100644 --- a/src/backends/cuda/sim_system.cpp +++ b/src/backends/cuda/sim_system.cpp @@ -33,6 +33,18 @@ void SimSystem::on_write_scene(std::function&& action) noexcept engine().m_on_write_scene.register_action(*this, std::move(action)); } +void SimSystem::on_before_collision_detection(std::function&& action) noexcept +{ + check_state(SimEngineState::BuildSystems, "on_before_collision_detection()"); + engine().m_on_before_collision_detection.register_action(*this, std::move(action)); +} + +void SimSystem::on_after_contact_assembly(std::function&& action) noexcept +{ + check_state(SimEngineState::BuildSystems, "on_after_contact_assembly()"); + engine().m_on_after_contact_assembly.register_action(*this, std::move(action)); +} + SimEngine& SimSystem::engine() noexcept { return static_cast(Base::engine()); diff --git a/src/backends/cuda/sim_system.h b/src/backends/cuda/sim_system.h index de58013f2..da7f4d62d 100644 --- a/src/backends/cuda/sim_system.h +++ b/src/backends/cuda/sim_system.h @@ -43,6 +43,10 @@ class SimSystem : public backend::SimSystem */ void on_write_scene(std::function&& action) noexcept; + void on_before_collision_detection(std::function&& action) noexcept; + + void on_after_contact_assembly(std::function&& action) noexcept; + WorldVisitor& world() noexcept; void check_state(SimEngineState state, std::string_view function_name) noexcept; diff --git a/src/constitution/embedded_collision_mesh.cpp b/src/constitution/embedded_collision_mesh.cpp new file mode 100644 index 000000000..c18270413 --- /dev/null +++ b/src/constitution/embedded_collision_mesh.cpp @@ -0,0 +1,145 @@ +#include +#include +#include +#include + +#include +#include + +namespace uipc::constitution +{ +// --------------------------------------------------------------- +// CPU helpers — point-in-tet and barycentric coordinates +// --------------------------------------------------------------- + +static Float scalar_triple(const Vector3& u, + const Vector3& v, + const Vector3& w) +{ + return u.dot(v.cross(w)); +} + +// Compute barycentric coords of p in tet {a,b,c,d}. +// Returns false if the tet is degenerate. +static bool bary_in_tet(const Vector3& p, + const Vector3& a, + const Vector3& b, + const Vector3& c, + const Vector3& d, + Vector4& bary) +{ + Float vol = scalar_triple(b - a, c - a, d - a); + if(std::abs(vol) < Float(1e-12)) + return false; + + Float inv = Float(1) / vol; + Float l0 = scalar_triple(b - p, c - p, d - p) * inv; + Float l1 = scalar_triple(a - p, d - p, c - p) * inv; + Float l2 = scalar_triple(a - p, b - p, d - p) * inv; + Float l3 = Float(1) - l0 - l1 - l2; + + bary = Vector4(l0, l1, l2, l3); + return true; +} + +static bool is_inside(const Vector4& bary, Float eps = Float(1e-5)) +{ + return bary[0] >= -eps && bary[1] >= -eps + && bary[2] >= -eps && bary[3] >= -eps; +} + +// --------------------------------------------------------------- +// EmbeddedCollisionMesh +// --------------------------------------------------------------- + +bool EmbeddedCollisionMesh::apply_to(geometry::SimplicialComplex& tet_sc, + geometry::SimplicialComplex& surface_sc) const +{ + auto tet_pos_view = tet_sc.positions().view(); + auto tets_attr = tet_sc.tetrahedra().find(builtin::topo); + if(!tets_attr) + return false; + auto tets_view = tets_attr->view(); + auto surf_pos_view = surface_sc.positions().view(); + + SizeT n_surface = surf_pos_view.size(); + SizeT n_tets = tets_view.size(); + + // Create or reset output attributes on the surface mesh. + auto tet_index_attr = surface_sc.vertices().find("ecm_tet_index"); + if(!tet_index_attr) + tet_index_attr = + surface_sc.vertices().create("ecm_tet_index", -1); + + auto bary_attr = surface_sc.vertices().find("ecm_bary"); + if(!bary_attr) + bary_attr = + surface_sc.vertices().create("ecm_bary", Vector4::Zero()); + + if(!surface_sc.meta().find("ecm_driven")) + surface_sc.meta().create("ecm_driven", 1); + + auto tet_index_view = geometry::view(*tet_index_attr); + auto bary_view = geometry::view(*bary_attr); + + bool all_embedded = true; + + for(SizeT si = 0; si < n_surface; ++si) + { + const Vector3& p = surf_pos_view[si]; + + // Brute-force O(n_surface * n_tets) — fine at init. + // Replace with BVH if tet count > ~10k. + bool found = false; + Float best_min_bary = -std::numeric_limits::max(); + IndexT best_tet = 0; + Vector4 best_bary = Vector4::Constant(Float(0.25)); + + for(SizeT ti = 0; ti < n_tets; ++ti) + { + const Vector4i& cell = tets_view[ti]; + const Vector3& a = tet_pos_view[cell[0]]; + const Vector3& b = tet_pos_view[cell[1]]; + const Vector3& c = tet_pos_view[cell[2]]; + const Vector3& d = tet_pos_view[cell[3]]; + + Vector4 bary; + if(!bary_in_tet(p, a, b, c, d, bary)) + continue; + + if(is_inside(bary)) + { + tet_index_view[si] = (IndexT)ti; + bary_view[si] = bary; + found = true; + break; + } + + Float min_b = bary.minCoeff(); + if(min_b > best_min_bary) + { + best_min_bary = min_b; + best_tet = (IndexT)ti; + best_bary = bary; + } + } + + if(!found) + { + // Fallback: clamp to tet boundary and renormalize. + best_bary = best_bary.cwiseMax(Float(0)); + Float sum = best_bary.sum(); + if(sum > Float(1e-12)) + best_bary /= sum; + else + best_bary = Vector4::Constant(Float(0.25)); + + tet_index_view[si] = best_tet; + bary_view[si] = best_bary; + all_embedded = false; + } + } + + return all_embedded; +} +}