Block-accelerated contact Hessian assembly via MeshFEMSparse

[zfergus]

Contact Hessian assembly spends almost none of its time on math. Across the eight benchmark scenes here, evaluating the local per-collision Hessians accounts for 2-8% of Potential::hessian(); everything else is building triplets, sorting them inside setFromTriplets, and then multiplying by the selection matrix twice in to_full_dof.

This PR replaces both halves of that. 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 went from ~600 ms to 49 ms. Values match the old path to floating point rounding.

What's in this PR

One commit per step, so this can be reviewed or reverted piecewise.

1. Benchmarks: a contact-scene fixture (8 scenes, 390 to 512k collisions) plus Catch2 benchmarks that separate local derivative evaluation from global assembly from to_full_dof. Every number below comes from these.
2. Full-DOF folding: Potential::gradient/hessian take a new in_full_dof flag. When the mesh's DOF map is a plain selection matrix, which is the case unless someone passed a custom displacement map, stencil vertex IDs get remapped to full-mesh IDs during assembly and the two to_full_dof products disappear. Non-selection maps fall back to the old behavior internally, so the flag is always safe to pass.
3. A HessianAssembler interface. Potential::assemble_hessian(collisions, mesh, X, assembler, ...) separates local derivative evaluation from global matrix construction. The old triplet code moves into TripletHessianAssembler unchanged.
4. MeshFEMHessianAssembler, the block-CSC backend. It builds a block sparsity pattern from the collision stencils, then scatters local Hessians into the value array with MeshFEM's sorted column-merge and per-column locks, which skips triplet construction and the setFromTriplets sort completely.
5. Pattern reuse. The assembler now lives across assemblies, one per Newton solve. It detects contact-set changes and rebuilds the pattern only when it has to, with a tolerance for stale blocks that mirrors MeshFEM's own sparsityPatternUpdateThreshold, plus an opt-in fast path for callers who know the set is unchanged.
6. Gradient assembly picks between a gather-based per-vertex reduction and the existing scatter-plus-reduce depending on problem shape. This also fixes full-DOF folding being a pessimization for gradients.
7. IPC_TOOLKIT_WITH_MESHFEM_SPARSE defaults to ON and Potential::hessian() routes through the block backend, with the triplet path kept as the fallback.
8. MeshFEMHessianAssembler::block_matrix() hands out the matrix in its native block-CSC form, and the assemblers are bound to Python.

Dependency notes

MeshFEMSparse and MeshFEMCore are fetched with CPM DOWNLOAD_ONLY (pinned SHAs, plus SHA256 hashes on the archives) and compiled into a small static target of five source files. That gets us the matrix data structures and the assembly routines without the sparse direct solver wrappers, so SuiteSparse never enters the picture. Transitive dependencies are Eigen and TBB, both of which we already build. MIT licensed.

Integrating it turned up three upstream bugs, which I fixed and submitted rather than worked around locally. MeshFEM/MeshFEMCore#1 is Eigen 5 support: Eigen 5 removed internal::make_coherent, which MeshFEMCore's autodiff header calls. MeshFEM/MeshFEMSparse#1 covers the other two, both of which come from the same assumption that every block column has a diagonal block. That holds for FE Hessians, where every node belongs to an element, and fails for contact Hessians, where only vertices currently in contact appear at all. The recipe pins fork commits carrying all three patches, marked TEMPORARY; once those PRs land it is a two-line change to point back at upstream.

Before approving this, note that it changes the default numerical output of hessian() within floating point reordering. Any downstream test that compares Hessians for exact equality will notice.

---

Performance report

Summary

End-to-end Hessian speedup against the current triplet baseline, measured within the same run (assembly plus DOF map, Eigen output).

scene (collisions)	fold to_full_dof	interface	block-CSC cold	pattern reuse	assume unchanged
two-cubes (390)	1.29x	=	3.0x	3.1x	3.5x
bunny (570)	1.47x	=	4.0x	10.5x	10.9x
cloth-ball (3.0k)	1.75x	=	9.9x	28.9x	28.7x
cloth-funnel (5.2k)	1.44x	=	2.5x	3.5x	3.3x
armadillo-rollers (15.5k)	1.83x	=	7.9x	11.3x	13.7x
n-body (16.0k)	1.66x	=	10.4x	16.7x	22.8x
rod-twist (327k)	1.48x	=	8.4x	7.9x	12.1x
puffer-ball (512k)	1.40x	=	11.0x	18.7x*	23.6x*

The interface column is deliberately a wash; it is a refactor and I checked that it costs nothing. Steps 3 and 4 subsume the folding from step 2, which stays useful for anyone building without MeshFEMSparse. "Assume unchanged" skips change detection when the caller asserts the collision set is identical.

Puffer-ball is thermally noisy across runs, so its baseline moves around; against the fastest baseline I ever measured for it those last two columns are 13.2x and 17.6x rather than 18.7x and 23.6x. Every other number is a same-run comparison.

Setup: Apple Silicon macOS (Darwin 25.5.0), AppleClang 21, Release, CUDA off. Harness is tests/src/tests/potential/benchmark_assembly.cpp, run as ./ipc_toolkit_tests "[assembly]" --benchmark-samples 10 --benchmark-no-analysis.

Each scene mesh is padded with an equal number of unreferenced interior vertices, so to_full_dof does a real surface-to-volume scatter rather than a permutation. That is how simulators actually embed collision meshes, and without it the DOF map looks artificially cheap.

scene	collisions	collision verts (DOF)	full verts (DOF)
two-cubes (two-cubes-close.ply, dhat=1e-1)	390	244 (732)	488 (1,464)
bunny (bunny.ply, dhat=1e-2)	570	2,503 (7,509)	5,006 (15,018)
cloth-ball (cloth_ball92.ply, dhat=1e-3)	3,012	46,598 (139,794)	93,196 (279,588)
cloth-funnel (cloth-funnel/227.ply, dhat=4e-3)	5,168	9,450 (28,350)	18,900 (56,700)
armadillo-rollers (armadillo-rollers/326.ply, dhat=1e-3)	15,533	12,135 (36,405)	24,270 (72,810)
n-body (n-body-simulation/balls16_18.ply, dhat=1.78e-2)	15,969	73,960 (221,880)	147,920 (443,760)
rod-twist (rod-twist/3036.ply, dhat=1.09e-3)	327,491	39,940 (119,820)	79,880 (239,640)
puffer-ball (puffer-ball/20.ply, dhat=1.44e-4)	512,190	532,110 (1,596,330)	1,064,220 (3,192,660)

The simulation-frame scenes use dhat = 1e-3 times the bbox diagonal. Be careful raising these: collision count grows superlinearly, and an order of magnitude more dhat on cloth-ball will exhaust memory on a 64 GB machine. There is an [assembly-probe] test that counts broad-phase candidates first, which is how these values were chosen.

Baseline (triplet assembly on main)

Hessian path, Catch2 means over 10 samples:

scene	local hessians	hessian (collision DOF)	to_full_dof	end to end
two-cubes	33 µs	345 µs	146 µs	~491 µs
bunny	32 µs	467 µs	558 µs	1.02 ms
cloth-ball	181 µs	4.78 ms	7.85 ms	12.1 ms
armadillo-rollers	844 µs	10.4 ms	8.92 ms	19.2 ms

Ten assemblies of an unchanged collision set cost exactly ten times one assembly (two-cubes 3.78 ms, bunny 10.2 ms, cloth-ball 121.3 ms, armadillo-rollers 192.9 ms). Nothing is amortized across Newton iterations today, which is the entire opportunity behind step 5.

Cost breakdown, median of 5, all eight scenes. "Assembly" is the full call minus the local-only pass; "full-DOF" is the to_full_dof products.

scene	collisions	local (ms)	assembly (ms)	to_full_dof (ms)	local %	assembly %	full-DOF %
two-cubes	390	0.034	0.269	0.142	7.6%	60.6%	31.9%
bunny	570	0.043	0.585	0.582	3.5%	48.3%	48.1%
cloth-ball	3,012	0.254	6.011	7.925	1.8%	42.4%	55.9%
cloth-funnel	5,168	0.238	2.041	1.361	6.5%	56.1%	37.4%
armadillo-rollers	15,533	0.855	10.695	9.272	4.1%	51.4%	44.5%
n-body	15,969	0.623	9.548	11.417	2.9%	44.2%	52.9%
rod-twist	327,491	15.168	163.824	89.645	5.6%	61.0%	33.4%
puffer-ball	512,190	21.331	362.784	197.430	3.7%	62.4%	33.9%

Peak RSS for the full eight-scene run is 4.4 GB.

Gradients, for comparison:

scene	local gradients	gradient (collision DOF)	to_full_dof (SpMV)
two-cubes	21.8 µs	25.1 µs	1.2 µs
bunny	24.4 µs	84.0 µs	11.7 µs
cloth-ball	53.6 µs	612 µs	204 µs
armadillo-rollers	164 µs	321 µs	52 µs

So 92-98% of a contact Hessian is bookkeeping. On puffer-ball that is 560 ms of every Newton iteration. Pure assembly is the larger share at scale (61-62% on the two biggest scenes), which is what block-CSC assembly removes, while to_full_dof is 30-56% everywhere and is what folding removes. The gradient path is in much better shape, though the tbb::combinable reduction still costs 2-11x the local work on big meshes: cloth-ball spends 612 µs to produce 54 µs worth of derivatives, almost all of it zeroing and combining thread-local VectorXd(ndof) at 139k DOF.

Folding the DOF map into assembly

The folded Hessian time is essentially the collision-DOF assembly time on its own. Remapping stencil IDs costs nothing measurable and the two sparse products vanish.

scene	before (ms)	folded (ms)	speedup
two-cubes	0.57	0.44	1.29x
bunny	1.06	0.72	1.47x
cloth-ball	14.7	8.4	1.75x
cloth-funnel	3.7	2.6	1.44x
armadillo-rollers	21.1	11.5	1.83x
n-body	22.6	13.6	1.66x
rod-twist	268	181	1.48x
puffer-ball	855	610	1.40x

(That puffer-ball baseline is the noisy one again; the ratios are within-run.)

Gradients were a different story, and this surprised me. Folding made the big scenes slower:

scene	gradient + to_full_dof	folded
two-cubes	38 µs	35 µs
bunny	458 µs	285 µs
cloth-ball	730 µs	955 µs
n-body	1.27 ms	1.62 ms
rod-twist	3.97 ms	4.19 ms
puffer-ball	14.6 ms	18.0 ms

The thread-local accumulators double in size when you assemble in full DOF, and zeroing and combining them costs more than the to_full_dof SpMV you saved, which was only microseconds to a couple of milliseconds. Step 6 fixes the underlying reduction and this inverts.

On correctness: test_full_dof_assembly.cpp checks folded against mapped to 1e-13 relative on scenes with genuine interior-vertex selection, with and without PSD projection, plus the non-selection fallback and empty-collision sizing. The comparison cannot be exact because thread partitioning makes the summation order vary between runs, typically at 1e-16.

The assembler interface

begin(ndof, dim, num_stencils, stencil), a thread-safe add_local_hessian(local_hess, vertex_ids), end(). The old triplet code moved into TripletHessianAssembler verbatim, including its low-memory serial fallback, and hessian() became a thin wrapper over assemble_hessian. Folding lives in the driver, so every backend gets it without doing anything.

Collision-DOF Hessian totals before and after, to confirm the virtual call per collision does not register:

scene	baseline	after folding	through the interface
bunny	0.63 ms	0.49 ms	0.60 ms
cloth-ball	6.27 ms	6.69 ms	5.61 ms
armadillo-rollers	11.6 ms	12.1 ms	11.1 ms
n-body	10.2 ms	10.8 ms	10.2 ms
rod-twist	179 ms	176 ms	179 ms
puffer-ball	384 ms	656 ms	512 ms

All within run-to-run noise. Puffer-ball is again the outlier for thermal reasons rather than algorithmic ones.

The MeshFEMSparse backend

begin() builds the block sparsity pattern from the stencils via SystemAssembler::blockSparsityPattern; add_local_hessian() scatters each local Hessian into the value array using MeshFEM's sorted column-merge with per-column spin locks.

I did not take it entirely as-is. I compile MeshFEMSparse myself from a DOWNLOAD_ONLY fetch rather than through its CMake, which sidesteps its solver wrappers, its PUBLIC -fvisibility=hidden, and its own fetching of Eigen, TBB and MeshFEMCore. I also wrote the block-CSC to Eigen conversion instead of calling toEigen: it symmetrizes directly in parallel rather than going through an intermediate upper-triangle scalar matrix, which is about 2x faster even before caching, and it can reuse a cached structure across assemblies, which toEigen cannot.

The Eigen version mismatch was the real integration cost. MeshFEM targets Eigen 3.4 and we pin 5.0.1, and Eigen 5 removed internal::make_coherent. Because the call is qualified, name lookup happens at template definition time, so merely including MeshFEMCore's autodiff header fails to compile, and SparseMatrices.hh includes it at the root of everything. That is now fixed upstream in MeshFEM/MeshFEMCore#1 and the local shim is gone. I also apply EIGEN_DONT_VECTORIZE=1 to the MeshFEM target when SIMD is on, since compiling the same Eigen templates two ways is an ODR violation with real alignment consequences.

The other two upstream bugs took longer to find, and they are the same mistake in two places. CSCMatrix::expandSparsityPattern detects diagonal blocks by reading Ai[Ap[j+1] - 1], which is out of bounds for an empty block column; that showed up as intermittent SIGSEGVs on the large scenes until AddressSanitizer pinned it in one run. visitDiagonalScalarEntries does the analogous thing with diagBlockScalarLoc(), which for an empty column points into the preceding column's storage, so trace() returned 1951.93 where the dense trace is 447.82 and addDiag() wrote to the wrong entries. I only noticed the second one because apply() on the same matrix matched Eigen exactly, which ruled out the values and left the offsets. Both are fixed in MeshFEM/MeshFEMSparse#1.

Worth knowing if you use block_matrix(): the guard that should have caught this, missingRequiredDiagonalBlocks(), excludes the storage layout we use, so assertSupportsAssembly() never fires. I left that alone upstream rather than widen it, because assembly itself is genuinely fine (the assembler only touches columns a stencil references, and those always have their diagonal block) and widening it would make setZero() throw for legitimately sparse patterns.

Results with the pattern rebuilt on every call, median of 5, in full DOF. The first speedup column is what an Eigen-consuming caller sees; the second is what a caller that can take the block format directly would see.

scene	collisions	baseline	to Eigen	speedup	block-CSC	speedup
two-cubes	390	0.41 ms	0.14 ms	3.0x	0.13 ms	3.2x
bunny	570	1.20 ms	0.30 ms	4.0x	0.11 ms	10.9x
cloth-ball	3,012	13.0 ms	1.32 ms	9.9x	1.05 ms	12.4x
cloth-funnel	5,168	3.6 ms	1.49 ms	2.5x	1.02 ms	3.6x
armadillo-rollers	15,533	21.9 ms	2.76 ms	7.9x	2.11 ms	10.4x
n-body	15,969	23.1 ms	2.21 ms	10.4x	1.54 ms	15.0x
rod-twist	327,491	257 ms	30.5 ms	8.4x	27.3 ms	9.4x
puffer-ball	512,190	571 ms	51.9 ms	11.0x	42.8 ms	13.3x

Because the conversion is cheap, Eigen consumers already get most of the win without any API change. Puffer-ball drops from 571 ms to 52 ms, and 21 ms of what remains is the local derivative evaluation itself. The tests compare against the triplet assembler at 1e-13 relative on every scene, in both PSD projection modes and both DOF spaces. Three consecutive full benchmark runs have been clean since the empty-column fix, where before I could not get through one reliably.

Reusing the pattern

begin() runs MeshFEM's detectChangedEntries against the cached pattern. It reuses the pattern unless the contact set gained a vertex pair, or lost more blocks than stale_block_tolerance() allows, which defaults to zero. Stale blocks assemble to explicit zeros, so the values stay correct either way.

The Eigen conversion structure gets cached the same way, so get_matrix() becomes a parallel value refill. The interface does not change at all. Holding one assembler across assemble_hessian calls is enough to get all of this.

scene	baseline	cold	reused	speedup
two-cubes	0.41 ms	0.151 ms	0.132 ms	3.1x
bunny	1.18 ms	0.182 ms	0.113 ms	10.5x
cloth-ball	14.3 ms	1.44 ms	0.493 ms	28.9x
cloth-funnel	3.9 ms	1.61 ms	1.11 ms	3.5x
armadillo-rollers	21.3 ms	2.55 ms	1.88 ms	11.3x
n-body	24.6 ms	2.93 ms	1.48 ms	16.7x
rod-twist	294 ms	37.1 ms	37.1 ms	7.9x
puffer-ball	804 ms	57.7 ms	43.1 ms	18.7x

Steady-state Hessians now cost 0.1 to 43 ms where the triplet path cost 0.4 to 800 ms. On cloth-ball the steady state (0.49 ms) is about twice the local derivative cost (0.23 ms), so assembly has more or less stopped being a line item. Tests cover the reuse decisions: identical set reuses, a shrunken set reuses only within tolerance and still produces exact values, and any new vertex pair forces a rebuild regardless of tolerance, because scattering into a missing entry would be undefined behavior.

Notice rod-twist gets nothing from reuse. At that size detectChangedEntries costs about as much as rebuilding, which is what motivated set_assume_unchanged_stencils(true): a caller who knows the set is identical, say because it is reassembling with a different PSD projection or stiffness, can skip detection. A differing stencil count disproves the assumption and falls back automatically, debug builds verify the claim, and a wrong assertion with matching counts is documented as undefined behavior.

scene	reused	assumed	speedup (assumed)
two-cubes	0.164 ms	0.141 ms	3.5x
bunny	0.141 ms	0.095 ms	10.9x
cloth-ball	0.505 ms	0.482 ms	28.7x
cloth-funnel	1.29 ms	1.08 ms	3.3x
armadillo-rollers	1.56 ms	1.42 ms	13.7x
n-body	1.15 ms	1.00 ms	22.8x
rod-twist	27.8 ms	23.3 ms	12.1x
puffer-ball	36.8 ms	32.5 ms	23.6x

The reused numbers for the two biggest scenes are better in this run than the previous one, so some of the detection cost is just noise. The fast path removes it either way.

Gradients

Two strategies with opposite scaling, picked by problem shape. The gather path writes local gradients to a flat buffer, builds a vertex-to-slot adjacency with a parallel counting sort, and lets each vertex sum its own contributions, so cost tracks the number of contributions. The scatter path is the existing thread-local accumulators, whose cost tracks ndof. Gather wins when out_ndof exceeds the number of contributions, which is the measured crossover.

scene	path	before	after
two-cubes	scatter	25 µs	28 µs
bunny	gather	84 µs	92 µs
cloth-ball	gather	512-612 µs	381 µs
cloth-funnel	gather	682 µs	554 µs
armadillo-rollers	scatter	324 µs	309 µs
rod-twist	scatter	3.78 ms	3.53 ms
puffer-ball	scatter	12.1 ms	11.1 ms

I tried gather everywhere first. It regressed rod-twist by 1.6x, which has 1.3M contributions over 120k DOF, so the buffer and adjacency traffic outweighs what thread-local accumulators cost there. Hence the split rather than a wholesale replacement.

This also removes the earlier caveat about folding gradients: cloth-ball goes from 595 µs to 444 µs folded, puffer-ball from 13.7 ms to 11.3 ms. Summation order is still not fixed on either path, so results vary by rounding between runs. Sorting each gather bucket would make that path reproducible if we ever want it.

Gradients were never the headline. These are 1.2-1.6x in the favorable regime and parity elsewhere, on a path that costs microseconds to milliseconds against the Hessian's milliseconds to hundreds of milliseconds. The remaining tbb::combinable sites (friction force, Gauss-Newton diagonal, smooth contact) could adopt the same helper later.

Getting at the block matrix, and Python

get_matrix() converts to Eigen::SparseMatrix, which is what most callers want. If you are already using MeshFEM, that conversion is wasted work, so block_matrix() returns the assembled BlockCSCHessianBase directly and you can hand it to MeshFEM's block SpMV or its Cholesky factorizers. On bunny that is 0.11 ms instead of 0.30 ms; on puffer-ball 42.8 instead of 51.9.

MeshFEM::BlockCSCHessianBase is forward declared in our header, so including meshfem_hessian_assembler.hpp still does not drag in any MeshFEM header. Callers who want the block matrix include <MeshFEMSparse/BlockCSCHessian.hh> themselves.

One caveat, since the pattern only covers vertices in contact: reading operations are fine (trace() skips the missing diagonals, apply() and the sparsity pattern are unaffected), but addDiag() and setDiag() have no entry to write and throw. Insert the missing diagonal blocks first if you need them.

assemble_hessian and all three assembler classes are bound to Python, so a Python caller can hold an assembler across iterations and get pattern reuse instead of being stuck on the cold path inside hessian():

asm = ipctk.MeshFEMHessianAssembler()
for _ in range(n_iters):
    barrier.assemble_hessian(collisions, mesh, X, asm, in_full_dof=True)
    H = asm.get_matrix()   # pattern reused from the second call on

The classes are py::is_final(). Subclassing them in Python would only be useful with a trampoline, and a Python-defined assembler would take the GIL once per collision, which is half a million times per assembly on puffer-ball.

Making it the default

IPC_TOOLKIT_WITH_MESHFEM_SPARSE now defaults to ON, and hessian() assembles through the block backend when it is compiled in, using a take_matrix() that moves the matrix out instead of copying it. With the option off you get the triplet path. A ColMajor derivative layout disables the option automatically with a warning, since the block layout needs per-vertex blocks.

scene	triplet (old default)	block-CSC (new default)
cloth-ball	5.0-6.7 ms	1.5 ms
armadillo-rollers	11-18 ms	2.2 ms
rod-twist	165-212 ms	29.5 ms
puffer-ball	375-1020 ms	48.9 ms

That shifts the bottleneck. to_full_dof is now 75-84% of end to end for anyone not passing in_full_dof=true, so that flag is the next thing worth wiring up in PolyFEM. The assume-unchanged path also beats cold assembly now on the two biggest scenes, 21.6 and 30.1 ms against 31.9 and 49.3 ms.

The full ctest suite (286 tests) passes with the option on and with it off. CI's existing matrix picks up the new default on all three platforms, which is the first time the backend gets built with MSVC, so that job is worth watching. The off path has no dedicated CI job.

Docs: the three assembler classes are in the C++ API docs, IPC_TOOLKIT_WITH_MESHFEM_SPARSE is in Doxygen's PREDEFINED so the guarded class renders, and MeshFEMSparse is in the optional-dependencies table. The Python bindings cover in_full_dof, is_selection_dof_map, assemble_hessian, and the assembler classes. MeshFEMHessianAssembler is inside the option's #ifdef, so a wheel built with the option off would not export it; if that matters we should define the name unconditionally and raise on construction instead.

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[codecov]

Codecov Report

❌ Patch coverage is 97.30640% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 96.65%. Comparing base (57344ee) to head (5bcea54).

Files with missing lines	Patch %	Lines
src/ipc/utils/hessian_assembler.cpp	93.54%	4 Missing ⚠️
src/ipc/utils/meshfem_hessian_assembler.cpp	97.84%	3 Missing ⚠️
src/ipc/potentials/potential.cpp	98.87%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #246      +/-   ##
==========================================
+ Coverage   96.58%   96.65%   +0.07%     
==========================================
  Files         163      168       +5     
  Lines       16673    16910     +237     
  Branches      922      954      +32     
==========================================
+ Hits        16103    16344     +241     
+ Misses        570      566       -4     

Flag	Coverage Δ
unittests	96.65% <97.30%> (+0.07%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.