ipc-sim · zfergus · Jul 30, 2026 · Jul 30, 2026 · Jul 30, 2026 · Jul 30, 2026
diff --git 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" ON)
 
 # Advanced options
 option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
@@ -151,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
@@ -263,6 +271,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
@@ -0,0 +1,97 @@
+# 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)
+
+# 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)
+CPMAddPackage(
+    NAME MeshFEMCore
+    URL "https://github.com/zfergus/MeshFEMCore/archive/8d0e84788189748d9e906cc7f807507a3cb4b2ef.zip"
+    URL_HASH SHA256=71fe52e49276a401ae64d692dc26aea4147b1baa636bc78082d3fd0061eb4ccb
+    DOWNLOAD_ONLY YES
+)
+CPMAddPackage(
+    NAME MeshFEMSparse
+    URL "https://github.com/zfergus/MeshFEMSparse/archive/15a92834189ade69ed9e00373adc3036a72cd40a.zip"
+    URL_HASH SHA256=882f3a126f59023b4d4aba6e1b96b89d5f82f2e80bfe086c1f3afd435f8e4c80
+    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
+)
+
+# 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/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
@@ -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).
+
 .. 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
@@ -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/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 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::
+
+    .. 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 agree with ``collision_mesh.to_full_dof(...)`` up to the order in which the local contributions are summed.
+
+.. note::
+    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
 ^^^^^^^^^^^^^^^^^^^^^^
 
@@ -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`.
+        - ``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`` 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`` accepts the assembler as a parameter, allowing a single instance to persist across the solve and retain its sparsity pattern:
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // A single assembler for the entire solve, rather than 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 to retain 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
+
+            # A single assembler for the entire solve, rather than 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 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.
+
+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 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::
+
+    .. 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`` 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`` (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.
diff --git 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/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,