From 7e79e39fa40e3eee073a016d3321314a5d92fa4c Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Wed, 29 Jul 2026 23:21:56 -0700
Subject: [PATCH 01/15] Add contact assembly benchmarks (Phase 0 of
 block-assembly plan)

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>
---
 tests/src/tests/potential/CMakeLists.txt      |   3 +
 tests/src/tests/potential/assembly_scene.cpp  | 106 ++++++
 tests/src/tests/potential/assembly_scene.hpp  |  84 +++++
 .../tests/potential/benchmark_assembly.cpp    | 347 ++++++++++++++++++
 4 files changed, 540 insertions(+)
 create mode 100644 tests/src/tests/potential/assembly_scene.cpp
 create mode 100644 tests/src/tests/potential/assembly_scene.hpp
 create mode 100644 tests/src/tests/potential/benchmark_assembly.cpp

diff --git a/tests/src/tests/potential/CMakeLists.txt b/tests/src/tests/potential/CMakeLists.txt
index f2ff672f6..4a40f0792 100644
--- a/tests/src/tests/potential/CMakeLists.txt
+++ b/tests/src/tests/potential/CMakeLists.txt
@@ -7,8 +7,11 @@ set(SOURCES
   test_distance_vector_methods.cpp
 
   # Benchmarks
+  benchmark_assembly.cpp
 
   # Utilities
+  assembly_scene.cpp
+  assembly_scene.hpp
 )
 
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
diff --git a/tests/src/tests/potential/assembly_scene.cpp b/tests/src/tests/potential/assembly_scene.cpp
new file mode 100644
index 000000000..f965017a7
--- /dev/null
+++ b/tests/src/tests/potential/assembly_scene.cpp
@@ -0,0 +1,106 @@
+#include "assembly_scene.hpp"
+
+#include <tests/utils.hpp>
+
+#include <spdlog/fmt/fmt.h>
+
+namespace ipc::tests {
+
+AssemblyScene::AssemblyScene(
+    std::string label,
+    const CollisionMesh& mesh,
+    Eigen::MatrixXd vertices,
+    NormalCollisions collisions,
+    const BarrierPotential& potential)
+    : m_label(std::move(label))
+    , m_mesh(mesh)
+    , m_vertices(std::move(vertices))
+    , m_collisions(std::move(collisions))
+    , m_potential(potential)
+{
+}
+
+std::string AssemblyScene::stats() const
+{
+    return fmt::format(
+        "{}: {} collisions | {} collision vertices ({} DOF) | "
+        "{} full vertices ({} DOF) | dim={}",
+        label(), num_collisions(), num_vertices(), ndof(), full_num_vertices(),
+        full_ndof(), m_mesh.dim());
+}
+
+std::optional<AssemblyScene> build_assembly_scene(const AssemblySceneSpec& spec)
+{
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    if (!load_mesh(spec.mesh_name, vertices, edges, faces)) {
+        return std::nullopt;
+    }
+
+    // Pad the full mesh with interior vertices that no edge or face
+    // references, so that the collision mesh is a strict subset of the full
+    // mesh (see AssemblySceneSpec::interior_vertex_ratio).
+    const Eigen::Index num_surface_vertices = vertices.rows();
+    const Eigen::Index num_interior_vertices = Eigen::Index(
+        std::llround(spec.interior_vertex_ratio * num_surface_vertices));
+
+    Eigen::MatrixXd full_vertices(
+        num_surface_vertices + num_interior_vertices, vertices.cols());
+    full_vertices.topRows(num_surface_vertices) = vertices;
+    if (num_interior_vertices > 0) {
+        // Position is irrelevant (these are excluded from the collision mesh),
+        // but the centroid keeps any bounding-box computation well behaved.
+        full_vertices.bottomRows(num_interior_vertices).rowwise() =
+            vertices.colwise().mean();
+    }
+
+    const CollisionMesh mesh(
+        CollisionMesh::construct_is_on_surface(full_vertices.rows(), edges),
+        std::vector<bool>(full_vertices.rows(), false), full_vertices, edges,
+        faces);
+
+    Eigen::MatrixXd collision_vertices = mesh.vertices(full_vertices);
+
+    NormalCollisions collisions;
+    collisions.build(mesh, collision_vertices, spec.dhat);
+    if (collisions.empty()) {
+        return std::nullopt;
+    }
+
+    // A stiffness of 1 keeps the numbers interpretable; assembly cost is
+    // independent of its value.
+    const BarrierPotential potential(spec.dhat, /*stiffness=*/1.0);
+
+    std::optional<AssemblyScene> scene;
+    scene.emplace(
+        spec.label.empty() ? spec.mesh_name : spec.label, mesh,
+        std::move(collision_vertices), std::move(collisions), potential);
+    return scene;
+}
+
+const std::vector<AssemblySceneSpec>& assembly_scene_specs()
+{
+    // dhat values were chosen so each scene has a non-trivial collision set;
+    // see the "Assembly scene statistics" test, which prints the actual counts.
+    // WARNING: increasing `dhat` grows the collision set superlinearly. The
+    // values below are measured (see the "Assembly scene statistics" test) and
+    // span 390 -> 512k collisions. Do not raise them without checking the
+    // resulting candidate/collision count first (use the "[assembly-probe]"
+    // test): `dhat` an order of magnitude larger on `cloth_ball92.ply`
+    // exhausts memory on a 64 GB host.
+    //
+    // The simulation-frame scenes use dhat = 1e-3 * bbox diagonal.
+    static const std::vector<AssemblySceneSpec> specs = {
+        { "two-cubes-close.ply", 1e-1, 1.0, "two-cubes" },
+        { "bunny.ply", 1e-2, 1.0, "bunny" },
+        { "cloth_ball92.ply", 1e-3, 1.0, "cloth-ball" },
+        { "cloth-funnel/227.ply", 4.0e-3, 1.0, "cloth-funnel" },
+        { "armadillo-rollers/326.ply", 1e-3, 1.0, "armadillo-rollers" },
+        { "n-body-simulation/balls16_18.ply", 1.78e-2, 1.0, "n-body" },
+        { "rod-twist/3036.ply", 1.09e-3, 1.0, "rod-twist" },
+        { "puffer-ball/20.ply", 1.44e-4, 1.0, "puffer-ball" },
+    };
+    return specs;
+}
+
+} // namespace ipc::tests
diff --git a/tests/src/tests/potential/assembly_scene.hpp b/tests/src/tests/potential/assembly_scene.hpp
new file mode 100644
index 000000000..7d1016205
--- /dev/null
+++ b/tests/src/tests/potential/assembly_scene.hpp
@@ -0,0 +1,84 @@
+#pragma once
+
+#include <ipc/collision_mesh.hpp>
+#include <ipc/collisions/normal/normal_collisions.hpp>
+#include <ipc/potentials/barrier_potential.hpp>
+
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace ipc::tests {
+
+/// @brief Description of a contact scene used to benchmark assembly.
+struct AssemblySceneSpec {
+    /// @brief Name of the mesh file relative to the test data directory.
+    std::string mesh_name;
+    /// @brief Barrier activation distance.
+    double dhat;
+    /// @brief Number of synthetic interior vertices per collision vertex.
+    ///
+    /// The collision meshes in the test data are surfaces, so every vertex ends
+    /// up in the collision mesh and the selection matrix used by
+    /// `CollisionMesh::to_full_dof` is (a permutation of) the identity. Real
+    /// users embed the surface in a volumetric FE mesh whose interior nodes are
+    /// absent from the collision mesh. We emulate that by padding the full mesh
+    /// with vertices that no edge or face references, which
+    /// `construct_is_on_surface` then excludes. A ratio of 1.0 means the full
+    /// mesh has twice as many vertices as the collision mesh.
+    double interior_vertex_ratio = 1.0;
+    /// @brief Short label used in benchmark names.
+    std::string label;
+};
+
+/// @brief A pre-built contact scene: mesh, positions, collisions, potential.
+///
+/// Building the collision set is deliberately *not* part of what the assembly
+/// benchmarks measure, so it happens once here.
+class AssemblyScene {
+public:
+    /// @note `CollisionMesh` and `BarrierPotential` declare destructors, which
+    /// suppresses their implicit move constructors, so they are taken by
+    /// const reference rather than by value.
+    AssemblyScene(
+        std::string label,
+        const CollisionMesh& mesh,
+        Eigen::MatrixXd vertices,
+        NormalCollisions collisions,
+        const BarrierPotential& potential);
+
+    const std::string& label() const { return m_label; }
+    const CollisionMesh& mesh() const { return m_mesh; }
+    const Eigen::MatrixXd& vertices() const { return m_vertices; }
+    const NormalCollisions& collisions() const { return m_collisions; }
+    const BarrierPotential& potential() const { return m_potential; }
+
+    size_t num_collisions() const { return m_collisions.size(); }
+    size_t num_vertices() const { return m_mesh.num_vertices(); }
+    size_t full_num_vertices() const { return m_mesh.full_num_vertices(); }
+    size_t ndof() const { return m_mesh.ndof(); }
+    size_t full_ndof() const { return m_mesh.full_ndof(); }
+
+    /// @brief A human-readable one-line summary of the scene's size.
+    std::string stats() const;
+
+private:
+    std::string m_label;
+    CollisionMesh m_mesh;
+    Eigen::MatrixXd m_vertices;
+    NormalCollisions m_collisions;
+    BarrierPotential m_potential;
+};
+
+/// @brief Build a scene, or return nullopt if its mesh is unavailable.
+///
+/// Some test meshes are private, so callers must handle absence by skipping.
+std::optional<AssemblyScene>
+build_assembly_scene(const AssemblySceneSpec& spec);
+
+/// @brief The standard set of scenes used by the assembly benchmarks.
+///
+/// Ordered by increasing collision count so a truncated run is still useful.
+const std::vector<AssemblySceneSpec>& assembly_scene_specs();
+
+} // namespace ipc::tests
diff --git a/tests/src/tests/potential/benchmark_assembly.cpp b/tests/src/tests/potential/benchmark_assembly.cpp
new file mode 100644
index 000000000..dedbde099
--- /dev/null
+++ b/tests/src/tests/potential/benchmark_assembly.cpp
@@ -0,0 +1,347 @@
+// Baseline measurements for the cost of assembling contact gradients and
+// Hessians. These exist to quantify, separately:
+//
+//   1. per-collision (local) derivative evaluation,
+//   2. global assembly (triplets + setFromTriplets),
+//   3. the reduced-DOF map (`CollisionMesh::to_full_dof`).
+//
+// (2) is not measured directly: it is the difference between the full call and
+// (1), because the two are interleaved inside a single `tbb::parallel_for`. The
+// "Assembly cost breakdown" test computes that difference and prints a table.
+//
+// Run with:
+//   ./ipc_toolkit_tests "[assembly]" --benchmark-samples 20
+//
+// The scenes are shared with the correctness tests via `assembly_scene.hpp`.
+
+#include "assembly_scene.hpp"
+
+#include <tests/utils.hpp>
+
+#include <catch2/benchmark/catch_benchmark.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+#include <catch2/generators/catch_generators_range.hpp>
+
+#include <ipc/candidates/candidates.hpp>
+#include <ipc/utils/eigen_ext.hpp>
+
+#include <tbb/blocked_range.h>
+#include <tbb/parallel_reduce.h>
+
+#include <spdlog/fmt/fmt.h>
+
+#include <algorithm>
+#include <chrono>
+#include <cstdio>
+#include <cstdlib>
+
+using namespace ipc;
+
+namespace {
+
+/// @brief Number of assemblies per "solve" in the amortized benchmark.
+///
+/// A Newton solve performs one Hessian assembly per iteration while the
+/// collision set is unchanged, which is exactly the situation a persistent
+/// sparsity pattern exploits. Assembling repeatedly in a single benchmark
+/// sample keeps that opportunity visible in the baseline numbers.
+constexpr int ASSEMBLIES_PER_SOLVE = 10;
+
+/// @brief Evaluate every local Hessian without assembling anything.
+///
+/// Mirrors the per-collision work inside `Potential::hessian` exactly: the DOF
+/// gather plus the local Hessian. Returns a value derived from every result so
+/// the computation cannot be optimized away.
+double local_hessians_only(
+    const ipc::tests::AssemblyScene& scene,
+    const PSDProjectionMethod psd = PSDProjectionMethod::NONE)
+{
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXi& edges = scene.mesh().edges();
+    const Eigen::MatrixXi& faces = scene.mesh().faces();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    return tbb::parallel_reduce(
+        tbb::blocked_range<size_t>(size_t(0), collisions.size()), 0.0,
+        [&](const tbb::blocked_range<size_t>& r, double partial) {
+            for (size_t i = r.begin(); i < r.end(); i++) {
+                const MatrixMax12d local_hess = potential.hessian(
+                    collisions[i], collisions[i].dof(X, edges, faces), psd);
+                partial += local_hess(0, 0);
+            }
+            return partial;
+        },
+        std::plus<double>());
+}
+
+/// @brief Evaluate every local gradient without assembling anything.
+double local_gradients_only(const ipc::tests::AssemblyScene& scene)
+{
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXi& edges = scene.mesh().edges();
+    const Eigen::MatrixXi& faces = scene.mesh().faces();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    return tbb::parallel_reduce(
+        tbb::blocked_range<size_t>(size_t(0), collisions.size()), 0.0,
+        [&](const tbb::blocked_range<size_t>& r, double partial) {
+            for (size_t i = r.begin(); i < r.end(); i++) {
+                const VectorMax12d local_grad = potential.gradient(
+                    collisions[i], collisions[i].dof(X, edges, faces));
+                partial += local_grad(0);
+            }
+            return partial;
+        },
+        std::plus<double>());
+}
+
+/// @brief Median wall-clock time of `f` over `num_samples` runs, in seconds.
+///
+/// Used by the breakdown table. Catch2's `BENCHMARK` reports richer statistics
+/// but does not expose its measurements programmatically, so the table below
+/// does its own (much simpler) timing. Treat the `BENCHMARK` output as
+/// authoritative and the table as a summary of the same ratios.
+template <typename F> double median_seconds(F&& f, const int num_samples = 5)
+{
+    std::vector<double> samples;
+    samples.reserve(num_samples);
+    for (int i = 0; i < num_samples; i++) {
+        const auto start = std::chrono::steady_clock::now();
+        f();
+        const auto end = std::chrono::steady_clock::now();
+        samples.push_back(std::chrono::duration<double>(end - start).count());
+    }
+    std::sort(samples.begin(), samples.end());
+    return samples[samples.size() / 2];
+}
+
+} // namespace
+
+TEST_CASE("Assembly scene statistics", "[!benchmark][assembly]")
+{
+    fmt::print("\n=== Assembly benchmark scenes ===\n");
+    for (const auto& spec : ipc::tests::assembly_scene_specs()) {
+        const std::optional<ipc::tests::AssemblyScene> scene =
+            ipc::tests::build_assembly_scene(spec);
+        if (!scene.has_value()) {
+            fmt::print(
+                "  {} (dhat={}): UNAVAILABLE (missing mesh or no collisions)\n",
+                spec.mesh_name, spec.dhat);
+        } else {
+            fmt::print("  {}\n", scene->stats());
+        }
+        // Flush eagerly: stdout is fully buffered when redirected, and an
+        // oversized scene can exhaust memory before the buffer is drained.
+        std::fflush(stdout);
+    }
+    fmt::print("\n");
+}
+
+// Safely probe a candidate scene for inclusion in `assembly_scene_specs()`.
+//
+// Building a collision set with a too-large dhat can exhaust host memory, so
+// this test (a) sizes dhat relative to the mesh's bounding-box diagonal and
+// (b) counts broad-phase candidates first, refusing to build the collision set
+// if there are too many. Run one scene per process:
+//
+//   IPC_ASSEMBLY_PROBE_MESH=puffer-ball/20.ply \
+//   IPC_ASSEMBLY_PROBE_DHAT_REL=1e-3 \
+//   ./ipc_toolkit_tests "[assembly-probe]"
+TEST_CASE("Assembly scene probe", "[.][assembly-probe]")
+{
+    const char* mesh_name = std::getenv("IPC_ASSEMBLY_PROBE_MESH");
+    if (mesh_name == nullptr) {
+        SKIP("Set IPC_ASSEMBLY_PROBE_MESH to probe a scene.");
+    }
+    const char* dhat_rel_str = std::getenv("IPC_ASSEMBLY_PROBE_DHAT_REL");
+    const double dhat_rel =
+        (dhat_rel_str != nullptr) ? std::atof(dhat_rel_str) : 1e-3;
+
+    // Above this many broad-phase candidates, do not attempt to build the
+    // collision set: candidate/collision storage grows superlinearly with dhat
+    // and has exhausted memory on a 64 GB host before.
+    constexpr size_t MAX_SAFE_CANDIDATES = 10'000'000;
+
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(ipc::tests::load_mesh(mesh_name, vertices, edges, faces));
+
+    const double bbox_diag =
+        (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).norm();
+    const double dhat = dhat_rel * bbox_diag;
+
+    const CollisionMesh mesh =
+        CollisionMesh::build_from_full_mesh(vertices, edges, faces);
+    vertices = mesh.vertices(vertices);
+
+    fmt::print(
+        "probe {}: V={} E={} F={} bbox_diag={:g} dhat={:g} (rel={:g})\n",
+        mesh_name, vertices.rows(), edges.rows(), faces.rows(), bbox_diag, dhat,
+        dhat_rel);
+    std::fflush(stdout);
+
+    Candidates candidates;
+    candidates.build(mesh, vertices, vertices, /*inflation_radius=*/dhat / 2);
+    fmt::print("probe {}: candidates={}\n", mesh_name, candidates.size());
+    std::fflush(stdout);
+
+    if (candidates.size() > MAX_SAFE_CANDIDATES) {
+        fmt::print(
+            "probe {}: REFUSING to build collisions (> {} candidates)\n",
+            mesh_name, MAX_SAFE_CANDIDATES);
+        return;
+    }
+
+    NormalCollisions collisions;
+    collisions.build(candidates, mesh, vertices, dhat);
+    fmt::print("probe {}: collisions={}\n", mesh_name, collisions.size());
+    std::fflush(stdout);
+}
+
+TEST_CASE("Benchmark contact Hessian assembly", "[!benchmark][assembly]")
+{
+    const auto spec = GENERATE(from_range(ipc::tests::assembly_scene_specs()));
+
+    const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+        ipc::tests::build_assembly_scene(spec);
+    if (!maybe_scene.has_value()) {
+        SKIP(
+            fmt::format(
+                "Scene '{}' is unavailable (missing mesh or no collisions).",
+                spec.mesh_name));
+    }
+    const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+
+    const CollisionMesh& mesh = scene.mesh();
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    fmt::print("\n{}\n", scene.stats());
+
+    // Precomputed so `to_full_dof` can be timed in isolation.
+    const Eigen::SparseMatrix<double> hess =
+        potential.hessian(collisions, mesh, X);
+
+    BENCHMARK(fmt::format("{}: local hessians", scene.label()))
+    {
+        return local_hessians_only(scene);
+    };
+
+    BENCHMARK(fmt::format("{}: hessian (collision DOF)", scene.label()))
+    {
+        return potential.hessian(collisions, mesh, X);
+    };
+
+    BENCHMARK(fmt::format("{}: to_full_dof (hessian)", scene.label()))
+    {
+        return mesh.to_full_dof(hess);
+    };
+
+    BENCHMARK(fmt::format("{}: hessian + to_full_dof", scene.label()))
+    {
+        return mesh.to_full_dof(potential.hessian(collisions, mesh, X));
+    };
+
+    BENCHMARK(
+        fmt::format(
+            "{}: {}x (hessian + to_full_dof)", scene.label(),
+            ASSEMBLIES_PER_SOLVE))
+    {
+        double checksum = 0;
+        for (int i = 0; i < ASSEMBLIES_PER_SOLVE; i++) {
+            checksum += mesh.to_full_dof(potential.hessian(collisions, mesh, X))
+                            .coeff(0, 0);
+        }
+        return checksum;
+    };
+}
+
+TEST_CASE("Benchmark contact gradient assembly", "[!benchmark][assembly]")
+{
+    const auto spec = GENERATE(from_range(ipc::tests::assembly_scene_specs()));
+
+    const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+        ipc::tests::build_assembly_scene(spec);
+    if (!maybe_scene.has_value()) {
+        SKIP(
+            fmt::format(
+                "Scene '{}' is unavailable (missing mesh or no collisions).",
+                spec.mesh_name));
+    }
+    const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+
+    const CollisionMesh& mesh = scene.mesh();
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    const Eigen::VectorXd grad = potential.gradient(collisions, mesh, X);
+
+    BENCHMARK(fmt::format("{}: local gradients", scene.label()))
+    {
+        return local_gradients_only(scene);
+    };
+
+    BENCHMARK(fmt::format("{}: gradient (collision DOF)", scene.label()))
+    {
+        return potential.gradient(collisions, mesh, X);
+    };
+
+    BENCHMARK(fmt::format("{}: to_full_dof (gradient)", scene.label()))
+    {
+        return mesh.to_full_dof(grad);
+    };
+}
+
+TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
+{
+    // Emits the table that goes into the performance report. Every row is
+    // measured on the same process/thread configuration so the shares are
+    // directly comparable.
+    fmt::print("\n=== Hessian assembly cost breakdown ===\n");
+    fmt::print(
+        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>8} {:>8} {:>8}\n", "scene",
+        "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)", "local%",
+        "asm%", "full%");
+
+    for (const auto& spec : ipc::tests::assembly_scene_specs()) {
+        const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+            ipc::tests::build_assembly_scene(spec);
+        if (!maybe_scene.has_value()) {
+            continue;
+        }
+        const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+
+        const CollisionMesh& mesh = scene.mesh();
+        const NormalCollisions& collisions = scene.collisions();
+        const BarrierPotential& potential = scene.potential();
+        const Eigen::MatrixXd& X = scene.vertices();
+
+        const Eigen::SparseMatrix<double> hess =
+            potential.hessian(collisions, mesh, X);
+
+        const double t_local =
+            median_seconds([&] { (void)local_hessians_only(scene); });
+        const double t_total = median_seconds(
+            [&] { (void)potential.hessian(collisions, mesh, X); });
+        const double t_full =
+            median_seconds([&] { (void)mesh.to_full_dof(hess); });
+
+        // Assembly is what the full call does beyond the local derivatives.
+        const double t_asm = t_total - t_local;
+        const double t_end_to_end = t_total + t_full;
+
+        constexpr double MS = 1e3;
+        fmt::print(
+            "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} "
+            "{:>7.1f}% {:>7.1f}% {:>7.1f}%\n",
+            scene.label(), scene.num_collisions(), t_local * MS, t_total * MS,
+            t_asm * MS, t_full * MS, 100.0 * t_local / t_end_to_end,
+            100.0 * t_asm / t_end_to_end, 100.0 * t_full / t_end_to_end);
+    }
+    fmt::print("\n");
+}

From fb7edda9215008248aa9fe1ec681c913044f34e9 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Wed, 29 Jul 2026 23:40:45 -0700
Subject: [PATCH 02/15] Assemble potential derivatives directly in full-mesh
 DOF (Phase 1)

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>
---
 python/src/collision_mesh.cpp                 |   7 +
 python/src/potentials/potential.hpp           |  17 +-
 src/ipc/collision_mesh.cpp                    |   3 +-
 src/ipc/collision_mesh.hpp                    |  11 ++
 src/ipc/potentials/potential.cpp              |  69 ++++++--
 src/ipc/potentials/potential.hpp              |  12 +-
 tests/src/tests/potential/CMakeLists.txt      |   1 +
 .../tests/potential/benchmark_assembly.cpp    |  38 ++++-
 .../potential/test_full_dof_assembly.cpp      | 151 ++++++++++++++++++
 9 files changed, 275 insertions(+), 34 deletions(-)
 create mode 100644 tests/src/tests/potential/test_full_dof_assembly.cpp

diff --git a/python/src/collision_mesh.cpp b/python/src/collision_mesh.cpp
index 13908f89b..088e8b304 100644
--- a/python/src/collision_mesh.cpp
+++ b/python/src/collision_mesh.cpp
@@ -410,6 +410,13 @@ void define_collision_mesh(py::module_& m)
                 Matrix quantity on the full mesh with size equal to full_ndof() × full_ndof().
             )ipc_Qu8mg5v7",
             "X"_a)
+        .def_property_readonly(
+            "is_selection_dof_map", &CollisionMesh::is_selection_dof_map,
+            R"ipc_Qu8mg5v7(
+            Whether the full ↔ collision DOF map is a pure selection matrix.
+
+            This is the case unless a (non-empty) displacement map was provided at construction. When true, to_full_dof() is equivalent to scattering entries from collision DOF i to full DOF dim * to_full_vertex_id(i // dim) + i % dim, so derivatives can be assembled directly in full-mesh DOFs instead of applying to_full_dof() after the fact.
+            )ipc_Qu8mg5v7")
         .def_property_readonly(
             "vertex_vertex_adjacencies",
             &CollisionMesh::vertex_vertex_adjacencies,
diff --git a/python/src/potentials/potential.hpp b/python/src/potentials/potential.hpp
index c8e676210..33ad4854d 100644
--- a/python/src/potentials/potential.hpp
+++ b/python/src/potentials/potential.hpp
@@ -36,7 +36,7 @@ void define_potential_methods(PyClass& potential)
             "gradient",
             py::overload_cast<
                 const TCollisions&, const CollisionMesh&,
-                Eigen::ConstRef<Eigen::MatrixXd>>(
+                Eigen::ConstRef<Eigen::MatrixXd>, const bool>(
                 &Potential<TCollisions>::gradient, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the gradient of the potential.
@@ -45,17 +45,18 @@ void define_potential_methods(PyClass& potential)
                 collisions: The set of collisions.
                 mesh: The collision mesh.
                 X: Degrees of freedom of the collision mesh (e.g., vertices or velocities).
+                in_full_dof: If true, return the gradient in full-mesh DOF (equivalent to mesh.to_full_dof(gradient), but assembled directly in full DOF when possible, avoiding the extra map).
 
             Returns:
-                The gradient of the potential w.r.t. X. This will have a size of X.size.
+                The gradient of the potential w.r.t. X. This will have a size of X.size (or mesh.full_ndof if in_full_dof).
             )ipc_Qu8mg5v7",
-            "collisions"_a, "mesh"_a, "X"_a)
+            "collisions"_a, "mesh"_a, "X"_a, "in_full_dof"_a = false)
         .def(
             "hessian",
             py::overload_cast<
                 const TCollisions&, const CollisionMesh&,
-                Eigen::ConstRef<Eigen::MatrixXd>, const PSDProjectionMethod>(
-                &Potential<TCollisions>::hessian, py::const_),
+                Eigen::ConstRef<Eigen::MatrixXd>, const PSDProjectionMethod,
+                const bool>(&Potential<TCollisions>::hessian, py::const_),
             R"ipc_Qu8mg5v7(
             Compute the hessian of the potential.
 
@@ -64,12 +65,14 @@ void define_potential_methods(PyClass& potential)
                 mesh: The collision mesh.
                 X: Degrees of freedom of the collision mesh (e.g., vertices or velocities).
                 project_hessian_to_psd: Make sure the hessian is positive semi-definite.
+                in_full_dof: If true, return the Hessian in full-mesh DOF (equivalent to mesh.to_full_dof(hessian), but assembled directly in full DOF when possible, avoiding the two sparse-matrix products).
 
             Returns:
-                The Hessian of the potential w.r.t. X. This will have a size of X.size by X.size.
+                The Hessian of the potential w.r.t. X. This will have a size of X.size by X.size (or mesh.full_ndof square if in_full_dof).
             )ipc_Qu8mg5v7",
             "collisions"_a, "mesh"_a, "X"_a,
-            "project_hessian_to_psd"_a = PSDProjectionMethod::NONE)
+            "project_hessian_to_psd"_a = PSDProjectionMethod::NONE,
+            "in_full_dof"_a = false)
         .def(
             "__call__",
             py::overload_cast<const TCollision&, Eigen::ConstRef<VectorMax12d>>(
diff --git a/src/ipc/collision_mesh.cpp b/src/ipc/collision_mesh.cpp
index 3105b1de2..3b2e15be1 100644
--- a/src/ipc/collision_mesh.cpp
+++ b/src/ipc/collision_mesh.cpp
@@ -70,7 +70,8 @@ CollisionMesh::CollisionMesh(
     // Initializes m_select_vertices and m_select_dof
     init_selection_matrices(dim);
 
-    if (displacement_map.size() == 0) {
+    m_is_selection_dof_map = displacement_map.size() == 0;
+    if (m_is_selection_dof_map) {
         m_displacement_map = m_select_vertices;
         m_displacement_dof_map = m_select_dof;
     } else {
diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp
index e867edf2b..fdf8a8d15 100644
--- a/src/ipc/collision_mesh.hpp
+++ b/src/ipc/collision_mesh.hpp
@@ -209,6 +209,14 @@ class CollisionMesh {
     Eigen::SparseMatrix<double>
     to_full_dof(const Eigen::SparseMatrix<double>& X) const;
 
+    /// @brief Whether the full ↔ collision DOF map is a pure selection matrix.
+    /// This is the case unless a (non-empty) displacement map was provided at
+    /// construction. When true, to_full_dof() is equivalent to scattering
+    /// entries from collision DOF i to full DOF `dim * to_full_vertex_id(i /
+    /// dim) + i % dim`, so derivatives can be assembled directly in full-mesh
+    /// DOFs instead of applying to_full_dof() after the fact.
+    bool is_selection_dof_map() const { return m_is_selection_dof_map; }
+
     // -----------------------------------------------------------------------
 
     /// @brief Get the vertex-vertex adjacency matrix.
@@ -410,6 +418,9 @@ class CollisionMesh {
     /// @brief Mapping from full displacements DOF to collision displacements DOF
     /// @note this is premultiplied by m_select_dof
     Eigen::SparseMatrix<double> m_displacement_dof_map;
+    /// @brief Whether the user-provided displacement map is the identity
+    /// (i.e., m_displacement_dof_map is a pure selection matrix).
+    bool m_is_selection_dof_map = true;
 
     /// @brief Vertices adjacent to vertices
     std::vector<std::vector<index_t>> m_vertex_vertex_adjacencies;
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index d4c08feae..8a305ca96 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -60,18 +60,28 @@ template <class TCollisions>
 Eigen::VectorXd Potential<TCollisions>::gradient(
     const TCollisions& collisions,
     const CollisionMesh& mesh,
-    Eigen::ConstRef<Eigen::MatrixXd> X) const
+    Eigen::ConstRef<Eigen::MatrixXd> X,
+    const bool in_full_dof) const
 {
     assert(X.rows() == mesh.num_vertices());
     IPC_TOOLKIT_PROFILE_BLOCK("Potential<T>::gradient()");
 
+    // Assemble directly in full-mesh DOF when the DOF map is a pure selection
+    // (remapping stencil vertex IDs is then equivalent to to_full_dof());
+    // otherwise assemble in collision DOF and apply the map at the end.
+    const bool fold_to_full = in_full_dof && mesh.is_selection_dof_map();
+    const bool map_to_full = in_full_dof && !fold_to_full;
+
+    const int out_ndof = fold_to_full ? mesh.full_ndof() : X.size();
+
     if (collisions.empty()) {
-        return Eigen::VectorXd::Zero(X.size());
+        Eigen::VectorXd grad = Eigen::VectorXd::Zero(out_ndof);
+        return map_to_full ? mesh.to_full_dof(grad) : grad;
     }
 
     const int dim = X.cols();
 
-    tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(X.size()));
+    tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(out_ndof));
 
     {
         IPC_TOOLKIT_PROFILE_BLOCK("Compute Local Gradients");
@@ -81,16 +91,26 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
             const VectorMaxNd local_grad = this->gradient(
                 collision, collision.dof(X, mesh.edges(), mesh.faces()));
 
+            auto ids = collision.vertex_ids(mesh.edges(), mesh.faces());
+            if (fold_to_full) {
+                for (auto& id : ids) {
+                    if (id >= 0) {
+                        id = mesh.to_full_vertex_id(id);
+                    }
+                }
+            }
+
             local_gradient_to_global_gradient(
-                local_grad, collision.vertex_ids(mesh.edges(), mesh.faces()),
-                dim, grad.local());
+                local_grad, ids, dim, grad.local());
         });
     }
 
     {
         IPC_TOOLKIT_PROFILE_BLOCK("Combine Local Gradients");
-        return grad.combine([](const Eigen::VectorXd& a,
-                               const Eigen::VectorXd& b) { return a + b; });
+        Eigen::VectorXd combined_grad = grad.combine(
+            [](const Eigen::VectorXd& a,
+               const Eigen::VectorXd& b) -> Eigen::VectorXd { return a + b; });
+        return map_to_full ? mesh.to_full_dof(combined_grad) : combined_grad;
     }
 }
 
@@ -99,20 +119,31 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const TCollisions& collisions,
     const CollisionMesh& mesh,
     Eigen::ConstRef<Eigen::MatrixXd> X,
-    const PSDProjectionMethod project_hessian_to_psd) const
+    const PSDProjectionMethod project_hessian_to_psd,
+    const bool in_full_dof) const
 {
     assert(X.rows() == mesh.num_vertices());
     IPC_TOOLKIT_PROFILE_BLOCK("Potential<T>::hessian()");
 
+    // Assemble directly in full-mesh DOF when the DOF map is a pure selection
+    // (remapping stencil vertex IDs is then equivalent to to_full_dof());
+    // otherwise assemble in collision DOF and apply the map at the end.
+    const bool fold_to_full = in_full_dof && mesh.is_selection_dof_map();
+    const bool map_to_full = in_full_dof && !fold_to_full;
+
+    const int n_verts = fold_to_full ? mesh.full_num_vertices() // NOLINT
+                                     : mesh.num_vertices();
+    const int ndof = fold_to_full ? mesh.full_ndof() : X.size(); // NOLINT
+
     if (collisions.empty()) {
-        return Eigen::SparseMatrix<double>(X.size(), X.size());
+        const Eigen::SparseMatrix<double> hess(ndof, ndof);
+        return map_to_full ? mesh.to_full_dof(hess) : hess;
     }
 
     const Eigen::MatrixXi& edges = mesh.edges();
     const Eigen::MatrixXi& faces = mesh.faces();
 
     const int dim = X.cols();
-    const int ndof = X.size();
 
     constexpr int MAX_TRIPLETS_SIZE = 10'000'000;
     const int buffer_size = std::min(MAX_TRIPLETS_SIZE, ndof);
@@ -138,14 +169,22 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
             {
                 IPC_TOOLKIT_PROFILE_BLOCK(
                     "Map Local Hessian to Global Triplets");
+                auto ids = collision.vertex_ids(edges, faces);
+                if (fold_to_full) {
+                    for (auto& id : ids) {
+                        if (id >= 0) {
+                            id = mesh.to_full_vertex_id(id);
+                        }
+                    }
+                }
                 local_hessian_to_global_triplets(
-                    local_hess, collision.vertex_ids(edges, faces), dim,
-                    *(hess_triplets.cache), mesh.num_vertices());
+                    local_hess, ids, dim, *(hess_triplets.cache), n_verts);
             }
         });
     }
     if (storage.empty()) {
-        return Eigen::SparseMatrix<double>();
+        const Eigen::SparseMatrix<double> hess(ndof, ndof);
+        return map_to_full ? mesh.to_full_dof(hess) : hess;
     }
 
     // Assemble the stiffness matrix by concatenating the tuples in each local
@@ -183,7 +222,7 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
             hess += local_storage.cache->get_matrix(false); // will also prune
         }
         hess.makeCompressed();
-        return hess;
+        return map_to_full ? mesh.to_full_dof(hess) : hess;
     }
 
     // Allocate triplets
@@ -218,7 +257,7 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
         hess.setFromTriplets(triplets.begin(), triplets.end());
     }
 
-    return hess;
+    return map_to_full ? mesh.to_full_dof(hess) : hess;
 }
 
 template class Potential<NormalCollisions>;
diff --git a/src/ipc/potentials/potential.hpp b/src/ipc/potentials/potential.hpp
index 2c6c42785..cc68426b9 100644
--- a/src/ipc/potentials/potential.hpp
+++ b/src/ipc/potentials/potential.hpp
@@ -40,24 +40,28 @@ template <class TCollisions> class Potential {
     /// @param collisions The set of collisions.
     /// @param mesh The collision mesh.
     /// @param X Degrees of freedom of the collision mesh (e.g., vertices or velocities).
-    /// @returns The gradient of the potential w.r.t. X. This will have a size of X.size().
+    /// @param in_full_dof If true, return the gradient in full-mesh DOF (equivalent to `mesh.to_full_dof(gradient)`, but assembled directly in full DOF when possible, avoiding the extra map).
+    /// @returns The gradient of the potential w.r.t. X. This will have a size of X.size() (or `mesh.full_ndof()` if in_full_dof).
     Eigen::VectorXd gradient(
         const TCollisions& collisions,
         const CollisionMesh& mesh,
-        Eigen::ConstRef<Eigen::MatrixXd> X) const;
+        Eigen::ConstRef<Eigen::MatrixXd> X,
+        const bool in_full_dof = false) const;
 
     /// @brief Compute the hessian of the potential.
     /// @param collisions The set of collisions.
     /// @param mesh The collision mesh.
     /// @param X Degrees of freedom of the collision mesh (e.g., vertices or velocities).
     /// @param project_hessian_to_psd Make sure the hessian is positive semi-definite.
-    /// @returns The Hessian of the potential w.r.t. X. This will have a size of X.size() by X.size().
+    /// @param in_full_dof If true, return the Hessian in full-mesh DOF (equivalent to `mesh.to_full_dof(hessian)`, but assembled directly in full DOF when possible, avoiding the two sparse-matrix products).
+    /// @returns The Hessian of the potential w.r.t. X. This will have a size of X.size() by X.size() (or `mesh.full_ndof()` square if in_full_dof).
     Eigen::SparseMatrix<double> hessian(
         const TCollisions& collisions,
         const CollisionMesh& mesh,
         Eigen::ConstRef<Eigen::MatrixXd> X,
         const PSDProjectionMethod project_hessian_to_psd =
-            PSDProjectionMethod::NONE) const;
+            PSDProjectionMethod::NONE,
+        const bool in_full_dof = false) const;
 
     // -- Single collision methods ---------------------------------------------
 
diff --git a/tests/src/tests/potential/CMakeLists.txt b/tests/src/tests/potential/CMakeLists.txt
index 4a40f0792..cc2ecde69 100644
--- a/tests/src/tests/potential/CMakeLists.txt
+++ b/tests/src/tests/potential/CMakeLists.txt
@@ -5,6 +5,7 @@ set(SOURCES
   test_smooth_potential.cpp
   test_friction_potential.cpp
   test_distance_vector_methods.cpp
+  test_full_dof_assembly.cpp
 
   # Benchmarks
   benchmark_assembly.cpp
diff --git a/tests/src/tests/potential/benchmark_assembly.cpp b/tests/src/tests/potential/benchmark_assembly.cpp
index dedbde099..864fc92d0 100644
--- a/tests/src/tests/potential/benchmark_assembly.cpp
+++ b/tests/src/tests/potential/benchmark_assembly.cpp
@@ -246,6 +246,15 @@ TEST_CASE("Benchmark contact Hessian assembly", "[!benchmark][assembly]")
         return mesh.to_full_dof(potential.hessian(collisions, mesh, X));
     };
 
+    // Phase 1: assemble directly in full DOF (folds to_full_dof into the
+    // triplet remap).
+    BENCHMARK(fmt::format("{}: hessian (full DOF, folded)", scene.label()))
+    {
+        return potential.hessian(
+            collisions, mesh, X, PSDProjectionMethod::NONE,
+            /*in_full_dof=*/true);
+    };
+
     BENCHMARK(
         fmt::format(
             "{}: {}x (hessian + to_full_dof)", scene.label(),
@@ -295,6 +304,12 @@ TEST_CASE("Benchmark contact gradient assembly", "[!benchmark][assembly]")
     {
         return mesh.to_full_dof(grad);
     };
+
+    // Phase 1: assemble directly in full DOF.
+    BENCHMARK(fmt::format("{}: gradient (full DOF, folded)", scene.label()))
+    {
+        return potential.gradient(collisions, mesh, X, /*in_full_dof=*/true);
+    };
 }
 
 TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
@@ -304,9 +319,10 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
     // directly comparable.
     fmt::print("\n=== Hessian assembly cost breakdown ===\n");
     fmt::print(
-        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>8} {:>8} {:>8}\n", "scene",
-        "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)", "local%",
-        "asm%", "full%");
+        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>8} {:>8} {:>8} "
+        "{:>8}\n",
+        "scene", "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)",
+        "fold(ms)", "local%", "asm%", "full%", "speedup");
 
     for (const auto& spec : ipc::tests::assembly_scene_specs()) {
         const std::optional<ipc::tests::AssemblyScene> maybe_scene =
@@ -330,6 +346,12 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
             [&] { (void)potential.hessian(collisions, mesh, X); });
         const double t_full =
             median_seconds([&] { (void)mesh.to_full_dof(hess); });
+        // Phase 1: fold to_full_dof into assembly.
+        const double t_folded = median_seconds([&] {
+            (void)potential.hessian(
+                collisions, mesh, X, PSDProjectionMethod::NONE,
+                /*in_full_dof=*/true);
+        });
 
         // Assembly is what the full call does beyond the local derivatives.
         const double t_asm = t_total - t_local;
@@ -337,11 +359,13 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
 
         constexpr double MS = 1e3;
         fmt::print(
-            "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} "
-            "{:>7.1f}% {:>7.1f}% {:>7.1f}%\n",
+            "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} "
+            "{:>7.1f}% {:>7.1f}% {:>7.1f}% {:>7.2f}x\n",
             scene.label(), scene.num_collisions(), t_local * MS, t_total * MS,
-            t_asm * MS, t_full * MS, 100.0 * t_local / t_end_to_end,
-            100.0 * t_asm / t_end_to_end, 100.0 * t_full / t_end_to_end);
+            t_asm * MS, t_full * MS, t_folded * MS,
+            100.0 * t_local / t_end_to_end, 100.0 * t_asm / t_end_to_end,
+            100.0 * t_full / t_end_to_end, t_end_to_end / t_folded);
+        std::fflush(stdout);
     }
     fmt::print("\n");
 }
diff --git a/tests/src/tests/potential/test_full_dof_assembly.cpp b/tests/src/tests/potential/test_full_dof_assembly.cpp
new file mode 100644
index 000000000..e69a3b1bb
--- /dev/null
+++ b/tests/src/tests/potential/test_full_dof_assembly.cpp
@@ -0,0 +1,151 @@
+// Tests for assembling potential derivatives directly in full-mesh DOF
+// (`in_full_dof=true`), which must match applying `CollisionMesh::to_full_dof`
+// to the collision-DOF result.
+
+#include "assembly_scene.hpp"
+
+#include <tests/config.hpp>
+#include <tests/utils.hpp>
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+
+#include <ipc/potentials/barrier_potential.hpp>
+
+using namespace ipc;
+
+TEST_CASE(
+    "Full-DOF assembly matches to_full_dof", "[potential][assembly][full_dof]")
+{
+    // Scenes are padded with interior vertices, so the selection matrix is a
+    // genuine subset selection (full_ndof == 2 * ndof), not a permutation.
+#ifdef NDEBUG
+    const size_t scene_index = GENERATE(size_t(0), size_t(1), size_t(3));
+#else
+    const size_t scene_index = 0; // two-cubes; larger scenes are slow in debug
+#endif
+    const auto& spec = ipc::tests::assembly_scene_specs().at(scene_index);
+
+    const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+        ipc::tests::build_assembly_scene(spec);
+    if (!maybe_scene.has_value()) {
+        SKIP(fmt::format("Scene '{}' is unavailable.", spec.mesh_name));
+    }
+    const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+    CAPTURE(scene.label());
+
+    const CollisionMesh& mesh = scene.mesh();
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    REQUIRE(mesh.is_selection_dof_map());
+    REQUIRE(mesh.full_ndof() > mesh.ndof());
+
+    SECTION("gradient")
+    {
+        const Eigen::VectorXd grad_folded =
+            potential.gradient(collisions, mesh, X, /*in_full_dof=*/true);
+        const Eigen::VectorXd grad_mapped =
+            mesh.to_full_dof(potential.gradient(collisions, mesh, X));
+
+        REQUIRE(grad_folded.size() == mesh.full_ndof());
+        // tbb::combinable partitions work nondeterministically, so the two
+        // calls may sum in different orders; compare with a tight tolerance.
+        const double scale = std::max(1.0, grad_mapped.norm());
+        CHECK((grad_folded - grad_mapped).norm() <= 1e-13 * scale);
+    }
+
+    SECTION("hessian")
+    {
+        const PSDProjectionMethod psd =
+            GENERATE(PSDProjectionMethod::NONE, PSDProjectionMethod::CLAMP);
+        CAPTURE(psd);
+
+        const Eigen::SparseMatrix<double> hess_folded =
+            potential.hessian(collisions, mesh, X, psd, /*in_full_dof=*/true);
+        const Eigen::SparseMatrix<double> hess_mapped =
+            mesh.to_full_dof(potential.hessian(collisions, mesh, X, psd));
+
+        REQUIRE(hess_folded.rows() == mesh.full_ndof());
+        REQUIRE(hess_folded.cols() == mesh.full_ndof());
+
+        const double scale = std::max(1.0, hess_mapped.norm());
+        CHECK((hess_folded - hess_mapped).norm() <= 1e-13 * scale);
+    }
+}
+
+TEST_CASE(
+    "Full-DOF assembly with a non-selection displacement map",
+    "[potential][assembly][full_dof]")
+{
+    // With a user-provided displacement map, folding by index remap is not
+    // possible; in_full_dof must fall back to applying to_full_dof and still
+    // produce the mapped result.
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(tests::load_mesh("two-cubes-close.ply", vertices, edges, faces));
+
+    // A diagonal (non-identity) displacement map: full vars are half-scale.
+    Eigen::SparseMatrix<double> displacement_map(
+        vertices.rows(), vertices.rows());
+    displacement_map.setIdentity();
+    displacement_map *= 0.5;
+
+    const CollisionMesh mesh(vertices, edges, faces, displacement_map);
+    REQUIRE(!mesh.is_selection_dof_map());
+
+    const double dhat = 1e-1;
+    NormalCollisions collisions;
+    collisions.build(mesh, vertices, dhat);
+    REQUIRE(!collisions.empty());
+
+    const BarrierPotential potential(dhat, /*stiffness=*/1.0);
+
+    // Both calls run the same collision-DOF assembly (whose tbb reduction is
+    // order-nondeterministic) followed by to_full_dof, so compare with a
+    // tight tolerance rather than exactly.
+    const Eigen::VectorXd grad_full =
+        potential.gradient(collisions, mesh, vertices, /*in_full_dof=*/true);
+    const Eigen::VectorXd grad_mapped =
+        mesh.to_full_dof(potential.gradient(collisions, mesh, vertices));
+    CHECK(
+        (grad_full - grad_mapped).norm()
+        <= 1e-13 * std::max(1.0, grad_mapped.norm()));
+
+    const Eigen::SparseMatrix<double> hess_full = potential.hessian(
+        collisions, mesh, vertices, PSDProjectionMethod::NONE,
+        /*in_full_dof=*/true);
+    const Eigen::SparseMatrix<double> hess_mapped =
+        mesh.to_full_dof(potential.hessian(collisions, mesh, vertices));
+    CHECK(
+        (hess_full - hess_mapped).norm()
+        <= 1e-13 * std::max(1.0, hess_mapped.norm()));
+}
+
+TEST_CASE(
+    "Full-DOF assembly with no collisions", "[potential][assembly][full_dof]")
+{
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(tests::load_mesh("two-cubes-far.ply", vertices, edges, faces));
+
+    const CollisionMesh mesh =
+        CollisionMesh::build_from_full_mesh(vertices, edges, faces);
+    vertices = mesh.vertices(vertices);
+
+    const NormalCollisions collisions; // empty
+    const BarrierPotential potential(1e-3, /*stiffness=*/1.0);
+
+    const Eigen::VectorXd grad =
+        potential.gradient(collisions, mesh, vertices, /*in_full_dof=*/true);
+    CHECK(grad.size() == mesh.full_ndof());
+    CHECK(grad.norm() == 0.0);
+
+    const Eigen::SparseMatrix<double> hess = potential.hessian(
+        collisions, mesh, vertices, PSDProjectionMethod::NONE,
+        /*in_full_dof=*/true);
+    CHECK(hess.rows() == mesh.full_ndof());
+    CHECK(hess.cols() == mesh.full_ndof());
+    CHECK(hess.nonZeros() == 0);
+}

From 9dd46ce2dd64510198073d4b502124ec1204a7ed Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Wed, 29 Jul 2026 23:52:40 -0700
Subject: [PATCH 03/15] Introduce HessianAssembler seam for pluggable assembly
 backends (Phase 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>
---
 src/ipc/potentials/potential.cpp    | 162 +++++++---------------------
 src/ipc/potentials/potential.hpp    |  24 +++++
 src/ipc/utils/CMakeLists.txt        |   2 +
 src/ipc/utils/hessian_assembler.cpp | 140 ++++++++++++++++++++++++
 src/ipc/utils/hessian_assembler.hpp |  80 ++++++++++++++
 5 files changed, 287 insertions(+), 121 deletions(-)
 create mode 100644 src/ipc/utils/hessian_assembler.cpp
 create mode 100644 src/ipc/utils/hessian_assembler.hpp

diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index 8a305ca96..a2c3ba542 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -2,37 +2,18 @@
 
 #include <ipc/collisions/normal/normal_collisions.hpp>
 #include <ipc/collisions/tangential/tangential_collisions.hpp>
+#include <ipc/utils/hessian_assembler.hpp>
 #include <ipc/utils/local_to_global.hpp>
+#include <ipc/utils/logger.hpp>
 #include <ipc/utils/profiler.hpp>
 
 #include <tbb/blocked_range.h>
 #include <tbb/combinable.h>
-#include <tbb/enumerable_thread_specific.h>
 #include <tbb/parallel_for.h>
-#include <tbb/parallel_for_each.h>
 #include <tbb/parallel_reduce.h>
 
 namespace ipc {
 
-namespace {
-    void set_triplets(
-        const Eigen::SparseMatrix<double>& M,
-        std::vector<Eigen::Triplet<double>>& triplets,
-        const size_t start_index)
-    {
-        using InnerIterator = Eigen::SparseMatrix<double>::InnerIterator;
-        assert(start_index + M.nonZeros() <= triplets.size());
-        int count = 0;
-        for (int k = 0; k < M.outerSize(); ++k) {
-            for (InnerIterator it(M, k); it; ++it) {
-                assert(count < M.nonZeros());
-                triplets[start_index + count++] =
-                    Eigen::Triplet<double>(it.row(), it.col(), it.value());
-            }
-        }
-    }
-} // namespace
-
 template <class TCollisions>
 double Potential<TCollisions>::operator()(
     const TCollisions& collisions,
@@ -122,7 +103,6 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const PSDProjectionMethod project_hessian_to_psd,
     const bool in_full_dof) const
 {
-    assert(X.rows() == mesh.num_vertices());
     IPC_TOOLKIT_PROFILE_BLOCK("Potential<T>::hessian()");
 
     // Assemble directly in full-mesh DOF when the DOF map is a pure selection
@@ -131,31 +111,49 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const bool fold_to_full = in_full_dof && mesh.is_selection_dof_map();
     const bool map_to_full = in_full_dof && !fold_to_full;
 
-    const int n_verts = fold_to_full ? mesh.full_num_vertices() // NOLINT
-                                     : mesh.num_vertices();
-    const int ndof = fold_to_full ? mesh.full_ndof() : X.size(); // NOLINT
+    TripletHessianAssembler assembler;
+    assemble_hessian(
+        collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full);
+    const Eigen::SparseMatrix<double> hess = assembler.get_matrix();
 
-    if (collisions.empty()) {
-        const Eigen::SparseMatrix<double> hess(ndof, ndof);
-        return map_to_full ? mesh.to_full_dof(hess) : hess;
+    return map_to_full ? mesh.to_full_dof(hess) : hess;
+}
+
+template <class TCollisions>
+void Potential<TCollisions>::assemble_hessian(
+    const TCollisions& collisions,
+    const CollisionMesh& mesh,
+    Eigen::ConstRef<Eigen::MatrixXd> X,
+    HessianAssembler& assembler,
+    const PSDProjectionMethod project_hessian_to_psd,
+    const bool in_full_dof) const
+{
+    assert(X.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("Potential<T>::assemble_hessian()");
+
+    // The HessianAssembler interface is sized for the universal collision
+    // stencil (≤ 4 vertices ⇒ ≤ 12 DOF).
+    static_assert(TCollision::STENCIL_SIZE == 4);
+    static_assert(std::is_same_v<MatrixMaxNd, MatrixMax12d>);
+
+    if (in_full_dof && !mesh.is_selection_dof_map()) {
+        log_and_throw_error(
+            "assemble_hessian: in_full_dof requires the mesh's DOF map to be "
+            "a pure selection (see CollisionMesh::is_selection_dof_map); "
+            "assemble in collision DOF and apply to_full_dof instead.");
     }
 
     const Eigen::MatrixXi& edges = mesh.edges();
     const Eigen::MatrixXi& faces = mesh.faces();
 
     const int dim = X.cols();
+    const int ndof = in_full_dof ? mesh.full_ndof() : X.size(); // NOLINT
 
-    constexpr int MAX_TRIPLETS_SIZE = 10'000'000;
-    const int buffer_size = std::min(MAX_TRIPLETS_SIZE, ndof);
-
-    tbb::enumerable_thread_specific<LocalThreadMatStorage> storage(
-        LocalThreadMatStorage(buffer_size, ndof, ndof));
+    assembler.begin(ndof, dim, collisions.size());
 
     {
-        IPC_TOOLKIT_PROFILE_BLOCK("compute local hessians and triplets");
+        IPC_TOOLKIT_PROFILE_BLOCK("compute and assemble local hessians");
         tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
-            auto& hess_triplets = storage.local();
-
             const TCollision& collision = collisions[i];
 
             MatrixMaxNd local_hess;
@@ -166,98 +164,20 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
                     project_hessian_to_psd);
             }
 
-            {
-                IPC_TOOLKIT_PROFILE_BLOCK(
-                    "Map Local Hessian to Global Triplets");
-                auto ids = collision.vertex_ids(edges, faces);
-                if (fold_to_full) {
-                    for (auto& id : ids) {
-                        if (id >= 0) {
-                            id = mesh.to_full_vertex_id(id);
-                        }
+            auto ids = collision.vertex_ids(edges, faces);
+            if (in_full_dof) {
+                for (auto& id : ids) {
+                    if (id >= 0) {
+                        id = mesh.to_full_vertex_id(id);
                     }
                 }
-                local_hessian_to_global_triplets(
-                    local_hess, ids, dim, *(hess_triplets.cache), n_verts);
             }
-        });
-    }
-    if (storage.empty()) {
-        const Eigen::SparseMatrix<double> hess(ndof, ndof);
-        return map_to_full ? mesh.to_full_dof(hess) : hess;
-    }
-
-    // Assemble the stiffness matrix by concatenating the tuples in each local
-    // storage
-
-    {
-        IPC_TOOLKIT_PROFILE_BLOCK("Prune Local Storages");
-        tbb::parallel_for_each(
-            storage.begin(), storage.end(),
-            [](const auto& local_storage) { local_storage.cache->prune(); });
-    }
-
-    // Prepares for parallel concatenation
-    std::vector<size_t> offsets(storage.size());
-
-    size_t index = 0;
-    size_t triplet_count = 0;
-    for (auto& local_storage : storage) {
-        offsets[index++] = triplet_count;
-        triplet_count += local_storage.cache->triplet_count();
-    }
 
-    std::vector<Eigen::Triplet<double>> triplets;
-
-    Eigen::SparseMatrix<double> hess(ndof, ndof);
-    if (triplet_count >= triplets.max_size()) {
-        // Serial fallback version in case the vector of triplets cannot be
-        // allocated
-        logger().warn(
-            "Unable to allocate sufficient memory for triplets. "
-            "Falling back to serial assembly, which may impact performance. "
-            "Consider reducing the problem size or optimizing memory usage.");
-        // Serially merge local storages
-        for (LocalThreadMatStorage& local_storage : storage) {
-            hess += local_storage.cache->get_matrix(false); // will also prune
-        }
-        hess.makeCompressed();
-        return map_to_full ? mesh.to_full_dof(hess) : hess;
-    }
-
-    // Allocate triplets
-    {
-        IPC_TOOLKIT_PROFILE_BLOCK("Allocate Triplets");
-        triplets.resize(triplet_count);
-    }
-
-    // Parallel copy into triplets
-    {
-        IPC_TOOLKIT_PROFILE_BLOCK("Parallel Copy into Triplets");
-        tbb::parallel_for(size_t(0), storage.size(), [&](size_t i) {
-            const SparseMatrixCache& cache =
-                dynamic_cast<const SparseMatrixCache&>(
-                    *((storage.begin() + i)->cache));
-            size_t offset = offsets[i];
-
-            std::copy(
-                cache.entries().begin(), cache.entries().end(),
-                triplets.begin() + offset);
-            offset += cache.entries().size();
-
-            if (cache.mat().nonZeros() > 0) {
-                set_triplets(cache.mat(), triplets, offset);
-            }
+            assembler.add_local_hessian(local_hess, ids);
         });
     }
 
-    // Sort and assemble
-    {
-        IPC_TOOLKIT_PROFILE_BLOCK("Assemble Hessian from Triplets");
-        hess.setFromTriplets(triplets.begin(), triplets.end());
-    }
-
-    return map_to_full ? mesh.to_full_dof(hess) : hess;
+    assembler.end();
 }
 
 template class Potential<NormalCollisions>;
diff --git a/src/ipc/potentials/potential.hpp b/src/ipc/potentials/potential.hpp
index cc68426b9..617297292 100644
--- a/src/ipc/potentials/potential.hpp
+++ b/src/ipc/potentials/potential.hpp
@@ -7,6 +7,8 @@
 
 namespace ipc {
 
+class HessianAssembler; // forward declaration (see utils/hessian_assembler.hpp)
+
 /// @brief Base class for potentials.
 /// @tparam TCollisions The type of the collisions.
 template <class TCollisions> class Potential {
@@ -63,6 +65,28 @@ template <class TCollisions> class Potential {
             PSDProjectionMethod::NONE,
         const bool in_full_dof = false) const;
 
+    /// @brief Assemble the Hessian of the potential using a custom assembler.
+    ///
+    /// Evaluates the local Hessian of every collision (in parallel) and feeds
+    /// each to `assembler` (see HessianAssembler). This decouples the local
+    /// derivative evaluation from the global matrix construction; hessian()
+    /// is a thin wrapper around this using a TripletHessianAssembler.
+    ///
+    /// @param collisions The set of collisions.
+    /// @param mesh The collision mesh.
+    /// @param X Degrees of freedom of the collision mesh (e.g., vertices or velocities).
+    /// @param assembler The assembler that accumulates the local Hessians.
+    /// @param project_hessian_to_psd Make sure the hessian is positive semi-definite.
+    /// @param in_full_dof If true, stencil vertex IDs are remapped to full-mesh vertex IDs (requires `mesh.is_selection_dof_map()`; throws otherwise).
+    void assemble_hessian(
+        const TCollisions& collisions,
+        const CollisionMesh& mesh,
+        Eigen::ConstRef<Eigen::MatrixXd> X,
+        HessianAssembler& assembler,
+        const PSDProjectionMethod project_hessian_to_psd =
+            PSDProjectionMethod::NONE,
+        const bool in_full_dof = false) const;
+
     // -- Single collision methods ---------------------------------------------
 
     /// @brief Compute the potential for a single collision.
diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt
index f0086228e..9845512d2 100644
--- a/src/ipc/utils/CMakeLists.txt
+++ b/src/ipc/utils/CMakeLists.txt
@@ -3,6 +3,8 @@ set(SOURCES
   default_init_allocator.hpp
   eigen_ext.hpp
   eigen_ext.tpp
+  hessian_assembler.cpp
+  hessian_assembler.hpp
   local_to_global.hpp
   logger.cpp
   logger.hpp
diff --git a/src/ipc/utils/hessian_assembler.cpp b/src/ipc/utils/hessian_assembler.cpp
new file mode 100644
index 000000000..e6de9746f
--- /dev/null
+++ b/src/ipc/utils/hessian_assembler.cpp
@@ -0,0 +1,140 @@
+#include "hessian_assembler.hpp"
+
+#include <ipc/utils/logger.hpp>
+#include <ipc/utils/profiler.hpp>
+
+#include <tbb/parallel_for.h>
+#include <tbb/parallel_for_each.h>
+
+#include <algorithm>
+#include <vector>
+
+namespace ipc {
+
+namespace {
+    void set_triplets(
+        const Eigen::SparseMatrix<double>& M,
+        std::vector<Eigen::Triplet<double>>& triplets,
+        const size_t start_index)
+    {
+        using InnerIterator = Eigen::SparseMatrix<double>::InnerIterator;
+        assert(start_index + M.nonZeros() <= triplets.size());
+        int count = 0;
+        for (int k = 0; k < M.outerSize(); ++k) {
+            for (InnerIterator it(M, k); it; ++it) {
+                assert(count < M.nonZeros());
+                triplets[start_index + count++] =
+                    Eigen::Triplet<double>(it.row(), it.col(), it.value());
+            }
+        }
+    }
+} // namespace
+
+void TripletHessianAssembler::begin(
+    const int ndof, const int dim, const size_t num_stencils)
+{
+    m_ndof = ndof;
+    m_dim = dim;
+
+    constexpr int MAX_TRIPLETS_SIZE = 10'000'000;
+    const int buffer_size = std::min(MAX_TRIPLETS_SIZE, ndof);
+
+    m_storage = std::make_unique<
+        tbb::enumerable_thread_specific<LocalThreadMatStorage>>(
+        LocalThreadMatStorage(buffer_size, ndof, ndof));
+}
+
+void TripletHessianAssembler::add_local_hessian(
+    const MatrixMax12d& local_hess, const std::array<index_t, 4>& vertex_ids)
+{
+    assert(m_storage != nullptr);
+    IPC_TOOLKIT_PROFILE_BLOCK("Map Local Hessian to Global Triplets");
+    local_hessian_to_global_triplets(
+        local_hess, vertex_ids, m_dim, *(m_storage->local().cache),
+        m_ndof / m_dim);
+}
+
+Eigen::SparseMatrix<double> TripletHessianAssembler::get_matrix()
+{
+    assert(m_storage != nullptr);
+    tbb::enumerable_thread_specific<LocalThreadMatStorage>& storage =
+        *m_storage;
+
+    Eigen::SparseMatrix<double> hess(m_ndof, m_ndof);
+    if (storage.empty()) {
+        return hess;
+    }
+
+    // Assemble the stiffness matrix by concatenating the triplets in each
+    // local storage.
+
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("Prune Local Storages");
+        tbb::parallel_for_each(
+            storage.begin(), storage.end(),
+            [](const auto& local_storage) { local_storage.cache->prune(); });
+    }
+
+    // Prepares for parallel concatenation
+    std::vector<size_t> offsets(storage.size());
+
+    size_t index = 0;
+    size_t triplet_count = 0;
+    for (auto& local_storage : storage) {
+        offsets[index++] = triplet_count;
+        triplet_count += local_storage.cache->triplet_count();
+    }
+
+    std::vector<Eigen::Triplet<double>> triplets;
+
+    if (triplet_count >= triplets.max_size()) {
+        // Serial fallback version in case the vector of triplets cannot be
+        // allocated
+        logger().warn(
+            "Unable to allocate sufficient memory for triplets. "
+            "Falling back to serial assembly, which may impact performance. "
+            "Consider reducing the problem size or optimizing memory usage.");
+        // Serially merge local storages
+        for (LocalThreadMatStorage& local_storage : storage) {
+            hess += local_storage.cache->get_matrix(false); // will also prune
+        }
+        hess.makeCompressed();
+        return hess;
+    }
+
+    // Allocate triplets
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("Allocate Triplets");
+        triplets.resize(triplet_count);
+    }
+
+    // Parallel copy into triplets
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("Parallel Copy into Triplets");
+        tbb::parallel_for(size_t(0), storage.size(), [&](size_t i) {
+            const SparseMatrixCache& cache =
+                dynamic_cast<const SparseMatrixCache&>(
+                    *((storage.begin() + i)->cache));
+            size_t offset = offsets[i];
+
+            std::copy(
+                cache.entries().begin(), cache.entries().end(),
+                triplets.begin() + offset);
+            offset += cache.entries().size();
+
+            if (cache.mat().nonZeros() > 0) {
+                set_triplets(cache.mat(), triplets, offset);
+            }
+        });
+    }
+
+    // Sort and assemble
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("Assemble Hessian from Triplets");
+        hess.setFromTriplets(triplets.begin(), triplets.end());
+    }
+
+    return hess;
+}
+
+} // namespace ipc
diff --git a/src/ipc/utils/hessian_assembler.hpp b/src/ipc/utils/hessian_assembler.hpp
new file mode 100644
index 000000000..c2f010ea3
--- /dev/null
+++ b/src/ipc/utils/hessian_assembler.hpp
@@ -0,0 +1,80 @@
+#pragma once
+
+#include <ipc/config.hpp>
+#include <ipc/utils/eigen_ext.hpp>
+#include <ipc/utils/local_to_global.hpp>
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <tbb/enumerable_thread_specific.h>
+
+#include <array>
+
+namespace ipc {
+
+/// @brief Abstract sink for assembling local (per-collision) Hessians into a
+/// global matrix.
+///
+/// The driver (e.g., Potential::hessian) evaluates one local Hessian per
+/// collision stencil in parallel and hands each to add_local_hessian(). A
+/// concrete assembler decides how the global matrix is stored and built (e.g.,
+/// triplets + setFromTriplets, or a persistent block-sparse pattern), and
+/// exposes its own accessors for the result.
+///
+/// Call sequence: begin(), then any number of add_local_hessian() calls
+/// (possibly concurrent), then end().
+class HessianAssembler {
+public:
+    virtual ~HessianAssembler() = default;
+
+    /// @brief Prepare for assembly. Called once before any add_local_hessian.
+    /// @param ndof Number of global scalar DOF (rows == cols of the result).
+    /// @param dim Spatial dimension (rows/cols per vertex block).
+    /// @param num_stencils Number of local Hessians that will be added.
+    virtual void begin(int ndof, int dim, size_t num_stencils) = 0;
+
+    /// @brief Add one local (stencil) Hessian to the global matrix.
+    ///
+    /// Must be safe to call concurrently from multiple threads between
+    /// begin() and end().
+    ///
+    /// @param local_hess Local Hessian of size (n·dim)×(n·dim), where
+    ///     n = local_hess.rows() / dim is the number of stencil vertices.
+    /// @param vertex_ids Global vertex IDs of the stencil; the first n entries
+    ///     are valid (remaining entries may be negative placeholders).
+    virtual void add_local_hessian(
+        const MatrixMax12d& local_hess,
+        const std::array<index_t, 4>& vertex_ids) = 0;
+
+    /// @brief Finish assembly. Called once after all add_local_hessian calls.
+    virtual void end() = 0;
+};
+
+/// @brief The default HessianAssembler: thread-local triplet caches merged
+/// into an Eigen::SparseMatrix via setFromTriplets.
+///
+/// This reproduces the historical behavior of Potential::hessian exactly.
+class TripletHessianAssembler final : public HessianAssembler {
+public:
+    TripletHessianAssembler() = default;
+
+    void begin(int ndof, int dim, size_t num_stencils) override;
+
+    void add_local_hessian(
+        const MatrixMax12d& local_hess,
+        const std::array<index_t, 4>& vertex_ids) override;
+
+    void end() override { } // All work happens in get_matrix().
+
+    /// @brief Merge the thread-local caches and build the global matrix.
+    /// Call once, after end(); the internal caches are consumed.
+    Eigen::SparseMatrix<double> get_matrix();
+
+private:
+    int m_ndof = 0;
+    int m_dim = 0;
+    std::unique_ptr<tbb::enumerable_thread_specific<LocalThreadMatStorage>>
+        m_storage;
+};
+
+} // namespace ipc

From 69fe4369d91a7e6df2f1431f63632ae593862563 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Thu, 30 Jul 2026 01:04:35 -0700
Subject: [PATCH 04/15] Add MeshFEMSparse block-CSC assembly backend (Phase 3)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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>
---
 CMakeLists.txt                                |   8 +
 cmake/recipes/meshfem_sparse.cmake            | 105 ++++++++
 src/ipc/config.hpp.in                         |   1 +
 src/ipc/potentials/potential.cpp              |  26 +-
 src/ipc/utils/CMakeLists.txt                  |   3 +
 src/ipc/utils/hessian_assembler.cpp           |   5 +-
 src/ipc/utils/hessian_assembler.hpp           |  21 +-
 src/ipc/utils/meshfem_eigen_compat.hpp        |  43 +++
 src/ipc/utils/meshfem_hessian_assembler.cpp   | 250 ++++++++++++++++++
 src/ipc/utils/meshfem_hessian_assembler.hpp   |  58 ++++
 tests/src/tests/potential/CMakeLists.txt      |   1 +
 .../tests/potential/benchmark_assembly.cpp    |  67 ++++-
 .../tests/potential/test_meshfem_assembly.cpp |  93 +++++++
 13 files changed, 660 insertions(+), 21 deletions(-)
 create mode 100644 cmake/recipes/meshfem_sparse.cmake
 create mode 100644 src/ipc/utils/meshfem_eigen_compat.hpp
 create mode 100644 src/ipc/utils/meshfem_hessian_assembler.cpp
 create mode 100644 src/ipc/utils/meshfem_hessian_assembler.hpp
 create mode 100644 tests/src/tests/potential/test_meshfem_assembly.cpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6b9f6329a..94ced74c3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -101,6 +101,7 @@ option(IPC_TOOLKIT_WITH_FILIB                 "Use filib for interval arithmetic
 option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD method of IPC"    OFF)
 option(IPC_TOOLKIT_WITH_PROFILER              "Enable performance profiler"                   OFF)
 option(IPC_TOOLKIT_WITH_TRACY                 "Enable Tracy frame profiler"                   OFF)
+option(IPC_TOOLKIT_WITH_MESHFEM_SPARSE        "Use MeshFEMSparse for block-accelerated assembly" OFF)
 
 # Advanced options
 option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
@@ -263,6 +264,13 @@ if(IPC_TOOLKIT_WITH_FILIB)
   target_link_libraries(ipc_toolkit PUBLIC filib::filib)
 endif()
 
+# Block-accelerated Hessian assembly
+if(IPC_TOOLKIT_WITH_MESHFEM_SPARSE)
+  include(meshfem_sparse)
+  # PUBLIC: the MeshFEMHessianAssembler header includes MeshFEMSparse headers.
+  target_link_libraries(ipc_toolkit PUBLIC MeshFEM::Sparse)
+endif()
+
 if(IPC_TOOLKIT_WITH_PROFILER)
   # Add nlohmann/json for the profiler
   include(json)
diff --git a/cmake/recipes/meshfem_sparse.cmake b/cmake/recipes/meshfem_sparse.cmake
new file mode 100644
index 000000000..ecc9805a1
--- /dev/null
+++ b/cmake/recipes/meshfem_sparse.cmake
@@ -0,0 +1,105 @@
+# MeshFEMSparse (https://github.com/MeshFEM/MeshFEMSparse)
+# License: MIT
+#
+# Block-CSC sparse matrix data structures and fast Hessian assembly routines
+# from the MeshFEM project (Mohammadian et al., "MeshFEM: A Block-accelerated
+# Solver for Nonlinear Finite Elements", SIGGRAPH 2026).
+#
+# Downloaded with DOWNLOAD_ONLY and compiled into our own minimal static
+# library rather than through upstream's CMake. We only need the matrix data
+# structures and assembly routines; upstream's target additionally compiles
+# sparse direct solver wrappers (which clash with Eigen 5's BLAS
+# declarations), declares -fvisibility=hidden PUBLIC, and fetches its own
+# Eigen/TBB/MeshFEMCore.
+if(TARGET MeshFEMSparse)
+    return()
+endif()
+
+message(STATUS "Third-party: creating target 'MeshFEMSparse'")
+
+include(eigen)
+include(onetbb)
+find_package(Threads REQUIRED)
+
+include(CPM)
+CPMAddPackage(
+    NAME MeshFEMCore
+    URL "https://github.com/MeshFEM/MeshFEMCore/archive/24e81c425e85eee3ed79af000e82ef1a75bbe696.zip"
+    URL_HASH SHA256=7d505f57af2a4b2fbc01668d57c1296d646c6aa0aa923c9d70e814b55987aff7
+    DOWNLOAD_ONLY YES
+)
+CPMAddPackage(
+    NAME MeshFEMSparse
+    URL "https://github.com/MeshFEM/MeshFEMSparse/archive/efe1af87cf6f7d04359628552e84c666b8612f4c.zip"
+    URL_HASH SHA256=e23cb168e9e709662083116e73876f97c23edd60413108bd70bde5a5a1684dd5
+    DOWNLOAD_ONLY YES
+)
+
+add_library(MeshFEMSparse STATIC
+    # Matrix data structures and assembly (no Solvers/)
+    "${MeshFEMSparse_SOURCE_DIR}/src/lib/MeshFEMSparse/BlockCSCHessian.cc"
+    "${MeshFEMSparse_SOURCE_DIR}/src/lib/MeshFEMSparse/BorderedSparseHessian.cc"
+    # MeshFEMCore support code (parallelism arenas, benchmark stubs, types)
+    "${MeshFEMCore_SOURCE_DIR}/src/lib/MeshFEMCore/GlobalBenchmark.cc"
+    "${MeshFEMCore_SOURCE_DIR}/src/lib/MeshFEMCore/Parallelism.cc"
+    "${MeshFEMCore_SOURCE_DIR}/src/lib/MeshFEMCore/Types.cc"
+)
+
+add_library(MeshFEM::Sparse ALIAS MeshFEMSparse)
+
+target_include_directories(MeshFEMSparse SYSTEM PUBLIC
+    "${MeshFEMCore_SOURCE_DIR}/src/lib"
+    "${MeshFEMSparse_SOURCE_DIR}/src/lib"
+    "${CMAKE_CURRENT_BINARY_DIR}/meshfem/exports"
+)
+
+# MeshFEMCore's headers include the CMake-generated <MeshFEM_export.h>. We
+# build a static library, so the export macros are empty.
+file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/meshfem/exports/MeshFEM_export.h" [[
+#pragma once
+#define MESHFEM_EXPORT
+#define MESHFEM_NO_EXPORT
+#define MESHFEM_DEPRECATED
+#define MESHFEM_DEPRECATED_EXPORT
+#define MESHFEM_DEPRECATED_NO_EXPORT
+]])
+
+target_link_libraries(MeshFEMSparse PUBLIC
+    Eigen3::Eigen
+    TBB::tbb
+    Threads::Threads
+)
+
+target_compile_features(MeshFEMSparse PUBLIC cxx_std_17)
+
+# Matches upstream MeshFEMCore's PUBLIC definitions (Parallelism.hh compiles
+# its TBB code paths only when MESHFEM_WITH_TBB is defined).
+target_compile_definitions(MeshFEMSparse PUBLIC
+    MESHFEM_WITH_TBB
+    NOMINMAX
+    _ENABLE_EXTENDED_ALIGNED_STORAGE
+    _USE_MATH_DEFINES
+)
+
+# MeshFEM targets Eigen 3.4, but Eigen 5 removed
+# Eigen::internal::make_coherent, which MeshFEMCore's
+# AutomaticDifferentiation.hh references. Force-include a shim reimplementing
+# it (see src/ipc/utils/meshfem_eigen_compat.hpp); TUs outside this target
+# that include MeshFEM headers must include the shim first themselves.
+set(MESHFEM_EIGEN_COMPAT_HEADER
+    "${PROJECT_SOURCE_DIR}/src/ipc/utils/meshfem_eigen_compat.hpp")
+if(MSVC)
+    target_compile_options(MeshFEMSparse PRIVATE "/FI${MESHFEM_EIGEN_COMPAT_HEADER}")
+else()
+    target_compile_options(MeshFEMSparse PRIVATE "-include" "${MESHFEM_EIGEN_COMPAT_HEADER}")
+endif()
+
+# ipc_toolkit is compiled with EIGEN_DONT_VECTORIZE=1 when SIMD is enabled;
+# compiling the same Eigen templates with different vectorization settings is
+# an ODR violation with real alignment/layout consequences.
+if(IPC_TOOLKIT_WITH_SIMD)
+    target_compile_definitions(MeshFEMSparse PRIVATE EIGEN_DONT_VECTORIZE=1)
+endif()
+
+# Folder name for IDE
+set_target_properties(MeshFEMSparse PROPERTIES FOLDER "ThirdParty")
diff --git a/src/ipc/config.hpp.in b/src/ipc/config.hpp.in
index 3c5c84e75..2107213b3 100644
--- a/src/ipc/config.hpp.in
+++ b/src/ipc/config.hpp.in
@@ -21,6 +21,7 @@
 #cmakedefine IPC_TOOLKIT_WITH_FILIB
 #cmakedefine IPC_TOOLKIT_WITH_PROFILER
 #cmakedefine IPC_TOOLKIT_WITH_TRACY
+#cmakedefine IPC_TOOLKIT_WITH_MESHFEM_SPARSE
 // #define IPC_TOOLKIT_DEBUG_AUTODIFF
 
 namespace ipc {
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index a2c3ba542..99cd1b3f6 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -149,7 +149,20 @@ void Potential<TCollisions>::assemble_hessian(
     const int dim = X.cols();
     const int ndof = in_full_dof ? mesh.full_ndof() : X.size(); // NOLINT
 
-    assembler.begin(ndof, dim, collisions.size());
+    // Stencil vertex IDs as assembled (remapped to full-mesh IDs if folding).
+    const auto stencil_vertex_ids = [&](const size_t i) {
+        auto ids = collisions[i].vertex_ids(edges, faces);
+        if (in_full_dof) {
+            for (auto& id : ids) {
+                if (id >= 0) {
+                    id = mesh.to_full_vertex_id(id);
+                }
+            }
+        }
+        return ids;
+    };
+
+    assembler.begin(ndof, dim, collisions.size(), stencil_vertex_ids);
 
     {
         IPC_TOOLKIT_PROFILE_BLOCK("compute and assemble local hessians");
@@ -164,16 +177,7 @@ void Potential<TCollisions>::assemble_hessian(
                     project_hessian_to_psd);
             }
 
-            auto ids = collision.vertex_ids(edges, faces);
-            if (in_full_dof) {
-                for (auto& id : ids) {
-                    if (id >= 0) {
-                        id = mesh.to_full_vertex_id(id);
-                    }
-                }
-            }
-
-            assembler.add_local_hessian(local_hess, ids);
+            assembler.add_local_hessian(local_hess, stencil_vertex_ids(i));
         });
     }
 
diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt
index 9845512d2..932d736a9 100644
--- a/src/ipc/utils/CMakeLists.txt
+++ b/src/ipc/utils/CMakeLists.txt
@@ -10,6 +10,9 @@ set(SOURCES
   logger.hpp
   matrix_cache.cpp
   matrix_cache.hpp
+  meshfem_eigen_compat.hpp
+  meshfem_hessian_assembler.cpp
+  meshfem_hessian_assembler.hpp
   merge_thread_local.hpp
   profiler.cpp
   profiler.hpp
diff --git a/src/ipc/utils/hessian_assembler.cpp b/src/ipc/utils/hessian_assembler.cpp
index e6de9746f..f8e142f16 100644
--- a/src/ipc/utils/hessian_assembler.cpp
+++ b/src/ipc/utils/hessian_assembler.cpp
@@ -31,7 +31,10 @@ namespace {
 } // namespace
 
 void TripletHessianAssembler::begin(
-    const int ndof, const int dim, const size_t num_stencils)
+    const int ndof,
+    const int dim,
+    const size_t num_stencils,
+    const StencilGetter& /*stencil*/)
 {
     m_ndof = ndof;
     m_dim = dim;
diff --git a/src/ipc/utils/hessian_assembler.hpp b/src/ipc/utils/hessian_assembler.hpp
index c2f010ea3..831fda478 100644
--- a/src/ipc/utils/hessian_assembler.hpp
+++ b/src/ipc/utils/hessian_assembler.hpp
@@ -9,6 +9,7 @@
 #include <tbb/enumerable_thread_specific.h>
 
 #include <array>
+#include <functional>
 
 namespace ipc {
 
@@ -25,13 +26,26 @@ namespace ipc {
 /// (possibly concurrent), then end().
 class HessianAssembler {
 public:
+    /// @brief Callable returning the global vertex IDs of stencil i.
+    /// The IDs match those later passed to add_local_hessian for the same
+    /// stencil index (invalid trailing entries are negative). Must be safe to
+    /// call concurrently.
+    using StencilGetter = std::function<std::array<index_t, 4>(size_t)>;
+
     virtual ~HessianAssembler() = default;
 
     /// @brief Prepare for assembly. Called once before any add_local_hessian.
     /// @param ndof Number of global scalar DOF (rows == cols of the result).
     /// @param dim Spatial dimension (rows/cols per vertex block).
     /// @param num_stencils Number of local Hessians that will be added.
-    virtual void begin(int ndof, int dim, size_t num_stencils) = 0;
+    /// @param stencil Enumerates the stencils' vertex IDs; pattern-based
+    ///     assemblers use this to build their sparsity pattern up front. Only
+    ///     valid for the duration of the begin() call.
+    virtual void begin(
+        int ndof,
+        int dim,
+        size_t num_stencils,
+        const StencilGetter& stencil) = 0;
 
     /// @brief Add one local (stencil) Hessian to the global matrix.
     ///
@@ -58,7 +72,10 @@ class TripletHessianAssembler final : public HessianAssembler {
 public:
     TripletHessianAssembler() = default;
 
-    void begin(int ndof, int dim, size_t num_stencils) override;
+    /// @note `stencil` is unused: the triplet path needs no sparsity pattern.
+    void
+    begin(int ndof, int dim, size_t num_stencils, const StencilGetter& stencil)
+        override;
 
     void add_local_hessian(
         const MatrixMax12d& local_hess,
diff --git a/src/ipc/utils/meshfem_eigen_compat.hpp b/src/ipc/utils/meshfem_eigen_compat.hpp
new file mode 100644
index 000000000..13599d9a3
--- /dev/null
+++ b/src/ipc/utils/meshfem_eigen_compat.hpp
@@ -0,0 +1,43 @@
+#pragma once
+
+// Eigen ≥ 5 removed Eigen::internal::make_coherent (previously provided by
+// unsupported/Eigen/AutoDiff), which MeshFEMCore's AutomaticDifferentiation.hh
+// still references. This shim reimplements it with the Eigen 3.4 semantics:
+// if exactly one of the two derivative vectors is empty, resize it to match
+// the other and zero it.
+//
+// It must be included (or force-included via -include//FI) before any
+// MeshFEMCore/MeshFEMSparse header. It intentionally has no includes of its
+// own so it can be force-included into the MeshFEMSparse build (see
+// cmake/recipes/meshfem_sparse.cmake).
+//
+// TODO: remove once MeshFEM supports Eigen 5 upstream.
+
+// NOTE: EIGEN_WORLD_VERSION is 3 forever; Eigen 5 moved to semver with
+// EIGEN_MAJOR_VERSION=5. If no Eigen header has been seen yet (the
+// force-include case), assume Eigen 5 — this repository pins Eigen 5.
+#if !defined(EIGEN_MAJOR_VERSION) || EIGEN_MAJOR_VERSION >= 5
+
+namespace Eigen {
+namespace internal {
+
+    template <typename DerTypeA, typename DerTypeB>
+    inline void make_coherent(const DerTypeA& a, const DerTypeB& b)
+    {
+        // Eigen 3.4's implementation const-casts too (the derivatives are
+        // semantically mutable scratch space of the AutoDiffScalar pair).
+        DerTypeA& a_ref = const_cast<DerTypeA&>(a); // NOLINT
+        DerTypeB& b_ref = const_cast<DerTypeB&>(b); // NOLINT
+        if (a_ref.size() == 0 && b_ref.size() != 0) {
+            a_ref.resize(b_ref.size());
+            a_ref.setZero();
+        } else if (b_ref.size() == 0 && a_ref.size() != 0) {
+            b_ref.resize(a_ref.size());
+            b_ref.setZero();
+        }
+    }
+
+} // namespace internal
+} // namespace Eigen
+
+#endif
diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
new file mode 100644
index 000000000..1696bcb5c
--- /dev/null
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -0,0 +1,250 @@
+#include "meshfem_hessian_assembler.hpp"
+
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+
+#include <ipc/utils/logger.hpp>
+#include <ipc/utils/meshfem_eigen_compat.hpp> // must precede MeshFEM headers
+#include <ipc/utils/profiler.hpp>
+
+#include <MeshFEMSparse/SystemAssembler.hh>
+#include <tbb/parallel_for.h>
+
+#include <cassert>
+#include <cstddef>
+#include <numeric>
+#include <vector>
+
+namespace ipc {
+
+struct MeshFEMHessianAssembler::ImplBase {
+    virtual ~ImplBase() = default;
+    virtual void
+    add(const MatrixMax12d& local_hess,
+        const std::array<index_t, 4>& vertex_ids) = 0;
+    virtual Eigen::SparseMatrix<double> to_eigen() const = 0;
+};
+
+/// Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM's own
+/// IPC integration (IPCWrapper.cc): one dim-sized block variable per vertex.
+template <int dim>
+struct MeshFEMHessianAssembler::Impl final
+    : public MeshFEMHessianAssembler::ImplBase {
+    using VarStructure = MeshFEM::OptimizationVarStructure<dim>;
+    /// Local (block) variables of a collision stencil: 1–4 vertices.
+    using Stencil = MeshFEM::ElementBlockVarsWithSizeRange<1, 4>;
+
+    Impl(
+        const size_t num_block_vars,
+        const size_t num_stencils,
+        const StencilGetter& stencil)
+        : m_assembler(num_block_vars)
+        , m_vars(num_block_vars)
+    {
+        {
+            IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM block sparsity pattern");
+            m_H = m_assembler.blockSparsityPattern(
+                num_stencils,
+                [&stencil](const size_t i) { return to_stencil(stencil(i)); });
+        }
+        m_H->setZero(); // Allocate (and zero) the value array.
+        m_locks.init(num_block_vars);
+    }
+
+    void
+    add(const MatrixMax12d& local_hess,
+        const std::array<index_t, 4>& vertex_ids) override
+    {
+        assert(local_hess.rows() == local_hess.cols());
+        assert(local_hess.rows() % dim == 0);
+
+        const Stencil evars = to_stencil(vertex_ids);
+        assert(evars.size() == size_t(local_hess.rows()) / dim);
+
+        // Per-vertex-pair dim×dim block of the local Hessian; (a, b) are
+        // scalar offsets computed by the assembler (multiples of dim).
+        const auto local_block = [&local_hess](
+                                     const size_t a, const size_t b,
+                                     const size_t /*block_rows*/,
+                                     const size_t /*block_cols*/) {
+            return local_hess.template block<dim, dim>(a, b);
+        };
+
+        // Sorted column-merge scatter into the block-CSC value array, with
+        // per-block-column spin locks for thread safety.
+        MeshFEM::ElementHessianContribAssembler<
+            /*UseBlockMergeAlgorithm=*/true>::
+            template run</*InParallel=*/true>(
+                m_H->Ax.data(), *m_H, local_block, evars, m_vars, m_locks);
+    }
+
+    // Convert the block-CSC upper-triangle matrix to a full symmetric Eigen
+    // matrix.
+    //
+    // NOTE: We deliberately do not use BlockCSCHessian::toEigen here: its
+    // toScalar step assumes every block column is non-empty (true for FE
+    // Hessians, where every node belongs to an element, but false for contact
+    // Hessians, where most vertices are collision-free) and reads out of
+    // bounds otherwise. This implementation also skips toEigen's intermediate
+    // upper-triangle scalar matrix, symmetrizing directly instead.
+    //
+    // Value layout (ContiguousBlocks=true, StoreFullDiagonalBlocks=true, the
+    // library default): block entry ii occupies Ax[N²·ii, N²·(ii+1)),
+    // column-major within the block; diagonal blocks are full N×N with a
+    // zeroed strict lower triangle.
+    Eigen::SparseMatrix<double> to_eigen() const override
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM block CSC to Eigen");
+        static constexpr int N = dim;
+
+        const auto& Ap = m_H->Ap; // block column pointers (size n_blocks + 1)
+        const auto& Ai = m_H->Ai; // block row indices
+        const auto& Ax = m_H->Ax; // scalar values (N² per block)
+
+        const std::ptrdiff_t n_blocks = m_H->n;
+        const std::ptrdiff_t n = N * n_blocks; // scalar dimension
+
+        // --- Symmetrized *block* pattern -------------------------------
+        // Every stored block (bi, bj) (bi ≤ bj) appears at (bi, bj) and, if
+        // off-diagonal, mirrored at (bj, bi) in the full matrix.
+        struct BlockEntry {
+            std::ptrdiff_t row; ///< Block row in the full (symmetrized) matrix
+            std::ptrdiff_t src; ///< Index of the stored block (into Ai/Ax)
+            bool transposed;    ///< Whether the stored block is mirrored
+        };
+
+        std::vector<std::ptrdiff_t> sym_col_start(n_blocks + 1, 0);
+        for (std::ptrdiff_t bj = 0; bj < n_blocks; bj++) {
+            for (std::ptrdiff_t ii = Ap[bj]; ii < Ap[bj + 1]; ii++) {
+                sym_col_start[bj + 1]++;
+                if (Ai[ii] != bj) {
+                    sym_col_start[Ai[ii] + 1]++;
+                }
+            }
+        }
+        std::partial_sum(
+            sym_col_start.begin(), sym_col_start.end(), sym_col_start.begin());
+        const std::ptrdiff_t num_sym_blocks = sym_col_start[n_blocks];
+
+        // Sweeping bj in ascending order appends rows to every column in
+        // ascending order: column c receives its upper entries (rows ≤ c) at
+        // step bj = c and its mirrored entries (rows bj > c) at later steps.
+        std::vector<BlockEntry> sym_entries(num_sym_blocks);
+        {
+            std::vector<std::ptrdiff_t> cursor(
+                sym_col_start.begin(), sym_col_start.end() - 1);
+            for (std::ptrdiff_t bj = 0; bj < n_blocks; bj++) {
+                for (std::ptrdiff_t ii = Ap[bj]; ii < Ap[bj + 1]; ii++) {
+                    const std::ptrdiff_t bi = Ai[ii];
+                    sym_entries[cursor[bj]++] = { bi, ii, false };
+                    if (bi != bj) {
+                        sym_entries[cursor[bi]++] = { bj, ii, true };
+                    }
+                }
+            }
+        }
+
+        // --- Expand to a scalar CSC matrix ------------------------------
+        Eigen::SparseMatrix<double> M(n, n);
+        M.makeCompressed();
+        M.resizeNonZeros(N * N * num_sym_blocks);
+
+        int* outer = M.outerIndexPtr();
+        int* inner = M.innerIndexPtr();
+        double* values = M.valuePtr();
+
+        for (std::ptrdiff_t bj = 0; bj <= n_blocks; bj++) {
+            const std::ptrdiff_t col_nnz = bj < n_blocks
+                ? N * (sym_col_start[bj + 1] - sym_col_start[bj])
+                : 0;
+            for (int c = 0; c < ((bj < n_blocks) ? N : 1); c++) {
+                outer[N * bj + c] =
+                    int(N * N * sym_col_start[bj] + c * col_nnz);
+            }
+        }
+
+        tbb::parallel_for(
+            std::ptrdiff_t(0), n_blocks, [&](const std::ptrdiff_t bj) {
+                for (int c = 0; c < N; c++) {
+                    std::ptrdiff_t out = outer[N * bj + c];
+                    for (std::ptrdiff_t ei = sym_col_start[bj];
+                         ei < sym_col_start[bj + 1]; ei++) {
+                        const BlockEntry& e = sym_entries[ei];
+                        const double* block = Ax.data() + N * N * e.src;
+                        const bool diagonal = e.row == bj;
+                        for (int r = 0; r < N; r++) {
+                            inner[out] = int(N * e.row + r);
+                            if (e.transposed || (diagonal && r > c)) {
+                                values[out] = block[r * N + c]; // transposed
+                            } else {
+                                values[out] = block[c * N + r];
+                            }
+                            out++;
+                        }
+                    }
+                }
+            });
+
+        return M;
+    }
+
+private:
+    /// Compact the (possibly -1-padded) vertex ID array into a MeshFEM
+    /// stencil of block variables.
+    static Stencil to_stencil(const std::array<index_t, 4>& vertex_ids)
+    {
+        Stencil bvars(0);
+        size_t back = 0;
+        for (const index_t id : vertex_ids) {
+            if (id >= 0) {
+                bvars[back++] = id;
+            }
+        }
+        bvars.resize(back);
+        return bvars;
+    }
+
+    MeshFEM::SystemAssembler<dim> m_assembler;
+    VarStructure m_vars;
+    std::unique_ptr<MeshFEM::BlockCSCHessian<VarStructure>> m_H;
+    MeshFEM::VarLocks m_locks;
+};
+
+MeshFEMHessianAssembler::MeshFEMHessianAssembler() = default;
+MeshFEMHessianAssembler::~MeshFEMHessianAssembler() = default;
+
+void MeshFEMHessianAssembler::begin(
+    const int ndof,
+    const int dim,
+    const size_t num_stencils,
+    const StencilGetter& stencil)
+{
+    assert(ndof % dim == 0);
+    const size_t num_block_vars = size_t(ndof) / dim;
+    if (dim == 2) {
+        m_impl =
+            std::make_unique<Impl<2>>(num_block_vars, num_stencils, stencil);
+    } else if (dim == 3) {
+        m_impl =
+            std::make_unique<Impl<3>>(num_block_vars, num_stencils, stencil);
+    } else {
+        log_and_throw_error(
+            "MeshFEMHessianAssembler: unsupported dimension {}!", dim);
+    }
+}
+
+void MeshFEMHessianAssembler::add_local_hessian(
+    const MatrixMax12d& local_hess, const std::array<index_t, 4>& vertex_ids)
+{
+    assert(m_impl != nullptr);
+    m_impl->add(local_hess, vertex_ids);
+}
+
+Eigen::SparseMatrix<double> MeshFEMHessianAssembler::get_matrix() const
+{
+    assert(m_impl != nullptr);
+    return m_impl->to_eigen();
+}
+
+} // namespace ipc
+
+#endif // IPC_TOOLKIT_WITH_MESHFEM_SPARSE
diff --git a/src/ipc/utils/meshfem_hessian_assembler.hpp b/src/ipc/utils/meshfem_hessian_assembler.hpp
new file mode 100644
index 000000000..7839e3f0d
--- /dev/null
+++ b/src/ipc/utils/meshfem_hessian_assembler.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+
+#include <ipc/utils/hessian_assembler.hpp>
+
+#include <memory>
+
+namespace ipc {
+
+// The MeshFEM backend scatters per-vertex d×d blocks, which requires
+// derivatives to be laid out [x0, y0, z0, x1, ...].
+static_assert(
+    VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor,
+    "The MeshFEMSparse assembly backend requires "
+    "IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor.");
+
+/// @brief HessianAssembler backed by MeshFEMSparse's block-CSC data
+/// structures (Mohammadian et al. 2026).
+///
+/// begin() builds a block sparsity pattern (one dim×dim block per interacting
+/// vertex pair) from the stencils; add_local_hessian() scatters each local
+/// Hessian directly into the pattern's value array using a sorted
+/// column-merge with per-column locks — no triplets, no setFromTriplets.
+///
+/// The assembled matrix is symmetric and stored upper-triangle-only in block
+/// CSC format; get_matrix() converts to a full symmetric Eigen matrix. This
+/// conversion costs a full copy — direct (zero-copy) access to the block-CSC
+/// format is planned as part of the persistent-pattern API (Phase 4).
+class MeshFEMHessianAssembler final : public HessianAssembler {
+public:
+    MeshFEMHessianAssembler();
+    ~MeshFEMHessianAssembler() override;
+
+    void
+    begin(int ndof, int dim, size_t num_stencils, const StencilGetter& stencil)
+        override;
+
+    void add_local_hessian(
+        const MatrixMax12d& local_hess,
+        const std::array<index_t, 4>& vertex_ids) override;
+
+    void end() override { } // Values are scattered in place; nothing to merge.
+
+    /// @brief Convert the assembled matrix to a full symmetric Eigen matrix.
+    Eigen::SparseMatrix<double> get_matrix() const;
+
+private:
+    struct ImplBase;
+    template <int dim> struct Impl;
+    std::unique_ptr<ImplBase> m_impl;
+};
+
+} // namespace ipc
+
+#endif // IPC_TOOLKIT_WITH_MESHFEM_SPARSE
diff --git a/tests/src/tests/potential/CMakeLists.txt b/tests/src/tests/potential/CMakeLists.txt
index cc2ecde69..2f4b3589f 100644
--- a/tests/src/tests/potential/CMakeLists.txt
+++ b/tests/src/tests/potential/CMakeLists.txt
@@ -6,6 +6,7 @@ set(SOURCES
   test_friction_potential.cpp
   test_distance_vector_methods.cpp
   test_full_dof_assembly.cpp
+  test_meshfem_assembly.cpp
 
   # Benchmarks
   benchmark_assembly.cpp
diff --git a/tests/src/tests/potential/benchmark_assembly.cpp b/tests/src/tests/potential/benchmark_assembly.cpp
index 864fc92d0..ecd2e2916 100644
--- a/tests/src/tests/potential/benchmark_assembly.cpp
+++ b/tests/src/tests/potential/benchmark_assembly.cpp
@@ -25,6 +25,9 @@
 
 #include <ipc/candidates/candidates.hpp>
 #include <ipc/utils/eigen_ext.hpp>
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+#include <ipc/utils/meshfem_hessian_assembler.hpp>
+#endif
 
 #include <tbb/blocked_range.h>
 #include <tbb/parallel_reduce.h>
@@ -35,6 +38,7 @@
 #include <chrono>
 #include <cstdio>
 #include <cstdlib>
+#include <limits>
 
 using namespace ipc;
 
@@ -255,6 +259,32 @@ TEST_CASE("Benchmark contact Hessian assembly", "[!benchmark][assembly]")
             /*in_full_dof=*/true);
     };
 
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+    // Phase 3: MeshFEMSparse block-CSC backend (pattern built per call; the
+    // pattern-reuse win is Phase 4).
+    BENCHMARK(fmt::format("{}: hessian (MeshFEM, full DOF)", scene.label()))
+    {
+        MeshFEMHessianAssembler assembler;
+        potential.assemble_hessian(
+            collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+            /*in_full_dof=*/true);
+        return assembler.get_matrix();
+    };
+
+    // Without the Eigen conversion: what a caller that consumes the block-CSC
+    // format directly (e.g., a block-aware solver) would pay.
+    BENCHMARK(
+        fmt::format("{}: hessian (MeshFEM, no Eigen conv.)", scene.label()))
+    {
+        // Pattern + scattered values only. The scatter writes to heap
+        // storage through virtual dispatch, so it cannot be optimized away.
+        MeshFEMHessianAssembler assembler;
+        potential.assemble_hessian(
+            collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+            /*in_full_dof=*/true);
+    };
+#endif
+
     BENCHMARK(
         fmt::format(
             "{}: {}x (hessian + to_full_dof)", scene.label(),
@@ -319,10 +349,11 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
     // directly comparable.
     fmt::print("\n=== Hessian assembly cost breakdown ===\n");
     fmt::print(
-        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>8} {:>8} {:>8} "
-        "{:>8}\n",
+        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>8} "
+        "{:>8} {:>8} {:>8}\n",
         "scene", "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)",
-        "fold(ms)", "local%", "asm%", "full%", "speedup");
+        "fold(ms)", "mfem(ms)", "mfblk(ms)", "local%", "asm%", "full%",
+        "speedup");
 
     for (const auto& spec : ipc::tests::assembly_scene_specs()) {
         const std::optional<ipc::tests::AssemblyScene> maybe_scene =
@@ -353,6 +384,27 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
                 /*in_full_dof=*/true);
         });
 
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+        // Phase 3: MeshFEM block-CSC backend (pattern rebuilt per call).
+        const double t_meshfem = median_seconds([&] {
+            MeshFEMHessianAssembler assembler;
+            potential.assemble_hessian(
+                collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+                /*in_full_dof=*/true);
+            (void)assembler.get_matrix();
+        });
+        // Same, but skipping the block-CSC -> Eigen conversion.
+        const double t_meshfem_blk = median_seconds([&] {
+            MeshFEMHessianAssembler assembler;
+            potential.assemble_hessian(
+                collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+                /*in_full_dof=*/true);
+        });
+#else
+        const double t_meshfem = std::numeric_limits<double>::quiet_NaN();
+        const double t_meshfem_blk = std::numeric_limits<double>::quiet_NaN();
+#endif
+
         // Assembly is what the full call does beyond the local derivatives.
         const double t_asm = t_total - t_local;
         const double t_end_to_end = t_total + t_full;
@@ -360,11 +412,12 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
         constexpr double MS = 1e3;
         fmt::print(
             "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} "
-            "{:>7.1f}% {:>7.1f}% {:>7.1f}% {:>7.2f}x\n",
+            "{:>10.3f} {:>10.3f} {:>7.1f}% {:>7.1f}% {:>7.1f}% {:>7.2f}x\n",
             scene.label(), scene.num_collisions(), t_local * MS, t_total * MS,
-            t_asm * MS, t_full * MS, t_folded * MS,
-            100.0 * t_local / t_end_to_end, 100.0 * t_asm / t_end_to_end,
-            100.0 * t_full / t_end_to_end, t_end_to_end / t_folded);
+            t_asm * MS, t_full * MS, t_folded * MS, t_meshfem * MS,
+            t_meshfem_blk * MS, 100.0 * t_local / t_end_to_end,
+            100.0 * t_asm / t_end_to_end, 100.0 * t_full / t_end_to_end,
+            t_end_to_end / t_folded);
         std::fflush(stdout);
     }
     fmt::print("\n");
diff --git a/tests/src/tests/potential/test_meshfem_assembly.cpp b/tests/src/tests/potential/test_meshfem_assembly.cpp
new file mode 100644
index 000000000..4b0234d78
--- /dev/null
+++ b/tests/src/tests/potential/test_meshfem_assembly.cpp
@@ -0,0 +1,93 @@
+// Tests for the MeshFEMSparse-backed Hessian assembler: it must produce the
+// same global matrix as the default triplet assembler.
+//
+// Only compiled to something when IPC_TOOLKIT_WITH_MESHFEM_SPARSE is enabled.
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+
+#include "assembly_scene.hpp"
+
+#include <tests/config.hpp>
+#include <tests/utils.hpp>
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators.hpp>
+
+#include <ipc/utils/hessian_assembler.hpp>
+#include <ipc/utils/meshfem_hessian_assembler.hpp>
+
+using namespace ipc;
+
+TEST_CASE(
+    "MeshFEM assembly matches triplet assembly",
+    "[potential][assembly][meshfem]")
+{
+#ifdef NDEBUG
+    const size_t scene_index = GENERATE(size_t(0), size_t(1), size_t(3));
+#else
+    const size_t scene_index = 0; // two-cubes; larger scenes are slow in debug
+#endif
+    const auto& spec = ipc::tests::assembly_scene_specs().at(scene_index);
+
+    const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+        ipc::tests::build_assembly_scene(spec);
+    if (!maybe_scene.has_value()) {
+        SKIP(fmt::format("Scene '{}' is unavailable.", spec.mesh_name));
+    }
+    const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+    CAPTURE(scene.label());
+
+    const CollisionMesh& mesh = scene.mesh();
+    const NormalCollisions& collisions = scene.collisions();
+    const BarrierPotential& potential = scene.potential();
+    const Eigen::MatrixXd& X = scene.vertices();
+
+    const PSDProjectionMethod psd =
+        GENERATE(PSDProjectionMethod::NONE, PSDProjectionMethod::CLAMP);
+    const bool in_full_dof = GENERATE(false, true);
+    CAPTURE(psd, in_full_dof);
+
+    TripletHessianAssembler triplet_assembler;
+    potential.assemble_hessian(
+        collisions, mesh, X, triplet_assembler, psd, in_full_dof);
+    const Eigen::SparseMatrix<double> expected = triplet_assembler.get_matrix();
+
+    MeshFEMHessianAssembler meshfem_assembler;
+    potential.assemble_hessian(
+        collisions, mesh, X, meshfem_assembler, psd, in_full_dof);
+    const Eigen::SparseMatrix<double> actual = meshfem_assembler.get_matrix();
+
+    REQUIRE(actual.rows() == expected.rows());
+    REQUIRE(actual.cols() == expected.cols());
+
+    // Same additions in a different order; compare with a tight tolerance.
+    const double scale = std::max(1.0, expected.norm());
+    CHECK((actual - expected).norm() <= 1e-13 * scale);
+}
+
+TEST_CASE(
+    "MeshFEM assembly with no collisions", "[potential][assembly][meshfem]")
+{
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(tests::load_mesh("two-cubes-far.ply", vertices, edges, faces));
+
+    const CollisionMesh mesh =
+        CollisionMesh::build_from_full_mesh(vertices, edges, faces);
+    vertices = mesh.vertices(vertices);
+
+    const NormalCollisions collisions; // empty
+    const BarrierPotential potential(1e-3, /*stiffness=*/1.0);
+
+    MeshFEMHessianAssembler assembler;
+    potential.assemble_hessian(collisions, mesh, vertices, assembler);
+    const Eigen::SparseMatrix<double> hess = assembler.get_matrix();
+
+    CHECK(hess.rows() == mesh.ndof());
+    CHECK(hess.cols() == mesh.ndof());
+    CHECK(hess.norm() == 0.0);
+}
+
+#endif // IPC_TOOLKIT_WITH_MESHFEM_SPARSE

From 2d866a7522723f4e78d8b58a87eb103e82d77018 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Thu, 30 Jul 2026 09:44:53 -0700
Subject: [PATCH 05/15] Reuse the block sparsity pattern across assemblies
 (Phase 4)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

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>
---
 src/ipc/utils/meshfem_hessian_assembler.cpp   | 242 +++++++++++++-----
 src/ipc/utils/meshfem_hessian_assembler.hpp   |  62 ++++-
 .../tests/potential/benchmark_assembly.cpp    |  93 ++++++-
 .../tests/potential/test_meshfem_assembly.cpp |  90 +++++++
 4 files changed, 400 insertions(+), 87 deletions(-)

diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
index 1696bcb5c..2bb1bd64b 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.cpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -18,10 +18,23 @@ namespace ipc {
 
 struct MeshFEMHessianAssembler::ImplBase {
     virtual ~ImplBase() = default;
+
+    virtual int dimension() const = 0;
+    virtual size_t num_block_vars() const = 0;
+
+    /// Reuse the cached pattern if the stencils still fit (returns true) or
+    /// rebuild it (returns false). Zeroes the values either way.
+    virtual bool update_pattern(
+        size_t num_stencils,
+        const StencilGetter& stencil,
+        size_t stale_block_tolerance,
+        bool assume_unchanged) = 0;
+
     virtual void
     add(const MatrixMax12d& local_hess,
         const std::array<index_t, 4>& vertex_ids) = 0;
-    virtual Eigen::SparseMatrix<double> to_eigen() const = 0;
+
+    virtual const Eigen::SparseMatrix<double>& to_eigen() const = 0;
 };
 
 /// Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM's own
@@ -33,27 +46,70 @@ struct MeshFEMHessianAssembler::Impl final
     /// Local (block) variables of a collision stencil: 1–4 vertices.
     using Stencil = MeshFEM::ElementBlockVarsWithSizeRange<1, 4>;
 
-    Impl(
-        const size_t num_block_vars,
-        const size_t num_stencils,
-        const StencilGetter& stencil)
+    explicit Impl(const size_t num_block_vars)
         : m_assembler(num_block_vars)
         , m_vars(num_block_vars)
     {
+        m_locks.init(num_block_vars);
+    }
+
+    int dimension() const override { return dim; }
+    size_t num_block_vars() const override { return m_vars.numBlocks(); }
+
+    bool update_pattern(
+        const size_t num_stencils,
+        const StencilGetter& stencil,
+        const size_t stale_block_tolerance,
+        const bool assume_unchanged) override
+    {
+        const auto block_vars_of_stencil = [&stencil](const size_t i) {
+            return to_stencil(stencil(i));
+        };
+
+        if (m_H != nullptr) {
+            // Caller-asserted fast path: skip change detection, which costs
+            // as much as a pattern rebuild on large scenes. A differing
+            // stencil count disproves the assumption, so fall through to
+            // detection in that case.
+            if (assume_unchanged && num_stencils == m_num_stencils) {
+                // Verify the caller's claim where it is cheap to do so.
+                assert(
+                    m_assembler.detectChangedEntries(
+                        *m_H, num_stencils, block_vars_of_stencil)
+                    != MeshFEM::SystemAssembler<dim>::NEW_ENTRIES);
+                m_H->setZero();
+                return true;
+            }
+
+            IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM detect pattern change");
+            const size_t changed = m_assembler.detectChangedEntries(
+                *m_H, num_stencils, block_vars_of_stencil);
+            // NEW_ENTRIES (size_t max) must always trigger a rebuild:
+            // scattering into a pattern that is missing an entry is undefined.
+            if (changed != MeshFEM::SystemAssembler<dim>::NEW_ENTRIES
+                && changed <= stale_block_tolerance) {
+                m_num_stencils = num_stencils;
+                m_H->setZero();
+                return true;
+            }
+        }
+
         {
             IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM block sparsity pattern");
             m_H = m_assembler.blockSparsityPattern(
-                num_stencils,
-                [&stencil](const size_t i) { return to_stencil(stencil(i)); });
+                num_stencils, block_vars_of_stencil);
         }
+        m_num_stencils = num_stencils;
         m_H->setZero(); // Allocate (and zero) the value array.
-        m_locks.init(num_block_vars);
+        m_eigen_structure_valid = false;
+        return false;
     }
 
     void
     add(const MatrixMax12d& local_hess,
         const std::array<index_t, 4>& vertex_ids) override
     {
+        assert(m_H != nullptr);
         assert(local_hess.rows() == local_hess.cols());
         assert(local_hess.rows() % dim == 0);
 
@@ -78,7 +134,9 @@ struct MeshFEMHessianAssembler::Impl final
     }
 
     // Convert the block-CSC upper-triangle matrix to a full symmetric Eigen
-    // matrix.
+    // matrix. The structure (symmetrized pattern + index arrays) is cached
+    // and reused while the block pattern is unchanged; only the values are
+    // recomputed per call.
     //
     // NOTE: We deliberately do not use BlockCSCHessian::toEigen here: its
     // toScalar step assumes every block column is non-empty (true for FE
@@ -91,74 +149,103 @@ struct MeshFEMHessianAssembler::Impl final
     // library default): block entry ii occupies Ax[N²·ii, N²·(ii+1)),
     // column-major within the block; diagonal blocks are full N×N with a
     // zeroed strict lower triangle.
-    Eigen::SparseMatrix<double> to_eigen() const override
+    const Eigen::SparseMatrix<double>& to_eigen() const override
     {
         IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM block CSC to Eigen");
-        static constexpr int N = dim;
+        assert(m_H != nullptr);
 
+        if (!m_eigen_structure_valid) {
+            build_eigen_structure();
+            m_eigen_structure_valid = true;
+        }
+        fill_eigen_values();
+        return m_M;
+    }
+
+private:
+    static constexpr int N = dim;
+
+    /// A block of the full (symmetrized) matrix and the stored block backing
+    /// it.
+    struct BlockEntry {
+        std::ptrdiff_t row; ///< Block row in the full (symmetrized) matrix
+        std::ptrdiff_t src; ///< Index of the stored block (into Ai/Ax)
+        bool transposed;    ///< Whether the stored block is mirrored
+    };
+
+    /// Compact the (possibly -1-padded) vertex ID array into a MeshFEM
+    /// stencil of block variables.
+    static Stencil to_stencil(const std::array<index_t, 4>& vertex_ids)
+    {
+        Stencil bvars(0);
+        size_t back = 0;
+        for (const index_t id : vertex_ids) {
+            if (id >= 0) {
+                bvars[back++] = id;
+            }
+        }
+        bvars.resize(back);
+        return bvars;
+    }
+
+    void build_eigen_structure() const
+    {
         const auto& Ap = m_H->Ap; // block column pointers (size n_blocks + 1)
         const auto& Ai = m_H->Ai; // block row indices
-        const auto& Ax = m_H->Ax; // scalar values (N² per block)
 
         const std::ptrdiff_t n_blocks = m_H->n;
         const std::ptrdiff_t n = N * n_blocks; // scalar dimension
 
-        // --- Symmetrized *block* pattern -------------------------------
+        // --- Symmetrized *block* pattern --------------------------------
         // Every stored block (bi, bj) (bi ≤ bj) appears at (bi, bj) and, if
         // off-diagonal, mirrored at (bj, bi) in the full matrix.
-        struct BlockEntry {
-            std::ptrdiff_t row; ///< Block row in the full (symmetrized) matrix
-            std::ptrdiff_t src; ///< Index of the stored block (into Ai/Ax)
-            bool transposed;    ///< Whether the stored block is mirrored
-        };
-
-        std::vector<std::ptrdiff_t> sym_col_start(n_blocks + 1, 0);
+        m_sym_col_start.assign(n_blocks + 1, 0);
         for (std::ptrdiff_t bj = 0; bj < n_blocks; bj++) {
             for (std::ptrdiff_t ii = Ap[bj]; ii < Ap[bj + 1]; ii++) {
-                sym_col_start[bj + 1]++;
+                m_sym_col_start[bj + 1]++;
                 if (Ai[ii] != bj) {
-                    sym_col_start[Ai[ii] + 1]++;
+                    m_sym_col_start[Ai[ii] + 1]++;
                 }
             }
         }
         std::partial_sum(
-            sym_col_start.begin(), sym_col_start.end(), sym_col_start.begin());
-        const std::ptrdiff_t num_sym_blocks = sym_col_start[n_blocks];
+            m_sym_col_start.begin(), m_sym_col_start.end(),
+            m_sym_col_start.begin());
+        const std::ptrdiff_t num_sym_blocks = m_sym_col_start[n_blocks];
 
         // Sweeping bj in ascending order appends rows to every column in
         // ascending order: column c receives its upper entries (rows ≤ c) at
         // step bj = c and its mirrored entries (rows bj > c) at later steps.
-        std::vector<BlockEntry> sym_entries(num_sym_blocks);
+        m_sym_entries.resize(num_sym_blocks);
         {
             std::vector<std::ptrdiff_t> cursor(
-                sym_col_start.begin(), sym_col_start.end() - 1);
+                m_sym_col_start.begin(), m_sym_col_start.end() - 1);
             for (std::ptrdiff_t bj = 0; bj < n_blocks; bj++) {
                 for (std::ptrdiff_t ii = Ap[bj]; ii < Ap[bj + 1]; ii++) {
                     const std::ptrdiff_t bi = Ai[ii];
-                    sym_entries[cursor[bj]++] = { bi, ii, false };
+                    m_sym_entries[cursor[bj]++] = { bi, ii, false };
                     if (bi != bj) {
-                        sym_entries[cursor[bi]++] = { bj, ii, true };
+                        m_sym_entries[cursor[bi]++] = { bj, ii, true };
                     }
                 }
             }
         }
 
-        // --- Expand to a scalar CSC matrix ------------------------------
-        Eigen::SparseMatrix<double> M(n, n);
-        M.makeCompressed();
-        M.resizeNonZeros(N * N * num_sym_blocks);
+        // --- Scalar CSC index arrays ------------------------------------
+        m_M.resize(n, n);
+        m_M.makeCompressed();
+        m_M.resizeNonZeros(N * N * num_sym_blocks);
 
-        int* outer = M.outerIndexPtr();
-        int* inner = M.innerIndexPtr();
-        double* values = M.valuePtr();
+        int* outer = m_M.outerIndexPtr();
+        int* inner = m_M.innerIndexPtr();
 
         for (std::ptrdiff_t bj = 0; bj <= n_blocks; bj++) {
             const std::ptrdiff_t col_nnz = bj < n_blocks
-                ? N * (sym_col_start[bj + 1] - sym_col_start[bj])
+                ? N * (m_sym_col_start[bj + 1] - m_sym_col_start[bj])
                 : 0;
             for (int c = 0; c < ((bj < n_blocks) ? N : 1); c++) {
                 outer[N * bj + c] =
-                    int(N * N * sym_col_start[bj] + c * col_nnz);
+                    int(N * N * m_sym_col_start[bj] + c * col_nnz);
             }
         }
 
@@ -166,13 +253,34 @@ struct MeshFEMHessianAssembler::Impl final
             std::ptrdiff_t(0), n_blocks, [&](const std::ptrdiff_t bj) {
                 for (int c = 0; c < N; c++) {
                     std::ptrdiff_t out = outer[N * bj + c];
-                    for (std::ptrdiff_t ei = sym_col_start[bj];
-                         ei < sym_col_start[bj + 1]; ei++) {
-                        const BlockEntry& e = sym_entries[ei];
+                    for (std::ptrdiff_t ei = m_sym_col_start[bj];
+                         ei < m_sym_col_start[bj + 1]; ei++) {
+                        for (int r = 0; r < N; r++) {
+                            inner[out++] = int(N * m_sym_entries[ei].row + r);
+                        }
+                    }
+                }
+            });
+    }
+
+    void fill_eigen_values() const
+    {
+        const auto& Ax = m_H->Ax; // scalar values (N² per block)
+        const std::ptrdiff_t n_blocks = m_H->n;
+
+        const int* outer = m_M.outerIndexPtr();
+        double* values = m_M.valuePtr();
+
+        tbb::parallel_for(
+            std::ptrdiff_t(0), n_blocks, [&](const std::ptrdiff_t bj) {
+                for (int c = 0; c < N; c++) {
+                    std::ptrdiff_t out = outer[N * bj + c];
+                    for (std::ptrdiff_t ei = m_sym_col_start[bj];
+                         ei < m_sym_col_start[bj + 1]; ei++) {
+                        const BlockEntry& e = m_sym_entries[ei];
                         const double* block = Ax.data() + N * N * e.src;
                         const bool diagonal = e.row == bj;
                         for (int r = 0; r < N; r++) {
-                            inner[out] = int(N * e.row + r);
                             if (e.transposed || (diagonal && r > c)) {
                                 values[out] = block[r * N + c]; // transposed
                             } else {
@@ -183,30 +291,20 @@ struct MeshFEMHessianAssembler::Impl final
                     }
                 }
             });
-
-        return M;
-    }
-
-private:
-    /// Compact the (possibly -1-padded) vertex ID array into a MeshFEM
-    /// stencil of block variables.
-    static Stencil to_stencil(const std::array<index_t, 4>& vertex_ids)
-    {
-        Stencil bvars(0);
-        size_t back = 0;
-        for (const index_t id : vertex_ids) {
-            if (id >= 0) {
-                bvars[back++] = id;
-            }
-        }
-        bvars.resize(back);
-        return bvars;
     }
 
     MeshFEM::SystemAssembler<dim> m_assembler;
     VarStructure m_vars;
     std::unique_ptr<MeshFEM::BlockCSCHessian<VarStructure>> m_H;
     MeshFEM::VarLocks m_locks;
+    /// Stencil count of the cached pattern (assume-unchanged sanity check).
+    size_t m_num_stencils = 0;
+
+    // Cached Eigen conversion structure; valid while the pattern is reused.
+    mutable bool m_eigen_structure_valid = false;
+    mutable std::vector<std::ptrdiff_t> m_sym_col_start;
+    mutable std::vector<BlockEntry> m_sym_entries;
+    mutable Eigen::SparseMatrix<double> m_M;
 };
 
 MeshFEMHessianAssembler::MeshFEMHessianAssembler() = default;
@@ -220,16 +318,22 @@ void MeshFEMHessianAssembler::begin(
 {
     assert(ndof % dim == 0);
     const size_t num_block_vars = size_t(ndof) / dim;
-    if (dim == 2) {
-        m_impl =
-            std::make_unique<Impl<2>>(num_block_vars, num_stencils, stencil);
-    } else if (dim == 3) {
-        m_impl =
-            std::make_unique<Impl<3>>(num_block_vars, num_stencils, stencil);
-    } else {
-        log_and_throw_error(
-            "MeshFEMHessianAssembler: unsupported dimension {}!", dim);
+
+    if (m_impl == nullptr || m_impl->dimension() != dim
+        || m_impl->num_block_vars() != num_block_vars) {
+        if (dim == 2) {
+            m_impl = std::make_unique<Impl<2>>(num_block_vars);
+        } else if (dim == 3) {
+            m_impl = std::make_unique<Impl<3>>(num_block_vars);
+        } else {
+            log_and_throw_error(
+                "MeshFEMHessianAssembler: unsupported dimension {}!", dim);
+        }
     }
+
+    m_reused_pattern = m_impl->update_pattern(
+        num_stencils, stencil, m_stale_block_tolerance,
+        m_assume_unchanged_stencils);
 }
 
 void MeshFEMHessianAssembler::add_local_hessian(
@@ -239,7 +343,7 @@ void MeshFEMHessianAssembler::add_local_hessian(
     m_impl->add(local_hess, vertex_ids);
 }
 
-Eigen::SparseMatrix<double> MeshFEMHessianAssembler::get_matrix() const
+const Eigen::SparseMatrix<double>& MeshFEMHessianAssembler::get_matrix() const
 {
     assert(m_impl != nullptr);
     return m_impl->to_eigen();
diff --git a/src/ipc/utils/meshfem_hessian_assembler.hpp b/src/ipc/utils/meshfem_hessian_assembler.hpp
index 7839e3f0d..1695bd25b 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.hpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.hpp
@@ -25,10 +25,17 @@ static_assert(
 /// Hessian directly into the pattern's value array using a sorted
 /// column-merge with per-column locks — no triplets, no setFromTriplets.
 ///
+/// The assembler is designed to be **reused across assemblies** (e.g., one
+/// instance per Newton solve): begin() compares the stencils against the
+/// cached sparsity pattern and rebuilds it only if the contact set gained new
+/// entries or lost more than stale_block_tolerance() blocks; otherwise the
+/// pattern (and the cached Eigen structure) are reused and only the values
+/// are recomputed. Stale blocks left in a reused pattern assemble to explicit
+/// zeros, which do not affect the matrix's value.
+///
 /// The assembled matrix is symmetric and stored upper-triangle-only in block
-/// CSC format; get_matrix() converts to a full symmetric Eigen matrix. This
-/// conversion costs a full copy — direct (zero-copy) access to the block-CSC
-/// format is planned as part of the persistent-pattern API (Phase 4).
+/// CSC format; get_matrix() converts to a full symmetric Eigen matrix,
+/// reusing the cached structure when the pattern is unchanged.
 class MeshFEMHessianAssembler final : public HessianAssembler {
 public:
     MeshFEMHessianAssembler();
@@ -45,12 +52,59 @@ class MeshFEMHessianAssembler final : public HessianAssembler {
     void end() override { } // Values are scattered in place; nothing to merge.
 
     /// @brief Convert the assembled matrix to a full symmetric Eigen matrix.
-    Eigen::SparseMatrix<double> get_matrix() const;
+    ///
+    /// The reference stays valid (and its values current) until the next
+    /// begin() call; copy it to keep a snapshot.
+    const Eigen::SparseMatrix<double>& get_matrix() const;
+
+    /// @brief Number of vanished blocks tolerated before a pattern rebuild.
+    ///
+    /// When a stencil introduces a vertex pair absent from the cached pattern,
+    /// the pattern is always rebuilt. When blocks merely disappear (e.g., a
+    /// contact separates), the pattern is reused as long as at most this many
+    /// blocks vanished. Mirrors MeshFEM's sparsityPatternUpdateThreshold.
+    size_t stale_block_tolerance() const { return m_stale_block_tolerance; }
+
+    /// @brief Set the number of vanished blocks tolerated before a rebuild.
+    void set_stale_block_tolerance(const size_t tolerance)
+    {
+        m_stale_block_tolerance = tolerance;
+    }
+
+    /// @brief Whether the last begin() call reused the cached pattern.
+    bool reused_pattern() const { return m_reused_pattern; }
+
+    /// @brief Whether begin() may skip change detection entirely.
+    bool assume_unchanged_stencils() const
+    {
+        return m_assume_unchanged_stencils;
+    }
+
+    /// @brief Allow begin() to skip change detection entirely.
+    ///
+    /// When enabled, begin() reuses the cached pattern without comparing the
+    /// stencils against it, as long as a pattern exists and the stencil count
+    /// is unchanged (a differing count falls back to normal detection). Use
+    /// this when the caller knows the collision set is identical to the
+    /// previous assembly (e.g., reassembling with a different PSD projection
+    /// or stiffness): on large scenes, change detection costs as much as a
+    /// pattern rebuild.
+    ///
+    /// @warning If the stencils did change (with an equal count), assembly
+    /// reads out of bounds. Debug builds verify the assumption.
+    void set_assume_unchanged_stencils(const bool assume)
+    {
+        m_assume_unchanged_stencils = assume;
+    }
 
 private:
     struct ImplBase;
     template <int dim> struct Impl;
     std::unique_ptr<ImplBase> m_impl;
+
+    size_t m_stale_block_tolerance = 0;
+    bool m_assume_unchanged_stencils = false;
+    bool m_reused_pattern = false;
 };
 
 } // namespace ipc
diff --git a/tests/src/tests/potential/benchmark_assembly.cpp b/tests/src/tests/potential/benchmark_assembly.cpp
index ecd2e2916..0d065ea09 100644
--- a/tests/src/tests/potential/benchmark_assembly.cpp
+++ b/tests/src/tests/potential/benchmark_assembly.cpp
@@ -260,21 +260,20 @@ TEST_CASE("Benchmark contact Hessian assembly", "[!benchmark][assembly]")
     };
 
 #ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
-    // Phase 3: MeshFEMSparse block-CSC backend (pattern built per call; the
-    // pattern-reuse win is Phase 4).
-    BENCHMARK(fmt::format("{}: hessian (MeshFEM, full DOF)", scene.label()))
+    // Phase 3: MeshFEMSparse block-CSC backend, cold (pattern built per
+    // call, e.g., a fresh assembler every iteration).
+    BENCHMARK(fmt::format("{}: hessian (MeshFEM, cold)", scene.label()))
     {
         MeshFEMHessianAssembler assembler;
         potential.assemble_hessian(
             collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
             /*in_full_dof=*/true);
-        return assembler.get_matrix();
+        return &assembler.get_matrix();
     };
 
     // Without the Eigen conversion: what a caller that consumes the block-CSC
     // format directly (e.g., a block-aware solver) would pay.
-    BENCHMARK(
-        fmt::format("{}: hessian (MeshFEM, no Eigen conv.)", scene.label()))
+    BENCHMARK(fmt::format("{}: hessian (MeshFEM, cold, block)", scene.label()))
     {
         // Pattern + scattered values only. The scatter writes to heap
         // storage through virtual dispatch, so it cannot be optimized away.
@@ -283,6 +282,48 @@ TEST_CASE("Benchmark contact Hessian assembly", "[!benchmark][assembly]")
             collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
             /*in_full_dof=*/true);
     };
+
+    // Phase 4: persistent assembler — the pattern (and the cached Eigen
+    // structure) are reused across assemblies, as in a Newton solve.
+    {
+        MeshFEMHessianAssembler assembler;
+        potential.assemble_hessian(
+            collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+            /*in_full_dof=*/true); // warm-up: builds the pattern
+
+        BENCHMARK(fmt::format("{}: hessian (MeshFEM, reused)", scene.label()))
+        {
+            potential.assemble_hessian(
+                collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+                /*in_full_dof=*/true);
+            return &assembler.get_matrix();
+        };
+
+        BENCHMARK(
+            fmt::format(
+                "{}: {}x (MeshFEM, reused)", scene.label(),
+                ASSEMBLIES_PER_SOLVE))
+        {
+            double checksum = 0;
+            for (int i = 0; i < ASSEMBLIES_PER_SOLVE; i++) {
+                potential.assemble_hessian(
+                    collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+                    /*in_full_dof=*/true);
+                checksum += assembler.get_matrix().coeff(0, 0);
+            }
+            return checksum;
+        };
+
+        // Caller-asserted unchanged stencils: change detection skipped too.
+        assembler.set_assume_unchanged_stencils(true);
+        BENCHMARK(fmt::format("{}: hessian (MeshFEM, assumed)", scene.label()))
+        {
+            potential.assemble_hessian(
+                collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
+                /*in_full_dof=*/true);
+            return &assembler.get_matrix();
+        };
+    }
 #endif
 
     BENCHMARK(
@@ -349,11 +390,11 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
     // directly comparable.
     fmt::print("\n=== Hessian assembly cost breakdown ===\n");
     fmt::print(
-        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>8} "
-        "{:>8} {:>8} {:>8}\n",
+        "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} "
+        "{:>10} {:>10} {:>8} {:>8} {:>8} {:>8}\n",
         "scene", "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)",
-        "fold(ms)", "mfem(ms)", "mfblk(ms)", "local%", "asm%", "full%",
-        "speedup");
+        "fold(ms)", "mfem(ms)", "mfblk(ms)", "mfr(ms)", "mfa(ms)", "local%",
+        "asm%", "full%", "speedup");
 
     for (const auto& spec : ipc::tests::assembly_scene_specs()) {
         const std::optional<ipc::tests::AssemblyScene> maybe_scene =
@@ -400,9 +441,32 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
                 collisions, mesh, X, assembler, PSDProjectionMethod::NONE,
                 /*in_full_dof=*/true);
         });
+        // Phase 4: persistent assembler (pattern + Eigen structure reused).
+        MeshFEMHessianAssembler persistent_assembler;
+        potential.assemble_hessian(
+            collisions, mesh, X, persistent_assembler,
+            PSDProjectionMethod::NONE, /*in_full_dof=*/true); // warm-up
+        const double t_meshfem_reused = median_seconds([&] {
+            potential.assemble_hessian(
+                collisions, mesh, X, persistent_assembler,
+                PSDProjectionMethod::NONE, /*in_full_dof=*/true);
+            (void)persistent_assembler.get_matrix();
+        });
+        // Same, with caller-asserted unchanged stencils (no detection).
+        persistent_assembler.set_assume_unchanged_stencils(true);
+        const double t_meshfem_assumed = median_seconds([&] {
+            potential.assemble_hessian(
+                collisions, mesh, X, persistent_assembler,
+                PSDProjectionMethod::NONE, /*in_full_dof=*/true);
+            (void)persistent_assembler.get_matrix();
+        });
 #else
         const double t_meshfem = std::numeric_limits<double>::quiet_NaN();
         const double t_meshfem_blk = std::numeric_limits<double>::quiet_NaN();
+        const double t_meshfem_reused =
+            std::numeric_limits<double>::quiet_NaN();
+        const double t_meshfem_assumed =
+            std::numeric_limits<double>::quiet_NaN();
 #endif
 
         // Assembly is what the full call does beyond the local derivatives.
@@ -412,12 +476,13 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
         constexpr double MS = 1e3;
         fmt::print(
             "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} "
-            "{:>10.3f} {:>10.3f} {:>7.1f}% {:>7.1f}% {:>7.1f}% {:>7.2f}x\n",
+            "{:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} {:>7.1f}% {:>7.1f}% "
+            "{:>7.1f}% {:>7.2f}x\n",
             scene.label(), scene.num_collisions(), t_local * MS, t_total * MS,
             t_asm * MS, t_full * MS, t_folded * MS, t_meshfem * MS,
-            t_meshfem_blk * MS, 100.0 * t_local / t_end_to_end,
-            100.0 * t_asm / t_end_to_end, 100.0 * t_full / t_end_to_end,
-            t_end_to_end / t_folded);
+            t_meshfem_blk * MS, t_meshfem_reused * MS, t_meshfem_assumed * MS,
+            100.0 * t_local / t_end_to_end, 100.0 * t_asm / t_end_to_end,
+            100.0 * t_full / t_end_to_end, t_end_to_end / t_folded);
         std::fflush(stdout);
     }
     fmt::print("\n");
diff --git a/tests/src/tests/potential/test_meshfem_assembly.cpp b/tests/src/tests/potential/test_meshfem_assembly.cpp
index 4b0234d78..f6ff1ad19 100644
--- a/tests/src/tests/potential/test_meshfem_assembly.cpp
+++ b/tests/src/tests/potential/test_meshfem_assembly.cpp
@@ -18,6 +18,8 @@
 #include <ipc/utils/hessian_assembler.hpp>
 #include <ipc/utils/meshfem_hessian_assembler.hpp>
 
+#include <limits>
+
 using namespace ipc;
 
 TEST_CASE(
@@ -67,6 +69,94 @@ TEST_CASE(
     CHECK((actual - expected).norm() <= 1e-13 * scale);
 }
 
+TEST_CASE("MeshFEM assembly pattern reuse", "[potential][assembly][meshfem]")
+{
+    // A persistent assembler must produce correct results across repeated
+    // assemblies, reusing the cached sparsity pattern when the contact set
+    // allows it.
+#ifndef NDEBUG
+    SKIP("Building the bunny's collision sets is too slow in debug mode.");
+#endif
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(tests::load_mesh("bunny.ply", vertices, edges, faces));
+
+    const CollisionMesh mesh =
+        CollisionMesh::build_from_full_mesh(vertices, edges, faces);
+    vertices = mesh.vertices(vertices);
+
+    // Two nested contact sets: the small-dhat set is a subset of the large-
+    // dhat set (fewer vertex pairs within the activation distance).
+    const double dhat_large = 1e-2, dhat_small = 5e-3;
+    NormalCollisions collisions_large, collisions_small;
+    collisions_large.build(mesh, vertices, dhat_large);
+    collisions_small.build(mesh, vertices, dhat_small);
+    REQUIRE(!collisions_small.empty());
+    REQUIRE(collisions_small.size() < collisions_large.size());
+
+    const BarrierPotential potential(dhat_large, /*stiffness=*/1.0);
+
+    const auto reference = [&](const NormalCollisions& collisions) {
+        TripletHessianAssembler triplet_assembler;
+        potential.assemble_hessian(
+            collisions, mesh, vertices, triplet_assembler);
+        return triplet_assembler.get_matrix();
+    };
+    const auto check_matches = [](const Eigen::SparseMatrix<double>& actual,
+                                  const Eigen::SparseMatrix<double>& expected) {
+        // Stale blocks assemble to explicit zeros, so compare values (the
+        // difference ignores pattern mismatches), not nonZeros().
+        const double scale = std::max(1.0, expected.norm());
+        CHECK((actual - expected).norm() <= 1e-13 * scale);
+    };
+
+    MeshFEMHessianAssembler assembler;
+
+    // 1. First assembly builds the pattern.
+    potential.assemble_hessian(collisions_large, mesh, vertices, assembler);
+    CHECK(!assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_large));
+
+    // 2. Same collision set: the pattern must be reused.
+    potential.assemble_hessian(collisions_large, mesh, vertices, assembler);
+    CHECK(assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_large));
+
+    // 3. Shrunken collision set with a permissive tolerance: reused pattern
+    //    with stale (explicitly zero) blocks; values must still be correct.
+    assembler.set_stale_block_tolerance(std::numeric_limits<int>::max());
+    potential.assemble_hessian(collisions_small, mesh, vertices, assembler);
+    CHECK(assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_small));
+
+    // 4. Shrunken collision set with zero tolerance: rebuild.
+    assembler.set_stale_block_tolerance(0);
+    potential.assemble_hessian(collisions_small, mesh, vertices, assembler);
+    CHECK(!assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_small));
+
+    // 5. Grown collision set: new entries always force a rebuild, regardless
+    //    of the tolerance.
+    assembler.set_stale_block_tolerance(std::numeric_limits<int>::max());
+    potential.assemble_hessian(collisions_large, mesh, vertices, assembler);
+    CHECK(!assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_large));
+
+    // 6. Assume-unchanged fast path: identical set, detection skipped.
+    assembler.set_stale_block_tolerance(0);
+    assembler.set_assume_unchanged_stencils(true);
+    potential.assemble_hessian(collisions_large, mesh, vertices, assembler);
+    CHECK(assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_large));
+
+    // 7. Assume-unchanged with a differing stencil count: the assumption is
+    //    disproven, so it falls back to detection (and rebuilds here, since
+    //    the tolerance is zero and blocks disappeared).
+    potential.assemble_hessian(collisions_small, mesh, vertices, assembler);
+    CHECK(!assembler.reused_pattern());
+    check_matches(assembler.get_matrix(), reference(collisions_small));
+}
+
 TEST_CASE(
     "MeshFEM assembly with no collisions", "[potential][assembly][meshfem]")
 {

From 69d4174cb76a08686fcf68fae40959fb338d85a9 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Thu, 30 Jul 2026 23:02:09 -0700
Subject: [PATCH 06/15] Hybrid gather/scatter gradient assembly (Phase 5)

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>
---
 src/ipc/potentials/potential.cpp              | 124 +++++++++++++++++-
 .../potential/test_full_dof_assembly.cpp      |  10 +-
 2 files changed, 131 insertions(+), 3 deletions(-)

diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index 99cd1b3f6..4db2a69f9 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -12,8 +12,85 @@
 #include <tbb/parallel_for.h>
 #include <tbb/parallel_reduce.h>
 
+#include <atomic>
+#include <vector>
+
 namespace ipc {
 
+namespace {
+
+    /// @brief Gather-based global gradient assembly.
+    ///
+    /// Sums per-stencil gradient contributions into the global gradient by
+    /// gathering per *vertex* instead of scattering per stencil: a
+    /// vertex→slot adjacency is built with a parallel counting sort, then
+    /// each vertex sums its contributions independently. This avoids the
+    /// dense per-thread accumulators of a scatter+reduce (whose zero+combine
+    /// cost grows with ndof, not with the number of collisions).
+    ///
+    /// @param out_ndof Size of the output gradient.
+    /// @param dim Spatial dimension.
+    /// @param local_grads Flat buffer of slot contributions; slot s occupies
+    ///     [dim·s, dim·(s+1)). Slots with an invalid vertex are never read.
+    /// @param slot_vertex Global vertex of each slot (negative = unused).
+    Eigen::VectorXd gather_global_gradient(
+        const int out_ndof,
+        const int dim,
+        Eigen::ConstRef<Eigen::VectorXd> local_grads,
+        const std::vector<index_t>& slot_vertex)
+    {
+        const size_t num_slots = slot_vertex.size();
+        const size_t out_num_verts = size_t(out_ndof) / dim;
+        assert(local_grads.size() == Eigen::Index(dim * num_slots));
+
+        // --- Vertex → slot adjacency (parallel counting sort) -----------
+        std::vector<std::atomic<std::ptrdiff_t>> cursor(out_num_verts);
+        tbb::parallel_for(size_t(0), out_num_verts, [&](const size_t v) {
+            cursor[v].store(0, std::memory_order_relaxed);
+        });
+        tbb::parallel_for(size_t(0), num_slots, [&](const size_t s) {
+            if (slot_vertex[s] >= 0) {
+                cursor[slot_vertex[s]].fetch_add(1, std::memory_order_relaxed);
+            }
+        });
+
+        std::vector<std::ptrdiff_t> bucket_start(out_num_verts + 1);
+        bucket_start[0] = 0;
+        for (size_t v = 0; v < out_num_verts; v++) {
+            bucket_start[v + 1] =
+                bucket_start[v] + cursor[v].load(std::memory_order_relaxed);
+            cursor[v].store(bucket_start[v], std::memory_order_relaxed);
+        }
+
+        std::vector<std::ptrdiff_t> bucket_slots(bucket_start.back());
+        tbb::parallel_for(size_t(0), num_slots, [&](const size_t s) {
+            if (slot_vertex[s] >= 0) {
+                bucket_slots[cursor[slot_vertex[s]].fetch_add(
+                    1, std::memory_order_relaxed)] = s;
+            }
+        });
+
+        // --- Per-vertex gather -------------------------------------------
+        Eigen::VectorXd grad(out_ndof);
+        tbb::parallel_for(size_t(0), out_num_verts, [&](const size_t v) {
+            VectorMax3d acc = VectorMax3d::Zero(dim);
+            for (auto bi = bucket_start[v]; bi < bucket_start[v + 1]; bi++) {
+                acc += local_grads.segment(dim * bucket_slots[bi], dim);
+            }
+
+            if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+                grad.segment(dim * v, dim) = acc;
+            } else {
+                for (int d = 0; d < dim; d++) {
+                    grad[out_num_verts * d + v] = acc[d];
+                }
+            }
+        });
+        return grad;
+    }
+
+} // namespace
+
 template <class TCollisions>
 double Potential<TCollisions>::operator()(
     const TCollisions& collisions,
@@ -61,12 +138,57 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
     }
 
     const int dim = X.cols();
+    const size_t num_collisions = collisions.size();
+    constexpr int STENCIL_SIZE = TCollision::STENCIL_SIZE;
+    const size_t num_slots = STENCIL_SIZE * num_collisions;
+
+    // Two assembly strategies with complementary scaling:
+    //  - gather: per-vertex reduction through a slot adjacency. Cost scales
+    //    with the number of contributions.
+    //  - scatter+reduce: thread-local dense accumulators. Cost scales with
+    //    out_ndof (zeroing + combining a gradient-sized vector per thread).
+    // The crossover is at out_ndof ≈ num_slots on the benchmarked scenes.
+    if (size_t(out_ndof) > num_slots) {
+        // Evaluate the local gradients into a flat buffer of per-vertex
+        // slots: vertex k of collision i is slot s = STENCIL_SIZE·i + k,
+        // with its dim gradient entries at local_grads[dim·s].
+        Eigen::VectorXd local_grads(dim * num_slots);
+        std::vector<index_t> slot_vertex(num_slots);
+
+        {
+            IPC_TOOLKIT_PROFILE_BLOCK("Compute Local Gradients");
+            tbb::parallel_for(size_t(0), num_collisions, [&](size_t i) {
+                const TCollision& collision = collisions[i];
+
+                const VectorMaxNd local_grad = this->gradient(
+                    collision, collision.dof(X, mesh.edges(), mesh.faces()));
+                local_grads.segment(STENCIL_SIZE * dim * i, local_grad.size()) =
+                    local_grad;
+
+                const auto ids =
+                    collision.vertex_ids(mesh.edges(), mesh.faces());
+                const int n_verts = local_grad.size() / dim;
+                for (int k = 0; k < STENCIL_SIZE; k++) {
+                    index_t id = k < n_verts ? ids[k] : index_t(-1);
+                    if (fold_to_full && id >= 0) {
+                        id = mesh.to_full_vertex_id(id);
+                    }
+                    slot_vertex[STENCIL_SIZE * i + k] = id;
+                }
+            });
+        }
+
+        IPC_TOOLKIT_PROFILE_BLOCK("Gather Local Gradients");
+        Eigen::VectorXd grad =
+            gather_global_gradient(out_ndof, dim, local_grads, slot_vertex);
+        return map_to_full ? mesh.to_full_dof(grad) : grad;
+    }
 
     tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(out_ndof));
 
     {
         IPC_TOOLKIT_PROFILE_BLOCK("Compute Local Gradients");
-        tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
+        tbb::parallel_for(size_t(0), num_collisions, [&](size_t i) {
             const TCollision& collision = collisions[i];
 
             const VectorMaxNd local_grad = this->gradient(
diff --git a/tests/src/tests/potential/test_full_dof_assembly.cpp b/tests/src/tests/potential/test_full_dof_assembly.cpp
index e69a3b1bb..655f0c582 100644
--- a/tests/src/tests/potential/test_full_dof_assembly.cpp
+++ b/tests/src/tests/potential/test_full_dof_assembly.cpp
@@ -50,10 +50,16 @@ TEST_CASE(
             mesh.to_full_dof(potential.gradient(collisions, mesh, X));
 
         REQUIRE(grad_folded.size() == mesh.full_ndof());
-        // tbb::combinable partitions work nondeterministically, so the two
-        // calls may sum in different orders; compare with a tight tolerance.
+        // The folded and mapped paths sum in different orders; compare with
+        // a tight tolerance.
         const double scale = std::max(1.0, grad_mapped.norm());
         CHECK((grad_folded - grad_mapped).norm() <= 1e-13 * scale);
+
+        // Repeated evaluations may differ by floating-point rounding (the
+        // parallel summation order is not fixed) but must agree tightly.
+        const Eigen::VectorXd grad_folded2 =
+            potential.gradient(collisions, mesh, X, /*in_full_dof=*/true);
+        CHECK((grad_folded - grad_folded2).norm() <= 1e-13 * scale);
     }
 
     SECTION("hessian")

From 49e1f35f5f06a34ae189ae81f0fbfa91053ed747 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 31 Jul 2026 10:55:58 -0700
Subject: [PATCH 07/15] Make MeshFEMSparse block assembly the default (Phase 6)

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 (https://github.com/MeshFEM/MeshFEMCore/pull/1 for
Eigen 5 support, https://github.com/MeshFEM/MeshFEMSparse/pull/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>
---
 CMakeLists.txt                              |  9 ++++-
 cmake/recipes/meshfem_sparse.cmake          | 26 +++++--------
 docs/source/Doxyfile                        |  2 +-
 docs/source/about/dependencies.rst          |  9 +++++
 docs/source/cpp-api/utils.rst               | 19 +++++++++
 src/ipc/potentials/potential.cpp            |  8 ++++
 src/ipc/utils/CMakeLists.txt                |  1 -
 src/ipc/utils/meshfem_eigen_compat.hpp      | 43 ---------------------
 src/ipc/utils/meshfem_hessian_assembler.cpp | 26 +++++++++----
 src/ipc/utils/meshfem_hessian_assembler.hpp |  7 ++++
 10 files changed, 80 insertions(+), 70 deletions(-)
 delete mode 100644 src/ipc/utils/meshfem_eigen_compat.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 94ced74c3..0f14e5922 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -101,7 +101,7 @@ option(IPC_TOOLKIT_WITH_FILIB                 "Use filib for interval arithmetic
 option(IPC_TOOLKIT_WITH_INEXACT_CCD           "Use the original inexact CCD method of IPC"    OFF)
 option(IPC_TOOLKIT_WITH_PROFILER              "Enable performance profiler"                   OFF)
 option(IPC_TOOLKIT_WITH_TRACY                 "Enable Tracy frame profiler"                   OFF)
-option(IPC_TOOLKIT_WITH_MESHFEM_SPARSE        "Use MeshFEMSparse for block-accelerated assembly" OFF)
+option(IPC_TOOLKIT_WITH_MESHFEM_SPARSE        "Use MeshFEMSparse for block-accelerated assembly" ON)
 
 # Advanced options
 option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
@@ -152,6 +152,13 @@ if(IPC_TOOLKIT_WITH_CUDA)
   enable_language(CUDA)
 endif()
 
+## MeshFEMSparse block assembly requires per-vertex blocks in the derivative
+## layout, i.e., the default RowMajor ([x0, y0, z0, x1, ...]) ordering.
+if(IPC_TOOLKIT_WITH_MESHFEM_SPARSE AND NOT IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT STREQUAL "RowMajor")
+  message(WARNING "MeshFEMSparse assembly requires IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor. Continuing without MeshFEMSparse.")
+  set(IPC_TOOLKIT_WITH_MESHFEM_SPARSE OFF CACHE BOOL "Use MeshFEMSparse for block-accelerated assembly" FORCE)
+endif()
+
 ## SIMD support
 if(IPC_TOOLKIT_WITH_SIMD)
   # Figure out SIMD support
diff --git a/cmake/recipes/meshfem_sparse.cmake b/cmake/recipes/meshfem_sparse.cmake
index ecc9805a1..b3d563b04 100644
--- a/cmake/recipes/meshfem_sparse.cmake
+++ b/cmake/recipes/meshfem_sparse.cmake
@@ -21,17 +21,22 @@ include(eigen)
 include(onetbb)
 find_package(Threads REQUIRED)
 
+# TEMPORARY: pinned to zfergus' forks, which carry two fixes submitted
+# upstream — Eigen 5 support (https://github.com/MeshFEM/MeshFEMCore/pull/1)
+# and an out-of-bounds read on empty block columns
+# (https://github.com/MeshFEM/MeshFEMSparse/pull/1). Repoint to
+# MeshFEM/MeshFEMCore and MeshFEM/MeshFEMSparse once the PRs merge.
 include(CPM)
 CPMAddPackage(
     NAME MeshFEMCore
-    URL "https://github.com/MeshFEM/MeshFEMCore/archive/24e81c425e85eee3ed79af000e82ef1a75bbe696.zip"
-    URL_HASH SHA256=7d505f57af2a4b2fbc01668d57c1296d646c6aa0aa923c9d70e814b55987aff7
+    URL "https://github.com/zfergus/MeshFEMCore/archive/8d0e84788189748d9e906cc7f807507a3cb4b2ef.zip"
+    URL_HASH SHA256=71fe52e49276a401ae64d692dc26aea4147b1baa636bc78082d3fd0061eb4ccb
     DOWNLOAD_ONLY YES
 )
 CPMAddPackage(
     NAME MeshFEMSparse
-    URL "https://github.com/MeshFEM/MeshFEMSparse/archive/efe1af87cf6f7d04359628552e84c666b8612f4c.zip"
-    URL_HASH SHA256=e23cb168e9e709662083116e73876f97c23edd60413108bd70bde5a5a1684dd5
+    URL "https://github.com/zfergus/MeshFEMSparse/archive/ade01a1775f1911c0bf373190001923f452ad6cc.zip"
+    URL_HASH SHA256=6545cb0aa7b8513a370dc8561ea5abdd5ec484d03ded7e64f48e5dbfc94c49b0
     DOWNLOAD_ONLY YES
 )
 
@@ -81,19 +86,6 @@ target_compile_definitions(MeshFEMSparse PUBLIC
     _USE_MATH_DEFINES
 )
 
-# MeshFEM targets Eigen 3.4, but Eigen 5 removed
-# Eigen::internal::make_coherent, which MeshFEMCore's
-# AutomaticDifferentiation.hh references. Force-include a shim reimplementing
-# it (see src/ipc/utils/meshfem_eigen_compat.hpp); TUs outside this target
-# that include MeshFEM headers must include the shim first themselves.
-set(MESHFEM_EIGEN_COMPAT_HEADER
-    "${PROJECT_SOURCE_DIR}/src/ipc/utils/meshfem_eigen_compat.hpp")
-if(MSVC)
-    target_compile_options(MeshFEMSparse PRIVATE "/FI${MESHFEM_EIGEN_COMPAT_HEADER}")
-else()
-    target_compile_options(MeshFEMSparse PRIVATE "-include" "${MESHFEM_EIGEN_COMPAT_HEADER}")
-endif()
-
 # ipc_toolkit is compiled with EIGEN_DONT_VECTORIZE=1 when SIMD is enabled;
 # compiling the same Eigen templates with different vectorization settings is
 # an ODR violation with real alignment/layout consequences.
diff --git a/docs/source/Doxyfile b/docs/source/Doxyfile
index af1193c20..6a7d7d1d2 100644
--- a/docs/source/Doxyfile
+++ b/docs/source/Doxyfile
@@ -2315,7 +2315,7 @@ INCLUDE_FILE_PATTERNS  =
 # recursively expanded use the := operator instead of the = operator.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-PREDEFINED             = IPC_TOOLKIT_WITH_CORRECT_CCD IPC_TOOLKIT_WITH_ROBIN_MAP IPC_TOOLKIT_WITH_ABSEIL IPC_TOOLKIT_WITH_FILIB
+PREDEFINED             = IPC_TOOLKIT_WITH_INEXACT_CCD IPC_TOOLKIT_WITH_ROBIN_MAP IPC_TOOLKIT_WITH_ABSEIL IPC_TOOLKIT_WITH_FILIB IPC_TOOLKIT_WITH_MESHFEM_SPARSE
 
 # If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
 # tag can be used to specify a list of macro names that should be expanded. The
diff --git a/docs/source/about/dependencies.rst b/docs/source/about/dependencies.rst
index 5ca6d6760..3a3a9ca09 100644
--- a/docs/source/about/dependencies.rst
+++ b/docs/source/about/dependencies.rst
@@ -87,6 +87,12 @@ Additionally, IPC Toolkit may optionally use the following libraries:
       - `github.com/zfergus/filib <https://github.com/zfergus/filib>`_
       - |:white_check_mark:|
       - ``IPC_TOOLKIT_WITH_FILIB``
+    * - MeshFEMSparse
+      - Block-CSC data structures for fast Hessian assembly (see :cpp:class:`ipc::MeshFEMHessianAssembler`)
+      - MIT
+      - `github.com/MeshFEM/MeshFEMSparse <https://github.com/MeshFEM/MeshFEMSparse>`_
+      - |:white_check_mark:|
+      - ``IPC_TOOLKIT_WITH_MESHFEM_SPARSE``
     * - nlohmann/json
       - JSON parsing for profiler and tests
       - MIT
@@ -114,6 +120,9 @@ Additionally, IPC Toolkit may optionally use the following libraries:
 
 Some of these libraries are enabled by default, and some are not. You can enable or disable them by passing the appropriate CMake option when you configure the IPC Toolkit build.
 
+.. note::
+    ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise), and it is currently pinned to a fork carrying two fixes submitted upstream (`MeshFEM/MeshFEMCore#1 <https://github.com/MeshFEM/MeshFEMCore/pull/1>`_, `MeshFEM/MeshFEMSparse#1 <https://github.com/MeshFEM/MeshFEMSparse/pull/1>`_).
+
 .. warning::
     ``filib`` is licensed under `LGPL-2.1 <https://github.com/zfergus/filib/blob/main/LICENSE>`_ and as such it is required to be dynamically linked. Doing so automatically is a challenge, so by default we use static linkage. Enabling dynamic linkage requires copying the ``.so``/``.dylib``/``.dll`` file to the binary directory or system path. To enable this, set the CMake option ``FILIB_BUILD_SHARED_LIBS`` to ``ON`` and add this CMake code to copy the shared library object to the binary directory:
 
diff --git a/docs/source/cpp-api/utils.rst b/docs/source/cpp-api/utils.rst
index 8d494bc4c..b2769c0f9 100644
--- a/docs/source/cpp-api/utils.rst
+++ b/docs/source/cpp-api/utils.rst
@@ -15,6 +15,25 @@ Positive Semi-Definite Projection
 
 .. doxygenenum:: ipc::PSDProjectionMethod
 
+Hessian Assembly
+----------------
+
+Pluggable backends for assembling per-collision Hessians into a global
+matrix (see :cpp:func:`ipc::Potential::assemble_hessian`).
+
+.. doxygenclass:: ipc::HessianAssembler
+.. doxygenclass:: ipc::TripletHessianAssembler
+
+The following backend is available when the toolkit is compiled with
+``IPC_TOOLKIT_WITH_MESHFEM_SPARSE`` (the default), in which case it is also
+what :cpp:func:`ipc::Potential::hessian` uses internally. It assembles into
+`MeshFEMSparse <https://github.com/MeshFEM/MeshFEMSparse>`_'s block-CSC data
+structures (no triplets, no ``setFromTriplets``) and reuses the sparsity
+pattern across assemblies, making repeated contact Hessians roughly an order
+of magnitude faster than the triplet path on large scenes.
+
+.. doxygenclass:: ipc::MeshFEMHessianAssembler
+
 Eigen Extensions
 ----------------
 
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index 4db2a69f9..43558028d 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -5,6 +5,7 @@
 #include <ipc/utils/hessian_assembler.hpp>
 #include <ipc/utils/local_to_global.hpp>
 #include <ipc/utils/logger.hpp>
+#include <ipc/utils/meshfem_hessian_assembler.hpp>
 #include <ipc/utils/profiler.hpp>
 
 #include <tbb/blocked_range.h>
@@ -233,10 +234,17 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const bool fold_to_full = in_full_dof && mesh.is_selection_dof_map();
     const bool map_to_full = in_full_dof && !fold_to_full;
 
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+    MeshFEMHessianAssembler assembler;
+    assemble_hessian(
+        collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full);
+    const Eigen::SparseMatrix<double> hess = assembler.take_matrix();
+#else
     TripletHessianAssembler assembler;
     assemble_hessian(
         collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full);
     const Eigen::SparseMatrix<double> hess = assembler.get_matrix();
+#endif
 
     return map_to_full ? mesh.to_full_dof(hess) : hess;
 }
diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt
index 932d736a9..25a2368f6 100644
--- a/src/ipc/utils/CMakeLists.txt
+++ b/src/ipc/utils/CMakeLists.txt
@@ -10,7 +10,6 @@ set(SOURCES
   logger.hpp
   matrix_cache.cpp
   matrix_cache.hpp
-  meshfem_eigen_compat.hpp
   meshfem_hessian_assembler.cpp
   meshfem_hessian_assembler.hpp
   merge_thread_local.hpp
diff --git a/src/ipc/utils/meshfem_eigen_compat.hpp b/src/ipc/utils/meshfem_eigen_compat.hpp
deleted file mode 100644
index 13599d9a3..000000000
--- a/src/ipc/utils/meshfem_eigen_compat.hpp
+++ /dev/null
@@ -1,43 +0,0 @@
-#pragma once
-
-// Eigen ≥ 5 removed Eigen::internal::make_coherent (previously provided by
-// unsupported/Eigen/AutoDiff), which MeshFEMCore's AutomaticDifferentiation.hh
-// still references. This shim reimplements it with the Eigen 3.4 semantics:
-// if exactly one of the two derivative vectors is empty, resize it to match
-// the other and zero it.
-//
-// It must be included (or force-included via -include//FI) before any
-// MeshFEMCore/MeshFEMSparse header. It intentionally has no includes of its
-// own so it can be force-included into the MeshFEMSparse build (see
-// cmake/recipes/meshfem_sparse.cmake).
-//
-// TODO: remove once MeshFEM supports Eigen 5 upstream.
-
-// NOTE: EIGEN_WORLD_VERSION is 3 forever; Eigen 5 moved to semver with
-// EIGEN_MAJOR_VERSION=5. If no Eigen header has been seen yet (the
-// force-include case), assume Eigen 5 — this repository pins Eigen 5.
-#if !defined(EIGEN_MAJOR_VERSION) || EIGEN_MAJOR_VERSION >= 5
-
-namespace Eigen {
-namespace internal {
-
-    template <typename DerTypeA, typename DerTypeB>
-    inline void make_coherent(const DerTypeA& a, const DerTypeB& b)
-    {
-        // Eigen 3.4's implementation const-casts too (the derivatives are
-        // semantically mutable scratch space of the AutoDiffScalar pair).
-        DerTypeA& a_ref = const_cast<DerTypeA&>(a); // NOLINT
-        DerTypeB& b_ref = const_cast<DerTypeB&>(b); // NOLINT
-        if (a_ref.size() == 0 && b_ref.size() != 0) {
-            a_ref.resize(b_ref.size());
-            a_ref.setZero();
-        } else if (b_ref.size() == 0 && a_ref.size() != 0) {
-            b_ref.resize(a_ref.size());
-            b_ref.setZero();
-        }
-    }
-
-} // namespace internal
-} // namespace Eigen
-
-#endif
diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
index 2bb1bd64b..4ad6e7e42 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.cpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -3,7 +3,6 @@
 #ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
 
 #include <ipc/utils/logger.hpp>
-#include <ipc/utils/meshfem_eigen_compat.hpp> // must precede MeshFEM headers
 #include <ipc/utils/profiler.hpp>
 
 #include <MeshFEMSparse/SystemAssembler.hh>
@@ -35,6 +34,7 @@ struct MeshFEMHessianAssembler::ImplBase {
         const std::array<index_t, 4>& vertex_ids) = 0;
 
     virtual const Eigen::SparseMatrix<double>& to_eigen() const = 0;
+    virtual Eigen::SparseMatrix<double> take_eigen() = 0;
 };
 
 /// Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM's own
@@ -138,12 +138,11 @@ struct MeshFEMHessianAssembler::Impl final
     // and reused while the block pattern is unchanged; only the values are
     // recomputed per call.
     //
-    // NOTE: We deliberately do not use BlockCSCHessian::toEigen here: its
-    // toScalar step assumes every block column is non-empty (true for FE
-    // Hessians, where every node belongs to an element, but false for contact
-    // Hessians, where most vertices are collision-free) and reads out of
-    // bounds otherwise. This implementation also skips toEigen's intermediate
-    // upper-triangle scalar matrix, symmetrizing directly instead.
+    // NOTE: We deliberately do not use BlockCSCHessian::toEigen here: it
+    // cannot reuse a cached structure across assemblies, and it goes through
+    // an intermediate upper-triangle scalar matrix (serially) where this
+    // implementation symmetrizes directly (in parallel) — roughly 2× faster
+    // even cold.
     //
     // Value layout (ContiguousBlocks=true, StoreFullDiagonalBlocks=true, the
     // library default): block entry ii occupies Ax[N²·ii, N²·(ii+1)),
@@ -162,6 +161,13 @@ struct MeshFEMHessianAssembler::Impl final
         return m_M;
     }
 
+    Eigen::SparseMatrix<double> take_eigen() override
+    {
+        to_eigen();
+        m_eigen_structure_valid = false; // m_M is about to be gutted
+        return std::move(m_M);
+    }
+
 private:
     static constexpr int N = dim;
 
@@ -349,6 +355,12 @@ const Eigen::SparseMatrix<double>& MeshFEMHessianAssembler::get_matrix() const
     return m_impl->to_eigen();
 }
 
+Eigen::SparseMatrix<double> MeshFEMHessianAssembler::take_matrix()
+{
+    assert(m_impl != nullptr);
+    return m_impl->take_eigen();
+}
+
 } // namespace ipc
 
 #endif // IPC_TOOLKIT_WITH_MESHFEM_SPARSE
diff --git a/src/ipc/utils/meshfem_hessian_assembler.hpp b/src/ipc/utils/meshfem_hessian_assembler.hpp
index 1695bd25b..fd123a1ea 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.hpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.hpp
@@ -57,6 +57,13 @@ class MeshFEMHessianAssembler final : public HessianAssembler {
     /// begin() call; copy it to keep a snapshot.
     const Eigen::SparseMatrix<double>& get_matrix() const;
 
+    /// @brief Like get_matrix(), but moves the matrix out of the assembler.
+    ///
+    /// Avoids a copy for one-shot use (e.g., Potential::hessian). The cached
+    /// Eigen structure is invalidated; the next get_matrix()/take_matrix()
+    /// rebuilds it.
+    Eigen::SparseMatrix<double> take_matrix();
+
     /// @brief Number of vanished blocks tolerated before a pattern rebuild.
     ///
     /// When a stencil introduces a vertex pair absent from the cached pattern,

From 7fc84528a2d063b39accc5bfc25c784c103dab8c Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 31 Jul 2026 11:24:08 -0700
Subject: [PATCH 08/15] Suppress clang-tidy identifier-naming on m_H/m_M matrix
 members

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>
---
 src/ipc/utils/meshfem_hessian_assembler.cpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
index 4ad6e7e42..7870ff560 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.cpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -301,7 +301,8 @@ struct MeshFEMHessianAssembler::Impl final
 
     MeshFEM::SystemAssembler<dim> m_assembler;
     VarStructure m_vars;
-    std::unique_ptr<MeshFEM::BlockCSCHessian<VarStructure>> m_H;
+    std::unique_ptr<MeshFEM::BlockCSCHessian<VarStructure>>
+        m_H; // NOLINT(readability-identifier-naming): H = Hessian matrix
     MeshFEM::VarLocks m_locks;
     /// Stencil count of the cached pattern (assume-unchanged sanity check).
     size_t m_num_stencils = 0;
@@ -310,7 +311,8 @@ struct MeshFEMHessianAssembler::Impl final
     mutable bool m_eigen_structure_valid = false;
     mutable std::vector<std::ptrdiff_t> m_sym_col_start;
     mutable std::vector<BlockEntry> m_sym_entries;
-    mutable Eigen::SparseMatrix<double> m_M;
+    mutable Eigen::SparseMatrix<double>
+        m_M; // NOLINT(readability-identifier-naming): M = matrix
 };
 
 MeshFEMHessianAssembler::MeshFEMHessianAssembler() = default;

From 9e07a0e1e77dea123558695c35a3e9ce7034e4cc Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 31 Jul 2026 13:10:53 -0700
Subject: [PATCH 09/15] Expose the block-CSC matrix and bind the assemblers to
 Python

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>
---
 cmake/recipes/meshfem_sparse.cmake            |  4 +-
 python/src/bindings.cpp                       |  1 +
 python/src/potentials/potential.hpp           | 20 +++++
 python/src/utils/CMakeLists.txt               |  1 +
 python/src/utils/bindings.hpp                 |  1 +
 python/src/utils/hessian_assembler.cpp        | 80 +++++++++++++++++++
 src/ipc/utils/meshfem_hessian_assembler.cpp   | 18 +++++
 src/ipc/utils/meshfem_hessian_assembler.hpp   | 30 +++++++
 .../tests/potential/test_meshfem_assembly.cpp | 56 +++++++++++++
 9 files changed, 209 insertions(+), 2 deletions(-)
 create mode 100644 python/src/utils/hessian_assembler.cpp

diff --git a/cmake/recipes/meshfem_sparse.cmake b/cmake/recipes/meshfem_sparse.cmake
index b3d563b04..7f4a99477 100644
--- a/cmake/recipes/meshfem_sparse.cmake
+++ b/cmake/recipes/meshfem_sparse.cmake
@@ -35,8 +35,8 @@ CPMAddPackage(
 )
 CPMAddPackage(
     NAME MeshFEMSparse
-    URL "https://github.com/zfergus/MeshFEMSparse/archive/ade01a1775f1911c0bf373190001923f452ad6cc.zip"
-    URL_HASH SHA256=6545cb0aa7b8513a370dc8561ea5abdd5ec484d03ded7e64f48e5dbfc94c49b0
+    URL "https://github.com/zfergus/MeshFEMSparse/archive/15a92834189ade69ed9e00373adc3036a72cd40a.zip"
+    URL_HASH SHA256=882f3a126f59023b4d4aba6e1b96b89d5f82f2e80bfe086c1f3afd435f8e4c80
     DOWNLOAD_ONLY YES
 )
 
diff --git a/python/src/bindings.cpp b/python/src/bindings.cpp
index d20970c2e..8ab40e974 100644
--- a/python/src/bindings.cpp
+++ b/python/src/bindings.cpp
@@ -117,6 +117,7 @@ PYBIND11_MODULE(ipctk, m)
     define_tangential_adhesion_potential(m);
 
     // utils
+    define_hessian_assembler(m);
     define_logger(m);
     define_profiler(m);
     define_thread_limiter(m);
diff --git a/python/src/potentials/potential.hpp b/python/src/potentials/potential.hpp
index 33ad4854d..04c4886b3 100644
--- a/python/src/potentials/potential.hpp
+++ b/python/src/potentials/potential.hpp
@@ -1,6 +1,7 @@
 #include <common.hpp>
 
 #include <ipc/potentials/potential.hpp>
+#include <ipc/utils/hessian_assembler.hpp>
 
 using namespace ipc;
 
@@ -73,6 +74,25 @@ void define_potential_methods(PyClass& potential)
             "collisions"_a, "mesh"_a, "X"_a,
             "project_hessian_to_psd"_a = PSDProjectionMethod::NONE,
             "in_full_dof"_a = false)
+        .def(
+            "assemble_hessian", &Potential<TCollisions>::assemble_hessian,
+            R"ipc_Qu8mg5v7(
+            Assemble the Hessian of the potential using a custom assembler.
+
+            Evaluates the local Hessian of every collision (in parallel) and feeds each to `assembler`, then leaves the
+            result in the assembler (use its get_matrix()). Reuse one assembler across calls to also reuse its sparsity pattern.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                X: Degrees of freedom of the collision mesh (e.g., vertices or velocities).
+                assembler: The assembler that accumulates the local Hessians.
+                project_hessian_to_psd: Make sure the hessian is positive semi-definite.
+                in_full_dof: If true, stencil vertex IDs are remapped to full-mesh vertex IDs (requires mesh.is_selection_dof_map).
+            )ipc_Qu8mg5v7",
+            "collisions"_a, "mesh"_a, "X"_a, "assembler"_a,
+            "project_hessian_to_psd"_a = PSDProjectionMethod::NONE,
+            "in_full_dof"_a = false)
         .def(
             "__call__",
             py::overload_cast<const TCollision&, Eigen::ConstRef<VectorMax12d>>(
diff --git a/python/src/utils/CMakeLists.txt b/python/src/utils/CMakeLists.txt
index 7a17ed16c..ea828fe0d 100644
--- a/python/src/utils/CMakeLists.txt
+++ b/python/src/utils/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(SOURCES
   eigen_ext.cpp
+  hessian_assembler.cpp
   logger.cpp
   profiler.cpp
   thread_limiter.cpp
diff --git a/python/src/utils/bindings.hpp b/python/src/utils/bindings.hpp
index 368962c93..db2de2056 100644
--- a/python/src/utils/bindings.hpp
+++ b/python/src/utils/bindings.hpp
@@ -3,6 +3,7 @@
 #include <pybind11/pybind11.h>
 
 void define_eigen_ext(py::module_& m);
+void define_hessian_assembler(py::module_& m);
 void define_logger(py::module_& m);
 void define_profiler(py::module_& m);
 void define_thread_limiter(py::module_& m);
diff --git a/python/src/utils/hessian_assembler.cpp b/python/src/utils/hessian_assembler.cpp
new file mode 100644
index 000000000..e9573b11a
--- /dev/null
+++ b/python/src/utils/hessian_assembler.cpp
@@ -0,0 +1,80 @@
+#include <common.hpp>
+
+#include <ipc/utils/hessian_assembler.hpp>
+#include <ipc/utils/meshfem_hessian_assembler.hpp>
+
+using namespace ipc;
+
+void define_hessian_assembler(py::module_& m)
+{
+    // NOTE: HessianAssembler is intentionally not extendable from Python. The
+    // driver calls add_local_hessian() once per collision, so a Python-defined
+    // assembler would acquire the GIL hundreds of thousands of times per
+    // assembly. Only the built-in backends are exposed.
+    py::class_<HessianAssembler>(
+        m, "HessianAssembler", py::is_final(), R"ipc_Qu8mg5v7(
+        Abstract sink for assembling local (per-collision) Hessians into a global matrix.
+
+        Cannot be constructed or subclassed from Python; use TripletHessianAssembler or MeshFEMHessianAssembler.
+        )ipc_Qu8mg5v7");
+
+    py::class_<TripletHessianAssembler, HessianAssembler>(
+        m, "TripletHessianAssembler", py::is_final(), R"ipc_Qu8mg5v7(
+        Assembles through thread-local triplet caches and Eigen's setFromTriplets.
+        )ipc_Qu8mg5v7")
+        .def(py::init())
+        .def(
+            "get_matrix", &TripletHessianAssembler::get_matrix,
+            R"ipc_Qu8mg5v7(
+            Merge the thread-local caches and build the global matrix.
+
+            Call once, after assembly; the internal caches are consumed.
+
+            Returns:
+                The assembled Hessian.
+            )ipc_Qu8mg5v7");
+
+#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
+    py::class_<MeshFEMHessianAssembler, HessianAssembler>(
+        m, "MeshFEMHessianAssembler", py::is_final(), R"ipc_Qu8mg5v7(
+        Assembles into MeshFEMSparse's block-CSC data structures.
+
+        Reuse one instance across assemblies (e.g., one per Newton solve) to also reuse the sparsity pattern, which is where most of the speedup over the triplet assembler comes from.
+        )ipc_Qu8mg5v7")
+        .def(py::init())
+        .def(
+            "get_matrix", &MeshFEMHessianAssembler::get_matrix,
+            R"ipc_Qu8mg5v7(
+            Convert the assembled matrix to a full symmetric sparse matrix.
+
+            Returns:
+                The assembled Hessian.
+            )ipc_Qu8mg5v7")
+        .def_property(
+            "stale_block_tolerance",
+            &MeshFEMHessianAssembler::stale_block_tolerance,
+            &MeshFEMHessianAssembler::set_stale_block_tolerance,
+            R"ipc_Qu8mg5v7(
+            Number of vanished blocks tolerated before the pattern is rebuilt.
+
+            A stencil that introduces a new vertex pair always forces a rebuild. When blocks merely disappear (e.g., a contact separates), the pattern is reused as long as at most this many blocks vanished. Stale blocks assemble to explicit zeros, so values stay correct. Defaults to 0.
+            )ipc_Qu8mg5v7")
+        .def_property(
+            "assume_unchanged_stencils",
+            &MeshFEMHessianAssembler::assume_unchanged_stencils,
+            &MeshFEMHessianAssembler::set_assume_unchanged_stencils,
+            R"ipc_Qu8mg5v7(
+            Whether to skip change detection entirely.
+
+            When enabled, the cached pattern is reused without comparing the stencils against it, as long as the stencil count is unchanged. Use this only when you know the collision set is identical to the previous assembly (e.g., reassembling with a different PSD projection or stiffness); on large scenes change detection costs about as much as rebuilding the pattern.
+
+            Warning:
+                If the stencils did change while the count stayed equal, assembly reads out of bounds.
+            )ipc_Qu8mg5v7")
+        .def_property_readonly(
+            "reused_pattern", &MeshFEMHessianAssembler::reused_pattern,
+            R"ipc_Qu8mg5v7(
+            Whether the last assembly reused the cached sparsity pattern.
+            )ipc_Qu8mg5v7");
+#endif
+}
diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
index 7870ff560..e36020a2d 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.cpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -35,6 +35,7 @@ struct MeshFEMHessianAssembler::ImplBase {
 
     virtual const Eigen::SparseMatrix<double>& to_eigen() const = 0;
     virtual Eigen::SparseMatrix<double> take_eigen() = 0;
+    virtual const MeshFEM::BlockCSCHessianBase& block_matrix() const = 0;
 };
 
 /// Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM's own
@@ -168,6 +169,12 @@ struct MeshFEMHessianAssembler::Impl final
         return std::move(m_M);
     }
 
+    const MeshFEM::BlockCSCHessianBase& block_matrix() const override
+    {
+        assert(m_H != nullptr);
+        return *m_H;
+    }
+
 private:
     static constexpr int N = dim;
 
@@ -363,6 +370,17 @@ Eigen::SparseMatrix<double> MeshFEMHessianAssembler::take_matrix()
     return m_impl->take_eigen();
 }
 
+const MeshFEM::BlockCSCHessianBase&
+MeshFEMHessianAssembler::block_matrix() const
+{
+    if (m_impl == nullptr) {
+        log_and_throw_error(
+            "MeshFEMHessianAssembler::block_matrix() called before the first "
+            "assembly!");
+    }
+    return m_impl->block_matrix();
+}
+
 } // namespace ipc
 
 #endif // IPC_TOOLKIT_WITH_MESHFEM_SPARSE
diff --git a/src/ipc/utils/meshfem_hessian_assembler.hpp b/src/ipc/utils/meshfem_hessian_assembler.hpp
index fd123a1ea..e5d7913bb 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.hpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.hpp
@@ -8,6 +8,13 @@
 
 #include <memory>
 
+namespace MeshFEM {
+/// Forward declaration so block_matrix() can be exposed without pulling
+/// MeshFEMSparse's headers into this one. Include
+/// <MeshFEMSparse/BlockCSCHessian.hh> to use the returned object.
+struct BlockCSCHessianBase;
+} // namespace MeshFEM
+
 namespace ipc {
 
 // The MeshFEM backend scatters per-vertex d×d blocks, which requires
@@ -64,6 +71,29 @@ class MeshFEMHessianAssembler final : public HessianAssembler {
     /// rebuilds it.
     Eigen::SparseMatrix<double> take_matrix();
 
+    /// @brief Access the assembled matrix in its native block-CSC form.
+    ///
+    /// Skips the conversion performed by get_matrix(), which is worth doing if
+    /// you can consume the block format directly (e.g., MeshFEM's block SpMV
+    /// or its Cholesky factorizers). Include
+    /// <MeshFEMSparse/BlockCSCHessian.hh> to use the result.
+    ///
+    /// The matrix is symmetric with only the upper triangle stored, and its
+    /// block sparsity pattern covers the interacting vertex pairs of the last
+    /// assembly (possibly with explicitly-zero stale blocks if the pattern was
+    /// reused). The reference stays valid until the next begin() call, which
+    /// may rebuild the pattern in place.
+    ///
+    /// @note A contact Hessian has no diagonal block for any vertex that is
+    /// not in contact, so most block columns are empty. Reading operations
+    /// handle that (trace() skips the missing diagonals; apply() and the
+    /// sparsity pattern are unaffected), but the diagonal-mutating operations
+    /// (addDiag(), setDiag()) have no entry to write and throw. Insert the
+    /// missing diagonal blocks first if you need them.
+    ///
+    /// @throws std::runtime_error if called before the first assembly.
+    const MeshFEM::BlockCSCHessianBase& block_matrix() const;
+
     /// @brief Number of vanished blocks tolerated before a pattern rebuild.
     ///
     /// When a stencil introduces a vertex pair absent from the cached pattern,
diff --git a/tests/src/tests/potential/test_meshfem_assembly.cpp b/tests/src/tests/potential/test_meshfem_assembly.cpp
index f6ff1ad19..ef14dd204 100644
--- a/tests/src/tests/potential/test_meshfem_assembly.cpp
+++ b/tests/src/tests/potential/test_meshfem_assembly.cpp
@@ -14,10 +14,16 @@
 
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/generators/catch_generators.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 
 #include <ipc/utils/hessian_assembler.hpp>
 #include <ipc/utils/meshfem_hessian_assembler.hpp>
 
+// Deliberately included here rather than by meshfem_hessian_assembler.hpp:
+// block_matrix() is declared against a forward declaration, so only callers
+// that want the block matrix pay for MeshFEM's headers.
+#include <MeshFEMSparse/BlockCSCHessian.hh>
+
 #include <limits>
 
 using namespace ipc;
@@ -157,6 +163,56 @@ TEST_CASE("MeshFEM assembly pattern reuse", "[potential][assembly][meshfem]")
     check_matches(assembler.get_matrix(), reference(collisions_small));
 }
 
+TEST_CASE(
+    "MeshFEM assembly exposes the block-CSC matrix",
+    "[potential][assembly][meshfem]")
+{
+    // A downstream user should be able to consume the native block-CSC matrix
+    // (e.g., to hand it to MeshFEM's solvers) instead of paying for the Eigen
+    // conversion. This also checks that the forward declaration in our header
+    // is enough: MeshFEMSparse's header is included by this test, not by
+    // meshfem_hessian_assembler.hpp.
+    const auto& spec = ipc::tests::assembly_scene_specs().at(0); // two-cubes
+
+    const std::optional<ipc::tests::AssemblyScene> maybe_scene =
+        ipc::tests::build_assembly_scene(spec);
+    if (!maybe_scene.has_value()) {
+        SKIP(fmt::format("Scene '{}' is unavailable.", spec.mesh_name));
+    }
+    const ipc::tests::AssemblyScene& scene = maybe_scene.value();
+
+    MeshFEMHessianAssembler assembler;
+    scene.potential().assemble_hessian(
+        scene.collisions(), scene.mesh(), scene.vertices(), assembler);
+
+    const MeshFEM::BlockCSCHessianBase& H = assembler.block_matrix();
+    const Eigen::SparseMatrix<double> H_eigen = assembler.get_matrix();
+
+    REQUIRE(H.numScalarCols() == size_t(H_eigen.cols()));
+    // Only the upper triangle is stored, so the block matrix holds fewer
+    // scalar entries than the symmetrized Eigen matrix.
+    CHECK(H.scalarNNZ() <= size_t(H_eigen.nonZeros()));
+
+    // trace() must skip the empty block columns of a contact pattern rather
+    // than reading the preceding column's storage.
+    CHECK_THAT(
+        H.trace(),
+        Catch::Matchers::WithinRel(Eigen::MatrixXd(H_eigen).trace(), 1e-12));
+
+    // addDiag()/setDiag() cannot work without every diagonal block present,
+    // so they must reject a contact pattern instead of corrupting it.
+    const Eigen::VectorXd ones =
+        Eigen::VectorXd::Ones(Eigen::Index(H.numScalarCols()));
+    CHECK_THROWS(const_cast<MeshFEM::BlockCSCHessianBase&>(H).addDiag(ones));
+
+    // Exercise MeshFEM's own API on the result: y = H x must match Eigen's.
+    const Eigen::VectorXd x =
+        Eigen::VectorXd::Random(Eigen::Index(H.numScalarCols()));
+    const Eigen::VectorXd y_expected = H_eigen * x;
+    const Eigen::VectorXd y = H.apply(x);
+    CHECK((y - y_expected).norm() <= 1e-12 * std::max(1.0, y_expected.norm()));
+}
+
 TEST_CASE(
     "MeshFEM assembly with no collisions", "[potential][assembly][meshfem]")
 {

From 54c43722af04beea0e0a09fa443435ebdb46de57 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 31 Jul 2026 13:48:31 -0700
Subject: [PATCH 10/15] Document in_full_dof and assembler reuse in the
 simulation tutorial

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>
---
 cmake/recipes/meshfem_sparse.cmake   |   6 +-
 docs/source/about/dependencies.rst   |   2 +-
 docs/source/tutorials/simulation.rst | 105 ++++++++++++++++++++++++++-
 3 files changed, 108 insertions(+), 5 deletions(-)

diff --git a/cmake/recipes/meshfem_sparse.cmake b/cmake/recipes/meshfem_sparse.cmake
index 7f4a99477..7bd0b2220 100644
--- a/cmake/recipes/meshfem_sparse.cmake
+++ b/cmake/recipes/meshfem_sparse.cmake
@@ -21,9 +21,9 @@ include(eigen)
 include(onetbb)
 find_package(Threads REQUIRED)
 
-# TEMPORARY: pinned to zfergus' forks, which carry two fixes submitted
-# upstream — Eigen 5 support (https://github.com/MeshFEM/MeshFEMCore/pull/1)
-# and an out-of-bounds read on empty block columns
+# TEMPORARY: pinned to zfergus' forks, which carry fixes submitted upstream:
+# Eigen 5 support (https://github.com/MeshFEM/MeshFEMCore/pull/1) and the
+# mishandling of empty block columns, which a contact pattern is full of
 # (https://github.com/MeshFEM/MeshFEMSparse/pull/1). Repoint to
 # MeshFEM/MeshFEMCore and MeshFEM/MeshFEMSparse once the PRs merge.
 include(CPM)
diff --git a/docs/source/about/dependencies.rst b/docs/source/about/dependencies.rst
index 3a3a9ca09..d738b1e3c 100644
--- a/docs/source/about/dependencies.rst
+++ b/docs/source/about/dependencies.rst
@@ -121,7 +121,7 @@ Additionally, IPC Toolkit may optionally use the following libraries:
 Some of these libraries are enabled by default, and some are not. You can enable or disable them by passing the appropriate CMake option when you configure the IPC Toolkit build.
 
 .. note::
-    ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise), and it is currently pinned to a fork carrying two fixes submitted upstream (`MeshFEM/MeshFEMCore#1 <https://github.com/MeshFEM/MeshFEMCore/pull/1>`_, `MeshFEM/MeshFEMSparse#1 <https://github.com/MeshFEM/MeshFEMSparse/pull/1>`_).
+    ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise), and it is currently pinned to a fork carrying fixes submitted upstream (`MeshFEM/MeshFEMCore#1 <https://github.com/MeshFEM/MeshFEMCore/pull/1>`_, `MeshFEM/MeshFEMSparse#1 <https://github.com/MeshFEM/MeshFEMSparse/pull/1>`_).
 
 .. warning::
     ``filib`` is licensed under `LGPL-2.1 <https://github.com/zfergus/filib/blob/main/LICENSE>`_ and as such it is required to be dynamically linked. Doing so automatically is a challenge, so by default we use static linkage. Enabling dynamic linkage requires copying the ``.so``/``.dylib``/``.dll`` file to the binary directory or system path. To enable this, set the CMake option ``FILIB_BUILD_SHARED_LIBS`` to ``ON`` and add this CMake code to copy the shared library object to the binary directory:
diff --git a/docs/source/tutorials/simulation.rst b/docs/source/tutorials/simulation.rst
index 7faf069be..7a13239ef 100644
--- a/docs/source/tutorials/simulation.rst
+++ b/docs/source/tutorials/simulation.rst
@@ -152,6 +152,36 @@ When computing the gradient and Hessian of the potentials, the derivatives will
             hess = B.hessian(collision, collision_mesh, vertices)
             hess_full = collision_mesh.to_full_dof(hess)
 
+If the full DOF derivatives are all you need, ask for them directly with ``in_full_dof`` instead of mapping afterwards. The stencil indices are remapped during assembly, so the gradient is scattered into full DOF as it is summed and the Hessian never goes through the two sparse matrix products that ``to_full_dof`` performs:
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            Eigen::VectorXd grad_full = B.gradient(
+                collisions, collision_mesh, vertices, /*in_full_dof=*/true);
+
+            Eigen::SparseMatrix<double> hess_full = B.hessian(
+                collisions, collision_mesh, vertices,
+                ipc::PSDProjectionMethod::NONE, /*in_full_dof=*/true);
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            grad_full = B.gradient(
+                collisions, collision_mesh, vertices, in_full_dof=True)
+
+            hess_full = B.hessian(
+                collisions, collision_mesh, vertices, in_full_dof=True)
+
+The results match ``collision_mesh.to_full_dof(...)`` up to the order in which the local contributions are summed.
+
+.. note::
+    Remapping indices only works when the map from collision to full DOF is a pure selection, which is the case for every collision mesh built without a displacement map. You can check with ``collision_mesh.is_selection_dof_map()``. When a displacement map is present, ``in_full_dof`` still returns the right answer, but it does so by falling back to ``to_full_dof``, so there is nothing to gain.
+
 Codimensional Vertices
 ^^^^^^^^^^^^^^^^^^^^^^
 
@@ -263,4 +293,77 @@ To remedy this, we can project the Hessian onto the positive semidefinite (PSD)
     .. md-tab-item:: Python
 
         - ``ProjectToPSD.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
-        - ``ProjectToPSD.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
\ No newline at end of file
+        - ``ProjectToPSD.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
+
+Reusing the Hessian Assembler
+-----------------------------
+
+Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anything but a small scene, evaluating the local Hessians is the cheap part; most of the time goes into figuring out where those numbers belong in the global matrix. A Newton solve pays that cost on every iteration even though the contact set usually barely moves between iterations.
+
+``Potential::assemble_hessian`` takes the assembler as an argument so you can keep it alive across the whole solve and let it remember its work:
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // One assembler for the solve, not one per iteration.
+            ipc::MeshFEMHessianAssembler assembler;
+
+            for (int i = 0; i < max_iterations; i++) {
+                // ... update vertices and rebuild the collision set ...
+
+                B.assemble_hessian(
+                    collisions, collision_mesh, vertices, assembler,
+                    ipc::PSDProjectionMethod::CLAMP, /*in_full_dof=*/true);
+
+                // Valid until the next assembly; copy it if you need it longer.
+                const Eigen::SparseMatrix<double>& hess = assembler.get_matrix();
+
+                // ... solve for the Newton direction, line search, etc. ...
+            }
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            # One assembler for the solve, not one per iteration.
+            assembler = ipctk.MeshFEMHessianAssembler()
+
+            for i in range(max_iterations):
+                # ... update vertices and rebuild the collision set ...
+
+                B.assemble_hessian(
+                    collisions, collision_mesh, vertices, assembler,
+                    ipctk.PSDProjectionMethod.CLAMP, in_full_dof=True)
+
+                hess = assembler.get_matrix()
+
+                # ... solve for the Newton direction, line search, etc. ...
+
+The first call walks the collision stencils and builds a block sparsity pattern: one :math:`d \times d` block per pair of interacting vertices, rather than one entry per scalar. Later calls check the new stencils against that pattern. A contact that appears adds a block the pattern does not have, which forces a rebuild. A contact that separates only leaves a block behind that assembles to zero, which costs a little memory but does not change the matrix, so the pattern is kept as long as no more than ``stale_block_tolerance`` blocks have gone stale.
+
+Internally the matrix is symmetric with only its upper triangle stored, and ``get_matrix()`` mirrors it into a full Eigen matrix, reusing the cached structure whenever the pattern survives. If your solver speaks block CSC natively, ``block_matrix()`` hands you MeshFEM's matrix and skips that conversion. Our header only forward declares the type, so include ``<MeshFEMSparse/BlockCSCHessian.hh>`` yourself when you want it.
+
+When you reassemble without touching the collision set (a new stiffness, say, or the same Hessian under a different PSD projection) you can skip the comparison as well:
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            assembler.set_assume_unchanged_stencils(true);
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            assembler.assume_unchanged_stencils = True
+
+.. warning::
+    ``assume_unchanged_stencils`` trusts you. If the stencils changed while their count stayed the same, assembly reads past the end of the pattern. Debug builds check the assumption; release builds do not.
+
+.. note::
+    ``MeshFEMHessianAssembler`` requires the toolkit to be built with ``IPC_TOOLKIT_WITH_MESHFEM_SPARSE`` (on by default, and also what ``Potential::hessian`` uses internally). ``TripletHessianAssembler`` is always available and does what ``hessian`` has always done, gathering thread-local triplets and handing them to ``setFromTriplets``, but it has no pattern to carry over so reusing it saves nothing.
\ No newline at end of file

From 5bcea54e3772f7e2ceb652641d836b7b0b5710a7 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 31 Jul 2026 13:53:54 -0700
Subject: [PATCH 11/15] Use a more formal register in the new tutorial prose

---
 docs/source/tutorials/simulation.rst | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/docs/source/tutorials/simulation.rst b/docs/source/tutorials/simulation.rst
index 7a13239ef..f2758966c 100644
--- a/docs/source/tutorials/simulation.rst
+++ b/docs/source/tutorials/simulation.rst
@@ -152,7 +152,7 @@ When computing the gradient and Hessian of the potentials, the derivatives will
             hess = B.hessian(collision, collision_mesh, vertices)
             hess_full = collision_mesh.to_full_dof(hess)
 
-If the full DOF derivatives are all you need, ask for them directly with ``in_full_dof`` instead of mapping afterwards. The stencil indices are remapped during assembly, so the gradient is scattered into full DOF as it is summed and the Hessian never goes through the two sparse matrix products that ``to_full_dof`` performs:
+If only the full DOF derivatives are required, the optional ``in_full_dof`` parameter requests them directly rather than mapping after the fact. The stencil indices are remapped during assembly, so the gradient is scattered into full DOF as it is accumulated and the Hessian avoids the two sparse matrix products performed by ``to_full_dof``:
 
 .. md-tab-set::
 
@@ -177,10 +177,10 @@ If the full DOF derivatives are all you need, ask for them directly with ``in_fu
             hess_full = B.hessian(
                 collisions, collision_mesh, vertices, in_full_dof=True)
 
-The results match ``collision_mesh.to_full_dof(...)`` up to the order in which the local contributions are summed.
+The results agree with ``collision_mesh.to_full_dof(...)`` up to the order in which the local contributions are summed.
 
 .. note::
-    Remapping indices only works when the map from collision to full DOF is a pure selection, which is the case for every collision mesh built without a displacement map. You can check with ``collision_mesh.is_selection_dof_map()``. When a displacement map is present, ``in_full_dof`` still returns the right answer, but it does so by falling back to ``to_full_dof``, so there is nothing to gain.
+    Remapping the indices is only valid when the map from collision to full DOF is a pure selection, which holds for any collision mesh constructed without a displacement map. Use ``collision_mesh.is_selection_dof_map()`` to query this. When a displacement map is present, ``in_full_dof`` still returns the correct result, but it does so by applying ``to_full_dof`` internally and therefore offers no advantage.
 
 Codimensional Vertices
 ^^^^^^^^^^^^^^^^^^^^^^
@@ -298,9 +298,9 @@ To remedy this, we can project the Hessian onto the positive semidefinite (PSD)
 Reusing the Hessian Assembler
 -----------------------------
 
-Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anything but a small scene, evaluating the local Hessians is the cheap part; most of the time goes into figuring out where those numbers belong in the global matrix. A Newton solve pays that cost on every iteration even though the contact set usually barely moves between iterations.
+Each call to ``Potential::hessian`` constructs a sparse matrix from scratch. Except on small scenes, evaluating the local Hessians accounts for a minority of the cost; the bulk is spent determining where each local contribution belongs in the global matrix. A Newton solve repeats that work every iteration, even though the contact set typically changes little between iterations.
 
-``Potential::assemble_hessian`` takes the assembler as an argument so you can keep it alive across the whole solve and let it remember its work:
+``Potential::assemble_hessian`` accepts the assembler as a parameter, allowing a single instance to persist across the solve and retain its sparsity pattern:
 
 .. md-tab-set::
 
@@ -308,7 +308,7 @@ Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anyt
 
         .. code-block:: c++
 
-            // One assembler for the solve, not one per iteration.
+            // A single assembler for the entire solve, rather than one per iteration.
             ipc::MeshFEMHessianAssembler assembler;
 
             for (int i = 0; i < max_iterations; i++) {
@@ -318,7 +318,7 @@ Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anyt
                     collisions, collision_mesh, vertices, assembler,
                     ipc::PSDProjectionMethod::CLAMP, /*in_full_dof=*/true);
 
-                // Valid until the next assembly; copy it if you need it longer.
+                // Valid until the next assembly; copy it to retain it longer.
                 const Eigen::SparseMatrix<double>& hess = assembler.get_matrix();
 
                 // ... solve for the Newton direction, line search, etc. ...
@@ -328,7 +328,7 @@ Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anyt
 
         .. code-block:: python
 
-            # One assembler for the solve, not one per iteration.
+            # A single assembler for the entire solve, rather than one per iteration.
             assembler = ipctk.MeshFEMHessianAssembler()
 
             for i in range(max_iterations):
@@ -342,11 +342,11 @@ Each call to ``Potential::hessian`` builds a sparse matrix from scratch. On anyt
 
                 # ... solve for the Newton direction, line search, etc. ...
 
-The first call walks the collision stencils and builds a block sparsity pattern: one :math:`d \times d` block per pair of interacting vertices, rather than one entry per scalar. Later calls check the new stencils against that pattern. A contact that appears adds a block the pattern does not have, which forces a rebuild. A contact that separates only leaves a block behind that assembles to zero, which costs a little memory but does not change the matrix, so the pattern is kept as long as no more than ``stale_block_tolerance`` blocks have gone stale.
+The first call traverses the collision stencils and builds a block sparsity pattern, allocating one :math:`d \times d` block per interacting vertex pair instead of one entry per scalar. Subsequent calls compare the new stencils against the cached pattern. A newly active contact introduces a block the pattern does not contain and therefore forces a rebuild. A separating contact merely leaves behind a block that assembles to zero, which consumes some memory but does not alter the matrix, so the pattern is retained provided no more than ``stale_block_tolerance`` blocks have become stale.
 
-Internally the matrix is symmetric with only its upper triangle stored, and ``get_matrix()`` mirrors it into a full Eigen matrix, reusing the cached structure whenever the pattern survives. If your solver speaks block CSC natively, ``block_matrix()`` hands you MeshFEM's matrix and skips that conversion. Our header only forward declares the type, so include ``<MeshFEMSparse/BlockCSCHessian.hh>`` yourself when you want it.
+The assembled matrix is symmetric and stored with only its upper triangle; ``get_matrix()`` mirrors it into a full Eigen matrix, reusing the cached structure whenever the pattern is unchanged. Solvers that consume block CSC directly can instead call ``block_matrix()`` to obtain MeshFEM's representation and avoid the conversion. Because our header only forward declares that type, such callers must include ``<MeshFEMSparse/BlockCSCHessian.hh>`` themselves.
 
-When you reassemble without touching the collision set (a new stiffness, say, or the same Hessian under a different PSD projection) you can skip the comparison as well:
+When reassembling without modifying the collision set, for example under a different stiffness or PSD projection, the comparison itself can also be skipped:
 
 .. md-tab-set::
 
@@ -363,7 +363,7 @@ When you reassemble without touching the collision set (a new stiffness, say, or
             assembler.assume_unchanged_stencils = True
 
 .. warning::
-    ``assume_unchanged_stencils`` trusts you. If the stencils changed while their count stayed the same, assembly reads past the end of the pattern. Debug builds check the assumption; release builds do not.
+    ``assume_unchanged_stencils`` is an unchecked assertion. If the stencils did change while their count remained equal, assembly reads past the end of the pattern. Debug builds verify the assumption; release builds do not.
 
 .. note::
-    ``MeshFEMHessianAssembler`` requires the toolkit to be built with ``IPC_TOOLKIT_WITH_MESHFEM_SPARSE`` (on by default, and also what ``Potential::hessian`` uses internally). ``TripletHessianAssembler`` is always available and does what ``hessian`` has always done, gathering thread-local triplets and handing them to ``setFromTriplets``, but it has no pattern to carry over so reusing it saves nothing.
\ No newline at end of file
+    ``MeshFEMHessianAssembler`` requires the toolkit to be built with ``IPC_TOOLKIT_WITH_MESHFEM_SPARSE`` (enabled by default, and the backend ``Potential::hessian`` uses internally). ``TripletHessianAssembler`` is always available and reproduces the historical behavior of ``hessian``, accumulating thread-local triplets and merging them with ``setFromTriplets``. It maintains no sparsity pattern, so reusing an instance of it confers no benefit.
\ No newline at end of file

From f18782ff474af82acdd53cc988660fa2dadf2d01 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Tue, 4 Aug 2026 12:48:51 -0500
Subject: [PATCH 12/15] Update dependencies.rst

---
 docs/source/about/dependencies.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/source/about/dependencies.rst b/docs/source/about/dependencies.rst
index d738b1e3c..82b6162d2 100644
--- a/docs/source/about/dependencies.rst
+++ b/docs/source/about/dependencies.rst
@@ -121,7 +121,7 @@ Additionally, IPC Toolkit may optionally use the following libraries:
 Some of these libraries are enabled by default, and some are not. You can enable or disable them by passing the appropriate CMake option when you configure the IPC Toolkit build.
 
 .. note::
-    ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise), and it is currently pinned to a fork carrying fixes submitted upstream (`MeshFEM/MeshFEMCore#1 <https://github.com/MeshFEM/MeshFEMCore/pull/1>`_, `MeshFEM/MeshFEMSparse#1 <https://github.com/MeshFEM/MeshFEMSparse/pull/1>`_).
+    ``MeshFEMSparse`` (and its transitive dependency ``MeshFEMCore``) is downloaded source-only and compiled into a minimal static library (matrix data structures and assembly routines; no sparse direct solvers). When enabled (the default), :cpp:func:`ipc::Potential::hessian` assembles through the block-CSC backend — several times faster than the triplet-based assembly, with identical results up to floating-point summation order — and a :cpp:class:`ipc::MeshFEMHessianAssembler` held across :cpp:func:`ipc::Potential::assemble_hessian` calls additionally reuses the sparsity pattern between assemblies. It requires ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor`` (the default; the option is automatically disabled otherwise).
 
 .. warning::
     ``filib`` is licensed under `LGPL-2.1 <https://github.com/zfergus/filib/blob/main/LICENSE>`_ and as such it is required to be dynamically linked. Doing so automatically is a challenge, so by default we use static linkage. Enabling dynamic linkage requires copying the ``.so``/``.dylib``/``.dll`` file to the binary directory or system path. To enable this, set the CMake option ``FILIB_BUILD_SHARED_LIBS`` to ``ON`` and add this CMake code to copy the shared library object to the binary directory:

From 385fd62074028d2ed76811289ecabecc1ca8e36f Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Thu, 6 Aug 2026 14:21:04 -0500
Subject: [PATCH 13/15] Measure the triplet baseline explicitly in the
 breakdown table

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>
---
 tests/src/tests/potential/benchmark_assembly.cpp | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/tests/src/tests/potential/benchmark_assembly.cpp b/tests/src/tests/potential/benchmark_assembly.cpp
index 0d065ea09..7e0b244cb 100644
--- a/tests/src/tests/potential/benchmark_assembly.cpp
+++ b/tests/src/tests/potential/benchmark_assembly.cpp
@@ -25,6 +25,7 @@
 
 #include <ipc/candidates/candidates.hpp>
 #include <ipc/utils/eigen_ext.hpp>
+#include <ipc/utils/hessian_assembler.hpp>
 #ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
 #include <ipc/utils/meshfem_hessian_assembler.hpp>
 #endif
@@ -392,7 +393,7 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
     fmt::print(
         "{:<20} {:>9} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} {:>10} "
         "{:>10} {:>10} {:>8} {:>8} {:>8} {:>8}\n",
-        "scene", "#collis", "local(ms)", "total(ms)", "asm(ms)", "full(ms)",
+        "scene", "#collis", "local(ms)", "triplet(ms)", "asm(ms)", "full(ms)",
         "fold(ms)", "mfem(ms)", "mfblk(ms)", "mfr(ms)", "mfa(ms)", "local%",
         "asm%", "full%", "speedup");
 
@@ -414,8 +415,13 @@ TEST_CASE("Assembly cost breakdown", "[!benchmark][assembly]")
 
         const double t_local =
             median_seconds([&] { (void)local_hessians_only(scene); });
-        const double t_total = median_seconds(
-            [&] { (void)potential.hessian(collisions, mesh, X); });
+        // The triplet baseline, measured explicitly: hessian() routes through
+        // whichever backend is compiled in, so it is no longer the baseline.
+        const double t_total = median_seconds([&] {
+            TripletHessianAssembler assembler;
+            potential.assemble_hessian(collisions, mesh, X, assembler);
+            (void)assembler.get_matrix();
+        });
         const double t_full =
             median_seconds([&] { (void)mesh.to_full_dof(hess); });
         // Phase 1: fold to_full_dof into assembly.

From 5a833ca35a68b7544d0ca63c76c61d6f7de64ec2 Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Fri, 7 Aug 2026 01:37:20 -0500
Subject: [PATCH 14/15] Harden the assembly entry points

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.
---
 cmake/recipes/meshfem_sparse.cmake          |  6 ++--
 src/ipc/potentials/potential.cpp            | 36 +++++++++++++++++----
 src/ipc/utils/hessian_assembler.hpp         |  1 +
 src/ipc/utils/meshfem_hessian_assembler.cpp | 10 +++---
 4 files changed, 41 insertions(+), 12 deletions(-)

diff --git a/cmake/recipes/meshfem_sparse.cmake b/cmake/recipes/meshfem_sparse.cmake
index 7bd0b2220..e0fd321aa 100644
--- a/cmake/recipes/meshfem_sparse.cmake
+++ b/cmake/recipes/meshfem_sparse.cmake
@@ -88,9 +88,11 @@ target_compile_definitions(MeshFEMSparse PUBLIC
 
 # ipc_toolkit is compiled with EIGEN_DONT_VECTORIZE=1 when SIMD is enabled;
 # compiling the same Eigen templates with different vectorization settings is
-# an ODR violation with real alignment/layout consequences.
+# an ODR violation with real alignment/layout consequences. PUBLIC so the
+# setting travels with the target: anything including MeshFEMSparse's headers
+# must agree with how its own translation units were compiled.
 if(IPC_TOOLKIT_WITH_SIMD)
-    target_compile_definitions(MeshFEMSparse PRIVATE EIGEN_DONT_VECTORIZE=1)
+    target_compile_definitions(MeshFEMSparse PUBLIC EIGEN_DONT_VECTORIZE=1)
 endif()
 
 # Folder name for IDE
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index 43558028d..7efdcf667 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -135,7 +135,10 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
 
     if (collisions.empty()) {
         Eigen::VectorXd grad = Eigen::VectorXd::Zero(out_ndof);
-        return map_to_full ? mesh.to_full_dof(grad) : grad;
+        if (map_to_full) {
+            return mesh.to_full_dof(grad);
+        }
+        return grad;
     }
 
     const int dim = X.cols();
@@ -182,7 +185,10 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
         IPC_TOOLKIT_PROFILE_BLOCK("Gather Local Gradients");
         Eigen::VectorXd grad =
             gather_global_gradient(out_ndof, dim, local_grads, slot_vertex);
-        return map_to_full ? mesh.to_full_dof(grad) : grad;
+        if (map_to_full) {
+            return mesh.to_full_dof(grad);
+        }
+        return grad; // implicitly moved
     }
 
     tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(out_ndof));
@@ -214,7 +220,10 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
         Eigen::VectorXd combined_grad = grad.combine(
             [](const Eigen::VectorXd& a,
                const Eigen::VectorXd& b) -> Eigen::VectorXd { return a + b; });
-        return map_to_full ? mesh.to_full_dof(combined_grad) : combined_grad;
+        if (map_to_full) {
+            return mesh.to_full_dof(combined_grad);
+        }
+        return combined_grad; // implicitly moved
     }
 }
 
@@ -234,19 +243,34 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const bool fold_to_full = in_full_dof && mesh.is_selection_dof_map();
     const bool map_to_full = in_full_dof && !fold_to_full;
 
+    if (collisions.empty()) {
+        // Short-circuit: building a sparsity pattern for an empty contact set
+        // costs O(ndof) (or more) work to produce an all-zero matrix.
+        assert(X.rows() == mesh.num_vertices());
+        const int out_ndof = fold_to_full ? mesh.full_ndof() : X.size();
+        Eigen::SparseMatrix<double> hess(out_ndof, out_ndof);
+        if (map_to_full) {
+            return mesh.to_full_dof(hess);
+        }
+        return hess;
+    }
+
 #ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE
     MeshFEMHessianAssembler assembler;
     assemble_hessian(
         collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full);
-    const Eigen::SparseMatrix<double> hess = assembler.take_matrix();
+    Eigen::SparseMatrix<double> hess = assembler.take_matrix();
 #else
     TripletHessianAssembler assembler;
     assemble_hessian(
         collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full);
-    const Eigen::SparseMatrix<double> hess = assembler.get_matrix();
+    Eigen::SparseMatrix<double> hess = assembler.get_matrix();
 #endif
 
-    return map_to_full ? mesh.to_full_dof(hess) : hess;
+    if (map_to_full) {
+        return mesh.to_full_dof(hess);
+    }
+    return hess; // implicitly moved
 }
 
 template <class TCollisions>
diff --git a/src/ipc/utils/hessian_assembler.hpp b/src/ipc/utils/hessian_assembler.hpp
index 831fda478..e14cc0b9b 100644
--- a/src/ipc/utils/hessian_assembler.hpp
+++ b/src/ipc/utils/hessian_assembler.hpp
@@ -10,6 +10,7 @@
 
 #include <array>
 #include <functional>
+#include <memory>
 
 namespace ipc {
 
diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp
index e36020a2d..db9899fbb 100644
--- a/src/ipc/utils/meshfem_hessian_assembler.cpp
+++ b/src/ipc/utils/meshfem_hessian_assembler.cpp
@@ -331,6 +331,11 @@ void MeshFEMHessianAssembler::begin(
     const size_t num_stencils,
     const StencilGetter& stencil)
 {
+    // Validate before dividing by dim: dim == 0 would trap.
+    if (dim != 2 && dim != 3) {
+        log_and_throw_error(
+            "MeshFEMHessianAssembler: unsupported dimension {}!", dim);
+    }
     assert(ndof % dim == 0);
     const size_t num_block_vars = size_t(ndof) / dim;
 
@@ -338,11 +343,8 @@ void MeshFEMHessianAssembler::begin(
         || m_impl->num_block_vars() != num_block_vars) {
         if (dim == 2) {
             m_impl = std::make_unique<Impl<2>>(num_block_vars);
-        } else if (dim == 3) {
-            m_impl = std::make_unique<Impl<3>>(num_block_vars);
         } else {
-            log_and_throw_error(
-                "MeshFEMHessianAssembler: unsupported dimension {}!", dim);
+            m_impl = std::make_unique<Impl<3>>(num_block_vars);
         }
     }
 

From 45f198340e3e9e251ddb54dde07b33dc99603a9d Mon Sep 17 00:00:00 2001
From: Zachary Ferguson <zy.fergus@gmail.com>
Date: Thu, 6 Aug 2026 15:24:24 -0500
Subject: [PATCH 15/15] Mark MeshFEMSparse as PRIVATE

---
 CMakeLists.txt       | 3 +--
 tests/CMakeLists.txt | 5 +++++
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0f14e5922..a48845d7d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -274,8 +274,7 @@ endif()
 # Block-accelerated Hessian assembly
 if(IPC_TOOLKIT_WITH_MESHFEM_SPARSE)
   include(meshfem_sparse)
-  # PUBLIC: the MeshFEMHessianAssembler header includes MeshFEMSparse headers.
-  target_link_libraries(ipc_toolkit PUBLIC MeshFEM::Sparse)
+  target_link_libraries(ipc_toolkit PRIVATE MeshFEM::Sparse)
 endif()
 
 if(IPC_TOOLKIT_WITH_PROFILER)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b534803f8..58a5174b4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -65,6 +65,11 @@ target_link_libraries(ipc_toolkit_tests PRIVATE finitediff::finitediff)
 include(json)
 target_link_libraries(ipc_toolkit_tests PRIVATE nlohmann_json::nlohmann_json)
 
+if(IPC_TOOLKIT_WITH_MESHFEM_SPARSE)
+  include(meshfem_sparse)
+  target_link_libraries(ipc_toolkit_tests PRIVATE MeshFEM::Sparse)
+endif()
+
 if(IPC_TOOLKIT_WITH_CUDA)
   # We need to explicitly state that we need all CUDA files in the particle
   # library to be built with -dc as the member functions could be called by