diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2dc83ba96..d52d55134 100644
--- a/CMakeLists.txt
+++ b/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)
diff --git a/docs/source/tutorials/getting_started.rst b/docs/source/tutorials/getting_started.rst
index 0f945cc37..caa518d61 100644
--- a/docs/source/tutorials/getting_started.rst
+++ b/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
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/docs/source/tutorials/misc.rst b/docs/source/tutorials/misc.rst
index 3abef8fa1..0b1e952a6 100644
--- a/docs/source/tutorials/misc.rst
+++ b/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()
\ No newline at end of file
+            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>`_.
\ No newline at end of file
diff --git a/src/ipc/collision_mesh.cpp b/src/ipc/collision_mesh.cpp
index e37133783..8526a0881 100644
--- a/src/ipc/collision_mesh.cpp
+++ b/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;
diff --git a/src/ipc/collision_mesh.hpp b/src/ipc/collision_mesh.hpp
index ae1a215a6..073193dda 100644
--- a/src/ipc/collision_mesh.hpp
+++ b/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).
diff --git a/src/ipc/config.hpp.in b/src/ipc/config.hpp.in
index 04711cf7b..2cfdba279 100644
--- a/src/ipc/config.hpp.in
+++ b/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
\ No newline at end of file
diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp
index 4e0c6f46f..94eebdda1 100644
--- a/src/ipc/potentials/normal_potential.cpp
+++ b/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
\ No newline at end of file
diff --git a/src/ipc/potentials/normal_potential.hpp b/src/ipc/potentials/normal_potential.hpp
index bae76a514..30bb4827f 100644
--- a/src/ipc/potentials/normal_potential.hpp
+++ b/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.
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index c4da5d2b2..44f1fbe91 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/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());
             }
         });
 
diff --git a/src/ipc/potentials/tangential_potential.cpp b/src/ipc/potentials/tangential_potential.cpp
index f594db53b..7020a3ee7 100644
--- a/src/ipc/potentials/tangential_potential.cpp
+++ b/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());
             }
         });
 
diff --git a/src/ipc/smooth_contact/smooth_contact_potential.cpp b/src/ipc/smooth_contact/smooth_contact_potential.cpp
index f3b34579d..e22b88b00 100644
--- a/src/ipc/smooth_contact/smooth_contact_potential.cpp
+++ b/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());
             }
         });
 
diff --git a/src/ipc/utils/local_to_global.hpp b/src/ipc/utils/local_to_global.hpp
index aa1f3c74f..ffd0f4360 100644
--- a/src/ipc/utils/local_to_global.hpp
+++ b/src/ipc/utils/local_to_global.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <ipc/config.hpp>
+#include <ipc/utils/logger.hpp>
 #include <ipc/utils/matrix_cache.hpp>
 
 #include <Eigen/Core>
@@ -9,70 +11,155 @@
 
 namespace ipc {
 
-template <typename IDContainer>
+/// @brief Add a local gradient to the global gradient, given the vertex ids and dimension.
+/// @tparam IDContainer A container type that supports operator[] and size() to access vertex ids. Can be a std::vector, Eigen::VectorXi, etc.
+/// @tparam GlobalOrder The storage order of the global gradient. Must be either Eigen::RowMajor or Eigen::ColMajor. If RowMajor, the global gradient is assumed to be ordered as [x0, y0, z0, x1, y1, z1, ...]. If ColMajor, the global gradient is assumed to be ordered as [x0, x1, ..., y0, y1, ..., z0, z1, ...].
+/// @param local_grad The local gradient to be added to the global gradient. Must be of size n_verts * dim, where n_verts is the number of vertices in the local element and dim is the dimension of the problem.
+/// @param ids The vertex ids corresponding to the local gradient. Must be of size at least n_verts, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param dim The dimension of the problem (e.g., 2 for 2D, 3 for 3D).
+/// @param grad The global gradient to which the local gradient will be added.
+template <typename IDContainer, int GlobalOrder = VERTEX_DERIVATIVE_LAYOUT>
 void local_gradient_to_global_gradient(
     Eigen::ConstRef<Eigen::VectorXd> local_grad,
     const IDContainer& ids,
     const int dim,
     Eigen::Ref<Eigen::VectorXd> grad)
 {
+    static_assert(
+        GlobalOrder == Eigen::ColMajor || GlobalOrder == Eigen::RowMajor);
     assert(local_grad.size() % dim == 0);
     const int n_verts = local_grad.size() / dim;
     assert(ids.size() >= n_verts); // Can be extra ids
+    const int n_total_verts = grad.size() / dim;
+    assert(grad.size() % dim == 0); // Ensure grad is properly sized
     for (int i = 0; i < n_verts; i++) {
-        grad.segment(dim * ids[i], dim) += local_grad.segment(dim * i, dim);
+        if constexpr (GlobalOrder == Eigen::RowMajor) {
+            grad.segment(dim * ids[i], dim) += local_grad.segment(dim * i, dim);
+        } else {
+            for (int d = 0; d < dim; d++) {
+                grad[n_total_verts * d + ids[i]] += local_grad[dim * i + d];
+            }
+        }
     }
 }
 
-template <typename IDContainer>
+/// @brief Add a local gradient to the global gradient, given the vertex ids and dimension.
+/// @tparam IDContainer A container type that supports operator[] and size() to access vertex ids. Can be a std::vector, Eigen::VectorXi, etc.
+/// @tparam GlobalOrder The storage order of the global gradient. Must be either Eigen::RowMajor or Eigen::ColMajor. If RowMajor, the global gradient is assumed to be ordered as [x0, y0, z0, x1, y1, z1, ...]. If ColMajor, the global gradient is assumed to be ordered as [x0, x1, ..., y0, y1, ..., z0, z1, ...].
+/// @param local_grad The local gradient to be added to the global gradient. Must be of size n_verts * dim, where n_verts is the number of vertices in the local element and dim is the dimension of the problem.
+/// @param ids The vertex ids corresponding to the local gradient. Must be of size at least n_verts, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param dim The dimension of the problem (e.g., 2 for 2D, 3 for 3D).
+/// @param grad The global gradient to which the local gradient will be added.
+template <typename IDContainer, int GlobalOrder = VERTEX_DERIVATIVE_LAYOUT>
 void local_gradient_to_global_gradient(
     Eigen::ConstRef<Eigen::VectorXd> local_grad,
     const IDContainer& ids,
     const int dim,
     Eigen::SparseVector<double>& grad)
 {
+    static_assert(
+        GlobalOrder == Eigen::ColMajor || GlobalOrder == Eigen::RowMajor);
     assert(local_grad.size() % dim == 0);
     const int n_verts = local_grad.size() / dim;
     assert(ids.size() >= n_verts); // Can be extra ids
+    const int n_total_verts = grad.size() / dim;
+    assert(grad.size() % dim == 0); // Ensure grad is properly sized
     for (int i = 0; i < n_verts; i++) {
         for (int d = 0; d < dim; d++) {
-            grad.coeffRef(dim * ids[i] + d) += local_grad(dim * i + d);
+            if constexpr (GlobalOrder == Eigen::RowMajor) {
+                grad.coeffRef(dim * ids[i] + d) += local_grad(dim * i + d);
+            } else {
+                grad.coeffRef(n_total_verts * d + ids[i]) +=
+                    local_grad(dim * i + d);
+            }
         }
     }
 }
 
-template <typename IDContainer>
+/// @brief Add a local Hessian to the global Hessian, given the vertex ids and dimension.
+///
+/// The local Hessian is added to the global Hessian in triplet form, so that it
+/// can be efficiently assembled into a sparse matrix later.
+///
+/// @tparam IDContainer A container type that supports operator[] and size() to access vertex ids. Can be a std::vector, Eigen::VectorXi, etc.
+/// @tparam GlobalOrder The storage order of the global Hessian. Must be either Eigen::RowMajor or Eigen::ColMajor. If RowMajor, the global Hessian is assumed to be ordered as [x0, y0, z0, x1, y1, z1, ...]. If ColMajor, the global Hessian is assumed to be ordered as [x0, x1, ..., y0, y1, ..., z0, z1, ...].
+/// @param local_hessian The local Hessian to be added to the global Hessian. Must be of size (n_verts * dim) x (n_verts * dim), where n_verts is the number of vertices in the local element and dim is the dimension of the problem.
+/// @param ids The vertex ids corresponding to the local gradient. Must be of size at least n_verts, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param dim The dimension of the problem (e.g., 2 for 2D, 3 for 3D).
+/// @param triplets The vector of triplets to which the local Hessian will be added.
+/// @param n_total_verts The total number of vertices in the global mesh. Only required when GlobalOrder is ColMajor. Ignored when GlobalOrder is RowMajor.
+template <typename IDContainer, int GlobalOrder = VERTEX_DERIVATIVE_LAYOUT>
 void local_hessian_to_global_triplets(
     Eigen::ConstRef<Eigen::MatrixXd> local_hessian,
     const IDContainer& ids,
     const int dim,
-    std::vector<Eigen::Triplet<double>>& triplets)
+    std::vector<Eigen::Triplet<double>>& triplets,
+    const int n_total_verts = -1)
 {
+    static_assert(
+        GlobalOrder == Eigen::ColMajor || GlobalOrder == Eigen::RowMajor);
     assert(local_hessian.rows() == local_hessian.cols());
     assert(local_hessian.rows() % dim == 0);
     const int n_verts = local_hessian.rows() / dim;
     assert(ids.size() >= n_verts); // Can be extra ids
+    if constexpr (GlobalOrder == Eigen::ColMajor) {
+        if (n_total_verts <= 0) {
+            log_and_throw_error(
+                "Total number of vertices must be provided for ColMajor storage order.");
+        }
+    }
     for (int i = 0; i < n_verts; i++) {
         for (int j = 0; j < n_verts; j++) {
             for (int k = 0; k < dim; k++) {
                 for (int l = 0; l < dim; l++) {
-                    triplets.emplace_back(
-                        dim * ids[i] + k, dim * ids[j] + l,
-                        local_hessian(dim * i + k, dim * j + l));
+                    if constexpr (GlobalOrder == Eigen::RowMajor) {
+                        triplets.emplace_back(
+                            dim * ids[i] + k, dim * ids[j] + l,
+                            local_hessian(dim * i + k, dim * j + l));
+                    } else {
+                        triplets.emplace_back(
+                            n_total_verts * k + ids[i],
+                            n_total_verts * l + ids[j],
+                            local_hessian(dim * i + k, dim * j + l));
+                    }
                 }
             }
         }
     }
 }
 
-template <typename Derived, typename IDContainer1, typename IDContainer2>
+/// @brief Add a local Jacobian to the global Jacobian, given the vertex ids and dimension.
+///
+/// The local Jacobian is added to the global Jacobian in triplet form, so that
+/// it can be efficiently assembled into a sparse matrix later.
+///
+/// @tparam Derived The derived type of the local Jacobian matrix. Must be a matrix expression that can be evaluated to an Eigen::MatrixXd. The local Jacobian must be of size (n_rows * dim) x (n_cols * dim), where n_rows and n_cols are the number of vertices in the local element for the row and column respectively, and dim is the dimension of the problem.
+/// @tparam IDContainer1 A container type that supports operator[] and size() to access vertex ids for the row indices. Can be a std::vector, Eigen::VectorXi, etc. Must be of size at least n_rows, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @tparam IDContainer2 A container type that supports operator[] and size() to access vertex ids for the column indices. Can be a std::vector, Eigen::VectorXi, etc. Must be of size at least n_cols, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @tparam GlobalOrder The storage order of the global Jacobian. Must be either Eigen::RowMajor or Eigen::ColMajor. If RowMajor, the global Jacobian is assumed to be ordered as [x0, y0, z0, x1, y1, z1, ...]. If ColMajor, the global Jacobian is assumed to be ordered as [x0, x1, ..., y0, y1, ..., z0, z1, ...].
+/// @param local_jacobian The local Jacobian to be added to the global Jacobian. Must be of size (n_rows * dim) x (n_cols * dim), where n_rows and n_cols are the number of vertices in the local element for the row and column respectively, and dim is the dimension of the problem.
+/// @param row_ids The vertex ids corresponding to the row indices of the local Jacobian. Must be of size at least n_rows, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param col_ids The vertex ids corresponding to the column indices of the local Jacobian. Must be of size at least n_cols, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param dim The dimension of the problem (e.g., 2 for 2D, 3 for 3D).
+/// @param triplets The vector of triplets to which the local Hessian will be added.
+/// @param n_total_rows The total number of row vertices in the global mesh. Only required when GlobalOrder is ColMajor. Ignored when GlobalOrder is RowMajor.
+/// @param n_total_cols The total number of column vertices in the global mesh. Only required when GlobalOrder is ColMajor. Ignored when GlobalOrder is RowMajor.
+template <
+    typename Derived,
+    typename IDContainer1,
+    typename IDContainer2,
+    int GlobalOrder = VERTEX_DERIVATIVE_LAYOUT>
 void local_jacobian_to_global_triplets(
     const Eigen::MatrixBase<Derived>& local_jacobian,
     const IDContainer1& row_ids,
     const IDContainer2& col_ids,
     int dim,
-    std::vector<Eigen::Triplet<double>>& triplets)
+    std::vector<Eigen::Triplet<double>>& triplets,
+    const int n_total_rows = -1,
+    const int n_total_cols = -1)
 {
+    static_assert(
+        GlobalOrder == Eigen::ColMajor || GlobalOrder == Eigen::RowMajor);
     // assert(local_jacobian.rows() == local_jacobian.cols());
     assert(local_jacobian.rows() % dim == 0);
     assert(local_jacobian.cols() % dim == 0);
@@ -80,21 +167,37 @@ void local_jacobian_to_global_triplets(
     const int n_cols = local_jacobian.cols() / dim;
     assert(row_ids.size() >= n_rows); // Can be extra ids
     assert(col_ids.size() >= n_cols); // Can be extra ids
+    if constexpr (GlobalOrder == Eigen::ColMajor) {
+        if (n_total_rows <= 0 || n_total_cols <= 0) {
+            log_and_throw_error(
+                "Total number of rows and columns must be provided for ColMajor storage order.");
+        }
+    }
     for (int i = 0; i < n_rows; i++) {
         for (int j = 0; j < n_cols; j++) {
             for (int k = 0; k < dim; k++) {
                 for (int l = 0; l < dim; l++) {
-                    triplets.emplace_back(
-                        dim * row_ids[i] + k, dim * col_ids[j] + l,
-                        local_jacobian(dim * i + k, dim * j + l));
+                    if constexpr (GlobalOrder == Eigen::RowMajor) {
+                        triplets.emplace_back(
+                            dim * row_ids[i] + k, dim * col_ids[j] + l,
+                            local_jacobian(dim * i + k, dim * j + l));
+                    } else {
+                        triplets.emplace_back(
+                            n_total_rows * k + row_ids[i],
+                            n_total_cols * l + col_ids[j],
+                            local_jacobian(dim * i + k, dim * j + l));
+                    }
                 }
             }
         }
     }
 }
 
-class LocalThreadMatStorage {
-public:
+/// @brief A helper class to store a local matrix cache for each thread.
+/// This is useful for parallelizing the assembly of the global Hessian, where
+/// each thread can have its own local cache to store the triplets before they
+/// are merged into the global cache.
+struct LocalThreadMatStorage {
     std::unique_ptr<MatrixCache> cache = nullptr;
 
     LocalThreadMatStorage() = delete;
@@ -138,25 +241,56 @@ class LocalThreadMatStorage {
     }
 };
 
-template <typename Derived, typename IDContainer>
+/// @brief Add a local Hessian to the global Hessian, given the vertex ids and dimension.
+///
+/// The local Hessian is added to the global Hessian in triplet form, so that it
+/// can be efficiently assembled into a sparse matrix later.
+///
+/// @tparam Derived The derived type of the local Hessian matrix. Must be a matrix expression that can be evaluated to an Eigen::MatrixXd. The local Hessian must be of size (n_verts * dim) x (n_verts * dim), where n_verts is the number of vertices in the local element and dim is the dimension of the problem.
+/// @tparam IDContainer A container type that supports operator[] and size() to access vertex ids. Can be a std::vector, Eigen::VectorXi, etc. Must be of size at least n_verts, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @tparam GlobalOrder The storage order of the global gradient. Must be either Eigen::RowMajor or Eigen::ColMajor. If RowMajor, the global gradient is assumed to be ordered as [x0, y0, z0, x1, y1, z1, ...]. If ColMajor, the global gradient is assumed to be ordered as [x0, x1, ..., y0, y1, ..., z0, z1, ...].
+/// @param local_hessian The local Hessian to be added to the global Hessian. Must be of size (n_verts * dim) x (n_verts * dim), where n_verts is the number of vertices in the local element and dim is the dimension of the problem.
+/// @param ids The vertex ids corresponding to the local gradient. Must be of size at least n_verts, but can be larger (e.g., if the local element has fewer vertices than the number of ids).
+/// @param dim The dimension of the problem (e.g., 2 for 2D, 3 for 3D).
+/// @param triplets The MatrixCache to which the local Hessian will be added. The triplets will be added to the first matrix in the cache (i.e., matrix index 0). The cache should be initialized with the correct number of rows and columns for the global Hessian, and should have enough reserved space for the number of triplets that will be added.
+/// @param n_total_verts The total number of vertices in the global mesh. Only required when GlobalOrder is ColMajor. Ignored when GlobalOrder is RowMajor.
+template <
+    typename Derived,
+    typename IDContainer,
+    int GlobalOrder = VERTEX_DERIVATIVE_LAYOUT>
 void local_hessian_to_global_triplets(
     const Eigen::MatrixBase<Derived>& local_hessian,
     const IDContainer& ids,
     int dim,
-    MatrixCache& triplets)
+    MatrixCache& triplets,
+    const int n_total_verts = -1)
 {
+    static_assert(
+        GlobalOrder == Eigen::ColMajor || GlobalOrder == Eigen::RowMajor);
     assert(local_hessian.rows() == local_hessian.cols());
     assert(local_hessian.rows() % dim == 0);
     const int n_verts = local_hessian.rows() / dim;
     assert(ids.size() >= n_verts); // Can be extra ids
+    if constexpr (GlobalOrder == Eigen::ColMajor) {
+        if (n_total_verts <= 0) {
+            log_and_throw_error(
+                "Total number of vertices must be provided for ColMajor storage order.");
+        }
+    }
     for (int i = 0; i < n_verts; i++) {
         for (int j = 0; j < n_verts; j++) {
             for (int k = 0; k < dim; k++) {
                 for (int l = 0; l < dim; l++) {
                     const auto& val = local_hessian(dim * i + k, dim * j + l);
                     if (val != 0) {
-                        triplets.add_value(
-                            0, dim * ids[i] + k, dim * ids[j] + l, val);
+                        if constexpr (GlobalOrder == Eigen::RowMajor) {
+                            triplets.add_value(
+                                0, dim * ids[i] + k, dim * ids[j] + l, val);
+                        } else {
+                            triplets.add_value(
+                                0, n_total_verts * k + ids[i],
+                                n_total_verts * l + ids[j], val);
+                        }
                     }
                 }
             }
diff --git a/tests/src/tests/collisions/test_normal_collisions.cpp b/tests/src/tests/collisions/test_normal_collisions.cpp
index 5f62a924e..c56dee4f9 100644
--- a/tests/src/tests/collisions/test_normal_collisions.cpp
+++ b/tests/src/tests/collisions/test_normal_collisions.cpp
@@ -3,6 +3,7 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
 
+#include <ipc/config.hpp>
 #include <ipc/collisions/normal/normal_collisions.hpp>
 #include <ipc/potentials/barrier_potential.hpp>
 
@@ -91,8 +92,14 @@ TEST_CASE("Codim. vertex-vertex collisions", "[collisions][codim]")
         CHECK(barrier_potential(collisions, mesh, vertices) > 0.0);
         const Eigen::VectorXd grad =
             barrier_potential.gradient(collisions, mesh, vertices);
-        for (int i = 0; i < vertices.rows(); i++) {
-            const Eigen::Vector3d f = -grad.segment<3>(3 * i);
+        const int n_verts = vertices.rows();
+        for (int i = 0; i < n_verts; i++) {
+            Eigen::Vector3d f;
+            if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+                f = -grad.segment<3>(3 * i);
+            } else {
+                f << -grad[i], -grad[n_verts + i], -grad[2 * n_verts + i];
+            }
             CHECK(f.normalized().isApprox(
                 vertices.row(i).normalized().transpose()));
         }
diff --git a/tests/src/tests/dof_layout.hpp b/tests/src/tests/dof_layout.hpp
new file mode 100644
index 000000000..8dbb26fa1
--- /dev/null
+++ b/tests/src/tests/dof_layout.hpp
@@ -0,0 +1,101 @@
+/// @file
+/// Utility helpers that make finite-difference tests work regardless of
+/// whether ipc::VERTEX_DERIVATIVE_LAYOUT is RowMajor or ColMajor.
+///
+/// fd::flatten / fd::unflatten always use **RowMajor** ordering:
+///   flatten   :  MatrixXd(n_verts, dim)  →  [x0 y0 z0 x1 y1 z1 …]
+///   unflatten :  [x0 y0 z0 x1 y1 z1 …]  →  MatrixXd(n_verts, dim)
+///
+/// When VERTEX_DERIVATIVE_LAYOUT == ColMajor the library returns gradient /
+/// hessian values laid out as
+///   [x0 x1 …  y0 y1 …  z0 z1 …]   (columns-of-vertices first)
+///
+/// The helpers below bridge the two worlds so that every test can be written as
+///
+///   fd::finite_gradient(dof_flatten(V), [&](auto& x){ … dof_unflatten(x, V) …
+///   }, fgrad);
+///   fgrad = dof_reorder_grad(fgrad, n_verts, dim);   // RowMajor → layout
+///
+/// and the comparison with the library result is valid.
+
+#pragma once
+
+#include <ipc/config.hpp>
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+#include <cassert>
+
+namespace ipc::tests {
+
+// ─── Flatten / Unflatten in VERTEX_DERIVATIVE_LAYOUT order ──────────────────
+
+/// Flatten a vertex matrix (n_verts × dim) to a DOF vector in the order
+/// dictated by VERTEX_DERIVATIVE_LAYOUT.
+///
+///   RowMajor → [x0 y0 z0  x1 y1 z1  …]   (same as fd::flatten)
+///   ColMajor → [x0 x1 …   y0 y1 …   z0 z1 …]
+inline auto flatten(const Eigen::MatrixXd& V)
+{
+    return V.reshaped<VERTEX_DERIVATIVE_LAYOUT>();
+}
+
+/// Unflatten a DOF vector back to a vertex matrix (n_verts × dim).
+/// The inverse of dof_flatten.
+inline auto unflatten(const Eigen::VectorXd& x, int dim)
+{
+    return x.reshaped<VERTEX_DERIVATIVE_LAYOUT>(x.size() / dim, dim);
+}
+
+// ─── Reorder FD results from RowMajor to VERTEX_DERIVATIVE_LAYOUT ──────────
+//
+// These are used when the FD callback is forced to stay in RowMajor order
+// (e.g. because it calls fd::flatten/unflatten internally) but the library
+// output is in VERTEX_DERIVATIVE_LAYOUT.
+
+/// Reorder a gradient vector computed with RowMajor flattening so that it
+/// matches VERTEX_DERIVATIVE_LAYOUT.
+///
+/// When the layout is already RowMajor this is a no-op (compiler will
+/// optimize the copy away or it can be elided).
+inline auto
+reorder_gradient(const Eigen::VectorXd& g_rowmajor, int n_verts, int dim)
+{
+    if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+        return g_rowmajor; // Already in the correct order, no copy needed.
+    } else {
+        return g_rowmajor.reshaped<Eigen::RowMajor>(n_verts, dim)
+            .reshaped<VERTEX_DERIVATIVE_LAYOUT>()
+            .eval();
+    }
+}
+
+/// Reorder a sparse Hessian from RowMajor DOF ordering to
+/// VERTEX_DERIVATIVE_LAYOUT DOF ordering.
+inline auto reorder_hessian(
+    const Eigen::SparseMatrix<double>& H_rowmajor, int n_verts, int dim)
+{
+    if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+        return H_rowmajor; // Already in the correct order, no copy needed.
+    } else {
+        const int n = n_verts * dim;
+        assert(H_rowmajor.rows() == n && H_rowmajor.cols() == n);
+        std::vector<Eigen::Triplet<double>> trips;
+        trips.reserve(H_rowmajor.nonZeros());
+        for (int k = 0; k < H_rowmajor.outerSize(); ++k) {
+            for (Eigen::SparseMatrix<double>::InnerIterator it(H_rowmajor, k);
+                 it; ++it) {
+                // row_major → col_major
+                trips.emplace_back(
+                    (it.row() % dim) * n_verts + it.row() / dim,
+                    (it.col() % dim) * n_verts + it.col() / dim, it.value());
+            }
+        }
+        Eigen::SparseMatrix<double> H(n, n);
+        H.setFromTriplets(trips.begin(), trips.end());
+        return H;
+    }
+}
+
+} // namespace ipc::tests
\ No newline at end of file
diff --git a/tests/src/tests/friction/test_force_jacobian.cpp b/tests/src/tests/friction/test_force_jacobian.cpp
index 048913619..591f6fc6e 100644
--- a/tests/src/tests/friction/test_force_jacobian.cpp
+++ b/tests/src/tests/friction/test_force_jacobian.cpp
@@ -1,4 +1,5 @@
 #include <tests/config.hpp>
+#include <tests/dof_layout.hpp>
 #include <tests/friction/friction_data_generator.hpp>
 #include <tests/utils.hpp>
 
@@ -60,11 +61,11 @@ void check_friction_force_jacobian(
     }
     auto PA_X = [&](const Eigen::VectorXd& x) {
         CollisionMesh fd_mesh(
-            fd::unflatten(x, X.cols()), mesh.edges(), mesh.faces());
+            tests::unflatten(x, X.cols()), mesh.edges(), mesh.faces());
         return fd_mesh.vertex_areas();
     };
     Eigen::MatrixXd fd_JPA_wrt_X;
-    fd::finite_jacobian(fd::flatten(X), PA_X, fd_JPA_wrt_X);
+    fd::finite_jacobian(tests::flatten(X), PA_X, fd_JPA_wrt_X);
 
     CHECK(fd::compare_jacobian(JPA_wrt_X, fd_JPA_wrt_X));
     if (!fd::compare_jacobian(JPA_wrt_X, fd_JPA_wrt_X)) {
@@ -79,11 +80,11 @@ void check_friction_force_jacobian(
     }
     auto EA_X = [&](const Eigen::VectorXd& x) {
         CollisionMesh fd_mesh(
-            fd::unflatten(x, X.cols()), mesh.edges(), mesh.faces());
+            tests::unflatten(x, X.cols()), mesh.edges(), mesh.faces());
         return fd_mesh.edge_areas();
     };
     Eigen::MatrixXd fd_JEA_wrt_X;
-    fd::finite_jacobian(fd::flatten(X), EA_X, fd_JEA_wrt_X);
+    fd::finite_jacobian(tests::flatten(X), EA_X, fd_JEA_wrt_X);
 
     CHECK(fd::compare_jacobian(JEA_wrt_X, fd_JEA_wrt_X));
     if (!fd::compare_jacobian(JEA_wrt_X, fd_JEA_wrt_X)) {
@@ -97,7 +98,7 @@ void check_friction_force_jacobian(
         FrictionPotential::DiffWRT::REST_POSITIONS);
 
     auto F_X = [&](const Eigen::VectorXd& x) {
-        Eigen::MatrixXd fd_X = fd::unflatten(x, X.cols());
+        Eigen::MatrixXd fd_X = tests::unflatten(x, X.cols());
 
         CollisionMesh fd_mesh(fd_X, mesh.edges(), mesh.faces());
         fd_mesh.init_area_jacobians();
@@ -122,7 +123,7 @@ void check_friction_force_jacobian(
             fd_friction_collisions, fd_mesh, fd_X, Ut, velocities, B);
     };
     Eigen::MatrixXd fd_JF_wrt_X;
-    fd::finite_jacobian(fd::flatten(X), F_X, fd_JF_wrt_X);
+    fd::finite_jacobian(tests::flatten(X), F_X, fd_JF_wrt_X);
 
     CHECK(fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X));
     if (!fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X)) {
@@ -135,7 +136,7 @@ void check_friction_force_jacobian(
         FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS);
 
     auto F_Ut = [&](const Eigen::VectorXd& ut) {
-        Eigen::MatrixXd fd_Ut = fd::unflatten(ut, Ut.cols());
+        Eigen::MatrixXd fd_Ut = tests::unflatten(ut, Ut.cols());
 
         TangentialCollisions fd_friction_collisions;
         if (recompute_collisions) {
@@ -156,7 +157,7 @@ void check_friction_force_jacobian(
         return D.force(tangential_collisions, mesh, X, fd_Ut, velocities, B);
     };
     Eigen::MatrixXd fd_JF_wrt_Ut;
-    fd::finite_jacobian(fd::flatten(Ut), F_Ut, fd_JF_wrt_Ut);
+    fd::finite_jacobian(tests::flatten(Ut), F_Ut, fd_JF_wrt_Ut);
 
     CHECK(fd::compare_jacobian(JF_wrt_Ut, fd_JF_wrt_Ut));
     if (!fd::compare_jacobian(JF_wrt_Ut, fd_JF_wrt_Ut)) {
@@ -172,10 +173,10 @@ void check_friction_force_jacobian(
     auto F_V = [&](const Eigen::VectorXd& v) {
         return D.force(
             tangential_collisions, mesh, X, Ut,
-            fd::unflatten(v, velocities.cols()), B);
+            tests::unflatten(v, velocities.cols()), B);
     };
     Eigen::MatrixXd fd_JF_wrt_V;
-    fd::finite_jacobian(fd::flatten(velocities), F_V, fd_JF_wrt_V);
+    fd::finite_jacobian(tests::flatten(velocities), F_V, fd_JF_wrt_V);
 
     CHECK(fd::compare_jacobian(JF_wrt_V, fd_JF_wrt_V));
     if (!fd::compare_jacobian(JF_wrt_V, fd_JF_wrt_V)) {
@@ -189,10 +190,11 @@ void check_friction_force_jacobian(
 
     auto grad = [&](const Eigen::VectorXd& v) {
         return D.gradient(
-            tangential_collisions, mesh, fd::unflatten(v, velocities.cols()));
+            tangential_collisions, mesh,
+            tests::unflatten(v, velocities.cols()));
     };
     Eigen::MatrixXd fd_hessian;
-    fd::finite_jacobian(fd::flatten(velocities), grad, fd_hessian);
+    fd::finite_jacobian(tests::flatten(velocities), grad, fd_hessian);
 
     CHECK(fd::compare_jacobian(hess_D, fd_hessian));
     if (!fd::compare_jacobian(hess_D, fd_hessian)) {
@@ -301,9 +303,9 @@ TEST_CASE(
         const auto dir = tests::DATA_DIR / "friction-force-jacobian" / scene;
         X = tests::loadMarketXd((dir / "X.mtx").string());
         Ut = tests::loadMarketXd((dir / "Ut.mtx").string());
-        Ut = fd::unflatten(Ut, X.cols());
+        Ut = fd::unflatten(Ut, X.cols()); // data files are always RowMajor
         U = tests::loadMarketXd((dir / "U.mtx").string());
-        U = fd::unflatten(U, X.cols());
+        U = fd::unflatten(U, X.cols()); // data files are always RowMajor
         if (is_2D) {
             E = tests::loadMarketXi((dir / "F.mtx").string());
         } else {
@@ -390,11 +392,11 @@ void check_smooth_friction_force_jacobian(
 
     auto grad = [&](const Eigen::VectorXd& v) {
         return D.gradient(
-            friction_collisions, mesh, fd::unflatten(v, velocities.cols()));
+            friction_collisions, mesh, tests::unflatten(v, velocities.cols()));
     };
     Eigen::MatrixXd fd_hessian;
     fd::finite_jacobian(
-        fd::flatten(velocities), grad, fd_hessian, fd::AccuracyOrder::FOURTH,
+        tests::flatten(velocities), grad, fd_hessian, fd::AccuracyOrder::FOURTH,
         1e-6 * dhat);
     // CHECK(fd::compare_jacobian(hess_D, fd_hessian));
     // if (!fd::compare_jacobian(hess_D, fd_hessian)) {
@@ -511,7 +513,7 @@ void check_smooth_friction_force_jacobian(
         FrictionPotential::DiffWRT::REST_POSITIONS);
 
     auto F_X = [&](const Eigen::VectorXd& x) {
-        Eigen::MatrixXd fd_X = fd::unflatten(x, X.cols());
+        Eigen::MatrixXd fd_X = tests::unflatten(x, X.cols());
         Eigen::MatrixXd fd_lagged_positions = fd_X + Ut;
 
         CollisionMesh fd_mesh(fd_X, mesh.edges(), mesh.faces());
@@ -530,7 +532,7 @@ void check_smooth_friction_force_jacobian(
     };
     Eigen::MatrixXd fd_JF_wrt_X;
     fd::finite_jacobian(
-        fd::flatten(X), F_X, fd_JF_wrt_X, fd::AccuracyOrder::FOURTH,
+        tests::flatten(X), F_X, fd_JF_wrt_X, fd::AccuracyOrder::FOURTH,
         1e-6 * dhat);
     // CHECK(fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X));
     // if (!fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X)) {
@@ -547,7 +549,7 @@ void check_smooth_friction_force_jacobian(
         FrictionPotential::DiffWRT::LAGGED_DISPLACEMENTS);
 
     auto F_Ut = [&](const Eigen::VectorXd& ut) {
-        Eigen::MatrixXd fd_Ut = fd::unflatten(ut, Ut.cols());
+        Eigen::MatrixXd fd_Ut = tests::unflatten(ut, Ut.cols());
         Eigen::MatrixXd fd_lagged_positions = X + fd_Ut;
 
         auto fd_collisions = create_smooth_collision(mesh, fd_lagged_positions);
@@ -563,7 +565,7 @@ void check_smooth_friction_force_jacobian(
     };
     Eigen::MatrixXd fd_JF_wrt_Ut;
     fd::finite_jacobian(
-        fd::flatten(Ut), F_Ut, fd_JF_wrt_Ut, fd::AccuracyOrder::FOURTH,
+        tests::flatten(Ut), F_Ut, fd_JF_wrt_Ut, fd::AccuracyOrder::FOURTH,
         1e-6 * dhat);
     // CHECK(fd::compare_jacobian(JF_wrt_Ut, fd_JF_wrt_Ut));
     // if (!fd::compare_jacobian(JF_wrt_Ut, fd_JF_wrt_Ut)) {
@@ -582,11 +584,11 @@ void check_smooth_friction_force_jacobian(
     auto F_V = [&](const Eigen::VectorXd& v) {
         return D.smooth_contact_force(
             friction_collisions, mesh, X, Ut,
-            fd::unflatten(v, velocities.cols()));
+            tests::unflatten(v, velocities.cols()));
     };
     Eigen::MatrixXd fd_JF_wrt_V;
     fd::finite_jacobian(
-        fd::flatten(velocities), F_V, fd_JF_wrt_V, fd::AccuracyOrder::FOURTH,
+        tests::flatten(velocities), F_V, fd_JF_wrt_V, fd::AccuracyOrder::FOURTH,
         1e-6 * dhat);
     CHECK(
         (JF_wrt_V.norm() == 0
diff --git a/tests/src/tests/friction/test_friction.cpp b/tests/src/tests/friction/test_friction.cpp
index 54926584d..24629f526 100644
--- a/tests/src/tests/friction/test_friction.cpp
+++ b/tests/src/tests/friction/test_friction.cpp
@@ -1,4 +1,5 @@
 #include <tests/config.hpp>
+#include <tests/dof_layout.hpp>
 #include <tests/utils.hpp>
 
 #include <catch2/catch_test_macros.hpp>
@@ -242,6 +243,15 @@ TEST_CASE(
 
     Eigen::VectorXd grad = D.gradient(tangential_collisions, mesh, velocity);
 
+    if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::ColMajor) {
+        // The saved expected data is in RowMajor DOF ordering; convert to the
+        // current VERTEX_DERIVATIVE_LAYOUT before comparison.
+        const int n_verts = V_start.rows();
+        const int dim = V_start.cols();
+        expected_grad = tests::reorder_gradient(expected_grad, n_verts, dim);
+        expected_hess = tests::reorder_hessian(expected_hess, n_verts, dim);
+    }
+
     CHECK(grad.isApprox(expected_grad));
 
     Eigen::SparseMatrix<double> hess =
diff --git a/tests/src/tests/potential/test_adhesion_potentials.cpp b/tests/src/tests/potential/test_adhesion_potentials.cpp
index 5dbf11435..1d4b06282 100644
--- a/tests/src/tests/potential/test_adhesion_potentials.cpp
+++ b/tests/src/tests/potential/test_adhesion_potentials.cpp
@@ -2,6 +2,7 @@
 #include <catch2/catch_approx.hpp>
 
 #include <tests/utils.hpp>
+#include <tests/dof_layout.hpp>
 
 #include <ipc/ipc.hpp>
 #include <ipc/adhesion/adhesion.hpp>
@@ -93,10 +94,10 @@ TEST_CASE("Normal adhesion potential", "[potential][adhesion]")
 
     Eigen::VectorXd fgrad;
     fd::finite_gradient(
-        fd::flatten(vertices),
+        tests::flatten(vertices),
         [&](const Eigen::VectorXd& x) {
             return potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         },
         fgrad);
 
@@ -111,10 +112,10 @@ TEST_CASE("Normal adhesion potential", "[potential][adhesion]")
         Eigen::MatrixXd fhess;
 
         fd::finite_jacobian(
-            fd::flatten(vertices),
+            tests::flatten(vertices),
             [&](const Eigen::VectorXd& x) {
                 return potential.gradient(
-                    collisions, mesh, fd::unflatten(x, vertices.cols()));
+                    collisions, mesh, tests::unflatten(x, vertices.cols()));
             },
             fhess);
 
diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp
index 55d344ac2..2bb56f509 100644
--- a/tests/src/tests/potential/test_barrier_potential.cpp
+++ b/tests/src/tests/potential/test_barrier_potential.cpp
@@ -1,4 +1,5 @@
 #include <tests/config.hpp>
+#include <tests/dof_layout.hpp>
 #include <tests/utils.hpp>
 
 #include <catch2/catch_test_macros.hpp>
@@ -86,9 +87,9 @@ TEST_CASE(
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return barrier_potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
-        fd::finite_gradient(fd::flatten(vertices), f, fgrad_b);
+        fd::finite_gradient(tests::flatten(vertices), f, fgrad_b);
     }
 
     REQUIRE(grad_b.squaredNorm() > 0);
@@ -106,9 +107,9 @@ TEST_CASE(
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return barrier_potential.gradient(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
-        fd::finite_jacobian(fd::flatten(vertices), f, fhess_b);
+        fd::finite_jacobian(tests::flatten(vertices), f, fhess_b);
     }
 
     REQUIRE(hess_b.squaredNorm() > 0);
@@ -228,7 +229,7 @@ TEST_CASE(
 
     // Compute the gradient using finite differences
     auto f = [&](const Eigen::VectorXd& x) {
-        const Eigen::MatrixXd fd_V = fd::unflatten(x, mesh.dim());
+        const Eigen::MatrixXd fd_V = tests::unflatten(x, mesh.dim());
 
         NormalCollisions fd_collisions;
         fd_collisions.set_use_area_weighting(use_area_weighting);
@@ -239,7 +240,7 @@ TEST_CASE(
         return barrier_potential(collisions, mesh, fd_V);
     };
     Eigen::VectorXd fgrad_b;
-    fd::finite_gradient(fd::flatten(vertices), f, fgrad_b);
+    fd::finite_gradient(tests::flatten(vertices), f, fgrad_b);
 
     CHECK(fd::compare_gradient(grad_b, fgrad_b));
 }
@@ -305,7 +306,8 @@ TEST_CASE(
         barrier_potential.shape_derivative(
             collisions[i], collisions[i].vertex_ids(edges, faces),
             collisions[i].dof(rest_positions, edges, faces),
-            collisions[i].dof(vertices, edges, faces), triplets);
+            collisions[i].dof(vertices, edges, faces), triplets,
+            vertices.rows());
         Eigen::SparseMatrix<double> JF_wrt_X_sparse(ndof, ndof);
         JF_wrt_X_sparse.setFromTriplets(triplets.begin(), triplets.end());
         const Eigen::MatrixXd JF_wrt_X = Eigen::MatrixXd(JF_wrt_X_sparse);
@@ -314,20 +316,19 @@ TEST_CASE(
             // TODO: Recompute weight based on x
             assert(use_area_weighting == false);
             // Recompute eps_x based on x
+            const Eigen::MatrixXd X =
+                tests::unflatten(x, rest_positions.cols());
             double prev_eps_x = -1;
             if (collisions.is_edge_edge(i)) {
                 EdgeEdgeNormalCollision& c =
                     dynamic_cast<EdgeEdgeNormalCollision&>(collisions[i]);
                 prev_eps_x = c.eps_x;
                 c.eps_x = edge_edge_mollifier_threshold(
-                    x.segment<3>(3 * edges(c.edge0_id, 0)),
-                    x.segment<3>(3 * edges(c.edge0_id, 1)),
-                    x.segment<3>(3 * edges(c.edge1_id, 0)),
-                    x.segment<3>(3 * edges(c.edge1_id, 1)));
+                    X.row(edges(c.edge0_id, 0)), X.row(edges(c.edge0_id, 1)),
+                    X.row(edges(c.edge1_id, 0)), X.row(edges(c.edge1_id, 1)));
             }
 
-            const Eigen::MatrixXd positions =
-                fd::unflatten(x, rest_positions.cols()) + displacements;
+            const Eigen::MatrixXd positions = X + displacements;
             const VectorMax12d dof = collisions[i].dof(positions, edges, faces);
 
             Eigen::VectorXd grad = Eigen::VectorXd::Zero(ndof);
@@ -346,7 +347,7 @@ TEST_CASE(
         };
 
         Eigen::MatrixXd fd_JF_wrt_X;
-        fd::finite_jacobian(fd::flatten(rest_positions), F_X, fd_JF_wrt_X);
+        fd::finite_jacobian(tests::flatten(rest_positions), F_X, fd_JF_wrt_X);
 
         CHECK(fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X));
         if (!fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X)) {
@@ -365,7 +366,8 @@ TEST_CASE(
         barrier_potential.shape_derivative(
             collisions[i], collisions[i].vertex_ids(edges, faces),
             collisions[i].dof(rest_positions, edges, faces),
-            collisions[i].dof(vertices, edges, faces), triplets);
+            collisions[i].dof(vertices, edges, faces), triplets,
+            vertices.rows());
         Eigen::SparseMatrix<double> JF_wrt_X_sparse(ndof, ndof);
         JF_wrt_X_sparse.setFromTriplets(triplets.begin(), triplets.end());
         sum += Eigen::MatrixXd(JF_wrt_X_sparse);
@@ -373,7 +375,7 @@ TEST_CASE(
     CHECK(fd::compare_jacobian(JF_wrt_X, sum));
 
     auto F_X = [&](const Eigen::VectorXd& x) {
-        const Eigen::MatrixXd fd_X = fd::unflatten(x, rest_positions.cols());
+        const Eigen::MatrixXd fd_X = tests::unflatten(x, rest_positions.cols());
         const Eigen::MatrixXd fd_V = fd_X + displacements;
 
         CollisionMesh fd_mesh(fd_X, mesh.edges(), mesh.faces());
@@ -388,7 +390,7 @@ TEST_CASE(
         return barrier_potential.gradient(collisions, fd_mesh, fd_V);
     };
     Eigen::MatrixXd fd_JF_wrt_X;
-    fd::finite_jacobian(fd::flatten(rest_positions), F_X, fd_JF_wrt_X);
+    fd::finite_jacobian(tests::flatten(rest_positions), F_X, fd_JF_wrt_X);
 
     CHECK(fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X));
     if (!fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X)) {
@@ -443,7 +445,7 @@ TEST_CASE(
         barrier_potential.shape_derivative(collisions, mesh, vertices);
 
     auto F_X = [&](const Eigen::VectorXd& x) {
-        const Eigen::MatrixXd fd_X = fd::unflatten(x, X.cols());
+        const Eigen::MatrixXd fd_X = tests::unflatten(x, X.cols());
         const Eigen::MatrixXd fd_V = fd_X + U;
 
         CollisionMesh fd_mesh(fd_X, mesh.edges(), mesh.faces());
@@ -456,7 +458,7 @@ TEST_CASE(
         return barrier_potential.gradient(fd_collisions, fd_mesh, fd_V);
     };
     Eigen::MatrixXd fd_JF_wrt_X;
-    fd::finite_jacobian(fd::flatten(X), F_X, fd_JF_wrt_X);
+    fd::finite_jacobian(tests::flatten(X), F_X, fd_JF_wrt_X);
 
     CHECK(fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X));
     if (!fd::compare_jacobian(JF_wrt_X, fd_JF_wrt_X)) {
diff --git a/tests/src/tests/potential/test_friction_potential.cpp b/tests/src/tests/potential/test_friction_potential.cpp
index f428adde4..a66495e0a 100644
--- a/tests/src/tests/potential/test_friction_potential.cpp
+++ b/tests/src/tests/potential/test_friction_potential.cpp
@@ -1,4 +1,5 @@
 #include <tests/config.hpp>
+#include <tests/dof_layout.hpp>
 #include <tests/friction/friction_data_generator.hpp>
 #include <tests/utils.hpp>
 
@@ -32,15 +33,16 @@ TEST_CASE("Friction gradient and hessian", "[friction][gradient][hessian]")
 
     // Compute the gradient using finite differences
     auto f = [&](const Eigen::VectorXd& x) {
-        const Eigen::MatrixXd fd_U = fd::unflatten(x, data.V1.cols()) - data.V0;
+        const Eigen::MatrixXd fd_U =
+            tests::unflatten(x, data.V1.cols()) - data.V0;
         return D(tangential_collisions, mesh, fd_U);
     };
     Eigen::VectorXd fgrad;
-    fd::finite_gradient(fd::flatten(V1), f, fgrad);
+    fd::finite_gradient(tests::flatten(V1), f, fgrad);
     CHECK(fd::compare_gradient(grad, fgrad));
 
     const Eigen::MatrixXd hess = D.hessian(tangential_collisions, mesh, U);
     Eigen::MatrixXd fhess;
-    fd::finite_hessian(fd::flatten(V1), f, fhess);
+    fd::finite_hessian(tests::flatten(V1), f, fhess);
     CHECK(fd::compare_hessian(hess, fhess, 1e-3));
 }
\ No newline at end of file
diff --git a/tests/src/tests/potential/test_smooth_potential.cpp b/tests/src/tests/potential/test_smooth_potential.cpp
index ca073a2c7..2e9af48e3 100644
--- a/tests/src/tests/potential/test_smooth_potential.cpp
+++ b/tests/src/tests/potential/test_smooth_potential.cpp
@@ -1,10 +1,12 @@
 #include <tests/config.hpp>
+#include <tests/dof_layout.hpp>
 #include <tests/utils.hpp>
 
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
 #include <catch2/benchmark/catch_benchmark.hpp>
 
+#include <ipc/ipc.hpp>
 #include <ipc/potentials/barrier_potential.hpp>
 #include <ipc/smooth_contact/smooth_contact_potential.hpp>
 #include <ipc/distance/line_line.hpp>
@@ -12,7 +14,6 @@
 #include <finitediff.hpp>
 #include <igl/edges.h>
 #include <igl/readCSV.h>
-#include <ipc/ipc.hpp>
 
 using namespace ipc;
 
@@ -65,10 +66,11 @@ TEST_CASE("Smooth barrier potential codim", "[smooth_potential]")
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_gradient(
-            fd::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     // REQUIRE(grad_b.squaredNorm() > 1e-8);
@@ -88,10 +90,11 @@ TEST_CASE("Smooth barrier potential codim", "[smooth_potential]")
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential.gradient(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_jacobian(
-            fd::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     REQUIRE(hess_b.squaredNorm() > 1e-8);
@@ -187,10 +190,11 @@ TEST_CASE("Smooth barrier potential full gradient and hessian 3D", tagsopt)
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_gradient(
-            fd::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     // REQUIRE(grad_b.squaredNorm() > 1e-8);
@@ -210,10 +214,11 @@ TEST_CASE("Smooth barrier potential full gradient and hessian 3D", tagsopt)
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential.gradient(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_jacobian(
-            fd::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     REQUIRE(hess_b.squaredNorm() > 1e-8);
@@ -279,10 +284,11 @@ TEST_CASE("Smooth barrier potential real sim 2D C^2", "[smooth_potential]")
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_gradient(
-            fd::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     REQUIRE(grad_b.squaredNorm() > 1e-8);
@@ -302,10 +308,11 @@ TEST_CASE("Smooth barrier potential real sim 2D C^2", "[smooth_potential]")
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential.gradient(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_jacobian(
-            fd::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fhess_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     REQUIRE(hess_b.squaredNorm() > 1e-3);
@@ -368,10 +375,11 @@ TEST_CASE("Smooth barrier potential real sim 2D C^1", "[smooth_potential]")
     {
         auto f = [&](const Eigen::VectorXd& x) {
             return potential(
-                collisions, mesh, fd::unflatten(x, vertices.cols()));
+                collisions, mesh, tests::unflatten(x, vertices.cols()));
         };
         fd::finite_gradient(
-            fd::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND, 1e-8);
+            tests::flatten(vertices), f, fgrad_b, fd::AccuracyOrder::SECOND,
+            1e-8);
     }
 
     REQUIRE(grad_b.squaredNorm() > 1e-8);
diff --git a/tests/src/tests/test_collision_mesh.cpp b/tests/src/tests/test_collision_mesh.cpp
index fa60d6328..2077cd74e 100644
--- a/tests/src/tests/test_collision_mesh.cpp
+++ b/tests/src/tests/test_collision_mesh.cpp
@@ -1,6 +1,8 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_approx.hpp>
+#include <catch2/generators/catch_generators.hpp>
 #include <ipc/collision_mesh.hpp>
+#include <ipc/config.hpp>
 
 using namespace ipc;
 
@@ -113,4 +115,73 @@ TEST_CASE("Codim points collision mesh", "[collision_mesh]")
     Eigen::VectorXi expected_codim_vertices(4);
     expected_codim_vertices << 0, 1, 2, 3;
     CHECK(mesh.codim_vertices() == expected_codim_vertices);
+}
+
+TEST_CASE(
+    "vertex_matrix_to_dof_matrix",
+    "[collision_mesh][vertex_matrix_to_dof_matrix]")
+{
+    // Build a small 2×3 vertex-level matrix:
+    //
+    //   M_V = [ 1  0  2 ]
+    //         [ 0  3  0 ]
+    //
+    std::vector<Eigen::Triplet<double>> triplets;
+    triplets.emplace_back(0, 0, 1.0);
+    triplets.emplace_back(0, 2, 2.0);
+    triplets.emplace_back(1, 1, 3.0);
+
+    Eigen::SparseMatrix<double> M_V(2, 3);
+    M_V.setFromTriplets(triplets.begin(), triplets.end());
+
+    const int n_rows = static_cast<int>(M_V.rows()); // 2
+    const int n_cols = static_cast<int>(M_V.cols()); // 3
+
+    const int dim = GENERATE(2, 3);
+    CAPTURE(dim);
+
+    Eigen::SparseMatrix<double> M_dof =
+        CollisionMesh::vertex_matrix_to_dof_matrix(M_V, dim);
+
+    CHECK(M_dof.rows() == n_rows * dim);
+    CHECK(M_dof.cols() == n_cols * dim);
+
+    // Check the sparsity pattern: for each non-zero (r, c, v) in M_V,
+    // it should appear at (global_row(r, d), global_col(c, d)) for
+    // every d in [0, dim).
+    Eigen::MatrixXd M_dense = Eigen::MatrixXd(M_dof);
+    Eigen::MatrixXd expected =
+        Eigen::MatrixXd::Zero(n_rows * dim, n_cols * dim);
+
+    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++) {
+                int gr, gc;
+                if constexpr (VERTEX_DERIVATIVE_LAYOUT == Eigen::RowMajor) {
+                    gr = dim * it.row() + d;
+                    gc = dim * it.col() + d;
+                } else {
+                    gr = n_rows * d + it.row();
+                    gc = n_cols * d + it.col();
+                }
+                expected(gr, gc) = it.value();
+            }
+        }
+    }
+
+    CHECK(M_dense == expected);
+
+    // Verify semantics: M_dof * flatten(V) == flatten(M_V * V)
+    // V is (n_cols × dim) because M_V has n_cols columns, one per vertex.
+    Eigen::MatrixXd V(n_cols, dim); // NOLINT(bugprone-argument-comment)
+    V.setRandom();
+
+    Eigen::VectorXd x_dof = V.reshaped<VERTEX_DERIVATIVE_LAYOUT>();
+    Eigen::VectorXd expected_dof =
+        (M_V * V).reshaped<VERTEX_DERIVATIVE_LAYOUT>();
+
+    Eigen::VectorXd actual_dof = M_dof * x_dof;
+
+    CHECK(actual_dof.isApprox(expected_dof));
 }
\ No newline at end of file
diff --git a/tests/src/tests/utils/CMakeLists.txt b/tests/src/tests/utils/CMakeLists.txt
index 2a547cbd7..e04072a47 100644
--- a/tests/src/tests/utils/CMakeLists.txt
+++ b/tests/src/tests/utils/CMakeLists.txt
@@ -2,6 +2,7 @@ set(SOURCES
   # Tests
   test_utils.cpp
   test_matrixcache.cpp
+  test_local_to_global.cpp
 
   # Benchmarks
 
diff --git a/tests/src/tests/utils/test_local_to_global.cpp b/tests/src/tests/utils/test_local_to_global.cpp
new file mode 100644
index 000000000..fb9098d37
--- /dev/null
+++ b/tests/src/tests/utils/test_local_to_global.cpp
@@ -0,0 +1,596 @@
+#include <ipc/utils/local_to_global.hpp>
+
+#include <catch2/catch_test_macros.hpp>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <vector>
+
+using namespace ipc;
+
+// ─── helpers ────────────────────────────────────────────────────────────────
+
+/// Convert a RowMajor-ordered dense vector to ColMajor ordering.
+///   RowMajor:  [x0 y0 z0  x1 y1 z1  …]
+///   ColMajor:  [x0 x1 …  y0 y1 …  z0 z1 …]
+static Eigen::VectorXd
+row_to_col(const Eigen::VectorXd& v, int n_verts, int dim)
+{
+    return v.reshaped<Eigen::RowMajor>(n_verts, dim)
+        .reshaped<Eigen::ColMajor>();
+}
+
+/// Convert a RowMajor-ordered sparse matrix to ColMajor DOF ordering.
+///   RowMajor rows/cols indexed as (dim*v + d)
+///   ColMajor rows/cols indexed as (n_verts*d + v)
+static Eigen::MatrixXd
+row_to_col_mat(const Eigen::MatrixXd& M, int n_verts, int dim)
+{
+    const int N = n_verts * dim;
+    assert(M.rows() == N && M.cols() == N);
+    Eigen::MatrixXd out = Eigen::MatrixXd::Zero(N, N);
+    for (int vi = 0; vi < n_verts; ++vi) {
+        for (int di = 0; di < dim; ++di) {
+            for (int vj = 0; vj < n_verts; ++vj) {
+                for (int dj = 0; dj < dim; ++dj) {
+                    out(n_verts * di + vi, n_verts * dj + vj) =
+                        M(dim * vi + di, dim * vj + dj);
+                }
+            }
+        }
+    }
+    return out;
+}
+
+/// Same as row_to_col_mat but for a non-square Jacobian where row and column
+/// vertex counts may differ.
+static Eigen::MatrixXd row_to_col_jac(
+    const Eigen::MatrixXd& M, int n_row_verts, int n_col_verts, int dim)
+{
+    const int Nr = n_row_verts * dim;
+    const int Nc = n_col_verts * dim;
+    assert(M.rows() == Nr && M.cols() == Nc);
+    Eigen::MatrixXd out = Eigen::MatrixXd::Zero(Nr, Nc);
+    for (int vi = 0; vi < n_row_verts; ++vi) {
+        for (int di = 0; di < dim; ++di) {
+            for (int vj = 0; vj < n_col_verts; ++vj) {
+                for (int dj = 0; dj < dim; ++dj) {
+                    out(n_row_verts * di + vi, n_col_verts * dj + vj) =
+                        M(dim * vi + di, dim * vj + dj);
+                }
+            }
+        }
+    }
+    return out;
+}
+
+/// Assemble a dense matrix from a vector of triplets.
+static Eigen::MatrixXd triplets_to_dense(
+    const std::vector<Eigen::Triplet<double>>& trips, int rows, int cols)
+{
+    Eigen::SparseMatrix<double> sp(rows, cols);
+    sp.setFromTriplets(trips.begin(), trips.end());
+    return Eigen::MatrixXd(sp);
+}
+
+// ─── Test parameters ────────────────────────────────────────────────────────
+
+// Global mesh: 5 vertices in 3D  → 15 DOFs.
+// Local element touches vertices {1, 3} (2 verts, 6 local DOFs).
+static constexpr int DIM = 3;
+static constexpr int N_TOTAL_VERTS = 5;
+static constexpr int N_DOF = N_TOTAL_VERTS * DIM;
+
+static const std::vector<long> ELEM_IDS = { { 1, 3 } };
+static constexpr int N_LOCAL_VERTS = 2;
+static constexpr int N_LOCAL_DOF = N_LOCAL_VERTS * DIM;
+
+/// Deterministic local gradient (6 entries).
+static Eigen::VectorXd make_local_grad()
+{
+    Eigen::VectorXd g(N_LOCAL_DOF);
+    // RowMajor local ordering: [x1 y1 z1 x3 y3 z3]
+    g << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0;
+    return g;
+}
+
+/// Deterministic local hessian (6×6).
+static Eigen::MatrixXd make_local_hessian()
+{
+    Eigen::MatrixXd H = Eigen::MatrixXd::Zero(N_LOCAL_DOF, N_LOCAL_DOF);
+    for (int i = 0; i < N_LOCAL_DOF; ++i) {
+        for (int j = 0; j < N_LOCAL_DOF; ++j) {
+            H(i, j) = 10.0 * (i + 1) + (j + 1); // nonzero everywhere
+        }
+    }
+    return H;
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  1. local_gradient_to_global_gradient – dense Eigen::VectorXd
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE(
+    "local_gradient_to_global_gradient dense RowMajor vs ColMajor",
+    "[utils][local_to_global][gradient][dense]")
+{
+    const Eigen::VectorXd local_grad = make_local_grad();
+
+    // --- RowMajor assembly ---
+    Eigen::VectorXd grad_row = Eigen::VectorXd::Zero(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_grad, ELEM_IDS, DIM, grad_row);
+
+    // Verify a few values manually.
+    // Vertex 1, dim 0 → index 3*1+0 = 3
+    CHECK(grad_row[3] == 1.0);
+    // Vertex 1, dim 1 → index 3*1+1 = 4
+    CHECK(grad_row[4] == 2.0);
+    // Vertex 3, dim 2 → index 3*3+2 = 11
+    CHECK(grad_row[11] == 6.0);
+
+    // --- ColMajor assembly ---
+    Eigen::VectorXd grad_col = Eigen::VectorXd::Zero(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_grad, ELEM_IDS, DIM, grad_col);
+
+    // ColMajor: vertex v, dim d → index N_TOTAL_VERTS*d + v
+    // Vertex 1, dim 0 → index 5*0+1 = 1
+    CHECK(grad_col[1] == 1.0);
+    // Vertex 1, dim 1 → index 5*1+1 = 6
+    CHECK(grad_col[6] == 2.0);
+    // Vertex 3, dim 2 → index 5*2+3 = 13
+    CHECK(grad_col[13] == 6.0);
+
+    // Structural check: both vectors should be the permutation of each other.
+    Eigen::VectorXd expected_col = row_to_col(grad_row, N_TOTAL_VERTS, DIM);
+    CHECK((grad_col.array() == expected_col.array()).all());
+}
+
+TEST_CASE(
+    "local_gradient_to_global_gradient dense accumulates from two elements",
+    "[utils][local_to_global][gradient][dense]")
+{
+    // Two elements that share vertex 3.
+    const std::vector<long> ids_a = { 0, 3 };
+    const std::vector<long> ids_b = { 3, 4 };
+
+    Eigen::VectorXd local_a(N_LOCAL_DOF), local_b(N_LOCAL_DOF);
+    local_a << 1, 2, 3, 10, 20, 30;
+    local_b << 100, 200, 300, 7, 8, 9;
+
+    Eigen::VectorXd grad_row = Eigen::VectorXd::Zero(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_a, ids_a, DIM, grad_row);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_b, ids_b, DIM, grad_row);
+
+    Eigen::VectorXd grad_col = Eigen::VectorXd::Zero(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_a, ids_a, DIM, grad_col);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_b, ids_b, DIM, grad_col);
+
+    // Vertex 3 receives contributions from both elements.
+    // RowMajor vertex 3 dim 0 → index 9: 10 + 100 = 110
+    CHECK(grad_row[9] == 110.0);
+    // ColMajor vertex 3 dim 0 → index 5*0+3 = 3: should also be 110
+    CHECK(grad_col[3] == 110.0);
+
+    Eigen::VectorXd expected_col = row_to_col(grad_row, N_TOTAL_VERTS, DIM);
+    CHECK((grad_col.array() == expected_col.array()).all());
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  2. local_gradient_to_global_gradient – sparse Eigen::SparseVector
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE(
+    "local_gradient_to_global_gradient sparse RowMajor vs ColMajor",
+    "[utils][local_to_global][gradient][sparse]")
+{
+    const Eigen::VectorXd local_grad = make_local_grad();
+
+    Eigen::SparseVector<double> grad_row(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_grad, ELEM_IDS, DIM, grad_row);
+
+    Eigen::SparseVector<double> grad_col(N_DOF);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_grad, ELEM_IDS, DIM, grad_col);
+
+    // Convert both to dense for easy comparison.
+    Eigen::VectorXd dr = Eigen::VectorXd(grad_row);
+    Eigen::VectorXd dc = Eigen::VectorXd(grad_col);
+
+    Eigen::VectorXd expected_col = row_to_col(dr, N_TOTAL_VERTS, DIM);
+    CHECK((dc.array() == expected_col.array()).all());
+
+    // Spot-check: only 2*DIM entries should be nonzero.
+    CHECK(grad_row.nonZeros() == N_LOCAL_DOF);
+    CHECK(grad_col.nonZeros() == N_LOCAL_DOF);
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  3. local_hessian_to_global_triplets – std::vector<Triplet>
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE(
+    "local_hessian_to_global_triplets (triplet vector) RowMajor vs ColMajor",
+    "[utils][local_to_global][hessian][triplet]")
+{
+    const Eigen::MatrixXd local_hess = make_local_hessian();
+
+    // --- RowMajor ---
+    std::vector<Eigen::Triplet<double>> trips_row;
+    local_hessian_to_global_triplets<std::vector<long>, Eigen::RowMajor>(
+        local_hess, ELEM_IDS, DIM, trips_row);
+
+    // --- ColMajor ---
+    std::vector<Eigen::Triplet<double>> trips_col;
+    local_hessian_to_global_triplets<std::vector<long>, Eigen::ColMajor>(
+        local_hess, ELEM_IDS, DIM, trips_col, N_TOTAL_VERTS);
+
+    // Both should produce N_LOCAL_DOF^2 triplets.
+    CHECK(trips_row.size() == N_LOCAL_DOF * N_LOCAL_DOF);
+    CHECK(trips_col.size() == N_LOCAL_DOF * N_LOCAL_DOF);
+
+    Eigen::MatrixXd H_row = triplets_to_dense(trips_row, N_DOF, N_DOF);
+    Eigen::MatrixXd H_col = triplets_to_dense(trips_col, N_DOF, N_DOF);
+
+    Eigen::MatrixXd expected_col = row_to_col_mat(H_row, N_TOTAL_VERTS, DIM);
+    CHECK((H_col.array() == expected_col.array()).all());
+}
+
+TEST_CASE(
+    "local_hessian_to_global_triplets (triplet vector) ColMajor spot-check",
+    "[utils][local_to_global][hessian][triplet]")
+{
+    // 2 verts, dim=2, ids={0,2}, n_total=4
+    const int dim = 2;
+    const int n_total = 4;
+    const std::vector<long> ids = { 0, 2 };
+
+    // 4×4 local hessian
+    Eigen::MatrixXd H(4, 4);
+    H << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;
+
+    std::vector<Eigen::Triplet<double>> trips;
+    local_hessian_to_global_triplets<std::vector<long>, Eigen::ColMajor>(
+        H, ids, dim, trips, n_total);
+
+    Eigen::MatrixXd M = triplets_to_dense(trips, n_total * dim, n_total * dim);
+
+    // RowMajor local (i=0,d=0) → global ColMajor row = n_total*0 + ids[0] = 0
+    // RowMajor local (i=0,d=1) → global ColMajor row = n_total*1 + ids[0] = 4
+    // RowMajor local (i=1,d=0) → global ColMajor row = n_total*0 + ids[1] = 2
+    // RowMajor local (i=1,d=1) → global ColMajor row = n_total*1 + ids[1] = 6
+
+    // H(0,0) = local(dim*0+0, dim*0+0) = H(0,0)=1, row=0 col=0
+    CHECK(M(0, 0) == 1.0);
+    // H(0,1) = local(0,1) = 2, row=0, col=n_total*1+0=4
+    CHECK(M(0, 4) == 2.0);
+    // H(1,0) = local(1,0) = 5, row=n_total*1+0=4, col=0
+    CHECK(M(4, 0) == 5.0);
+    // H(2,3) = local(dim*1+0, dim*1+1) = H(2,3) = 12
+    //   row = n_total*0+ids[1] = 2, col = n_total*1+ids[1] = 6
+    CHECK(M(2, 6) == 12.0);
+    // H(3,2) = local(dim*1+1, dim*1+0) = H(3,2) = 15
+    //   row = n_total*1+ids[1] = 6, col = n_total*0+ids[1] = 2
+    CHECK(M(6, 2) == 15.0);
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  4. local_hessian_to_global_triplets – MatrixCache
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE(
+    "local_hessian_to_global_triplets (MatrixCache) RowMajor vs ColMajor",
+    "[utils][local_to_global][hessian][cache]")
+{
+    const Eigen::MatrixXd local_hess = make_local_hessian();
+
+    // --- RowMajor via MatrixCache ---
+    SparseMatrixCache cache_row(N_DOF, N_DOF);
+    cache_row.reserve(N_LOCAL_DOF * N_LOCAL_DOF);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::RowMajor>(
+        local_hess, ELEM_IDS, DIM, cache_row);
+
+    Eigen::MatrixXd H_row = Eigen::MatrixXd(cache_row.get_matrix());
+
+    // --- ColMajor via MatrixCache ---
+    SparseMatrixCache cache_col(N_DOF, N_DOF);
+    cache_col.reserve(N_LOCAL_DOF * N_LOCAL_DOF);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::ColMajor>(
+        local_hess, ELEM_IDS, DIM, cache_col, N_TOTAL_VERTS);
+
+    Eigen::MatrixXd H_col = Eigen::MatrixXd(cache_col.get_matrix());
+
+    Eigen::MatrixXd expected_col = row_to_col_mat(H_row, N_TOTAL_VERTS, DIM);
+    CHECK((H_col.array() == expected_col.array()).all());
+}
+
+TEST_CASE(
+    "local_hessian_to_global_triplets (MatrixCache) ColMajor with two elements",
+    "[utils][local_to_global][hessian][cache]")
+{
+    // Two elements that share vertex 2.
+    const int dim = 2;
+    const int n_total = 4;
+    const int n_dof = n_total * dim;
+    const std::vector<long> ids_a = { 0, 2 };
+    const std::vector<long> ids_b = { 2, 3 };
+
+    Eigen::MatrixXd Ha = Eigen::MatrixXd::Ones(4, 4) * 2.0;
+    Eigen::MatrixXd Hb = Eigen::MatrixXd::Ones(4, 4) * 3.0;
+
+    // RowMajor
+    SparseMatrixCache cr(n_dof, n_dof);
+    cr.reserve(32);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::RowMajor>(
+        Ha, ids_a, dim, cr);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::RowMajor>(
+        Hb, ids_b, dim, cr);
+    Eigen::MatrixXd M_row = Eigen::MatrixXd(cr.get_matrix());
+
+    // ColMajor
+    SparseMatrixCache cc(n_dof, n_dof);
+    cc.reserve(32);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::ColMajor>(
+        Ha, ids_a, dim, cc, n_total);
+    local_hessian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, Eigen::ColMajor>(
+        Hb, ids_b, dim, cc, n_total);
+    Eigen::MatrixXd M_col = Eigen::MatrixXd(cc.get_matrix());
+
+    Eigen::MatrixXd expected = row_to_col_mat(M_row, n_total, dim);
+    CHECK((M_col.array() == expected.array()).all());
+
+    // The (2,2) block in RowMajor should accumulate 2+3=5 for the
+    // diagonal dim entries of vertex 2.
+    CHECK(M_row(dim * 2, dim * 2) == 5.0);
+    CHECK(M_col(n_total * 0 + 2, n_total * 0 + 2) == 5.0);
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  5. local_jacobian_to_global_triplets
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE(
+    "local_jacobian_to_global_triplets RowMajor vs ColMajor square",
+    "[utils][local_to_global][jacobian]")
+{
+    // Square case (same as hessian but using the jacobian API).
+    const Eigen::MatrixXd local_jac = make_local_hessian();
+
+    std::vector<Eigen::Triplet<double>> trips_row;
+    local_jacobian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, std::vector<long>, Eigen::RowMajor>(
+        local_jac, ELEM_IDS, ELEM_IDS, DIM, trips_row);
+
+    std::vector<Eigen::Triplet<double>> trips_col;
+    local_jacobian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, std::vector<long>, Eigen::ColMajor>(
+        local_jac, ELEM_IDS, ELEM_IDS, DIM, trips_col, N_TOTAL_VERTS,
+        N_TOTAL_VERTS);
+
+    Eigen::MatrixXd J_row = triplets_to_dense(trips_row, N_DOF, N_DOF);
+    Eigen::MatrixXd J_col = triplets_to_dense(trips_col, N_DOF, N_DOF);
+
+    Eigen::MatrixXd expected = row_to_col_mat(J_row, N_TOTAL_VERTS, DIM);
+    CHECK((J_col.array() == expected.array()).all());
+}
+
+TEST_CASE(
+    "local_jacobian_to_global_triplets ColMajor non-square",
+    "[utils][local_to_global][jacobian]")
+{
+    // Non-square: 2 row-verts (ids {0,3}), 3 col-verts (ids {1,2,4}), dim=2.
+    const int dim = 2;
+    const int n_row_total = 5; // global row vertices
+    const int n_col_total = 5; // global col vertices
+    const std::vector<long> row_ids = { 0, 3 };
+    const std::vector<long> col_ids = { 1, 2, 4 };
+    const int n_row_local = 2;
+    const int n_col_local = 3;
+    const int lr = n_row_local * dim; // 4
+    const int lc = n_col_local * dim; // 6
+
+    Eigen::MatrixXd J(lr, lc);
+    for (int i = 0; i < lr; ++i) {
+        for (int j = 0; j < lc; ++j) {
+            J(i, j) = 100.0 * (i + 1) + (j + 1);
+        }
+    }
+
+    // RowMajor
+    std::vector<Eigen::Triplet<double>> trips_row;
+    local_jacobian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, std::vector<long>, Eigen::RowMajor>(
+        J, row_ids, col_ids, dim, trips_row);
+
+    // ColMajor
+    std::vector<Eigen::Triplet<double>> trips_col;
+    local_jacobian_to_global_triplets<
+        Eigen::MatrixXd, std::vector<long>, std::vector<long>, Eigen::ColMajor>(
+        J, row_ids, col_ids, dim, trips_col, n_row_total, n_col_total);
+
+    const int Nr = n_row_total * dim;
+    const int Nc = n_col_total * dim;
+    Eigen::MatrixXd J_row = triplets_to_dense(trips_row, Nr, Nc);
+    Eigen::MatrixXd J_col = triplets_to_dense(trips_col, Nr, Nc);
+
+    Eigen::MatrixXd expected =
+        row_to_col_jac(J_row, n_row_total, n_col_total, dim);
+    CHECK((J_col.array() == expected.array()).all());
+
+    // Spot-check a single entry.
+    // J(0,0) = 101.  RowMajor: row=dim*row_ids[0]+0=0, col=dim*col_ids[0]+0=2
+    CHECK(J_row(0, 2) == 101.0);
+    // ColMajor: row=n_row_total*0+row_ids[0]=0, col=n_col_total*0+col_ids[0]=1
+    CHECK(J_col(0, 1) == 101.0);
+}
+
+// ═══════════════════════════════════════════════════════════════════════════
+//  6. Edge cases
+// ═══════════════════════════════════════════════════════════════════════════
+
+TEST_CASE("ColMajor gradient 2D", "[utils][local_to_global][gradient][2D]")
+{
+    // 2D problem: 4 vertices, dim=2.
+    const int dim = 2;
+    const int n_total = 4;
+    const int n_dof = n_total * dim; // 8
+    const std::vector<long> ids = { 1, 3 };
+    Eigen::VectorXd local_grad(4);
+    local_grad << 10, 20, 30, 40; // [x1 y1 x3 y3]
+
+    // RowMajor
+    Eigen::VectorXd gr = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_grad, ids, dim, gr);
+    // Expected RowMajor: index 2*1+0=2 →10, 2*1+1=3 →20, 2*3+0=6 →30, 2*3+1=7
+    // →40
+    CHECK(gr[2] == 10.0);
+    CHECK(gr[3] == 20.0);
+    CHECK(gr[6] == 30.0);
+    CHECK(gr[7] == 40.0);
+
+    // ColMajor
+    Eigen::VectorXd gc = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_grad, ids, dim, gc);
+    // Expected ColMajor: n_total*0+1=1 →10, n_total*1+1=5 →20,
+    //                    n_total*0+3=3 →30, n_total*1+3=7 →40
+    CHECK(gc[1] == 10.0);
+    CHECK(gc[5] == 20.0);
+    CHECK(gc[3] == 30.0);
+    CHECK(gc[7] == 40.0);
+
+    Eigen::VectorXd expected = row_to_col(gr, n_total, dim);
+    CHECK((gc.array() == expected.array()).all());
+}
+
+TEST_CASE(
+    "ColMajor single vertex element",
+    "[utils][local_to_global][gradient][single]")
+{
+    // A single vertex local element.
+    const int dim = 3;
+    const int n_total = 3;
+    const int n_dof = n_total * dim;
+    const std::vector<long> ids = { 2 };
+    Eigen::VectorXd local_grad(3);
+    local_grad << 7, 8, 9;
+
+    Eigen::VectorXd gr = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_grad, ids, dim, gr);
+
+    Eigen::VectorXd gc = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_grad, ids, dim, gc);
+
+    Eigen::VectorXd expected = row_to_col(gr, n_total, dim);
+    CHECK((gc.array() == expected.array()).all());
+
+    // RowMajor: vertex 2, dim 0 → 3*2+0=6
+    CHECK(gr[6] == 7.0);
+    // ColMajor: vertex 2, dim 0 → 3*0+2=2
+    CHECK(gc[2] == 7.0);
+}
+
+TEST_CASE(
+    "ColMajor hessian identity-like local matrix",
+    "[utils][local_to_global][hessian][identity]")
+{
+    // Local identity hessian to verify diagonal placement.
+    const int dim = 2;
+    const int n_total = 4;
+    const int n_dof = n_total * dim;
+    const std::vector<long> ids = { 0, 3 };
+
+    Eigen::MatrixXd H = Eigen::MatrixXd::Identity(4, 4);
+
+    std::vector<Eigen::Triplet<double>> trips_row;
+    local_hessian_to_global_triplets<std::vector<long>, Eigen::RowMajor>(
+        H, ids, dim, trips_row);
+
+    std::vector<Eigen::Triplet<double>> trips_col;
+    local_hessian_to_global_triplets<std::vector<long>, Eigen::ColMajor>(
+        H, ids, dim, trips_col, n_total);
+
+    Eigen::MatrixXd M_row = triplets_to_dense(trips_row, n_dof, n_dof);
+    Eigen::MatrixXd M_col = triplets_to_dense(trips_col, n_dof, n_dof);
+
+    // RowMajor diagonal entries at (0,0),(1,1),(6,6),(7,7).
+    CHECK(M_row(0, 0) == 1.0);
+    CHECK(M_row(1, 1) == 1.0);
+    CHECK(M_row(6, 6) == 1.0);
+    CHECK(M_row(7, 7) == 1.0);
+
+    // ColMajor diagonal entries at (0,0),(4,4),(3,3),(7,7).
+    CHECK(M_col(0, 0) == 1.0);
+    CHECK(M_col(4, 4) == 1.0);
+    CHECK(M_col(3, 3) == 1.0);
+    CHECK(M_col(7, 7) == 1.0);
+
+    Eigen::MatrixXd expected = row_to_col_mat(M_row, n_total, dim);
+    CHECK((M_col.array() == expected.array()).all());
+}
+
+TEST_CASE(
+    "ColMajor with extra ids in container",
+    "[utils][local_to_global][gradient][extra_ids]")
+{
+    // ids container is larger than needed – only first n_verts are used.
+    const int dim = 2;
+    const int n_total = 5;
+    const int n_dof = n_total * dim;
+    const std::vector<long> ids = { 1, 4, 99 }; // 3 ids, but local only needs 2
+
+    Eigen::VectorXd local_grad(4); // 2 verts * 2 dim
+    local_grad << 5, 6, 7, 8;
+
+    Eigen::VectorXd gr = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::RowMajor>(
+        local_grad, ids, dim, gr);
+
+    Eigen::VectorXd gc = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<std::vector<long>, Eigen::ColMajor>(
+        local_grad, ids, dim, gc);
+
+    Eigen::VectorXd expected = row_to_col(gr, n_total, dim);
+    CHECK((gc.array() == expected.array()).all());
+}
+
+TEST_CASE(
+    "ColMajor with Eigen::VectorXi ids",
+    "[utils][local_to_global][gradient][eigen_ids]")
+{
+    // Ensure IDContainer works with Eigen types.
+    const int dim = 3;
+    const int n_total = 4;
+    const int n_dof = n_total * dim;
+    Eigen::VectorXi ids(2);
+    ids << 0, 3;
+
+    Eigen::VectorXd local_grad(6);
+    local_grad << 1, 2, 3, 4, 5, 6;
+
+    Eigen::VectorXd gr = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<Eigen::VectorXi, Eigen::RowMajor>(
+        local_grad, ids, dim, gr);
+
+    Eigen::VectorXd gc = Eigen::VectorXd::Zero(n_dof);
+    local_gradient_to_global_gradient<Eigen::VectorXi, Eigen::ColMajor>(
+        local_grad, ids, dim, gc);
+
+    Eigen::VectorXd expected = row_to_col(gr, n_total, dim);
+    CHECK((gc.array() == expected.array()).all());
+}