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 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/zfergus/MeshFEMCore/archive/8d0e84788189748d9e906cc7f807507a3cb4b2ef.zip"
+    URL_HASH SHA256=71fe52e49276a401ae64d692dc26aea4147b1baa636bc78082d3fd0061eb4ccb
+    DOWNLOAD_ONLY YES
+)
+CPMAddPackage(
+    NAME MeshFEMSparse
+    URL "https://github.com/zfergus/MeshFEMSparse/archive/ade01a1775f1911c0bf373190001923f452ad6cc.zip"
+    URL_HASH SHA256=6545cb0aa7b8513a370dc8561ea5abdd5ec484d03ded7e64f48e5dbfc94c49b0
+    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), 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
@@ -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/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
@@ -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
@@ -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
@@ -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/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 {