diff --git a/CMakeLists.txt b/CMakeLists.txt index 6b9f6329a..0f14e5922 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" 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 new file mode 100644 index 000000000..7bd0b2220 --- /dev/null +++ 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 . 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 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..82b6162d2 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 `_ - |: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 `_ + - |: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 `_ 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 `_'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 index 7faf069be..f2758966c 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 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 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`. \ 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`` 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& 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 ```` 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. \ No newline at end of file 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/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..04c4886b3 100644 --- a/python/src/potentials/potential.hpp +++ b/python/src/potentials/potential.hpp @@ -1,6 +1,7 @@ #include #include +#include using namespace ipc; @@ -36,7 +37,7 @@ void define_potential_methods(PyClass& potential) "gradient", py::overload_cast< const TCollisions&, const CollisionMesh&, - Eigen::ConstRef>( + Eigen::ConstRef, const bool>( &Potential::gradient, py::const_), R"ipc_Qu8mg5v7( Compute the gradient of the potential. @@ -45,17 +46,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, const PSDProjectionMethod>( - &Potential::hessian, py::const_), + Eigen::ConstRef, const PSDProjectionMethod, + const bool>(&Potential::hessian, py::const_), R"ipc_Qu8mg5v7( Compute the hessian of the potential. @@ -64,12 +66,33 @@ 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( + "assemble_hessian", &Potential::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>( 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 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 + +#include +#include + +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_( + 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_( + 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_( + 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/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 to_full_dof(const Eigen::SparseMatrix& 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 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> m_vertex_vertex_adjacencies; 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 d4c08feae..43558028d 100644 --- a/src/ipc/potentials/potential.cpp +++ b/src/ipc/potentials/potential.cpp @@ -2,35 +2,94 @@ #include #include +#include #include +#include +#include #include #include #include -#include #include -#include #include +#include +#include + namespace ipc { namespace { - void set_triplets( - const Eigen::SparseMatrix& M, - std::vector>& triplets, - const size_t start_index) + + /// @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 local_grads, + const std::vector& slot_vertex) { - using InnerIterator = Eigen::SparseMatrix::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(it.row(), it.col(), it.value()); + 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> 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 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 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 @@ -60,37 +119,102 @@ template Eigen::VectorXd Potential::gradient( const TCollisions& collisions, const CollisionMesh& mesh, - Eigen::ConstRef X) const + Eigen::ConstRef X, + const bool in_full_dof) const { assert(X.rows() == mesh.num_vertices()); IPC_TOOLKIT_PROFILE_BLOCK("Potential::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(); + 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 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 grad(Eigen::VectorXd::Zero(X.size())); + tbb::combinable 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( 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,32 +223,80 @@ Eigen::SparseMatrix Potential::hessian( const TCollisions& collisions, const CollisionMesh& mesh, Eigen::ConstRef 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::hessian()"); - if (collisions.empty()) { - return Eigen::SparseMatrix(X.size(), X.size()); + // 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; + +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE + MeshFEMHessianAssembler assembler; + assemble_hessian( + collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full); + const Eigen::SparseMatrix hess = assembler.take_matrix(); +#else + TripletHessianAssembler assembler; + assemble_hessian( + collisions, mesh, X, assembler, project_hessian_to_psd, fold_to_full); + const Eigen::SparseMatrix hess = assembler.get_matrix(); +#endif + + return map_to_full ? mesh.to_full_dof(hess) : hess; +} + +template +void Potential::assemble_hessian( + const TCollisions& collisions, + const CollisionMesh& mesh, + Eigen::ConstRef 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::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); + + 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 = X.size(); - - constexpr int MAX_TRIPLETS_SIZE = 10'000'000; - const int buffer_size = std::min(MAX_TRIPLETS_SIZE, ndof); + const int ndof = in_full_dof ? mesh.full_ndof() : X.size(); // NOLINT + + // 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; + }; - tbb::enumerable_thread_specific storage( - LocalThreadMatStorage(buffer_size, ndof, ndof)); + assembler.begin(ndof, dim, collisions.size(), stencil_vertex_ids); { - 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; @@ -135,90 +307,11 @@ Eigen::SparseMatrix Potential::hessian( project_hessian_to_psd); } - { - IPC_TOOLKIT_PROFILE_BLOCK( - "Map Local Hessian to Global Triplets"); - local_hessian_to_global_triplets( - local_hess, collision.vertex_ids(edges, faces), dim, - *(hess_triplets.cache), mesh.num_vertices()); - } + assembler.add_local_hessian(local_hess, stencil_vertex_ids(i)); }); } - if (storage.empty()) { - return Eigen::SparseMatrix(); - } - - // 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 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> triplets; - - Eigen::SparseMatrix 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 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( - *((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; + assembler.end(); } template class Potential; diff --git a/src/ipc/potentials/potential.hpp b/src/ipc/potentials/potential.hpp index 2c6c42785..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 Potential { @@ -40,24 +42,50 @@ template 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 X) const; + Eigen::ConstRef 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 hessian( const TCollisions& collisions, const CollisionMesh& mesh, Eigen::ConstRef X, const PSDProjectionMethod project_hessian_to_psd = - PSDProjectionMethod::NONE) const; + 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 X, + HessianAssembler& assembler, + const PSDProjectionMethod project_hessian_to_psd = + PSDProjectionMethod::NONE, + const bool in_full_dof = false) const; // -- Single collision methods --------------------------------------------- diff --git a/src/ipc/utils/CMakeLists.txt b/src/ipc/utils/CMakeLists.txt index f0086228e..25a2368f6 100644 --- a/src/ipc/utils/CMakeLists.txt +++ b/src/ipc/utils/CMakeLists.txt @@ -3,11 +3,15 @@ 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 matrix_cache.cpp matrix_cache.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 new file mode 100644 index 000000000..f8e142f16 --- /dev/null +++ b/src/ipc/utils/hessian_assembler.cpp @@ -0,0 +1,143 @@ +#include "hessian_assembler.hpp" + +#include +#include + +#include +#include + +#include +#include + +namespace ipc { + +namespace { + void set_triplets( + const Eigen::SparseMatrix& M, + std::vector>& triplets, + const size_t start_index) + { + using InnerIterator = Eigen::SparseMatrix::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(it.row(), it.col(), it.value()); + } + } + } +} // namespace + +void TripletHessianAssembler::begin( + const int ndof, + const int dim, + const size_t num_stencils, + const StencilGetter& /*stencil*/) +{ + 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(buffer_size, ndof, ndof)); +} + +void TripletHessianAssembler::add_local_hessian( + const MatrixMax12d& local_hess, const std::array& 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 TripletHessianAssembler::get_matrix() +{ + assert(m_storage != nullptr); + tbb::enumerable_thread_specific& storage = + *m_storage; + + Eigen::SparseMatrix 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 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> 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( + *((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..831fda478 --- /dev/null +++ b/src/ipc/utils/hessian_assembler.hpp @@ -0,0 +1,97 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +#include +#include + +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: + /// @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(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. + /// @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. + /// + /// 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& 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; + + /// @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, + const std::array& 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 get_matrix(); + +private: + int m_ndof = 0; + int m_dim = 0; + std::unique_ptr> + m_storage; +}; + +} // namespace ipc diff --git a/src/ipc/utils/meshfem_hessian_assembler.cpp b/src/ipc/utils/meshfem_hessian_assembler.cpp new file mode 100644 index 000000000..e36020a2d --- /dev/null +++ b/src/ipc/utils/meshfem_hessian_assembler.cpp @@ -0,0 +1,386 @@ +#include "meshfem_hessian_assembler.hpp" + +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE + +#include +#include + +#include +#include + +#include +#include +#include +#include + +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& vertex_ids) = 0; + + virtual const Eigen::SparseMatrix& to_eigen() const = 0; + virtual Eigen::SparseMatrix take_eigen() = 0; + virtual const MeshFEM::BlockCSCHessianBase& block_matrix() const = 0; +}; + +/// Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM's own +/// IPC integration (IPCWrapper.cc): one dim-sized block variable per vertex. +template +struct MeshFEMHessianAssembler::Impl final + : public MeshFEMHessianAssembler::ImplBase { + using VarStructure = MeshFEM::OptimizationVarStructure; + /// Local (block) variables of a collision stencil: 1–4 vertices. + using Stencil = MeshFEM::ElementBlockVarsWithSizeRange<1, 4>; + + 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::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::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, block_vars_of_stencil); + } + m_num_stencils = num_stencils; + m_H->setZero(); // Allocate (and zero) the value array. + m_eigen_structure_valid = false; + return false; + } + + void + add(const MatrixMax12d& local_hess, + const std::array& vertex_ids) override + { + assert(m_H != nullptr); + 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(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( + 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. 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: 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)), + // column-major within the block; diagonal blocks are full N×N with a + // zeroed strict lower triangle. + const Eigen::SparseMatrix& to_eigen() const override + { + IPC_TOOLKIT_PROFILE_BLOCK("MeshFEM block CSC to Eigen"); + assert(m_H != nullptr); + + if (!m_eigen_structure_valid) { + build_eigen_structure(); + m_eigen_structure_valid = true; + } + fill_eigen_values(); + return m_M; + } + + Eigen::SparseMatrix take_eigen() override + { + to_eigen(); + m_eigen_structure_valid = false; // m_M is about to be gutted + 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; + + /// 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& 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 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. + 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++) { + m_sym_col_start[bj + 1]++; + if (Ai[ii] != bj) { + m_sym_col_start[Ai[ii] + 1]++; + } + } + } + std::partial_sum( + 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. + m_sym_entries.resize(num_sym_blocks); + { + std::vector cursor( + 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]; + m_sym_entries[cursor[bj]++] = { bi, ii, false }; + if (bi != bj) { + m_sym_entries[cursor[bi]++] = { bj, ii, true }; + } + } + } + } + + // --- Scalar CSC index arrays ------------------------------------ + m_M.resize(n, n); + m_M.makeCompressed(); + m_M.resizeNonZeros(N * N * num_sym_blocks); + + 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 * (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 * m_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 = 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++) { + if (e.transposed || (diagonal && r > c)) { + values[out] = block[r * N + c]; // transposed + } else { + values[out] = block[c * N + r]; + } + out++; + } + } + } + }); + } + + MeshFEM::SystemAssembler m_assembler; + VarStructure m_vars; + std::unique_ptr> + 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; + + // Cached Eigen conversion structure; valid while the pattern is reused. + mutable bool m_eigen_structure_valid = false; + mutable std::vector m_sym_col_start; + mutable std::vector m_sym_entries; + mutable Eigen::SparseMatrix + m_M; // NOLINT(readability-identifier-naming): M = matrix +}; + +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 (m_impl == nullptr || m_impl->dimension() != dim + || m_impl->num_block_vars() != num_block_vars) { + if (dim == 2) { + m_impl = std::make_unique>(num_block_vars); + } else if (dim == 3) { + m_impl = std::make_unique>(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( + const MatrixMax12d& local_hess, const std::array& vertex_ids) +{ + assert(m_impl != nullptr); + m_impl->add(local_hess, vertex_ids); +} + +const Eigen::SparseMatrix& MeshFEMHessianAssembler::get_matrix() const +{ + assert(m_impl != nullptr); + return m_impl->to_eigen(); +} + +Eigen::SparseMatrix MeshFEMHessianAssembler::take_matrix() +{ + assert(m_impl != nullptr); + 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 new file mode 100644 index 000000000..e5d7913bb --- /dev/null +++ b/src/ipc/utils/meshfem_hessian_assembler.hpp @@ -0,0 +1,149 @@ +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE + +#include + +#include + +namespace MeshFEM { +/// Forward declaration so block_matrix() can be exposed without pulling +/// MeshFEMSparse's headers into this one. Include +/// to use the returned object. +struct BlockCSCHessianBase; +} // namespace MeshFEM + +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 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, +/// reusing the cached structure when the pattern is unchanged. +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& 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. + /// + /// The reference stays valid (and its values current) until the next + /// begin() call; copy it to keep a snapshot. + const Eigen::SparseMatrix& 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 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 + /// 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, + /// 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 struct Impl; + std::unique_ptr m_impl; + + size_t m_stale_block_tolerance = 0; + bool m_assume_unchanged_stencils = false; + bool m_reused_pattern = false; +}; + +} // 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 f2ff672f6..2f4b3589f 100644 --- a/tests/src/tests/potential/CMakeLists.txt +++ b/tests/src/tests/potential/CMakeLists.txt @@ -5,10 +5,15 @@ set(SOURCES test_smooth_potential.cpp test_friction_potential.cpp test_distance_vector_methods.cpp + test_full_dof_assembly.cpp + test_meshfem_assembly.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 + +#include + +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 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(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 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& 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 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 +#include +#include + +#include +#include +#include + +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 +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& 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..0d065ea09 --- /dev/null +++ b/tests/src/tests/potential/benchmark_assembly.cpp @@ -0,0 +1,489 @@ +// 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 + +#include +#include +#include +#include + +#include +#include +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE +#include +#endif + +#include +#include + +#include + +#include +#include +#include +#include +#include + +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(0), collisions.size()), 0.0, + [&](const tbb::blocked_range& 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()); +} + +/// @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(0), collisions.size()), 0.0, + [&](const tbb::blocked_range& 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()); +} + +/// @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 double median_seconds(F&& f, const int num_samples = 5) +{ + std::vector 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(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 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 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 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)); + }; + + // 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); + }; + +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE + // 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(); + }; + + // 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, cold, block)", 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); + }; + + // 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( + 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 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); + }; + + // 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]") +{ + // 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} {:>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)", "mfr(ms)", "mfa(ms)", "local%", + "asm%", "full%", "speedup"); + + for (const auto& spec : ipc::tests::assembly_scene_specs()) { + const std::optional 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 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); }); + // 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); + }); + +#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); + }); + // 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::quiet_NaN(); + const double t_meshfem_blk = std::numeric_limits::quiet_NaN(); + const double t_meshfem_reused = + std::numeric_limits::quiet_NaN(); + const double t_meshfem_assumed = + std::numeric_limits::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; + + constexpr double MS = 1e3; + fmt::print( + "{:<20} {:>9} {:>10.3f} {:>10.3f} {:>10.3f} {:>10.3f} {:>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, t_folded * MS, t_meshfem * MS, + 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_full_dof_assembly.cpp b/tests/src/tests/potential/test_full_dof_assembly.cpp new file mode 100644 index 000000000..655f0c582 --- /dev/null +++ b/tests/src/tests/potential/test_full_dof_assembly.cpp @@ -0,0 +1,157 @@ +// 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 +#include + +#include +#include + +#include + +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 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()); + // 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") + { + const PSDProjectionMethod psd = + GENERATE(PSDProjectionMethod::NONE, PSDProjectionMethod::CLAMP); + CAPTURE(psd); + + const Eigen::SparseMatrix hess_folded = + potential.hessian(collisions, mesh, X, psd, /*in_full_dof=*/true); + const Eigen::SparseMatrix 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 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 hess_full = potential.hessian( + collisions, mesh, vertices, PSDProjectionMethod::NONE, + /*in_full_dof=*/true); + const Eigen::SparseMatrix 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 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); +} 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..ef14dd204 --- /dev/null +++ b/tests/src/tests/potential/test_meshfem_assembly.cpp @@ -0,0 +1,239 @@ +// 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 + +#ifdef IPC_TOOLKIT_WITH_MESHFEM_SPARSE + +#include "assembly_scene.hpp" + +#include +#include + +#include +#include +#include + +#include +#include + +// 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 + +#include + +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 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 expected = triplet_assembler.get_matrix(); + + MeshFEMHessianAssembler meshfem_assembler; + potential.assemble_hessian( + collisions, mesh, X, meshfem_assembler, psd, in_full_dof); + const Eigen::SparseMatrix 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 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& actual, + const Eigen::SparseMatrix& 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::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::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 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 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 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(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]") +{ + 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 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