Configurable derivative layout

CMakeLists.txt

@@ -85,9 +85,12 @@ option(IPC_TOOLKIT_WITH_PROFILER              "Enable performance profiler"
 
 # Advanced options
 option(IPC_TOOLKIT_WITH_CODE_COVERAGE         "Enable coverage reporting"                     OFF)
-
 mark_as_advanced(IPC_TOOLKIT_WITH_CODE_COVERAGE) # This is used in GitHub Actions
 
+set(IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT "RowMajor" CACHE STRING "Layout of the vertex derivatives")
+set_property(CACHE IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT PROPERTY STRINGS "ColMajor" "RowMajor")
+mark_as_advanced(IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT) # This is an advanced option that most users won't care about.
+
 # Set default minimum C++ standard
 if(IPC_TOOLKIT_TOPLEVEL_PROJECT)
   set(CMAKE_CXX_STANDARD 17)

docs/source/tutorials/getting_started.rst

@@ -242,6 +242,10 @@ and the Hessian of size :math:`|V|d \times |V|d` with the order
     \frac{\partial^2 B}{\partial z_n^2}
     \end{bmatrix}.
 
+.. note::
+
+    This row-major derivative ordering is the default. If your simulation framework uses a column-major DOF layout, the IPC Toolkit can be configured to use that instead. See :ref:`Vertex Derivative Layout <tutorials/misc:Vertex Derivative Layout>` for details.
+
 Adaptive Barrier Stiffness
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 

docs/source/tutorials/misc.rst

@@ -142,4 +142,45 @@ You can also get the current maximum number of threads as follows:
 
         .. code-block:: python
 
-            nthreads = ipctk.get_num_threads()
+            nthreads = ipctk.get_num_threads()
+
+Vertex Derivative Layout
+------------------------
+
+By default, the IPC Toolkit orders derivative vectors (gradients and Hessians) with respect to vertex DOFs in a **row-major** layout:
+
+.. math::
+
+    \mathbf{x} = [x_0, y_0, z_0, x_1, y_1, z_1, \ldots]
+
+Some simulation frameworks instead use a **column-major** layout, where all coordinates of the same dimension are grouped together:
+
+.. math::
+
+    \mathbf{x} = [x_0, x_1, \ldots, y_0, y_1, \ldots, z_0, z_1, \ldots]
+
+The IPC Toolkit provides a compile-time CMake option to select which layout is used for all gradient and Hessian assembly. This is controlled by the ``IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT`` cache variable, which accepts either ``RowMajor`` (default) or ``ColMajor``.
+
+To build with column-major derivative ordering, pass the option when configuring CMake:
+
+.. code-block:: bash
+
+    cmake -DIPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=ColMajor -B build -S .
+
+Or to explicitly use the default row-major layout:
+
+.. code-block:: bash
+
+    cmake -DIPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT=RowMajor -B build -S .
+
+.. note::
+
+    This is an advanced option marked as such in CMake. Most users will not need to change it. It is primarily useful when integrating the IPC Toolkit into a simulation framework that stores DOFs in a column-major layout, avoiding the need to permute vectors and matrices at the interface boundary.
+
+This setting affects all functions that assemble local per-element gradients and Hessians into global vectors and matrices, including ``local_gradient_to_global_gradient``, ``local_hessian_to_global_triplets``, and ``local_jacobian_to_global_triplets``. It also affects ``CollisionMesh::vertex_matrix_to_dof_matrix``, which converts vertex-indexed sparse matrices to DOF-indexed sparse matrices.
+
+The chosen layout is stored as the compile-time constant ``ipc::VERTEX_DERIVATIVE_LAYOUT`` (equal to ``Eigen::RowMajor`` or ``Eigen::ColMajor``) in the generated ``config.hpp`` header.
+
+.. warning::
+
+    This feature is experimental. If you encounter any bugs, please report them on `GitHub <https://github.com/ipc-sim/ipc-toolkit/issues>`_.

src/ipc/collision_mesh.cpp

@@ -212,18 +212,27 @@ void CollisionMesh::init_selection_matrices(const int dim)
 Eigen::SparseMatrix<double> CollisionMesh::vertex_matrix_to_dof_matrix(
     const Eigen::SparseMatrix<double>& M_V, int dim)
 {
+    const int n_rows = M_V.rows();
+    const int n_cols = M_V.cols();
+
     std::vector<Eigen::Triplet<double>> triplets;
     using InnerIterator = Eigen::SparseMatrix<double>::InnerIterator;
     for (int k = 0; k < M_V.outerSize(); ++k) {
         for (InnerIterator it(M_V, k); it; ++it) {
             for (int d = 0; d < dim; d++) {
-                triplets.emplace_back(
-                    dim * it.row() + d, dim * it.col() + d, it.value());
+                if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+                    triplets.emplace_back(
+                        dim * it.row() + d, dim * it.col() + d, it.value());
+                } else {
+                    triplets.emplace_back(
+                        n_rows * d + it.row(), n_cols * d + it.col(),
+                        it.value());
+                }
             }
         }
     }
 
-    Eigen::SparseMatrix<double> M_dof(M_V.rows() * dim, M_V.cols() * dim);
+    Eigen::SparseMatrix<double> M_dof(n_rows * dim, n_cols * dim);
     M_dof.setFromTriplets(triplets.begin(), triplets.end());
     M_dof.makeCompressed();
     return M_dof;

src/ipc/collision_mesh.hpp

@@ -322,6 +322,10 @@ class CollisionMesh {
         Eigen::ConstRef<Eigen::MatrixXi> faces,
         Eigen::ConstRef<Eigen::MatrixXi> edges);
 
+    /// @brief Convert a matrix meant for M_V * vertices to M_dof * x by duplicating the entries dim times.
+    static Eigen::SparseMatrix<double> vertex_matrix_to_dof_matrix(
+        const Eigen::SparseMatrix<double>& M_V, int dim);
+
     /// A function that takes two vertex IDs and returns true if the vertices
     /// (and faces or edges containing the vertices) can collide. By default all
     /// primitives can collide with all other primitives.
@@ -351,10 +355,6 @@ class CollisionMesh {
     /// @brief Initialize vertex and edge areas.
     void init_areas();
 
-    /// @brief Convert a matrix meant for M_V * vertices to M_dof * x by duplicating the entries dim times.
-    static Eigen::SparseMatrix<double> vertex_matrix_to_dof_matrix(
-        const Eigen::SparseMatrix<double>& M_V, int dim);
-
     // -----------------------------------------------------------------------
 
     /// @brief The full vertex positions at rest (|FV| × dim).

src/ipc/config.hpp.in

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <Eigen/Core>
+
 #include <cstdint>
 
 // WARNING: Do not modify config.hpp directly. Instead, modify config.hpp.in.
@@ -37,4 +39,9 @@ using index_t = int32_t;
 // using scalar_t = float;
 // #endif
 
+// clang-format off
+/// @brief The layout of derivatives of vertex positions. The default is `Eigen::RowMajor`.
+inline constexpr int VERTEX_DERIVATIVE_LAYOUT = Eigen::@IPC_TOOLKIT_VERTEX_DERIVATIVE_LAYOUT@;
+// clang-format on
+
 } // namespace ipc

src/ipc/potentials/normal_potential.cpp

@@ -38,7 +38,8 @@ Eigen::SparseMatrix<double> NormalPotential::shape_derivative(
                 this->shape_derivative(
                     collisions[i], collisions[i].vertex_ids(edges, faces),
                     collisions[i].dof(rest_positions, edges, faces),
-                    collisions[i].dof(vertices, edges, faces), local_triplets);
+                    collisions[i].dof(vertices, edges, faces), local_triplets,
+                    mesh.num_vertices());
             }
         });
 
@@ -150,7 +151,8 @@ void NormalPotential::shape_derivative(
     const std::array<index_t, 4>& vertex_ids,
     Eigen::ConstRef<VectorMax12d> rest_positions, // = x̄
     Eigen::ConstRef<VectorMax12d> positions,      // = x̄ + u
-    std::vector<Eigen::Triplet<double>>& out) const
+    std::vector<Eigen::Triplet<double>>& out,
+    const int n_total_verts) const
 {
     assert(rest_positions.size() == positions.size());
 
@@ -174,11 +176,17 @@ void NormalPotential::shape_derivative(
 
         for (int i = 0; i < collision.num_vertices(); i++) {
             for (int d = 0; d < dim; d++) {
+                int row_idx;
+                if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+                    row_idx = vertex_ids[i] * dim + d;
+                } else {
+                    assert(n_total_verts > 0);
+                    row_idx = n_total_verts * d + vertex_ids[i];
+                }
                 using Itr = Eigen::SparseVector<double>::InnerIterator;
                 for (Itr j(collision.weight_gradient); j; ++j) {
                     out.emplace_back(
-                        vertex_ids[i] * dim + d, j.index(),
-                        grad_f[dim * i + d] * j.value());
+                        row_idx, j.index(), grad_f[dim * i + d] * j.value());
                 }
             }
         }
@@ -231,7 +239,8 @@ void NormalPotential::shape_derivative(
             + gradu_f * gradu_m.transpose() + m * hessu_f;
     }
 
-    local_hessian_to_global_triplets(local_hess, vertex_ids, dim, out);
+    local_hessian_to_global_triplets(
+        local_hess, vertex_ids, dim, out, n_total_verts);
 }
 
 } // namespace ipc

src/ipc/potentials/normal_potential.hpp

@@ -67,12 +67,14 @@ class NormalPotential : public Potential<NormalCollisions> {
     /// @param[in] rest_positions The collision stencil's rest positions.
     /// @param[in] positions The collision stencil's positions.
     /// @param[in,out] out Store the triplets of the shape derivative here.
+    /// @param[in] n_total_verts The total number of vertices in the mesh, used for computing global indices in the triplets. See also `local_hessian_to_global_triplets`.
     void shape_derivative(
         const NormalCollision& collision,
         const std::array<index_t, 4>& vertex_ids,
         Eigen::ConstRef<VectorMax12d> rest_positions,
         Eigen::ConstRef<VectorMax12d> positions,
-        std::vector<Eigen::Triplet<double>>& out) const;
+        std::vector<Eigen::Triplet<double>>& out,
+        const int n_total_verts = -1) const;
 
     /// @brief Compute the force magnitude for a collision.
     /// @param distance_squared The squared distance between elements.

src/ipc/potentials/potential.cpp

@@ -134,7 +134,8 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
                     collision.vertex_ids(edges, faces);
 
                 local_hessian_to_global_triplets(
-                    local_hess, vids, dim, *(hess_triplets.cache));
+                    local_hess, vids, dim, *(hess_triplets.cache),
+                    mesh.num_vertices());
             }
         });
 

src/ipc/potentials/tangential_potential.cpp

@@ -96,7 +96,8 @@ Eigen::SparseMatrix<double> TangentialPotential::force_jacobian(
                     collision.vertex_ids(mesh.edges(), mesh.faces());
 
                 local_hessian_to_global_triplets(
-                    local_force_jacobian, vis, dim, jac_triplets);
+                    local_force_jacobian, vis, dim, jac_triplets,
+                    mesh.num_vertices());
             }
         });
 
@@ -585,7 +586,8 @@ Eigen::SparseMatrix<double> TangentialPotential::smooth_contact_force_jacobian(
                     collision.vertex_ids(mesh.edges(), mesh.faces());
 
                 local_hessian_to_global_triplets(
-                    local_force_jacobian, vis, dim, jac_triplets);
+                    local_force_jacobian, vis, dim, jac_triplets,
+                    mesh.num_vertices());
 
                 if (wrt == DiffWRT::VELOCITIES) {
                     continue;
@@ -612,7 +614,8 @@ Eigen::SparseMatrix<double> TangentialPotential::smooth_contact_force_jacobian(
 
                 local_jacobian_to_global_triplets(
                     local_force * normal_force_grad.transpose(), vis,
-                    cc_vert_ids, dim, jac_triplets);
+                    cc_vert_ids, dim, jac_triplets, mesh.num_vertices(),
+                    mesh.num_vertices());
             }
         });
 

src/ipc/smooth_contact/smooth_contact_potential.cpp

@@ -109,7 +109,7 @@ Eigen::SparseMatrix<double> SmoothContactPotential::hessian(
 
                 local_hessian_to_global_triplets(
                     local_hess, collision.vertex_ids(), dim,
-                    *(hess_triplets.cache));
+                    *(hess_triplets.cache), mesh.num_vertices());
             }
         });