Block-accelerated contact Hessian assembly via MeshFEMSparse - #246
Draft
zfergus wants to merge 18 commits into
Draft
Block-accelerated contact Hessian assembly via MeshFEMSparse#246zfergus wants to merge 18 commits into
zfergus wants to merge 18 commits into
Conversation
Adds a reusable contact-scene fixture (8 scenes spanning 390 to 512k collisions, each padded with interior vertices so to_full_dof performs a genuine surface-to-volume scatter) and Catch2 benchmarks that isolate the three costs of contact Hessian/gradient assembly: 1. per-collision (local) derivative evaluation, 2. global assembly (triplets + setFromTriplets), 3. the reduced-DOF map (CollisionMesh::to_full_dof). Baseline findings: local derivative evaluation is only 1.8-7.6% of Hessian cost; the rest is assembly bookkeeping (42-62%) and to_full_dof SpGEMMs (30-56%). On the largest scene (puffer-ball, 512k collisions) bookkeeping costs ~560 ms per Newton iteration vs 21 ms of derivative math. Also adds a memory-guarded scene probe ([assembly-probe], hidden) that counts broad-phase candidates before building the collision set, since an oversized dhat can exhaust host memory. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds an in_full_dof parameter to Potential<T>::gradient/hessian. When the mesh's DOF map is a pure selection matrix (the default; tracked by the new CollisionMesh::is_selection_dof_map()), stencil vertex IDs are remapped to full-mesh IDs during triplet generation, producing the full-DOF result directly instead of applying to_full_dof afterwards. This eliminates the two serial SpGEMMs (S^T H S), which were 30-56% of end-to-end Hessian cost. With a user-provided displacement map, in_full_dof falls back to to_full_dof internally, so the flag is always safe. Measured end-to-end Hessian speedups: 1.29-1.83x across 8 scenes (390-512k collisions). Gradient folding is not beneficial on large scenes (the thread-local accumulators grow to full_ndof while the SpMV saved is cheap) and is left off by default; documented in the benchmark. Note: the defensive storage-empty path in hessian() now returns a correctly-sized (ndof x ndof) empty matrix instead of 0x0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…e 2) Extracts the global-matrix construction out of Potential<T>::hessian into an abstract HessianAssembler interface (begin / thread-safe add_local_hessian / end). The historical triplet + setFromTriplets path moves verbatim into TripletHessianAssembler, and hessian() becomes a thin wrapper over the new public Potential<T>::assemble_hessian driver, which also owns the Phase 1 full-DOF stencil remap so every future backend gets it for free. No behavior change; benchmarks confirm collision-DOF assembly times are within run-to-run noise of the previous implementation on all 8 scenes. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds MeshFEMHessianAssembler, a HessianAssembler backed by MeshFEMSparse's block-CSC data structures (Mohammadian et al., SIGGRAPH 2026): begin() builds a block sparsity pattern from the collision stencils and add_local_hessian() scatters each local Hessian directly into the value array via MeshFEM's sorted column-merge with per-column spin locks — no triplets, no setFromTriplets. Guarded by IPC_TOOLKIT_WITH_MESHFEM_SPARSE (default OFF). The HessianAssembler seam gains a StencilGetter argument to begin() so pattern-based backends can see stencils up front. The dependency is fetched with CPM DOWNLOAD_ONLY (pinned SHAs + SHA256 archive hashes) and compiled into a minimal static target (matrix data structures and assembly only, no sparse direct solvers), avoiding upstream's PUBLIC -fvisibility=hidden, its solver sources (which clash with Eigen 5's BLAS declarations), and its transitive dependency fetching. Compatibility notes: - MeshFEM targets Eigen 3.4; Eigen 5 removed internal::make_coherent, which MeshFEMCore/AutomaticDifferentiation.hh references (included by SparseMatrices.hh at the root of the header chain). A force-included shim (meshfem_eigen_compat.hpp) reimplements the Eigen 3.4 semantics. - BlockCSCHessian::toEigen/toScalar read out of bounds on empty block columns (impossible for FE Hessians, ubiquitous for contact Hessians: most vertices are collision-free), causing intermittent segfaults. Replaced with a custom direct block-CSC -> symmetric Eigen conversion, which is also ~2x faster than upstream's two-step expansion. Measured on 8 scenes (390-512k collisions), full-DOF Hessian, pattern rebuilt every call: 2.5-11x end-to-end vs the triplet path to an Eigen matrix, 3-15x to the block-CSC format. Matches the triplet assembler to <= 1e-13 relative across scenes x PSD projection x DOF space. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
MeshFEMHessianAssembler is now designed to live across assemblies (e.g., one instance per Newton solve). begin() compares the stencils against the cached block pattern via MeshFEM's detectChangedEntries and reuses it (values-only reset + scatter) unless the contact set gained a new vertex pair or lost more than stale_block_tolerance() blocks; stale blocks assemble to explicit zeros. The Eigen conversion structure (symmetrized pattern + index arrays) is cached the same way, so get_matrix() — now returning a const reference valid until the next begin() — reduces to a parallel value refill while the pattern holds. For callers that know the collision set is identical to the previous assembly (change detection costs a sizable fraction of a rebuild on large scenes), set_assume_unchanged_stencils(true) skips detection entirely; a differing stencil count falls back to detection automatically and debug builds verify the assumption. Amortization is automatic through the existing assemble_hessian seam — no API changes beyond the new accessors. Steady-state contact Hessians (Eigen output included) reach 3.3-29x over the triplet baseline across the 8 benchmark scenes (e.g., cloth-ball 14 ms -> 0.48 ms, puffer-ball ~600 ms -> 32 ms), with reuse semantics covered by new tests (identical/shrunken/grown sets, tolerance behavior, assume-unchanged fallback). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Potential<T>::gradient now selects between two assembly strategies based on problem shape (no API change, no new dependency): - gather (new): local gradients are written to a flat per-slot buffer, a vertex->slot adjacency is built with a parallel counting sort, and each vertex sums its contributions independently. Cost scales with the number of contributions rather than ndof. - scatter+reduce (previous behavior): thread-local dense accumulators whose zero+combine cost scales with ndof. Gather is selected when out_ndof > 4 * num_collisions, the empirical crossover on the benchmark scenes: contact-sparse large meshes get gather (cloth-ball 512-612 -> 381 us, n-body 917 -> 695 us), while collision-dense scenes (rod-twist: 1.3M contributions on 120k DOF, where gather's buffer + adjacency traffic measured 1.6x worse) keep the scatter path. This also removes the Phase 1 caveat that in_full_dof gradients could be slower: with gather the accumulators no longer grow with full_ndof (cloth-ball 595 -> 444 us, puffer-ball 13.7 -> 11.3 ms folded). Summation order remains floating-point nondeterministic on both paths; sorting each gather bucket would make that path reproducible if ever needed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
IPC_TOOLKIT_WITH_MESHFEM_SPARSE now defaults to ON (auto-disabled for IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=ColMajor, which the block layout does not support), and Potential<T>::hessian() assembles through the block-CSC backend when compiled in, via a new zero-copy MeshFEMHessianAssembler::take_matrix(). The triplet path remains as the fallback when the option is off. Every existing hessian() caller gets the speedup with no code change: cloth-ball 5-6.7 -> 1.5 ms, armadillo-rollers 11-18 -> 2.2 ms, rod-twist 165-212 -> 29.5 ms, puffer-ball 375-1020 -> 48.9 ms (identical results up to floating-point summation order; full 286-test suite passes in both configurations). The dependency is now pinned to fork commits carrying the two fixes submitted upstream (MeshFEM/MeshFEMCore#1 for Eigen 5 support, MeshFEM/MeshFEMSparse#1 for an out-of-bounds read on empty block columns) -- marked TEMPORARY in the recipe; repoint to upstream SHAs once merged. This allowed deleting the force-included make_coherent compatibility shim entirely. Also: document the HessianAssembler classes in the C++ API docs (with IPC_TOOLKIT_WITH_MESHFEM_SPARSE added to Doxygen's PREDEFINED so the guarded class renders) and add MeshFEMSparse to the optional-dependency docs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Single-capital-letter names for matrices (H = Hessian, M = matrix) are the codebase's mathematical convention; NOLINT the readability-identifier-naming check on them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #246 +/- ##
==========================================
+ Coverage 96.58% 96.63% +0.05%
==========================================
Files 163 168 +5
Lines 16668 16919 +251
Branches 921 959 +38
==========================================
+ Hits 16099 16350 +251
Misses 569 569
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Adds MeshFEMHessianAssembler::block_matrix(), which returns the assembled matrix in MeshFEMSparse's native block-CSC form so a downstream user can feed it to MeshFEM's block SpMV or Cholesky factorizers instead of paying for the Eigen conversion (0.11 vs 0.30 ms on bunny, 42.8 vs 51.9 ms on puffer-ball). MeshFEM::BlockCSCHessianBase is forward declared, so our header still does not pull in MeshFEMSparse's; callers that want the block matrix include <MeshFEMSparse/BlockCSCHessian.hh> themselves and everyone else pays nothing. Binds assemble_hessian, HessianAssembler, TripletHessianAssembler, and MeshFEMHessianAssembler to Python, so Python callers can now hold an assembler across iterations and get pattern reuse (previously they were limited to the cold path inside hessian()). All three classes are py::is_final(): a Python-defined assembler would take the GIL once per collision, which is hundreds of thousands of times per assembly on the larger scenes. Exercising block_matrix() turned up a third instance of the empty-block- column assumption upstream, in visitDiagonalScalarEntries, which made trace() read the preceding column's storage (1951.93 against a dense trace of 447.82) and addDiag()/setDiag() write to the wrong entries. Fixed in the pinned fork commit alongside the other two (MeshFEM/MeshFEMSparse#1); the tests now cover trace() agreement and that addDiag() rejects a pattern with missing diagonal blocks. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The tutorial still showed to_full_dof as the only way to get full-mesh derivatives, and said nothing about holding an assembler across a Newton solve, which is where most of the speedup lives. Adds an in_full_dof example next to the existing to_full_dof one (with a note on the pure-selection requirement and the silent fallback when a displacement map is present), and a section on reusing a MeshFEMHessianAssembler: what the cached pattern covers, when it is rebuilt, block_matrix() for solvers that speak block CSC, and the assume_unchanged_stencils escape hatch and its caveat. Also drops the now-wrong "two fixes" count for the pinned forks; the MeshFEMSparse PR carries two empty-block-column fixes of its own. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
zfergus
marked this pull request as draft
August 6, 2026 20:11
hessian() now routes through whichever backend is compiled in, so the table's local%/asm%/full% columns were comparing the MeshFEM path against itself. Time the triplet assembler directly instead, and rename the column to match. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Short-circuit an empty collision set in hessian(): building a sparsity pattern to produce an all-zero matrix costs O(ndof) for nothing. Validate dim in MeshFEMHessianAssembler::begin() before dividing by it, so an unsupported dimension throws instead of trapping. Make EIGEN_DONT_VECTORIZE PUBLIC on the MeshFEMSparse target: the setting has to travel with the target, since anything including its headers must agree with how its own translation units were compiled. Return by value rather than through the ternary so the returns are implicitly moved, and include <memory> where unique_ptr is used.
zfergus
force-pushed
the
feature/meshfem-assembly
branch
from
August 7, 2026 06:37
0f14389 to
45f1983
Compare
Guards against being included twice (via MeshFEM::Sparse or MeshFEMSparse targets) and aliases MeshFEMSparse as MeshFEM::Sparse for consistent namespaced usage. Also switches CPMAddPackage calls to the gh: URI shorthand instead of manual URL/URL_HASH pairs.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Contact Hessian assembly spends almost none of its time on math. Across the eight benchmark scenes here, evaluating the local per-collision Hessians is 2-8% of
Potential::hessian(). The rest is building triplets, sorting them insidesetFromTriplets, and multiplying by the selection matrix twice into_full_dof.This replaces both halves. Local Hessians scatter straight into MeshFEMSparse's block-CSC storage (Mohammadian et al., MeshFEM: A Block-accelerated Solver for Nonlinear Finite Elements, SIGGRAPH 2026), and the full-mesh DOF map is folded into that scatter instead of applied afterwards. Callers do not have to change anything to benefit: on puffer-ball (512k collisions) a Hessian goes from 918 ms to 46 ms.
Where the time goes
to_full_dofis absent from every MeshFEM bar because folding removes it outright, and local evaluation goes from a sliver to roughly half the bar. On the largest scenes the arithmetic is now the majority of the cost, which is where a contact Hessian should sit.What each change contributes
Four steps, each measured against the triplet baseline in the same run:
to_full_dofinto assembly (1.2-1.8x). When the mesh's DOF map is a plain selection matrix, which holds unless a custom displacement map was supplied, stencil vertex IDs are remapped during assembly and the two sparse products disappear. Non-selection maps fall back internally, so the newin_full_dofflag is always safe to pass.setFromTripletssort.Most of the win lands before any reuse. Reuse pays where pattern construction dominates and adds almost nothing on rod-twist, where detection costs about what a rebuild does, which is exactly what the opt-in fast path in step 4 exists for.
Gradients
Gradient assembly picks between a gather-based per-vertex reduction and the existing scatter-plus-reduce by problem shape. Note the axis: this path costs microseconds to milliseconds against the Hessian's milliseconds to hundreds, so it was never the headline.
API
Potential::gradientandhessiantake a newin_full_dofflag. A newHessianAssemblerinterface separates local derivative evaluation from global matrix construction; the old triplet code lives on unchanged asTripletHessianAssemblerand remains the fallback when the option is off. Hold aMeshFEMHessianAssembleracrossassemble_hessiancalls to get pattern reuse, or callblock_matrix()for the native block-CSC matrix if your solver can consume it. All of it is bound to Python.Dependency status
MeshFEMSparse and MeshFEMCore are fetched with CPM
DOWNLOAD_ONLY(pinned SHAs, SHA256 hashes on the archives) and compiled into a five-file static target: matrix data structures and assembly routines, no sparse direct solver wrappers, so SuiteSparse never enters the picture. Transitively Eigen and TBB, both already built. MIT licensed.Blocker: the recipe pins my forks, which carry three fixes submitted upstream. MeshFEMCore#1 is Eigen 5 support. MeshFEMSparse#1 covers two out-of-bounds reads that share a root cause: both assume every block column has a diagonal block, which holds for FE Hessians and fails for contact Hessians, where only vertices currently in contact appear at all. Pointing back at upstream is a two-line change once those land.
Measured on Apple Silicon macOS, AppleClang 21, Release, CUDA off, via
tests/src/tests/potential/benchmark_assembly.cpp. Every figure above is a median of three runs; single runs on this machine drift 3-21%.