Utils

Logger

spdlog::logger& ipc::logger();

Retrieves the current logger.

Returns:

A const reference to the logger object.

void ipc::set_logger(std::shared_ptr<spdlog::logger> logger);

Setup a logger object.

Calling this function with other function is not thread-safe.

Parameters:
std::shared_ptr<spdlog::logger> logger

New logger object to be used.

Positive Semi-Definite Projection

template <typename _Scalar, int _Rows, int _Cols, int _Options,
    int _MaxRows, int _MaxCols>
Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
ipc::project_to_psd(const Eigen::Matrix<_Scalar, _Rows, _Cols,
                        _Options, _MaxRows, _MaxCols>& A,
    const PSDProjectionMethod method = PSDProjectionMethod::CLAMP);

Matrix projection onto positive semi-definite cone.

Parameters:
const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> &A

Symmetric matrix to project

const PSDProjectionMethod method = PSDProjectionMethod::CLAMP

PSD projection method

Returns:

Projected matrix

template <typename _Scalar, int _Rows, int _Cols, int _Options,
    int _MaxRows, int _MaxCols>
Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
ipc::project_to_pd(const Eigen::Matrix<_Scalar, _Rows, _Cols,
                       _Options, _MaxRows, _MaxCols>& A,
    double eps = 1e-8);

Matrix projection onto positive definite cone.

Parameters:
const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> &A

Symmetric matrix to project

double eps = 1e-8

Minimum eigenvalue threshold

Returns:

Projected matrix

enum class ipc::PSDProjectionMethod : uint8_t;

Enumeration of implemented PSD projection methods.

Values:

enumerator NONE;

No PSD projection.

enumerator CLAMP;

Clamp negative eigenvalues to zero.

enumerator ABS;

Flip negative eigenvalues to positive.

Hessian Assembly

Pluggable backends for assembling per-collision Hessians into a global matrix (see ipc::Potential::assemble_hessian()).

class HessianAssembler;

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().

Subclassed by ipc::MeshFEMHessianAssembler, ipc::TripletHessianAssembler

Public Types

using StencilGetter = std::function<std::array<index_t, 4>(size_t)>;

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.

Public Functions

virtual ~HessianAssembler() = default;
virtual void begin(int ndof, int dim, size_t num_stencils,
    const StencilGetter& stencil)
    = 0;

Prepare for assembly.

Called once before any add_local_hessian.

Parameters:
int ndof

Number of global scalar DOF (rows == cols of the result).

int dim

Spatial dimension (rows/cols per vertex block).

size_t num_stencils

Number of local Hessians that will be added.

const StencilGetter &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 add_local_hessian(const MatrixMax12d& local_hess,
    const std::array<index_t, 4>& vertex_ids)
    = 0;

Add one local (stencil) Hessian to the global matrix.

Must be safe to call concurrently from multiple threads between begin() and end().

Parameters:
const MatrixMax12d &local_hess

Local Hessian of size (n·dim)×(n·dim), where n = local_hess.rows() / dim is the number of stencil vertices.

const std::array<index_t, 4> &vertex_ids

Global vertex IDs of the stencil; the first n entries are valid (remaining entries may be negative placeholders).

virtual void end() = 0;

Finish assembly. Called once after all add_local_hessian calls.

class TripletHessianAssembler : public ipc::HessianAssembler;

The default HessianAssembler: thread-local triplet caches merged into an Eigen::SparseMatrix via setFromTriplets.

This reproduces the historical behavior of Potential::hessian exactly.

Public Functions

TripletHessianAssembler() = default;
virtual void begin(int ndof, int dim, size_t num_stencils,
    const StencilGetter& stencil) override;

Note

stencil is unused: the triplet path needs no sparsity pattern.

virtual void add_local_hessian(const MatrixMax12d& local_hess,
    const std::array<index_t, 4>& vertex_ids) override;

Add one local (stencil) Hessian to the global matrix.

Must be safe to call concurrently from multiple threads between begin() and end().

Parameters:
const MatrixMax12d &local_hess

Local Hessian of size (n·dim)×(n·dim), where n = local_hess.rows() / dim is the number of stencil vertices.

const std::array<index_t, 4> &vertex_ids

Global vertex IDs of the stencil; the first n entries are valid (remaining entries may be negative placeholders).

inline virtual void end() override;

Finish assembly. Called once after all add_local_hessian calls.

Eigen::SparseMatrix<double> get_matrix();

Merge the thread-local caches and build the global matrix.

Call once, after end(); the internal caches are consumed.

Throws:

std::runtime_error – if called before the first assembly.

Private Members

int m_ndof = 0;
int m_dim = 0;
std::unique_ptr<
    tbb::enumerable_thread_specific<LocalThreadMatStorage>>
m_storage;

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 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.

class MeshFEMHessianAssembler : public ipc::HessianAssembler;

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.

Public Functions

MeshFEMHessianAssembler();
~MeshFEMHessianAssembler() override;
virtual void begin(int ndof, int dim, size_t num_stencils,
    const StencilGetter& stencil) override;

Prepare for assembly.

Called once before any add_local_hessian.

Parameters:
int ndof

Number of global scalar DOF (rows == cols of the result).

int dim

Spatial dimension (rows/cols per vertex block).

size_t num_stencils

Number of local Hessians that will be added.

const StencilGetter &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 add_local_hessian(const MatrixMax12d& local_hess,
    const std::array<index_t, 4>& vertex_ids) override;

Add one local (stencil) Hessian to the global matrix.

Must be safe to call concurrently from multiple threads between begin() and end().

Parameters:
const MatrixMax12d &local_hess

Local Hessian of size (n·dim)×(n·dim), where n = local_hess.rows() / dim is the number of stencil vertices.

const std::array<index_t, 4> &vertex_ids

Global vertex IDs of the stencil; the first n entries are valid (remaining entries may be negative placeholders).

inline virtual void end() override;

Finish assembly. Called once after all add_local_hessian calls.

const Eigen::SparseMatrix<double>& get_matrix() const;

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.

Throws:

std::runtime_error – if called before the first assembly.

Eigen::SparseMatrix<double> take_matrix();

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.

Throws:

std::runtime_error – if called before the first assembly.

const MeshFEM::BlockCSCHessianBase& block_matrix() const;

Access the assembled matrix in its native block-CSC form.

Skips the conversion performed by get_matrix(), which is worth doing if you can consume the block format directly (e.g., MeshFEM’s block SpMV or its Cholesky factorizers). Include <MeshFEMSparse/BlockCSCHessian.hh> to use the result.

The matrix is symmetric with only the upper triangle stored, and its block sparsity pattern covers the interacting vertex pairs of the last assembly (possibly with explicitly-zero stale blocks if the pattern was reused). The reference stays valid until the next begin() call, which may rebuild the pattern in place.

Note

A contact Hessian has no diagonal block for any vertex that is not in contact, so most block columns are empty. Reading operations handle that (trace() skips the missing diagonals; apply() and the sparsity pattern are unaffected), but the diagonal-mutating operations (addDiag(), setDiag()) have no entry to write and throw. Insert the missing diagonal blocks first if you need them.

Throws:

std::runtime_error – if called before the first assembly.

inline size_t stale_block_tolerance() const;

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.

inline void set_stale_block_tolerance(const size_t tolerance);

Set the number of vanished blocks tolerated before a rebuild.

inline bool reused_pattern() const;

Whether the last begin() call reused the cached pattern.

inline bool assume_unchanged_stencils() const;

Whether begin() may skip change detection entirely.

inline void set_assume_unchanged_stencils(const bool assume);

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.

Private Members

std::unique_ptr<ImplBase> m_impl;
size_t m_stale_block_tolerance = 0;
bool m_assume_unchanged_stencils = false;
bool m_reused_pattern = false;
template <int dim>
struct Impl : public ipc::MeshFEMHessianAssembler::ImplBase;

Dimension-specific implementation (dim ∈ {2, 3}), mirroring MeshFEM’s own IPC integration (IPCWrapper.cc): one dim-sized block variable per vertex.

Public Types

using VarStructure = MeshFEM::OptimizationVarStructure<dim>;
using Stencil = MeshFEM::ElementBlockVarsWithSizeRange<1, 4>;

Local (block) variables of a collision stencil: 1–4 vertices.

Public Functions

inline explicit Impl(const size_t num_block_vars);
inline virtual int dimension() const override;
inline virtual size_t num_block_vars() const override;
inline virtual bool update_pattern(const size_t num_stencils,
    const StencilGetter& stencil,
    const size_t stale_block_tolerance,
    const bool assume_unchanged) override;

Reuse the cached pattern if the stencils still fit (returns true) or rebuild it (returns false).

Zeroes the values either way.

inline virtual void add(const MatrixMax12d& local_hess,
    const std::array<index_t, 4>& vertex_ids) override;
inline virtual const Eigen::SparseMatrix<double>&
to_eigen() const override;
inline virtual Eigen::SparseMatrix<double> take_eigen() override;
inline virtual const MeshFEM::BlockCSCHessianBase&
block_matrix() const override;

Private Functions

inline void build_eigen_structure() const;
inline void fill_eigen_values() const;

Private Members

MeshFEM::SystemAssembler<dim> m_assembler;
VarStructure m_vars;
std::unique_ptr<MeshFEM::BlockCSCHessian<VarStructure>> m_H;
MeshFEM::VarLocks m_locks;
size_t m_num_stencils = 0;

Stencil count of the cached pattern (assume-unchanged sanity check).

mutable bool m_eigen_structure_valid = false;
mutable std::vector<std::ptrdiff_t> m_sym_col_start;
mutable std::vector<BlockEntry> m_sym_entries;
mutable Eigen::SparseMatrix<double> m_M;

Private Static Functions

static inline bool has_distinct_ids(
    const std::array<index_t, 4>& vertex_ids);

Whether the valid (non-negative) entries of vertex_ids are pairwise distinct, i.e., whether add() may take the column-merge fast path.

static inline Stencil to_stencil(
    const std::array<index_t, 4>& vertex_ids);

Compact the (possibly -1-padded) vertex ID array into a MeshFEM stencil of block variables.

Note

Repeated IDs are kept as-is: the stencil’s size must match the local Hessian’s block count, and both scatter routines are indexed by local block. add() picks the routine that handles duplicates.

Private Static Attributes

static constexpr int N = dim;
struct BlockEntry;

A block of the full (symmetrized) matrix and the stored block backing it.

Public Members

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.

struct ImplBase;

Subclassed by ipc::MeshFEMHessianAssembler::Impl< dim >

Public Functions

virtual ~ImplBase() = default;
virtual int dimension() const = 0;
virtual size_t num_block_vars() const = 0;
virtual bool update_pattern(size_t num_stencils,
    const StencilGetter& stencil, size_t stale_block_tolerance,
    bool assume_unchanged)
    = 0;

Reuse the cached pattern if the stencils still fit (returns true) or rebuild it (returns false).

Zeroes the values either way.

virtual void add(const MatrixMax12d& local_hess,
    const std::array<index_t, 4>& vertex_ids)
    = 0;
virtual const Eigen::SparseMatrix<double>& to_eigen() const = 0;
virtual Eigen::SparseMatrix<double> take_eigen() = 0;
virtual const MeshFEM::BlockCSCHessianBase& block_matrix() const
    = 0;

Eigen Extensions

using ArrayXb = Eigen::Array<bool, Eigen::Dynamic, 1>;

An array of boolean scalars.

using VectorXb = Eigen::Matrix<bool, Eigen::Dynamic, 1>;

A Vector of boolean scalars.

using Vector3b = Eigen::Matrix<bool, 3, 1>;

A Vector of boolean scalars with a fixed size of 3x1.

using MatrixXb
    = Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>;

A dynamic size matrix of boolean scalars.

template <typename T, int max_dim>
using VectorMax = Eigen::Matrix<T, Eigen::Dynamic, 1,
    Eigen::ColMajor, max_dim, 1>;

A dynamic size vector with a fixed maximum size.

Template Parameters:
typename T

The type of the vector elements.

int max_dim

The maximum size of the vector.

template <typename T, int max_dim>
using RowVectorMax = Eigen::Matrix<T, 1, Eigen::Dynamic,
    Eigen::RowMajor, 1, max_dim>;

A dynamic size row vector with a fixed maximum size.

Template Parameters:
typename T

The type of the vector elements.

int max_dim

The maximum size of the vector.

using Vector1f = Eigen::Vector<float, 1>;

A static size matrix of size of 1×1.

using Vector1d = Eigen::Vector<double, 1>;

A static size matrix of size of 1×1.

using Vector6f = Eigen::Vector<float, 6>;

A static size matrix of size of 6×1.

using Vector6d = Eigen::Vector<double, 6>;

A static size matrix of size of 6×1.

using Vector9f = Eigen::Vector<float, 9>;

A static size matrix of size of 9×1.

using Vector9d = Eigen::Vector<double, 9>;

A static size matrix of size of 9×1.

using Vector12f = Eigen::Vector<float, 12>;

A static size matrix of size of 12×1.

using Vector12d = Eigen::Vector<double, 12>;

A static size matrix of size of 12×1.

using Vector15f = Eigen::Vector<float, 15>;

A static size matrix of size of 15×1.

using Vector15d = Eigen::Vector<double, 15>;

A static size matrix of size of 15×1.

using Matrix6f = Eigen::Matrix<float, 6, 6>;

A static size matrix of size of 6×6.

using Matrix6d = Eigen::Matrix<double, 6, 6>;

A static size matrix of size of 6×6.

using Matrix9f = Eigen::Matrix<float, 9, 9>;

A static size matrix of size of 9×9.

using Matrix9d = Eigen::Matrix<double, 9, 9>;

A static size matrix of size of 9×9.

using Matrix12f = Eigen::Matrix<float, 12, 12>;

A static size matrix of size of 12×12.

using Matrix12d = Eigen::Matrix<double, 12, 12>;

A static size matrix of size of 12×12.

using Matrix15f = Eigen::Matrix<float, 15, 15>;

A static size matrix of size of 15×15.

using Matrix15d = Eigen::Matrix<double, 15, 15>;

A static size matrix of size of 15×15.

template <typename T> using VectorMax2 = VectorMax<T, 2>;

A dynamic size matrix with a fixed maximum size of 3×1.

template <typename T> using VectorMax3 = VectorMax<T, 3>;

A dynamic size matrix with a fixed maximum size of 3×1.

template <typename T> using VectorMax4 = VectorMax<T, 4>;

A dynamic size matrix with a fixed maximum size of 4×1.

template <typename T> using VectorMax6 = VectorMax<T, 6>;

A dynamic size matrix with a fixed maximum size of 6×1.

template <typename T> using VectorMax9 = VectorMax<T, 9>;

A dynamic size matrix with a fixed maximum size of 9×1.

template <typename T> using VectorMax12 = VectorMax<T, 12>;

A dynamic size matrix with a fixed maximum size of 12×1.

using VectorMax2f = VectorMax2<float>;

A dynamic size matrix with a fixed maximum size of 2×1.

using VectorMax2d = VectorMax2<double>;

A dynamic size matrix with a fixed maximum size of 2×1.

using VectorMax3f = VectorMax3<float>;

A dynamic size matrix with a fixed maximum size of 3×1.

using VectorMax3d = VectorMax3<double>;

A dynamic size matrix with a fixed maximum size of 3×1.

using VectorMax3i = VectorMax3<int>;

A dynamic size matrix with a fixed maximum size of 3×1.

using VectorMax4f = VectorMax4<float>;

A dynamic size matrix with a fixed maximum size of 4×1.

using VectorMax4d = VectorMax4<double>;

A dynamic size matrix with a fixed maximum size of 4×1.

using VectorMax4i = VectorMax4<int>;

A dynamic size matrix with a fixed maximum size of 4×1.

using VectorMax6f = VectorMax6<float>;

A dynamic size matrix with a fixed maximum size of 6×1.

using VectorMax6d = VectorMax6<double>;

A dynamic size matrix with a fixed maximum size of 6×1.

using VectorMax6b = VectorMax6<bool>;

A dynamic size matrix with a fixed maximum size of 6×1.

using VectorMax9f = VectorMax9<float>;

A dynamic size matrix with a fixed maximum size of 9×1.

using VectorMax9d = VectorMax9<double>;

A dynamic size matrix with a fixed maximum size of 9×1.

using VectorMax12f = VectorMax12<float>;

A dynamic size matrix with a fixed maximum size of 12×1.

using VectorMax12d = VectorMax12<double>;

A dynamic size matrix with a fixed maximum size of 12×1.

template <typename T> using RowVectorMax2 = RowVectorMax<T, 2>;

A dynamic size matrix with a fixed maximum size of 1×2.

template <typename T> using RowVectorMax3 = RowVectorMax<T, 3>;

A dynamic size matrix with a fixed maximum size of 1×3.

using RowVectorMax2d = RowVectorMax2<double>;

A dynamic size matrix with a fixed maximum size of 1×2.

using RowVectorMax3d = RowVectorMax3<double>;

A dynamic size matrix with a fixed maximum size of 1×3.

using RowVectorMax6d = RowVectorMax<double, 6>;

A dynamic size matrix with a fixed maximum size of 6×1.

using RowVectorMax9d = RowVectorMax<double, 9>;

A dynamic size matrix with a fixed maximum size of 9×1.

using RowVectorMax12d = RowVectorMax<double, 12>;

A dynamic size matrix with a fixed maximum size of 12×1.

template <typename T, int max_rows, int max_cols>
using MatrixMax = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
    Eigen::ColMajor, max_rows, max_cols>;
template <typename T> using MatrixMax2 = MatrixMax<T, 2, 2>;

A dynamic size matrix with a fixed maximum size of 3×3.

template <typename T> using MatrixMax3 = MatrixMax<T, 3, 3>;

A dynamic size matrix with a fixed maximum size of 3×3.

template <typename T> using MatrixMax6 = MatrixMax<T, 6, 6>;

A dynamic size matrix with a fixed maximum size of 6×6.

template <typename T> using MatrixMax9 = MatrixMax<T, 9, 9>;

A dynamic size matrix with a fixed maximum size of 9×9.

template <typename T> using MatrixMax12 = MatrixMax<T, 12, 12>;

A dynamic size matrix with a fixed maximum size of 12×12.

using MatrixMax2f = MatrixMax2<float>;

A dynamic size matrix with a fixed maximum size of 2×2.

using MatrixMax2d = MatrixMax2<double>;

A dynamic size matrix with a fixed maximum size of 2×2.

using MatrixMax3f = MatrixMax3<float>;

A dynamic size matrix with a fixed maximum size of 3×3.

using MatrixMax3d = MatrixMax3<double>;

A dynamic size matrix with a fixed maximum size of 3×3.

using MatrixMax6f = MatrixMax6<float>;

A dynamic size matrix with a fixed maximum size of 6×6.

using MatrixMax6d = MatrixMax6<double>;

A dynamic size matrix with a fixed maximum size of 6×6.

using MatrixMax9f = MatrixMax9<float>;

A dynamic size matrix with a fixed maximum size of 9×9.

using MatrixMax9d = MatrixMax9<double>;

A dynamic size matrix with a fixed maximum size of 9×9.

using MatrixMax12f = MatrixMax12<float>;

A dynamic size matrix with a fixed maximum size of 12×12.

using MatrixMax12d = MatrixMax12<double>;

A dynamic size matrix with a fixed maximum size of 12×12.

using DiagonalMatrixXd
    = Eigen::DiagonalMatrix<double, Eigen::Dynamic>;

A dynamic size diagonal matrix.

using DiagonalMatrixMax6d
    = Eigen::DiagonalMatrix<double, Eigen::Dynamic, 6>;

A dynamic size diagonal matrix with a fixed maximum size of 6×6.

template <typename T>
using ArrayMax2
    = Eigen::Array<T, Eigen::Dynamic, 1, Eigen::ColMajor, 2, 1>;

A dynamic size array with a fixed maximum size of 2×1.

template <typename T>
using ArrayMax3
    = Eigen::Array<T, Eigen::Dynamic, 1, Eigen::ColMajor, 3, 1>;

A dynamic size array with a fixed maximum size of 2×1.

template <typename T>
using ArrayMax4
    = Eigen::Array<T, Eigen::Dynamic, 1, Eigen::ColMajor, 4, 1>;

A dynamic size array with a fixed maximum size of 4×1.

using ArrayMax2d = ArrayMax2<double>;

A dynamic size array with a fixed maximum size of 2×1.

using ArrayMax3d = ArrayMax3<double>;

A dynamic size array with a fixed maximum size of 3×1.

using ArrayMax3i = ArrayMax3<int>;

A dynamic size array with a fixed maximum size of 3×1.

using ArrayMax4d = ArrayMax4<double>;

A dynamic size array with a fixed maximum size of 4×1.

using ArrayMax4i = ArrayMax4<int>;

A dynamic size array with a fixed maximum size of 4×1.

template <int dim>
using GradientType = std::tuple<double, Eigen::Vector<double, dim>>;
template <int dim>
using HessianType = std::tuple<double, Eigen::Vector<double, dim>,
    Eigen::Matrix<double, dim, dim>>;