diff --git a/src/ipc/collisions/tangential/tangential_collisions.cpp b/src/ipc/collisions/tangential/tangential_collisions.cpp
index 6eb77115c..2034cc174 100644
--- a/src/ipc/collisions/tangential/tangential_collisions.cpp
+++ b/src/ipc/collisions/tangential/tangential_collisions.cpp
@@ -116,7 +116,7 @@ void TangentialCollisions::build(
 
 void TangentialCollisions::build(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const SmoothCollisions& collisions,
     const SmoothContactParameters& params,
     const double normal_stiffness,
diff --git a/src/ipc/collisions/tangential/tangential_collisions.hpp b/src/ipc/collisions/tangential/tangential_collisions.hpp
index b117fc9f4..9f6c13196 100644
--- a/src/ipc/collisions/tangential/tangential_collisions.hpp
+++ b/src/ipc/collisions/tangential/tangential_collisions.hpp
@@ -97,7 +97,7 @@ class TangentialCollisions {
     /// @param blend_mu Function to blend vertex-based coefficients of friction. Defaults to average.
     void build(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const SmoothCollisions& collisions,
         const SmoothContactParameters& params,
         const double normal_stiffness,
diff --git a/src/ipc/math/math.cpp b/src/ipc/math/math.cpp
index 3d7954bbf..175ba657b 100644
--- a/src/ipc/math/math.cpp
+++ b/src/ipc/math/math.cpp
@@ -49,7 +49,7 @@ double opposite_direction_penalty(
     return Math<double>::smooth_heaviside(d.dot(t) / t.norm(), alpha, beta);
 }
 
-GradType<6> opposite_direction_penalty_grad(
+GradientType<6> opposite_direction_penalty_grad(
     Eigen::ConstRef<Eigen::Vector3d> t,
     Eigen::ConstRef<Eigen::Vector3d> d,
     const double alpha,
@@ -109,7 +109,7 @@ double negative_orientation_penalty(
     return opposite_direction_penalty(n, d, alpha, beta);
 }
 
-GradType<9> negative_orientation_penalty_grad(
+GradientType<9> negative_orientation_penalty_grad(
     Eigen::ConstRef<Eigen::Vector3d> t1,
     Eigen::ConstRef<Eigen::Vector3d> t2,
     Eigen::ConstRef<Eigen::Vector3d> d,
diff --git a/src/ipc/math/math.hpp b/src/ipc/math/math.hpp
index 9ff0ce769..0811e825b 100644
--- a/src/ipc/math/math.hpp
+++ b/src/ipc/math/math.hpp
@@ -119,7 +119,7 @@ double opposite_direction_penalty(
     const double alpha,
     const double beta);
 
-GradType<6> opposite_direction_penalty_grad(
+GradientType<6> opposite_direction_penalty_grad(
     Eigen::ConstRef<Eigen::Vector3d> t,
     Eigen::ConstRef<Eigen::Vector3d> d,
     const double alpha,
@@ -139,7 +139,7 @@ double negative_orientation_penalty(
     const double alpha,
     const double beta);
 
-GradType<9> negative_orientation_penalty_grad(
+GradientType<9> negative_orientation_penalty_grad(
     Eigen::ConstRef<Eigen::Vector3d> t1,
     Eigen::ConstRef<Eigen::Vector3d> t2,
     Eigen::ConstRef<Eigen::Vector3d> d,
diff --git a/src/ipc/potentials/potential.hpp b/src/ipc/potentials/potential.hpp
index 3b9d815a9..baf842653 100644
--- a/src/ipc/potentials/potential.hpp
+++ b/src/ipc/potentials/potential.hpp
@@ -12,7 +12,7 @@ template <class TCollisions> class Potential {
     using TCollision = typename TCollisions::value_type;
     /// @brief Maximum degrees of freedom per collision
     static constexpr int STENCIL_NDOF = 3 * TCollision::STENCIL_SIZE;
-    using VectorMaxNd = Vector<double, Eigen::Dynamic, STENCIL_NDOF>;
+    using VectorMaxNd = VectorMax<double, STENCIL_NDOF>;
     using MatrixMaxNd = MatrixMax<double, STENCIL_NDOF, STENCIL_NDOF>;
 
 public:
diff --git a/src/ipc/potentials/tangential_potential.cpp b/src/ipc/potentials/tangential_potential.cpp
index 798a9de28..f1fef39c3 100644
--- a/src/ipc/potentials/tangential_potential.cpp
+++ b/src/ipc/potentials/tangential_potential.cpp
@@ -683,7 +683,7 @@ TangentialPotential::VectorMaxNd TangentialPotential::smooth_contact_force(
     assert(rest_positions.size() == velocities.size());
 
     // const VectorMaxNd x = dof(rest_positions, edges, faces);
-    // const Vector<double, -1, STENCIL_NDOF> u =
+    // const VectorMax<double, STENCIL_NDOF> u =
     //     dof(lagged_displacements, edges, faces);
     // const VectorMaxNd v = dof(velocities, edges, faces);
     const VectorMaxNd lagged_positions =
diff --git a/src/ipc/smooth_contact/collisions/smooth_collision.cpp b/src/ipc/smooth_contact/collisions/smooth_collision.cpp
index 668cad2d7..5f569f09e 100644
--- a/src/ipc/smooth_contact/collisions/smooth_collision.cpp
+++ b/src/ipc/smooth_contact/collisions/smooth_collision.cpp
@@ -42,9 +42,9 @@ Eigen::VectorXd SmoothCollision::dof(Eigen::ConstRef<Eigen::MatrixXd> X) const
 
 template <typename PrimitiveA, typename PrimitiveB>
 auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::get_core_indices() const
-    -> Vector<int, N_CORE_DOFS>
+    -> Eigen::Vector<int, N_CORE_DOFS>
 {
-    Vector<int, N_CORE_DOFS> core_indices;
+    Eigen::Vector<int, N_CORE_DOFS> core_indices;
     core_indices << Eigen::VectorXi::LinSpaced(
         N_CORE_DOFS_A, 0, N_CORE_DOFS_A - 1),
         Eigen::VectorXi::LinSpaced(
@@ -61,7 +61,7 @@ SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::SmoothCollisionTemplate(
     const CollisionMesh& mesh,
     const SmoothContactParameters& params,
     const double _dhat,
-    const Eigen::MatrixXd& V)
+    Eigen::ConstRef<Eigen::MatrixXd> V)
     : SmoothCollision(_primitive0, _primitive1, _dhat, mesh)
 {
     VectorMax3d d =
@@ -104,16 +104,16 @@ SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::SmoothCollisionTemplate(
 
 template <typename PrimitiveA, typename PrimitiveB>
 double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::operator()(
-    Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+    Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
     const SmoothContactParameters& params) const
 {
-    Vector<double, N_CORE_POINTS * DIM> x;
+    Eigen::Vector<double, N_CORE_POINTS * DIM> x;
     x << positions.head(PrimitiveA::N_CORE_POINTS * DIM),
         positions.segment(
             primitive_a->n_dofs(), PrimitiveB::N_CORE_POINTS * DIM);
 
     // grad of "d" wrt. points
-    const Vector<double, DIM> closest_direction =
+    const Eigen::Vector<double, DIM> closest_direction =
         PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, double>::
             compute_closest_direction(x, DTYPE::AUTO);
     const double dist = closest_direction.norm();
@@ -143,19 +143,19 @@ double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::operator()(
 
 template <typename PrimitiveA, typename PrimitiveB>
 auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
-    Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+    Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
     const SmoothContactParameters& params) const
-    -> Vector<double, -1, ELEMENT_SIZE>
+    -> VectorMax<double, ELEMENT_SIZE>
 {
     const auto core_indices = get_core_indices();
 
-    Vector<double, N_CORE_DOFS> x;
+    Eigen::Vector<double, N_CORE_DOFS> x;
     x = positions(core_indices);
 
     const auto dtype =
         PrimitiveDistance<PrimitiveA, PrimitiveB>::compute_distance_type(x);
 
-    Vector<double, DIM> closest_direction;
+    Eigen::Vector<double, DIM> closest_direction;
     Eigen::Matrix<double, DIM, N_CORE_DOFS> closest_direction_grad;
     std::tie(closest_direction, closest_direction_grad) = PrimitiveDistance<
         PrimitiveA, PrimitiveB>::compute_closest_direction_gradient(x, dtype);
@@ -171,15 +171,16 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
 
     // gradient of barrier potential
     double barrier = 0;
-    Vector<double, N_CORE_DOFS> gBarrier = Vector<double, N_CORE_DOFS>::Zero();
+    Eigen::Vector<double, N_CORE_DOFS> gBarrier =
+        Eigen::Vector<double, N_CORE_DOFS>::Zero();
     {
         barrier = Math<double>::inv_barrier(dist / dhat(), params.r);
 
-        const Vector<double, DIM> closest_direction_normalized =
+        const Eigen::Vector<double, DIM> closest_direction_normalized =
             closest_direction / dist;
         const double barrier_1st_deriv =
             Math<double>::inv_barrier_grad(dist / dhat(), params.r) / dhat();
-        const Vector<double, DIM> gBarrier_wrt_d =
+        const Eigen::Vector<double, DIM> gBarrier_wrt_d =
             barrier_1st_deriv * closest_direction_normalized;
         gBarrier = closest_direction_grad.transpose() * gBarrier_wrt_d;
     }
@@ -187,14 +188,16 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
     // gradient of mollifier
     {
         double mollifier = 0;
-        Vector<double, N_CORE_DOFS> gMollifier =
-            Vector<double, N_CORE_DOFS>::Zero();
+        Eigen::Vector<double, N_CORE_DOFS> gMollifier =
+            Eigen::Vector<double, N_CORE_DOFS>::Zero();
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADGrad<N_CORE_DOFS>;
-        Vector<T, N_CORE_DOFS> xAD = slice_positions<T, N_CORE_DOFS, 1>(x);
-        Vector<T, DIM> closest_direction_autodiff = PrimitiveDistanceTemplate<
-            PrimitiveA, PrimitiveB, T>::compute_closest_direction(xAD, dtype);
+        Eigen::Vector<T, N_CORE_DOFS> xAD =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
+        Eigen::Vector<T, DIM> closest_direction_autodiff =
+            PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::
+                compute_closest_direction(xAD, dtype);
         const auto dist_sqr_AD = closest_direction_autodiff.squaredNorm();
         auto mollifier_autodiff =
             PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
@@ -202,11 +205,11 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
         mollifier = mollifier_autodiff.val;
         gMollifier = mollifier_autodiff.grad;
 #else
-        Vector<double, N_CORE_DOFS + 1> mollifier_grad;
+        Eigen::Vector<double, N_CORE_DOFS + 1> mollifier_grad;
         std::tie(mollifier, mollifier_grad) = PrimitiveDistance<
             PrimitiveA, PrimitiveB>::compute_mollifier_gradient(x, dist * dist);
 
-        const Vector<double, N_CORE_DOFS> dist_sqr_grad =
+        const Eigen::Vector<double, N_CORE_DOFS> dist_sqr_grad =
             2 * closest_direction_grad.transpose() * closest_direction;
         mollifier_grad.head(N_CORE_DOFS) +=
             mollifier_grad(N_CORE_DOFS) * dist_sqr_grad;
@@ -219,11 +222,11 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
 
     // grad of tangent/normal terms
     double orient = 0;
-    Vector<double, -1, ELEMENT_SIZE> gOrient;
+    VectorMax<double, ELEMENT_SIZE> gOrient;
     {
-        Vector<double, -1, ELEMENT_SIZE>
-            gA = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs()),
-            gB = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs());
+        VectorMax<double, ELEMENT_SIZE>
+            gA = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs()),
+            gB = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs());
         {
             gA(core_indices) =
                 closest_direction_grad.transpose() * gA_reduced.head(DIM);
@@ -254,19 +257,19 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::gradient(
 
 template <typename PrimitiveA, typename PrimitiveB>
 auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
-    Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+    Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
     const SmoothContactParameters& params) const
     -> MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>
 {
     const auto core_indices = get_core_indices();
 
-    Vector<double, N_CORE_DOFS> x;
+    Eigen::Vector<double, N_CORE_DOFS> x;
     x = positions(core_indices);
 
     const auto dtype =
         PrimitiveDistance<PrimitiveA, PrimitiveB>::compute_distance_type(x);
 
-    Vector<double, DIM> closest_direction;
+    Eigen::Vector<double, DIM> closest_direction;
     Eigen::Matrix<double, DIM, N_CORE_DOFS> closest_direction_grad;
     std::array<Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>, DIM>
         closest_direction_hess;
@@ -289,17 +292,18 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
 
     // hessian of barrier potential
     double barrier = 0;
-    Vector<double, N_CORE_DOFS> gBarrier = Vector<double, N_CORE_DOFS>::Zero();
+    Eigen::Vector<double, N_CORE_DOFS> gBarrier =
+        Eigen::Vector<double, N_CORE_DOFS>::Zero();
     Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS> hBarrier =
         Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>::Zero();
     {
         barrier = Math<double>::inv_barrier(dist / dhat(), params.r);
 
-        const Vector<double, DIM> closest_direction_normalized =
+        const Eigen::Vector<double, DIM> closest_direction_normalized =
             closest_direction / dist;
         const double barrier_1st_deriv =
             Math<double>::inv_barrier_grad(dist / dhat(), params.r) / dhat();
-        const Vector<double, DIM> gBarrier_wrt_d =
+        const Eigen::Vector<double, DIM> gBarrier_wrt_d =
             barrier_1st_deriv * closest_direction_normalized;
         gBarrier = closest_direction_grad.transpose() * gBarrier_wrt_d;
 
@@ -322,16 +326,18 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
     // hessian of mollifier
     {
         double mollifier = 0;
-        Vector<double, N_CORE_DOFS> gMollifier =
-            Vector<double, N_CORE_DOFS>::Zero();
+        Eigen::Vector<double, N_CORE_DOFS> gMollifier =
+            Eigen::Vector<double, N_CORE_DOFS>::Zero();
         Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS> hMollifier =
             Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>::Zero();
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADHessian<N_CORE_DOFS>;
-        Vector<T, N_CORE_DOFS> xAD = slice_positions<T, N_CORE_DOFS, 1>(x);
-        Vector<T, DIM> closest_direction_autodiff = PrimitiveDistanceTemplate<
-            PrimitiveA, PrimitiveB, T>::compute_closest_direction(xAD, dtype);
+        Eigen::Vector<T, N_CORE_DOFS> xAD =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
+        Eigen::Vector<T, DIM> closest_direction_autodiff =
+            PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::
+                compute_closest_direction(xAD, dtype);
         const auto dist_sqr_AD = closest_direction_autodiff.squaredNorm();
         auto mollifier_autodiff =
             PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
@@ -341,12 +347,12 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
         gMollifier = mollifier_autodiff.grad;
         hMollifier = mollifier_autodiff.Hess;
 #else
-        Vector<double, N_CORE_DOFS + 1> mollifier_grad;
+        Eigen::Vector<double, N_CORE_DOFS + 1> mollifier_grad;
         Eigen::Matrix<double, N_CORE_DOFS + 1, N_CORE_DOFS + 1> mollifier_hess;
         std::tie(mollifier, mollifier_grad, mollifier_hess) = PrimitiveDistance<
             PrimitiveA, PrimitiveB>::compute_mollifier_hessian(x, dist * dist);
 
-        const Vector<double, N_CORE_DOFS> dist_sqr_grad =
+        const Eigen::Vector<double, N_CORE_DOFS> dist_sqr_grad =
             2 * closest_direction_grad.transpose() * closest_direction;
         mollifier_grad.head(N_CORE_DOFS) +=
             mollifier_grad(N_CORE_DOFS) * dist_sqr_grad;
@@ -378,12 +384,12 @@ auto SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::hessian(
 
     // grad of tangent/normal terms
     double orient = 0;
-    Vector<double, -1, ELEMENT_SIZE> gOrient;
+    VectorMax<double, ELEMENT_SIZE> gOrient;
     MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE> hOrient;
     {
-        Vector<double, -1, ELEMENT_SIZE>
-            gA = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs()),
-            gB = Vector<double, -1, ELEMENT_SIZE>::Zero(n_dofs());
+        VectorMax<double, ELEMENT_SIZE>
+            gA = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs()),
+            gB = VectorMax<double, ELEMENT_SIZE>::Zero(n_dofs());
         MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>
             hA = MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE>::Zero(
                 n_dofs(), n_dofs()),
@@ -467,9 +473,9 @@ template <typename PrimitiveA, typename PrimitiveB>
 double SmoothCollisionTemplate<PrimitiveA, PrimitiveB>::compute_distance(
     Eigen::ConstRef<Eigen::MatrixXd> vertices) const
 {
-    Vector<double, -1, ELEMENT_SIZE> positions = dof(vertices);
+    VectorMax<double, ELEMENT_SIZE> positions = dof(vertices);
 
-    Vector<double, N_CORE_POINTS * DIM> x;
+    Eigen::Vector<double, N_CORE_POINTS * DIM> x;
     x << positions.head(PrimitiveA::N_CORE_POINTS * DIM),
         positions.segment(
             primitive_a->n_dofs(), PrimitiveB::N_CORE_POINTS * DIM);
diff --git a/src/ipc/smooth_contact/collisions/smooth_collision.hpp b/src/ipc/smooth_contact/collisions/smooth_collision.hpp
index b9dd5e6ea..21f15fc14 100644
--- a/src/ipc/smooth_contact/collisions/smooth_collision.hpp
+++ b/src/ipc/smooth_contact/collisions/smooth_collision.hpp
@@ -81,17 +81,17 @@ class SmoothCollision {
 
     /// @brief Compute the value of the GCP potential
     virtual double operator()(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const = 0;
 
     /// @brief Compute the gradient of the GCP potential wrt. vertices involved
-    virtual Vector<double, -1, ELEMENT_SIZE> gradient(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+    virtual VectorMax<double, ELEMENT_SIZE> gradient(
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const = 0;
 
     /// @brief Compute the Hessian of the GCP potential wrt. vertices involved
     virtual MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE> hessian(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const = 0;
 
     bool operator==(const SmoothCollision& other) const
@@ -153,7 +153,7 @@ class SmoothCollisionTemplate : public SmoothCollision {
         const CollisionMesh& mesh,
         const SmoothContactParameters& params,
         const double dhat,
-        const Eigen::MatrixXd& V);
+        Eigen::ConstRef<Eigen::MatrixXd> V);
 
     virtual ~SmoothCollisionTemplate() = default;
 
@@ -165,7 +165,7 @@ class SmoothCollisionTemplate : public SmoothCollision {
     }
     CollisionType type() const override;
 
-    Vector<int, N_CORE_DOFS> get_core_indices() const;
+    Eigen::Vector<int, N_CORE_DOFS> get_core_indices() const;
     std::array<index_t, N_CORE_DOFS> core_vertex_ids() const;
 
     int num_vertices() const override
@@ -173,8 +173,8 @@ class SmoothCollisionTemplate : public SmoothCollision {
         return primitive_a->n_vertices() + primitive_b->n_vertices();
     }
 
-    template <typename T>
-    Vector<T, N_CORE_DOFS> core_dof(const Eigen::MatrixX<T>& X) const
+    Eigen::Vector<double, N_CORE_DOFS>
+    core_dof(Eigen::ConstRef<Eigen::MatrixXd> X) const
     {
         return this->dof(X)(get_core_indices());
     }
@@ -186,15 +186,15 @@ class SmoothCollisionTemplate : public SmoothCollision {
     /// @param params GCP parameters
     /// @return GCP potential value
     double operator()(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const override;
 
     /// @brief Compute the potential gradient wrt. positions
     /// @param positions Vertex positions
     /// @param params GCP parameters
     /// @return GCP potential gradient
-    Vector<double, -1, ELEMENT_SIZE> gradient(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+    VectorMax<double, ELEMENT_SIZE> gradient(
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const override;
 
     /// @brief Compute the potential Hessian wrt. positions
@@ -202,7 +202,7 @@ class SmoothCollisionTemplate : public SmoothCollision {
     /// @param params GCP parameters
     /// @return GCP potential Hessian
     MatrixMax<double, ELEMENT_SIZE, ELEMENT_SIZE> hessian(
-        Eigen::ConstRef<Vector<double, -1, ELEMENT_SIZE>> positions,
+        Eigen::ConstRef<VectorMax<double, ELEMENT_SIZE>> positions,
         const SmoothContactParameters& params) const override;
 
     // ---- distance ----
diff --git a/src/ipc/smooth_contact/distance/mollifier.cpp b/src/ipc/smooth_contact/distance/mollifier.cpp
index 9ad370c1d..d6c62d53b 100644
--- a/src/ipc/smooth_contact/distance/mollifier.cpp
+++ b/src/ipc/smooth_contact/distance/mollifier.cpp
@@ -16,15 +16,15 @@ namespace {
     //     return Math<double>::mollifier((a - c) / c / eps);
     // }
 
-    Vector<int, 9> get_indices(int i, int j, int k)
+    Eigen::Vector<int, 9> get_indices(int i, int j, int k)
     {
-        Vector<int, 9> out;
+        Eigen::Vector<int, 9> out;
         out << 3 * i + 0, 3 * i + 1, 3 * i + 2, 3 * j + 0, 3 * j + 1, 3 * j + 2,
             3 * k + 0, 3 * k + 1, 3 * k + 2;
         return out;
     }
 
-    GradType<3> func_aux_grad(
+    GradientType<3> func_aux_grad(
         const double a, const double b, const double c, const double eps)
     {
         const double val = (a - c) / c / eps;
@@ -54,7 +54,7 @@ namespace {
 } // namespace
 
 /// @brief Compute the gradient of the mollifier function wrt. 4 edge points and the distance squared
-GradType<13> edge_edge_mollifier_gradient(
+GradientType<13> edge_edge_mollifier_gradient(
     Eigen::ConstRef<Eigen::Vector3d> ea0,
     Eigen::ConstRef<Eigen::Vector3d> ea1,
     Eigen::ConstRef<Eigen::Vector3d> eb0,
@@ -62,7 +62,7 @@ GradType<13> edge_edge_mollifier_gradient(
     const std::array<HeavisideType, 4>& mtypes,
     const double dist_sqr)
 {
-    Vector<double, 13> input;
+    Eigen::Vector<double, 13> input;
     input << ea0, ea1, eb0, eb1, dist_sqr;
 
     Eigen::Matrix<int, 4, 3> vert_indices;
@@ -133,11 +133,11 @@ GradType<13> edge_edge_mollifier_gradient(
 
     // derivatives of mollifier products, input order : [dist_sqr, edge_lengths
     // (1 x 2), point_edge_dists (1 x 4)]
-    Vector<double, 7> product_grad =
+    Eigen::Vector<double, 7> product_grad =
         mollifier_grad.transpose() * partial_products_1;
 
     // derivatives wrt. input : [ea0, ea1, eb0, eb1, dist_sqr]
-    Vector<double, 13> grad;
+    Eigen::Vector<double, 13> grad;
     grad << edge_lengths_grad.transpose() * product_grad.segment<2>(1)
             + point_edge_dists_grad.transpose() * product_grad.segment<4>(3),
         product_grad(0);
@@ -154,7 +154,7 @@ HessianType<13> edge_edge_mollifier_hessian(
     const std::array<HeavisideType, 4>& mtypes,
     const double dist_sqr)
 {
-    Vector<double, 13> input;
+    Eigen::Vector<double, 13> input;
     input << ea0, ea1, eb0, eb1, dist_sqr;
 
     Eigen::Matrix<int, 4, 3> vert_indices;
@@ -248,7 +248,7 @@ HessianType<13> edge_edge_mollifier_hessian(
 
     // derivatives of mollifier products, input order : [dist_sqr, edge_lengths
     // (1 x 2), point_edge_dists (1 x 4)]
-    Vector<double, 7> product_grad =
+    Eigen::Vector<double, 7> product_grad =
         mollifier_grad.transpose() * partial_products_1;
     Eigen::Matrix<double, 7, 7> product_hess =
         mollifier_grad.transpose() * partial_products_2 * mollifier_grad;
@@ -257,7 +257,7 @@ HessianType<13> edge_edge_mollifier_hessian(
     }
 
     // derivatives wrt. input : [ea0, ea1, eb0, eb1, dist_sqr]
-    Vector<double, 13> grad = Vector<double, 13>::Zero();
+    Eigen::Vector<double, 13> grad = Eigen::Vector<double, 13>::Zero();
     Eigen::Matrix<double, 13, 13> hess = Eigen::Matrix<double, 13, 13>::Zero();
     {
         Eigen::Matrix<double, 6, 12> grads;
@@ -316,7 +316,7 @@ std::array<HeavisideType, 4> edge_edge_mollifier_type(
 }
 
 /// @brief Compute the gradient of the mollifier function wrt. 4 edge points and the distance squared
-GradType<13> point_face_mollifier_gradient(
+GradientType<13> point_face_mollifier_gradient(
     Eigen::ConstRef<Eigen::Vector3d> p,
     Eigen::ConstRef<Eigen::Vector3d> e0,
     Eigen::ConstRef<Eigen::Vector3d> e1,
@@ -333,7 +333,7 @@ GradType<13> point_face_mollifier_gradient(
         const int ej = ((i + 1) % 3) * 3 + 3;
         const int ek = ((i + 2) % 3) * 3 + 3;
 
-        Vector<int, 9> ind;
+        Eigen::Vector<int, 9> ind;
 
         Eigen::Matrix<double, 2, 13> dist_grad =
             Eigen::Matrix<double, 2, 13>::Zero();
@@ -358,7 +358,7 @@ GradType<13> point_face_mollifier_gradient(
         grads(i, 12) += tmp_grad(2);
     }
 
-    Vector<double, 13> grad = (vals(0) * vals(1)) * grads.row(2)
+    Eigen::Vector<double, 13> grad = (vals(0) * vals(1)) * grads.row(2)
         + (vals(0) * vals(2)) * grads.row(1)
         + (vals(1) * vals(2)) * grads.row(0);
 
@@ -384,7 +384,7 @@ HessianType<13> point_face_mollifier_hessian(
         const int ej = ((i + 1) % 3) * 3 + 3;
         const int ek = ((i + 2) % 3) * 3 + 3;
 
-        Vector<int, 9> ind;
+        Eigen::Vector<int, 9> ind;
 
         Eigen::Matrix<double, 2, 13> dist_grad =
             Eigen::Matrix<double, 2, 13>::Zero();
@@ -426,7 +426,7 @@ HessianType<13> point_face_mollifier_hessian(
         hesses[i].row(12) += tmp_hess.block<1, 2>(2, 0) * dist_grad;
     }
 
-    Vector<double, 13> grad = (vals(0) * vals(1)) * grads.row(2)
+    Eigen::Vector<double, 13> grad = (vals(0) * vals(1)) * grads.row(2)
         + (vals(0) * vals(2)) * grads.row(1)
         + (vals(1) * vals(2)) * grads.row(0);
     Eigen::Matrix<double, 13, 13> hess = (vals(0) * vals(1)) * hesses[2]
diff --git a/src/ipc/smooth_contact/distance/mollifier.hpp b/src/ipc/smooth_contact/distance/mollifier.hpp
index 662d1cfea..c3da5768d 100644
--- a/src/ipc/smooth_contact/distance/mollifier.hpp
+++ b/src/ipc/smooth_contact/distance/mollifier.hpp
@@ -6,9 +6,9 @@
 namespace ipc {
 template <typename scalar, int dim>
 scalar point_edge_mollifier(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1,
     const scalar& dist_sqr);
 
 std::array<HeavisideType, 4> edge_edge_mollifier_type(
@@ -28,7 +28,7 @@ scalar edge_edge_mollifier(
     const scalar& dist_sqr);
 
 /// @brief Compute the gradient of the mollifier function wrt. 4 edge points and the distance squared
-GradType<13> edge_edge_mollifier_gradient(
+GradientType<13> edge_edge_mollifier_gradient(
     Eigen::ConstRef<Eigen::Vector3d> ea0,
     Eigen::ConstRef<Eigen::Vector3d> ea1,
     Eigen::ConstRef<Eigen::Vector3d> eb0,
@@ -54,7 +54,7 @@ scalar point_face_mollifier(
     const scalar& dist_sqr);
 
 /// @brief Compute the gradient of the mollifier function wrt. 4 edge points and the distance squared
-GradType<13> point_face_mollifier_gradient(
+GradientType<13> point_face_mollifier_gradient(
     Eigen::ConstRef<Eigen::Vector3d> p,
     Eigen::ConstRef<Eigen::Vector3d> e0,
     Eigen::ConstRef<Eigen::Vector3d> e1,
diff --git a/src/ipc/smooth_contact/distance/mollifier.tpp b/src/ipc/smooth_contact/distance/mollifier.tpp
index 9718c8d72..701f1de12 100644
--- a/src/ipc/smooth_contact/distance/mollifier.tpp
+++ b/src/ipc/smooth_contact/distance/mollifier.tpp
@@ -5,9 +5,9 @@
 namespace ipc {
 template <typename scalar, int dim>
 scalar point_edge_mollifier(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1,
     const scalar& dist_sqr)
 {
     const scalar denominator = dist_sqr * MOLLIFIER_THRESHOLD_EPS;
diff --git a/src/ipc/smooth_contact/distance/point_edge.cpp b/src/ipc/smooth_contact/distance/point_edge.cpp
index 303c98d53..beeddc148 100644
--- a/src/ipc/smooth_contact/distance/point_edge.cpp
+++ b/src/ipc/smooth_contact/distance/point_edge.cpp
@@ -6,17 +6,17 @@
 namespace ipc {
 template <typename scalar, int dim>
 scalar PointEdgeDistance<scalar, dim>::point_point_sqr_distance(
-    Eigen::ConstRef<Vector<scalar, dim>> a,
-    Eigen::ConstRef<Vector<scalar, dim>> b)
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> a,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> b)
 {
     return (a - b).squaredNorm();
 }
 
 template <typename scalar, int dim>
 scalar PointEdgeDistance<scalar, dim>::point_line_sqr_distance(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1)
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1)
 {
     if constexpr (dim == 2) {
         return Math<scalar>::sqr(Math<scalar>::cross2(e0 - p, e1 - p))
@@ -28,9 +28,9 @@ scalar PointEdgeDistance<scalar, dim>::point_line_sqr_distance(
 
 template <typename scalar, int dim>
 scalar PointEdgeDistance<scalar, dim>::point_edge_sqr_distance(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1,
     const PointEdgeDistanceType dtype)
 {
     switch (dtype) {
@@ -42,31 +42,31 @@ scalar PointEdgeDistance<scalar, dim>::point_edge_sqr_distance(
         return point_point_sqr_distance(p, e1);
     case PointEdgeDistanceType::AUTO:
     default:
-        const Vector<scalar, dim> t = e1 - e0;
-        const Vector<scalar, dim> pos = p - e0;
+        const Eigen::Vector<scalar, dim> t = e1 - e0;
+        const Eigen::Vector<scalar, dim> pos = p - e0;
         const scalar s = pos.dot(t) / t.squaredNorm();
         return (pos - Math<scalar>::l_ns(s) * t).squaredNorm();
     }
 }
 
 template <typename scalar, int dim>
-Vector<scalar, dim>
+Eigen::Vector<scalar, dim>
 PointEdgeDistance<scalar, dim>::point_line_closest_point_direction(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1)
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1)
 {
-    const Vector<scalar, dim> d = p - e0;
-    const Vector<scalar, dim> t = e1 - e0;
+    const Eigen::Vector<scalar, dim> d = p - e0;
+    const Eigen::Vector<scalar, dim> t = e1 - e0;
     return d - (d.dot(t) / t.squaredNorm()) * t;
 }
 
 template <typename scalar, int dim>
-Vector<scalar, dim>
+Eigen::Vector<scalar, dim>
 PointEdgeDistance<scalar, dim>::point_edge_closest_point_direction(
-    Eigen::ConstRef<Vector<scalar, dim>> p,
-    Eigen::ConstRef<Vector<scalar, dim>> e0,
-    Eigen::ConstRef<Vector<scalar, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<scalar, dim>> e1,
     const PointEdgeDistanceType dtype)
 {
     switch (dtype) {
@@ -78,30 +78,30 @@ PointEdgeDistance<scalar, dim>::point_edge_closest_point_direction(
         return p - e1;
     case PointEdgeDistanceType::AUTO:
     default:
-        Vector<scalar, dim> t = e1 - e0;
-        const Vector<scalar, dim> pos = p - e0;
+        Eigen::Vector<scalar, dim> t = e1 - e0;
+        const Eigen::Vector<scalar, dim> pos = p - e0;
         const scalar s = pos.dot(t) / t.squaredNorm();
         return pos - Math<scalar>::l_ns(s) * t;
     }
 }
 
 template <int dim>
-std::tuple<Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
+std::tuple<Eigen::Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
 PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_grad(
-    Eigen::ConstRef<Vector<double, dim>> p,
-    Eigen::ConstRef<Vector<double, dim>> e0,
-    Eigen::ConstRef<Vector<double, dim>> e1)
+    Eigen::ConstRef<Eigen::Vector<double, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e1)
 {
-    Vector<double, dim> val;
+    Eigen::Vector<double, dim> val;
     Eigen::Matrix<double, dim, 3 * dim> grad;
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADGrad<3 * dim>;
     ScalarBase::setVariableCount(3 * dim);
-    const Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
-    const Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
-    const Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
+    const Eigen::Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
+    const Eigen::Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
+    const Eigen::Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
 
-    const Vector<T, dim> out =
+    const Eigen::Vector<T, dim> out =
         PointEdgeDistance<T, dim>::point_line_closest_point_direction(
             pT, e0T, e1T);
     for (int i = 0; i < dim; i++) {
@@ -110,7 +110,7 @@ PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_grad(
     }
 #else
     const double uv = point_edge_closest_point(p, e0, e1);
-    const Vector<double, 3 * dim> g =
+    const Eigen::Vector<double, 3 * dim> g =
         point_edge_closest_point_jacobian(p, e0, e1);
 
     val = (p - e0) - uv * (e1 - e0);
@@ -127,26 +127,26 @@ PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_grad(
 
 template <int dim>
 std::tuple<
-    Vector<double, dim>,
+    Eigen::Vector<double, dim>,
     Eigen::Matrix<double, dim, 3 * dim>,
     std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim>>
 PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_hessian(
-    Eigen::ConstRef<Vector<double, dim>> p,
-    Eigen::ConstRef<Vector<double, dim>> e0,
-    Eigen::ConstRef<Vector<double, dim>> e1)
+    Eigen::ConstRef<Eigen::Vector<double, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e1)
 {
-    Vector<double, dim> val;
+    Eigen::Vector<double, dim> val;
     Eigen::Matrix<double, dim, 3 * dim> grad;
     std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim> hess;
 
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADHessian<3 * dim>;
     ScalarBase::setVariableCount(3 * dim);
-    Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
-    Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
-    Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
+    Eigen::Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
+    Eigen::Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
+    Eigen::Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
 
-    Vector<T, dim> out =
+    Eigen::Vector<T, dim> out =
         PointEdgeDistance<T, dim>::point_line_closest_point_direction(
             pT, e0T, e1T);
     for (int i = 0; i < dim; i++) {
@@ -156,7 +156,7 @@ PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_hessian(
     }
 #else
     const double uv = point_edge_closest_point(p, e0, e1);
-    const Vector<double, 3 * dim> g =
+    const Eigen::Vector<double, 3 * dim> g =
         point_edge_closest_point_jacobian(p, e0, e1);
     const Eigen::Matrix<double, 3 * dim, 3 * dim> h =
         point_edge_closest_point_hessian(p, e0, e1);
@@ -190,24 +190,24 @@ PointEdgeDistanceDerivatives<dim>::point_line_closest_point_direction_hessian(
 }
 
 template <int dim>
-std::tuple<Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
+std::tuple<Eigen::Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
 PointEdgeDistanceDerivatives<dim>::point_edge_closest_point_direction_grad(
-    Eigen::ConstRef<Vector<double, dim>> p,
-    Eigen::ConstRef<Vector<double, dim>> e0,
-    Eigen::ConstRef<Vector<double, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e1,
     const PointEdgeDistanceType dtype)
 {
-    Vector<double, dim> vec;
+    Eigen::Vector<double, dim> vec;
     Eigen::Matrix<double, dim, 3 * dim> grad =
         Eigen::Matrix<double, dim, 3 * dim>::Zero();
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADGrad<3 * dim>;
     ScalarBase::setVariableCount(3 * dim);
-    Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
-    Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
-    Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
+    Eigen::Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
+    Eigen::Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
+    Eigen::Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
 
-    Vector<T, dim> out =
+    Eigen::Vector<T, dim> out =
         PointEdgeDistance<T, dim>::point_edge_closest_point_direction(
             pT, e0T, e1T, dtype);
     for (int i = 0; i < dim; i++) {
@@ -240,16 +240,16 @@ PointEdgeDistanceDerivatives<dim>::point_edge_closest_point_direction_grad(
 
 template <int dim>
 std::tuple<
-    Vector<double, dim>,
+    Eigen::Vector<double, dim>,
     Eigen::Matrix<double, dim, 3 * dim>,
     std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim>>
 PointEdgeDistanceDerivatives<dim>::point_edge_closest_point_direction_hessian(
-    Eigen::ConstRef<Vector<double, dim>> p,
-    Eigen::ConstRef<Vector<double, dim>> e0,
-    Eigen::ConstRef<Vector<double, dim>> e1,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> p,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e0,
+    Eigen::ConstRef<Eigen::Vector<double, dim>> e1,
     const PointEdgeDistanceType dtype)
 {
-    Vector<double, dim> vec;
+    Eigen::Vector<double, dim> vec;
     Eigen::Matrix<double, dim, 3 * dim> grad =
         Eigen::Matrix<double, dim, 3 * dim>::Zero();
     std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim> hess;
@@ -259,11 +259,11 @@ PointEdgeDistanceDerivatives<dim>::point_edge_closest_point_direction_hessian(
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADHessian<3 * dim>;
     ScalarBase::setVariableCount(3 * dim);
-    Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
-    Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
-    Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
+    Eigen::Vector<T, dim> pT = slice_positions<T, 1, dim>(p);
+    Eigen::Vector<T, dim> e0T = slice_positions<T, 1, dim>(e0, dim);
+    Eigen::Vector<T, dim> e1T = slice_positions<T, 1, dim>(e1, 2 * dim);
 
-    Vector<T, dim> out =
+    Eigen::Vector<T, dim> out =
         PointEdgeDistance<T, dim>::point_edge_closest_point_direction(
             pT, e0T, e1T);
 
diff --git a/src/ipc/smooth_contact/distance/point_edge.hpp b/src/ipc/smooth_contact/distance/point_edge.hpp
index 9749234a1..11bade859 100644
--- a/src/ipc/smooth_contact/distance/point_edge.hpp
+++ b/src/ipc/smooth_contact/distance/point_edge.hpp
@@ -9,76 +9,75 @@
 #include <iostream>
 
 namespace ipc {
-template <typename scalar, int dim> class PointEdgeDistance {
+template <typename T, int dim> class PointEdgeDistance {
 public:
+    using VectorNT = Eigen::Vector<T, dim>;
+
     PointEdgeDistance() = delete;
     PointEdgeDistance(const PointEdgeDistance&) = delete;
     PointEdgeDistance& operator=(const PointEdgeDistance&) = delete;
 
-    static scalar point_point_sqr_distance(
-        Eigen::ConstRef<Vector<scalar, dim>> a,
-        Eigen::ConstRef<Vector<scalar, dim>> b);
+    static T point_point_sqr_distance(
+        Eigen::ConstRef<VectorNT> a, Eigen::ConstRef<VectorNT> b);
 
-    static scalar point_line_sqr_distance(
-        Eigen::ConstRef<Vector<scalar, dim>> p,
-        Eigen::ConstRef<Vector<scalar, dim>> e0,
-        Eigen::ConstRef<Vector<scalar, dim>> e1);
+    static T point_line_sqr_distance(
+        Eigen::ConstRef<VectorNT> p,
+        Eigen::ConstRef<VectorNT> e0,
+        Eigen::ConstRef<VectorNT> e1);
 
-    static scalar point_edge_sqr_distance(
-        Eigen::ConstRef<Vector<scalar, dim>> p,
-        Eigen::ConstRef<Vector<scalar, dim>> e0,
-        Eigen::ConstRef<Vector<scalar, dim>> e1,
+    static T point_edge_sqr_distance(
+        Eigen::ConstRef<VectorNT> p,
+        Eigen::ConstRef<VectorNT> e0,
+        Eigen::ConstRef<VectorNT> e1,
         const PointEdgeDistanceType dtype = PointEdgeDistanceType::AUTO);
 
-    static Vector<scalar, dim> point_line_closest_point_direction(
-        Eigen::ConstRef<Vector<scalar, dim>> p,
-        Eigen::ConstRef<Vector<scalar, dim>> e0,
-        Eigen::ConstRef<Vector<scalar, dim>> e1);
+    static VectorNT point_line_closest_point_direction(
+        Eigen::ConstRef<VectorNT> p,
+        Eigen::ConstRef<VectorNT> e0,
+        Eigen::ConstRef<VectorNT> e1);
 
-    static Vector<scalar, dim> point_edge_closest_point_direction(
-        Eigen::ConstRef<Vector<scalar, dim>> p,
-        Eigen::ConstRef<Vector<scalar, dim>> e0,
-        Eigen::ConstRef<Vector<scalar, dim>> e1,
+    static VectorNT point_edge_closest_point_direction(
+        Eigen::ConstRef<VectorNT> p,
+        Eigen::ConstRef<VectorNT> e0,
+        Eigen::ConstRef<VectorNT> e1,
         const PointEdgeDistanceType dtype = PointEdgeDistanceType::AUTO);
 };
 
 template <int dim> class PointEdgeDistanceDerivatives {
 public:
+    using VectorNd = Eigen::Vector<double, dim>;
+    using JacobianType =
+        std::tuple<VectorNd, Eigen::Matrix<double, dim, 3 * dim>>;
+    using HessianType = std::tuple<
+        VectorNd,
+        Eigen::Matrix<double, dim, 3 * dim>,
+        std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim>>;
+
     PointEdgeDistanceDerivatives() = delete;
     PointEdgeDistanceDerivatives(const PointEdgeDistanceDerivatives&) = delete;
     PointEdgeDistanceDerivatives&
     operator=(const PointEdgeDistanceDerivatives&) = delete;
 
-    static std::tuple<Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
-    point_line_closest_point_direction_grad(
-        Eigen::ConstRef<Vector<double, dim>> p,
-        Eigen::ConstRef<Vector<double, dim>> e0,
-        Eigen::ConstRef<Vector<double, dim>> e1);
+    static JacobianType point_line_closest_point_direction_grad(
+        Eigen::ConstRef<VectorNd> p,
+        Eigen::ConstRef<VectorNd> e0,
+        Eigen::ConstRef<VectorNd> e1);
 
-    static std::tuple<
-        Vector<double, dim>,
-        Eigen::Matrix<double, dim, 3 * dim>,
-        std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim>>
-    point_line_closest_point_direction_hessian(
-        Eigen::ConstRef<Vector<double, dim>> p,
-        Eigen::ConstRef<Vector<double, dim>> e0,
-        Eigen::ConstRef<Vector<double, dim>> e1);
+    static HessianType point_line_closest_point_direction_hessian(
+        Eigen::ConstRef<VectorNd> p,
+        Eigen::ConstRef<VectorNd> e0,
+        Eigen::ConstRef<VectorNd> e1);
 
-    static std::tuple<Vector<double, dim>, Eigen::Matrix<double, dim, 3 * dim>>
-    point_edge_closest_point_direction_grad(
-        Eigen::ConstRef<Vector<double, dim>> p,
-        Eigen::ConstRef<Vector<double, dim>> e0,
-        Eigen::ConstRef<Vector<double, dim>> e1,
+    static JacobianType point_edge_closest_point_direction_grad(
+        Eigen::ConstRef<VectorNd> p,
+        Eigen::ConstRef<VectorNd> e0,
+        Eigen::ConstRef<VectorNd> e1,
         const PointEdgeDistanceType dtype = PointEdgeDistanceType::AUTO);
 
-    static std::tuple<
-        Vector<double, dim>,
-        Eigen::Matrix<double, dim, 3 * dim>,
-        std::array<Eigen::Matrix<double, 3 * dim, 3 * dim>, dim>>
-    point_edge_closest_point_direction_hessian(
-        Eigen::ConstRef<Vector<double, dim>> p,
-        Eigen::ConstRef<Vector<double, dim>> e0,
-        Eigen::ConstRef<Vector<double, dim>> e1,
+    static HessianType point_edge_closest_point_direction_hessian(
+        Eigen::ConstRef<VectorNd> p,
+        Eigen::ConstRef<VectorNd> e0,
+        Eigen::ConstRef<VectorNd> e1,
         const PointEdgeDistanceType dtype = PointEdgeDistanceType::AUTO);
 };
 } // namespace ipc
diff --git a/src/ipc/smooth_contact/distance/primitive_distance.cpp b/src/ipc/smooth_contact/distance/primitive_distance.cpp
index 0d5454003..b27ff7919 100644
--- a/src/ipc/smooth_contact/distance/primitive_distance.cpp
+++ b/src/ipc/smooth_contact/distance/primitive_distance.cpp
@@ -14,7 +14,7 @@ namespace ipc {
 template <>
 typename PrimitiveDistType<Face, Point3>::type
 PrimitiveDistance<Face, Point3>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return point_triangle_distance_type(
         x.tail(3) /* point */, x.head(3), x.segment(3, 3),
@@ -24,7 +24,7 @@ PrimitiveDistance<Face, Point3>::compute_distance_type(
 template <>
 typename PrimitiveDistType<Edge3, Edge3>::type
 PrimitiveDistance<Edge3, Edge3>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return edge_edge_distance_type(
         x.head(3) /* edge 0 */, x.segment(3, 3) /* edge 0 */,
@@ -34,7 +34,7 @@ PrimitiveDistance<Edge3, Edge3>::compute_distance_type(
 template <>
 typename PrimitiveDistType<Edge3, Point3>::type
 PrimitiveDistance<Edge3, Point3>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return point_edge_distance_type(
         x.tail(3) /* point */, x.head(3) /* edge */,
@@ -44,7 +44,7 @@ PrimitiveDistance<Edge3, Point3>::compute_distance_type(
 template <>
 typename PrimitiveDistType<Edge2, Point2>::type
 PrimitiveDistance<Edge2, Point2>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return PointEdgeDistanceType::AUTO;
 }
@@ -52,7 +52,7 @@ PrimitiveDistance<Edge2, Point2>::compute_distance_type(
 template <>
 typename PrimitiveDistType<Point2, Point2>::type
 PrimitiveDistance<Point2, Point2>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return PointPointDistanceType::P_P;
 }
@@ -60,16 +60,16 @@ PrimitiveDistance<Point2, Point2>::compute_distance_type(
 template <>
 typename PrimitiveDistType<Point3, Point3>::type
 PrimitiveDistance<Point3, Point3>::compute_distance_type(
-    const Vector<double, N_CORE_DOFS>& x)
+    const Eigen::Vector<double, N_CORE_DOFS>& x)
 {
     return PointPointDistanceType::P_P;
 }
 
 template <>
-Vector<double, PrimitiveDistance<Face, Point3>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Face, Point3>::DIM>
 PrimitiveDistance<Face, Point3>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Face, Point3>::type dtype)
@@ -80,10 +80,10 @@ PrimitiveDistance<Face, Point3>::compute_closest_direction(
 }
 
 template <>
-Vector<double, PrimitiveDistance<Edge3, Edge3>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Edge3, Edge3>::DIM>
 PrimitiveDistance<Edge3, Edge3>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Edge3, Edge3>::type dtype)
@@ -94,10 +94,10 @@ PrimitiveDistance<Edge3, Edge3>::compute_closest_direction(
 }
 
 template <>
-Vector<double, PrimitiveDistance<Edge3, Point3>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Edge3, Point3>::DIM>
 PrimitiveDistance<Edge3, Point3>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Edge3, Point3>::type dtype)
@@ -107,10 +107,10 @@ PrimitiveDistance<Edge3, Point3>::compute_closest_direction(
 }
 
 template <>
-Vector<double, PrimitiveDistance<Edge2, Point2>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Edge2, Point2>::DIM>
 PrimitiveDistance<Edge2, Point2>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Edge2, Point2>::type dtype)
@@ -120,10 +120,10 @@ PrimitiveDistance<Edge2, Point2>::compute_closest_direction(
 }
 
 template <>
-Vector<double, PrimitiveDistance<Point2, Point2>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Point2, Point2>::DIM>
 PrimitiveDistance<Point2, Point2>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Point2, Point2>::type dtype)
@@ -132,10 +132,10 @@ PrimitiveDistance<Point2, Point2>::compute_closest_direction(
 }
 
 template <>
-Vector<double, PrimitiveDistance<Point3, Point3>::DIM>
+Eigen::Vector<double, PrimitiveDistance<Point3, Point3>::DIM>
 PrimitiveDistance<Point3, Point3>::compute_closest_direction(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& V,
+    Eigen::ConstRef<Eigen::MatrixXd> V,
     const index_t a,
     const index_t b,
     typename PrimitiveDistType<Point3, Point3>::type dtype)
@@ -147,7 +147,7 @@ PrimitiveDistance<Point3, Point3>::compute_closest_direction(
 
 template <>
 std::tuple<
-    Vector<double, PrimitiveDistance<Edge3, Edge3>::DIM>,
+    Eigen::Vector<double, PrimitiveDistance<Edge3, Edge3>::DIM>,
     Eigen::Matrix<
         double,
         PrimitiveDistance<Edge3, Edge3>::DIM,
@@ -159,7 +159,7 @@ std::tuple<
             PrimitiveDistance<Edge3, Edge3>::N_CORE_DOFS>,
         PrimitiveDistance<Edge3, Edge3>::DIM>>
 PrimitiveDistance<Edge3, Edge3>::compute_closest_direction_hessian(
-    const Vector<double, N_CORE_DOFS>& x,
+    const Eigen::Vector<double, N_CORE_DOFS>& x,
     typename PrimitiveDistType<Edge3, Edge3>::type dtype)
 {
     assert(dtype == EdgeEdgeDistanceType::EA_EB);
@@ -170,7 +170,7 @@ PrimitiveDistance<Edge3, Edge3>::compute_closest_direction_hessian(
 
 template <>
 std::tuple<
-    Vector<double, PrimitiveDistance<Point2, Point2>::DIM>,
+    Eigen::Vector<double, PrimitiveDistance<Point2, Point2>::DIM>,
     Eigen::Matrix<
         double,
         PrimitiveDistance<Point2, Point2>::DIM,
@@ -182,7 +182,7 @@ std::tuple<
             PrimitiveDistance<Point2, Point2>::N_CORE_DOFS>,
         PrimitiveDistance<Point2, Point2>::DIM>>
 PrimitiveDistance<Point2, Point2>::compute_closest_direction_hessian(
-    const Vector<double, N_CORE_DOFS>& x,
+    const Eigen::Vector<double, N_CORE_DOFS>& x,
     typename PrimitiveDistType<Point2, Point2>::type dtype)
 {
     Eigen::Vector2d out = x.tail(2) - x.head(2);
@@ -198,7 +198,7 @@ PrimitiveDistance<Point2, Point2>::compute_closest_direction_hessian(
 
 template <>
 std::tuple<
-    Vector<double, PrimitiveDistance<Point3, Point3>::DIM>,
+    Eigen::Vector<double, PrimitiveDistance<Point3, Point3>::DIM>,
     Eigen::Matrix<
         double,
         PrimitiveDistance<Point3, Point3>::DIM,
@@ -210,7 +210,7 @@ std::tuple<
             PrimitiveDistance<Point3, Point3>::N_CORE_DOFS>,
         PrimitiveDistance<Point3, Point3>::DIM>>
 PrimitiveDistance<Point3, Point3>::compute_closest_direction_hessian(
-    const Vector<double, N_CORE_DOFS>& x,
+    const Eigen::Vector<double, N_CORE_DOFS>& x,
     typename PrimitiveDistType<Point3, Point3>::type dtype)
 {
     Eigen::Vector3d out = x.tail(3) - x.head(3);
@@ -225,9 +225,9 @@ PrimitiveDistance<Point3, Point3>::compute_closest_direction_hessian(
 }
 
 template <>
-GradType<PrimitiveDistance<Edge3, Edge3>::N_CORE_DOFS + 1>
+GradientType<PrimitiveDistance<Edge3, Edge3>::N_CORE_DOFS + 1>
 PrimitiveDistance<Edge3, Edge3>::compute_mollifier_gradient(
-    const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+    const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
 {
     const auto otypes = edge_edge_mollifier_type(
         x.head(3), x.segment(3, 3), x.segment(6, 3), x.tail(3), dist_sqr);
@@ -240,7 +240,7 @@ PrimitiveDistance<Edge3, Edge3>::compute_mollifier_gradient(
 template <>
 HessianType<PrimitiveDistance<Edge3, Edge3>::N_CORE_DOFS + 1>
 PrimitiveDistance<Edge3, Edge3>::compute_mollifier_hessian(
-    const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+    const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
 {
     const auto otypes = edge_edge_mollifier_type(
         x.head<3>(), x.segment<3>(3), x.segment<3>(6), x.tail<3>(), dist_sqr);
@@ -251,13 +251,13 @@ PrimitiveDistance<Edge3, Edge3>::compute_mollifier_hessian(
 }
 
 template <>
-GradType<PrimitiveDistance<Face, Point3>::N_CORE_DOFS + 1>
+GradientType<PrimitiveDistance<Face, Point3>::N_CORE_DOFS + 1>
 PrimitiveDistance<Face, Point3>::compute_mollifier_gradient(
-    const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+    const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
 {
     const auto [val, grad] = point_face_mollifier_gradient(
         x.tail<3>(), x.head<3>(), x.segment<3>(3), x.segment<3>(6), dist_sqr);
-    Vector<int, 13> indices;
+    Eigen::Vector<int, 13> indices;
     indices << 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 12;
     return std::make_tuple(val, grad(indices));
 }
@@ -265,11 +265,11 @@ PrimitiveDistance<Face, Point3>::compute_mollifier_gradient(
 template <>
 HessianType<PrimitiveDistance<Face, Point3>::N_CORE_DOFS + 1>
 PrimitiveDistance<Face, Point3>::compute_mollifier_hessian(
-    const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+    const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
 {
     const auto [val, grad, hess] = point_face_mollifier_hessian(
         x.tail<3>(), x.head<3>(), x.segment<3>(3), x.segment<3>(6), dist_sqr);
-    Vector<int, 13> indices;
+    Eigen::Vector<int, 13> indices;
     indices << 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 12;
     return std::make_tuple(val, grad(indices), hess(indices, indices));
 }
diff --git a/src/ipc/smooth_contact/distance/primitive_distance.hpp b/src/ipc/smooth_contact/distance/primitive_distance.hpp
index e77a79f72..9debba70d 100644
--- a/src/ipc/smooth_contact/distance/primitive_distance.hpp
+++ b/src/ipc/smooth_contact/distance/primitive_distance.hpp
@@ -53,17 +53,22 @@ class PrimitiveDistanceTemplate {
         PrimitiveA::N_CORE_POINTS * PrimitiveA::DIM
         + PrimitiveB::N_CORE_POINTS * PrimitiveB::DIM;
 
+    using VectorNT = Eigen::Vector<T, N_CORE_DOFS>;
+
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        Eigen::ConstRef<VectorNT> x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype);
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        Eigen::ConstRef<VectorNT> x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype);
+
     static Eigen::Matrix<T, DIM, 2> compute_closest_point_pairs(
-        const Vector<T, N_CORE_DOFS>& x,
+        Eigen::ConstRef<VectorNT> x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype);
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr);
+
+    static T mollifier(Eigen::ConstRef<VectorNT> x, const T& dist_sqr);
 };
 
 template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
@@ -78,22 +83,23 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
         + PrimitiveB::N_CORE_POINTS * PrimitiveB::DIM;
 
     static typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type
-    compute_distance_type(const Vector<double, N_CORE_DOFS>& x);
+    compute_distance_type(const Eigen::Vector<double, N_CORE_DOFS>& x);
 
     static double compute_distance(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& V,
+        Eigen::ConstRef<Eigen::MatrixXd> V,
         const index_t a,
         const index_t b,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype);
 
-    static GradType<N_CORE_DOFS> compute_distance_gradient(
-        const Vector<double, N_CORE_DOFS>& x,
+    static GradientType<N_CORE_DOFS> compute_distance_gradient(
+        const Eigen::Vector<double, N_CORE_DOFS>& x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADGrad<N_CORE_DOFS>;
-        const Vector<T, N_CORE_DOFS> X = slice_positions<T, N_CORE_DOFS, 1>(x);
+        const Eigen::Vector<T, N_CORE_DOFS> X =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
         const T d = PrimitiveDistanceTemplate<
             PrimitiveA, PrimitiveB, T>::compute_distance(X, dtype);
 
@@ -101,12 +107,13 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
     }
 
     static HessianType<N_CORE_DOFS> compute_distance_hessian(
-        const Vector<double, N_CORE_DOFS>& x,
+        const Eigen::Vector<double, N_CORE_DOFS>& x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADHessian<N_CORE_DOFS>;
-        const Vector<T, N_CORE_DOFS> X = slice_positions<T, N_CORE_DOFS, 1>(x);
+        const Eigen::Vector<T, N_CORE_DOFS> X =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
         const T d = PrimitiveDistanceTemplate<
             PrimitiveA, PrimitiveB, T>::compute_distance(X, dtype);
 
@@ -114,26 +121,28 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
     }
 
     // points from primitiveA to primitiveB
-    static Vector<double, DIM> compute_closest_direction(
+    static Eigen::Vector<double, DIM> compute_closest_direction(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& V,
+        Eigen::ConstRef<Eigen::MatrixXd> V,
         const index_t a,
         const index_t b,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype);
 
-    static std::
-        tuple<Vector<double, DIM>, Eigen::Matrix<double, DIM, N_CORE_DOFS>>
-        compute_closest_direction_gradient(
-            const Vector<double, N_CORE_DOFS>& x,
-            typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype)
+    static std::tuple<
+        Eigen::Vector<double, DIM>,
+        Eigen::Matrix<double, DIM, N_CORE_DOFS>>
+    compute_closest_direction_gradient(
+        const Eigen::Vector<double, N_CORE_DOFS>& x,
+        typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADGrad<N_CORE_DOFS>;
-        const Vector<T, N_CORE_DOFS> X = slice_positions<T, N_CORE_DOFS, 1>(x);
-        const Vector<T, DIM> d = PrimitiveDistanceTemplate<
+        const Eigen::Vector<T, N_CORE_DOFS> X =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
+        const Eigen::Vector<T, DIM> d = PrimitiveDistanceTemplate<
             PrimitiveA, PrimitiveB, T>::compute_closest_direction(X, dtype);
 
-        Vector<double, DIM> out;
+        Eigen::Vector<double, DIM> out;
         Eigen::Matrix<double, DIM, N_CORE_DOFS> J =
             Eigen::Matrix<double, DIM, N_CORE_DOFS>::Zero();
         for (int i = 0; i < DIM; i++) {
@@ -144,20 +153,21 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
     }
 
     static std::tuple<
-        Vector<double, DIM>,
+        Eigen::Vector<double, DIM>,
         Eigen::Matrix<double, DIM, N_CORE_DOFS>,
         std::array<Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>, DIM>>
     compute_closest_direction_hessian(
-        const Vector<double, N_CORE_DOFS>& x,
+        const Eigen::Vector<double, N_CORE_DOFS>& x,
         typename PrimitiveDistType<PrimitiveA, PrimitiveB>::type dtype)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS);
         using T = ADHessian<N_CORE_DOFS>;
-        const Vector<T, N_CORE_DOFS> X = slice_positions<T, N_CORE_DOFS, 1>(x);
-        const Vector<T, DIM> d = PrimitiveDistanceTemplate<
+        const Eigen::Vector<T, N_CORE_DOFS> X =
+            slice_positions<T, N_CORE_DOFS, 1>(x);
+        const Eigen::Vector<T, DIM> d = PrimitiveDistanceTemplate<
             PrimitiveA, PrimitiveB, T>::compute_closest_direction(X, dtype);
 
-        Vector<double, DIM> out;
+        Eigen::Vector<double, DIM> out;
         Eigen::Matrix<double, DIM, N_CORE_DOFS> J =
             Eigen::Matrix<double, DIM, N_CORE_DOFS>::Zero();
         std::array<Eigen::Matrix<double, N_CORE_DOFS, N_CORE_DOFS>, DIM> H;
@@ -169,14 +179,15 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
         return std::make_tuple(out, J, H);
     }
 
-    static GradType<N_CORE_DOFS + 1> compute_mollifier_gradient(
-        const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+    static GradientType<N_CORE_DOFS + 1> compute_mollifier_gradient(
+        const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS + 1);
         using T = ADGrad<N_CORE_DOFS + 1>;
-        const Vector<T, N_CORE_DOFS + 1> X =
+        const Eigen::Vector<T, N_CORE_DOFS + 1> X =
             slice_positions<T, N_CORE_DOFS + 1, 1>(
-                (Vector<double, N_CORE_DOFS + 1>() << x, dist_sqr).finished());
+                (Eigen::Vector<double, N_CORE_DOFS + 1>() << x, dist_sqr)
+                    .finished());
         const T out =
             PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
                 X.head(N_CORE_DOFS), X(N_CORE_DOFS));
@@ -185,13 +196,14 @@ template <typename PrimitiveA, typename PrimitiveB> class PrimitiveDistance {
     }
 
     static HessianType<N_CORE_DOFS + 1> compute_mollifier_hessian(
-        const Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
+        const Eigen::Vector<double, N_CORE_DOFS>& x, const double dist_sqr)
     {
         ScalarBase::setVariableCount(N_CORE_DOFS + 1);
         using T = ADHessian<N_CORE_DOFS + 1>;
-        const Vector<T, N_CORE_DOFS + 1> X =
+        const Eigen::Vector<T, N_CORE_DOFS + 1> X =
             slice_positions<T, N_CORE_DOFS + 1, 1>(
-                (Vector<double, N_CORE_DOFS + 1>() << x, dist_sqr).finished());
+                (Eigen::Vector<double, N_CORE_DOFS + 1>() << x, dist_sqr)
+                    .finished());
         const T out =
             PrimitiveDistanceTemplate<PrimitiveA, PrimitiveB, T>::mollifier(
                 X.head(N_CORE_DOFS), X(N_CORE_DOFS));
diff --git a/src/ipc/smooth_contact/distance/primitive_distance.tpp b/src/ipc/smooth_contact/distance/primitive_distance.tpp
index 98ec97427..0fe0f1604 100644
--- a/src/ipc/smooth_contact/distance/primitive_distance.tpp
+++ b/src/ipc/smooth_contact/distance/primitive_distance.tpp
@@ -13,7 +13,7 @@ template <typename T> class PrimitiveDistanceTemplate<Face, Point3, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Face, Point3>::type dtype)
     {
         return point_triangle_sqr_distance<T>(
@@ -22,8 +22,8 @@ public:
             dtype);
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Face, Point3>::type dtype)
     {
         return point_triangle_closest_point_direction<T>(
@@ -32,7 +32,8 @@ public:
             dtype);
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         return point_face_mollifier<T>(
             x.template tail<3>() /* point */, x.template head<3>(),
@@ -48,7 +49,7 @@ template <typename T> class PrimitiveDistanceTemplate<Edge3, Edge3, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge3, Edge3>::type dtype)
     {
         return edge_edge_sqr_distance<T>(
@@ -58,8 +59,8 @@ public:
             x.template tail<3>() /* edge 1 */, dtype);
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge3, Edge3>::type dtype)
     {
         return edge_edge_closest_point_direction<T>(
@@ -69,7 +70,8 @@ public:
             x.template tail<3>() /* edge 1 */, dtype);
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         std::array<HeavisideType, 4> types {};
         types.fill(HeavisideType::VARIANT);
@@ -90,7 +92,7 @@ template <typename T> class PrimitiveDistanceTemplate<Edge2, Point2, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge2, Point2>::type dtype)
     {
         return PointEdgeDistance<T, DIM>::point_edge_sqr_distance(
@@ -98,8 +100,8 @@ public:
             x.template segment<2>(2) /* edge */, dtype);
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge2, Point2>::type dtype)
     {
         return PointEdgeDistance<T, DIM>::point_edge_closest_point_direction(
@@ -107,7 +109,8 @@ public:
             x.template segment<2>(2) /* edge */, dtype);
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         return point_edge_mollifier<T, 2>(
             x.template tail<2>() /* point */,
@@ -125,7 +128,7 @@ template <typename T> class PrimitiveDistanceTemplate<Edge3, Point3, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge3, Point3>::type dtype)
     {
         return PointEdgeDistance<T, DIM>::point_edge_sqr_distance(
@@ -133,8 +136,8 @@ public:
             x.template segment<3>(3) /* edge */, dtype);
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Edge3, Point3>::type dtype)
     {
         return PointEdgeDistance<T, DIM>::point_edge_closest_point_direction(
@@ -142,7 +145,8 @@ public:
             x.template segment<3>(3) /* edge */, dtype);
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         return point_edge_mollifier<T, 3>(
             x.template tail<3>() /* point */,
@@ -158,20 +162,21 @@ template <typename T> class PrimitiveDistanceTemplate<Point2, Point2, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Point2, Point2>::type dtype)
     {
         return (x.template tail<2>() - x.template head<2>()).squaredNorm();
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Point2, Point2>::type dtype)
     {
         return x.template tail<2>() - x.template head<2>();
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         return T(1.);
     }
@@ -184,20 +189,21 @@ template <typename T> class PrimitiveDistanceTemplate<Point3, Point3, T> {
 
 public:
     static T compute_distance(
-        const Vector<T, N_CORE_DOFS>& x,
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Point3, Point3>::type dtype)
     {
         return (x.template tail<3>() - x.template head<3>()).squaredNorm();
     }
 
-    static Vector<T, DIM> compute_closest_direction(
-        const Vector<T, N_CORE_DOFS>& x,
+    static Eigen::Vector<T, DIM> compute_closest_direction(
+        const Eigen::Vector<T, N_CORE_DOFS>& x,
         typename PrimitiveDistType<Point3, Point3>::type dtype)
     {
         return x.template tail<3>() - x.template head<3>();
     }
 
-    static T mollifier(const Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
+    static T
+    mollifier(const Eigen::Vector<T, N_CORE_DOFS>& x, const T& dist_sqr)
     {
         return T(1.);
     }
diff --git a/src/ipc/smooth_contact/primitives/edge.cpp b/src/ipc/smooth_contact/primitives/edge.cpp
index 8bf800454..15d50d354 100644
--- a/src/ipc/smooth_contact/primitives/edge.cpp
+++ b/src/ipc/smooth_contact/primitives/edge.cpp
@@ -8,7 +8,7 @@ template <int DIM>
 Edge<DIM>::Edge(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
diff --git a/src/ipc/smooth_contact/primitives/edge.hpp b/src/ipc/smooth_contact/primitives/edge.hpp
index d2878829d..e0711c7ed 100644
--- a/src/ipc/smooth_contact/primitives/edge.hpp
+++ b/src/ipc/smooth_contact/primitives/edge.hpp
@@ -9,7 +9,7 @@ template <int DIM> class Edge : public Primitive {
     static constexpr int N_CORE_POINTS = 2;
     using DVector = Eigen::Vector<double, DIM>;
     using XVector = Eigen::Vector<double, DIM == 2 ? 4 : 12>;
-    using GradType = Eigen::Vector<double, DIM == 2 ? 6 : 15>;
+    using GradientType = Eigen::Vector<double, DIM == 2 ? 6 : 15>;
     static constexpr int HESSIAN_ROWS = DIM == 2 ? 6 : 15;
     static constexpr int HESSIAN_COLS = HESSIAN_ROWS;
     using HessianType = Eigen::Matrix<double, HESSIAN_ROWS, HESSIAN_COLS>;
@@ -19,7 +19,7 @@ template <int DIM> class Edge : public Primitive {
     Edge(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
@@ -37,7 +37,8 @@ template <int DIM> class Edge : public Primitive {
     /// @param d Vector from closest point on the edge to the point outside of the edge
     /// @param x Positions of the two vertices of this edge
     /// @return Gradient of the potential wrt. d and x
-    GradType grad(Eigen::ConstRef<DVector> d, Eigen::ConstRef<XVector> x) const;
+    GradientType
+    grad(Eigen::ConstRef<DVector> d, Eigen::ConstRef<XVector> x) const;
 
     /// @brief Compute the Hessian of potential wrt. d and x
     /// @param d Vector from closest point on the edge to the point outside of the edge
diff --git a/src/ipc/smooth_contact/primitives/edge2.cpp b/src/ipc/smooth_contact/primitives/edge2.cpp
index 83fbc19a1..c009fdd56 100644
--- a/src/ipc/smooth_contact/primitives/edge2.cpp
+++ b/src/ipc/smooth_contact/primitives/edge2.cpp
@@ -6,7 +6,7 @@ namespace ipc {
 Edge2::Edge2(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
diff --git a/src/ipc/smooth_contact/primitives/edge2.hpp b/src/ipc/smooth_contact/primitives/edge2.hpp
index aa8dc217a..b52c13b43 100644
--- a/src/ipc/smooth_contact/primitives/edge2.hpp
+++ b/src/ipc/smooth_contact/primitives/edge2.hpp
@@ -14,7 +14,7 @@ class Edge2 : public Primitive {
     Edge2(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
diff --git a/src/ipc/smooth_contact/primitives/edge3.cpp b/src/ipc/smooth_contact/primitives/edge3.cpp
index 8377526f1..67ec74c7a 100644
--- a/src/ipc/smooth_contact/primitives/edge3.cpp
+++ b/src/ipc/smooth_contact/primitives/edge3.cpp
@@ -85,7 +85,7 @@ namespace {
         return tangent_term;
     }
 
-    GradType<15> smooth_edge3_tangent_term_gradient(
+    GradientType<15> smooth_edge3_tangent_term_gradient(
         Eigen::ConstRef<Eigen::Vector3d> dn,
         Eigen::ConstRef<Eigen::Vector3d> e0,
         Eigen::ConstRef<Eigen::Vector3d> e1,
@@ -97,7 +97,7 @@ namespace {
     {
         Eigen::Vector2d vals;
         vals << 1., 1.;
-        std::array<Vector<double, 15>, 2> grads;
+        std::array<Eigen::Vector<double, 15>, 2> grads;
         for (auto& g : grads) {
             g.setZero();
         }
@@ -112,11 +112,11 @@ namespace {
 
                 vals[d] = tmp_val;
 
-                Vector<double, 12> gradient_tmp;
+                Eigen::Vector<double, 12> gradient_tmp;
                 gradient_tmp << -tmp_grad.tail<3>(),
                     g.transpose() * tmp_grad.head<3>();
 
-                Vector<int, 12> indices;
+                Eigen::Vector<int, 12> indices;
                 indices << 0, 1, 2, 9, 10, 11, 3, 4, 5, 6, 7, 8;
                 if (d == 1) {
                     indices.segment<3>(3).array() += 3;
@@ -142,7 +142,7 @@ namespace {
     {
         Eigen::Vector2d vals;
         vals << 1., 1.;
-        std::array<Vector<double, 15>, 2> grads;
+        std::array<Eigen::Vector<double, 15>, 2> grads;
         std::array<Eigen::Matrix<double, 15, 15>, 2> hesses;
         for (auto& g : grads) {
             g.setZero();
@@ -164,7 +164,7 @@ namespace {
 
                 vals[d] = tmp_val;
 
-                Vector<double, 12> gradient_tmp;
+                Eigen::Vector<double, 12> gradient_tmp;
                 gradient_tmp << tmp_grad.tail<3>(),
                     g.transpose() * tmp_grad.head<3>();
 
@@ -182,7 +182,7 @@ namespace {
                 hessian_tmp.block<9, 3>(3, 0) =
                     g.transpose() * tmp_hess.block<3, 3>(0, 3);
 
-                Vector<int, 12> indices;
+                Eigen::Vector<int, 12> indices;
                 indices << 0, 1, 2, 9, 10, 11, 3, 4, 5, 6, 7, 8;
                 if (d == 1) {
                     indices.segment<3>(3).array() += 3;
@@ -221,7 +221,7 @@ namespace {
         return (e1 - e0).squaredNorm() * tangent_term * normal_term;
     }
 
-    GradType<15> smooth_edge3_term_gradient(
+    GradientType<15> smooth_edge3_term_gradient(
         Eigen::ConstRef<Eigen::Vector3d> direc,
         Eigen::ConstRef<Eigen::Vector3d> e0,
         Eigen::ConstRef<Eigen::Vector3d> e1,
@@ -383,7 +383,7 @@ namespace {
 Edge3::Edge3(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
@@ -550,7 +550,7 @@ double smooth_edge3_normal_term(
         (d - t0).cross(d - t1).dot(edge), alpha, beta);
 }
 
-GradType<15> smooth_edge3_normal_term_gradient(
+GradientType<15> smooth_edge3_normal_term_gradient(
     Eigen::ConstRef<Eigen::Vector3d> dn,
     Eigen::ConstRef<Eigen::Vector3d> e0,
     Eigen::ConstRef<Eigen::Vector3d> e1,
diff --git a/src/ipc/smooth_contact/primitives/edge3.hpp b/src/ipc/smooth_contact/primitives/edge3.hpp
index 6fbff3b25..32ce6f307 100644
--- a/src/ipc/smooth_contact/primitives/edge3.hpp
+++ b/src/ipc/smooth_contact/primitives/edge3.hpp
@@ -14,7 +14,7 @@ class Edge3 : public Primitive {
     Edge3(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
@@ -48,7 +48,7 @@ double smooth_edge3_normal_term(
     const double beta,
     const OrientationTypes& otypes);
 
-GradType<15> smooth_edge3_normal_term_gradient(
+GradientType<15> smooth_edge3_normal_term_gradient(
     Eigen::ConstRef<Eigen::Vector3d> dn,
     Eigen::ConstRef<Eigen::Vector3d> e0,
     Eigen::ConstRef<Eigen::Vector3d> e1,
diff --git a/src/ipc/smooth_contact/primitives/face.cpp b/src/ipc/smooth_contact/primitives/face.cpp
index 54b4e76fc..2d1254e85 100644
--- a/src/ipc/smooth_contact/primitives/face.cpp
+++ b/src/ipc/smooth_contact/primitives/face.cpp
@@ -23,7 +23,7 @@ namespace ipc {
 Face::Face(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
diff --git a/src/ipc/smooth_contact/primitives/face.hpp b/src/ipc/smooth_contact/primitives/face.hpp
index 4372f1f9a..962457acd 100644
--- a/src/ipc/smooth_contact/primitives/face.hpp
+++ b/src/ipc/smooth_contact/primitives/face.hpp
@@ -14,7 +14,7 @@ class Face : public Primitive {
     Face(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
diff --git a/src/ipc/smooth_contact/primitives/point2.cpp b/src/ipc/smooth_contact/primitives/point2.cpp
index 07635d232..52f415dfa 100644
--- a/src/ipc/smooth_contact/primitives/point2.cpp
+++ b/src/ipc/smooth_contact/primitives/point2.cpp
@@ -106,7 +106,7 @@ namespace {
 Point2::Point2(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
@@ -142,7 +142,8 @@ Point2::Point2(
 int Point2::n_vertices() const { return m_vertex_ids.size(); }
 
 double Point2::potential(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
     if (has_neighbor_1 && has_neighbor_2) {
         return smooth_point2_term<double>(
@@ -155,13 +156,14 @@ double Point2::potential(
         return 1.;
     }
 }
-Vector<double, -1, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
     if (has_neighbor_1 && has_neighbor_2) {
         ScalarBase::setVariableCount(4 * DIM);
         using T = ADGrad<4 * DIM>;
-        Vector<double, 4 * DIM> tmp;
+        Eigen::Vector<double, 4 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         return smooth_point2_term<T>(
@@ -170,14 +172,14 @@ Vector<double, -1, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
         using T = ADGrad<3 * DIM>;
-        Vector<double, 3 * DIM> tmp;
+        Eigen::Vector<double, 3 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
         return smooth_point2_term_one_side<T>(
                    X.row(1), X.row(0), X.row(2), m_params)
             .grad;
     } else {
-        return Vector<double, -1, Point2::MAX_SIZE + Point2::DIM>::Zero(
+        return VectorMax<double, Point2::MAX_SIZE + Point2::DIM>::Zero(
             x.size() + d.size());
     }
 }
@@ -186,12 +188,13 @@ MatrixMax<
     Point2::MAX_SIZE + Point2::DIM,
     Point2::MAX_SIZE + Point2::DIM>
 Point2::hessian(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
     if (has_neighbor_1 && has_neighbor_2) {
         ScalarBase::setVariableCount(4 * DIM);
         using T = ADHessian<4 * DIM>;
-        Vector<double, 4 * DIM> tmp;
+        Eigen::Vector<double, 4 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         return smooth_point2_term<T>(
@@ -200,7 +203,7 @@ Point2::hessian(
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
         using T = ADHessian<3 * DIM>;
-        Vector<double, 3 * DIM> tmp;
+        Eigen::Vector<double, 3 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
         return smooth_point2_term_one_side<T>(
diff --git a/src/ipc/smooth_contact/primitives/point2.hpp b/src/ipc/smooth_contact/primitives/point2.hpp
index 125bb4f4f..120468a73 100644
--- a/src/ipc/smooth_contact/primitives/point2.hpp
+++ b/src/ipc/smooth_contact/primitives/point2.hpp
@@ -14,29 +14,29 @@ class Point2 : public Primitive {
     Point2(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
     Point2(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices);
+        Eigen::ConstRef<Eigen::MatrixXd> vertices);
 
     int n_vertices() const override;
     int n_dofs() const override { return n_vertices() * DIM; }
 
     // assume the following functions are only called if active
     double potential(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
     // derivatives including wrt. d (the closest direction) in front
-    Vector<double, -1, MAX_SIZE + DIM> grad(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+    VectorMax<double, MAX_SIZE + DIM> grad(
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
     MatrixMax<double, MAX_SIZE + DIM, MAX_SIZE + DIM> hessian(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
 
 private:
     bool has_neighbor_1, has_neighbor_2;
diff --git a/src/ipc/smooth_contact/primitives/point3.cpp b/src/ipc/smooth_contact/primitives/point3.cpp
index e9a6099bc..e82fd3a65 100644
--- a/src/ipc/smooth_contact/primitives/point3.cpp
+++ b/src/ipc/smooth_contact/primitives/point3.cpp
@@ -9,7 +9,7 @@ namespace ipc {
 Point3::Point3(
     const index_t id,
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const VectorMax3d& d,
     const SmoothContactParameters& params)
     : Primitive(id, params)
@@ -74,15 +74,17 @@ Point3::Point3(
 int Point3::n_vertices() const { return local_to_global_vids.size(); }
 
 double Point3::potential(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
     const Eigen::Matrix<double, -1, DIM> X =
         slice_positions<double, -1, DIM>(x);
     return smooth_point3_term<double, -1>(X, d);
 }
 
-Vector<double, -1, Point3::MAX_SIZE + Point3::DIM> Point3::grad(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+VectorMax<double, Point3::MAX_SIZE + Point3::DIM> Point3::grad(
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADGrad<-1>;
@@ -104,7 +106,8 @@ MatrixMax<
     Point3::MAX_SIZE + Point3::DIM,
     Point3::MAX_SIZE + Point3::DIM>
 Point3::hessian(
-    const Vector<double, DIM>& d, const Vector<double, -1, MAX_SIZE>& x) const
+    const Eigen::Vector<double, DIM>& d,
+    const VectorMax<double, MAX_SIZE>& x) const
 {
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
     using T = ADHessian<-1>;
@@ -120,7 +123,7 @@ Point3::hessian(
 #endif
 }
 
-GradType<-1> Point3::smooth_point3_term_tangent_gradient(
+GradientType<-1> Point3::smooth_point3_term_tangent_gradient(
     Eigen::ConstRef<Eigen::RowVector3d> direc,
     Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> tangents,
     const double alpha,
@@ -231,7 +234,7 @@ HessianType<-1> Point3::smooth_point3_term_tangent_hessian(
     return std::make_tuple(values.prod(), tangent_grad, tangent_hess);
 }
 
-GradType<-1> Point3::smooth_point3_term_normal_gradient(
+GradientType<-1> Point3::smooth_point3_term_normal_gradient(
     Eigen::ConstRef<Eigen::RowVector3d> direc,
     Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> tangents,
     const double alpha,
@@ -400,7 +403,7 @@ bool Point3::smooth_point3_term_type(
     return normal_term > 0;
 }
 
-GradType<-1> Point3::smooth_point3_term_gradient(
+GradientType<-1> Point3::smooth_point3_term_gradient(
     Eigen::ConstRef<Eigen::RowVector3d> direc,
     Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> X,
     const SmoothContactParameters& params) const
diff --git a/src/ipc/smooth_contact/primitives/point3.hpp b/src/ipc/smooth_contact/primitives/point3.hpp
index 424cbb960..32407d7cd 100644
--- a/src/ipc/smooth_contact/primitives/point3.hpp
+++ b/src/ipc/smooth_contact/primitives/point3.hpp
@@ -14,14 +14,14 @@ class Point3 : public Primitive {
     Point3(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const VectorMax3d& d,
         const SmoothContactParameters& params);
 
     Point3(
         const index_t id,
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices);
+        Eigen::ConstRef<Eigen::MatrixXd> vertices);
     virtual ~Point3() = default;
 
     int n_vertices() const override;
@@ -29,15 +29,15 @@ class Point3 : public Primitive {
 
     // assume the following functions are only called if active
     double potential(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
     // derivatives including wrt. d (the closest direction) in front
-    Vector<double, -1, MAX_SIZE + DIM> grad(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+    VectorMax<double, MAX_SIZE + DIM> grad(
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
     MatrixMax<double, MAX_SIZE + DIM, MAX_SIZE + DIM> hessian(
-        const Vector<double, DIM>& d,
-        const Vector<double, -1, MAX_SIZE>& x) const;
+        const Eigen::Vector<double, DIM>& d,
+        const VectorMax<double, MAX_SIZE>& x) const;
 
     /// @brief
     /// @tparam scalar
@@ -52,7 +52,7 @@ class Point3 : public Primitive {
         const Eigen::Matrix<scalar, n_verts, 3>& X,
         Eigen::ConstRef<Eigen::RowVector3<scalar>> direc) const;
 
-    GradType<-1> smooth_point3_term_gradient(
+    GradientType<-1> smooth_point3_term_gradient(
         Eigen::ConstRef<Eigen::RowVector3d> direc,
         Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> X,
         const SmoothContactParameters& params) const;
@@ -62,7 +62,7 @@ class Point3 : public Primitive {
         Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> X,
         const SmoothContactParameters& params) const;
 
-    GradType<-1> smooth_point3_term_tangent_gradient(
+    GradientType<-1> smooth_point3_term_tangent_gradient(
         Eigen::ConstRef<Eigen::RowVector3d> direc,
         Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> tangents,
         const double alpha,
@@ -74,7 +74,7 @@ class Point3 : public Primitive {
         const double alpha,
         const double beta) const;
 
-    GradType<-1> smooth_point3_term_normal_gradient(
+    GradientType<-1> smooth_point3_term_normal_gradient(
         Eigen::ConstRef<Eigen::RowVector3d> direc,
         Eigen::ConstRef<Eigen::Matrix<double, -1, 3>> tangents,
         const double alpha,
diff --git a/src/ipc/smooth_contact/smooth_collisions_builder.cpp b/src/ipc/smooth_contact/smooth_collisions_builder.cpp
index 531b58fd2..40701c969 100644
--- a/src/ipc/smooth_contact/smooth_collisions_builder.cpp
+++ b/src/ipc/smooth_contact/smooth_collisions_builder.cpp
@@ -38,7 +38,7 @@ namespace {
 
 void SmoothCollisionsBuilder<2>::add_edge_vertex_collisions(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const std::vector<EdgeVertexCandidate>& candidates,
     const SmoothContactParameters& params,
     const std::function<double(const index_t)>& vert_dhat,
@@ -80,7 +80,7 @@ void SmoothCollisionsBuilder<2>::add_edge_vertex_collisions(
 
 void SmoothCollisionsBuilder<3>::add_edge_edge_collisions(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const std::vector<EdgeEdgeCandidate>& candidates,
     const SmoothContactParameters& params,
     const std::function<double(const index_t)>& vert_dhat,
@@ -115,7 +115,7 @@ void SmoothCollisionsBuilder<3>::add_edge_edge_collisions(
 
 void SmoothCollisionsBuilder<3>::add_face_vertex_collisions(
     const CollisionMesh& mesh,
-    const Eigen::MatrixXd& vertices,
+    Eigen::ConstRef<Eigen::MatrixXd> vertices,
     const std::vector<FaceVertexCandidate>& candidates,
     const SmoothContactParameters& params,
     const std::function<double(const index_t)>& vert_dhat,
diff --git a/src/ipc/smooth_contact/smooth_collisions_builder.hpp b/src/ipc/smooth_contact/smooth_collisions_builder.hpp
index ef1804057..c03afe544 100644
--- a/src/ipc/smooth_contact/smooth_collisions_builder.hpp
+++ b/src/ipc/smooth_contact/smooth_collisions_builder.hpp
@@ -17,7 +17,7 @@ template <> class SmoothCollisionsBuilder<2> {
 
     void add_edge_vertex_collisions(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const std::vector<EdgeVertexCandidate>& candidates,
         const SmoothContactParameters& params,
         const std::function<double(const index_t)>& vert_dhat,
@@ -53,7 +53,7 @@ template <> class SmoothCollisionsBuilder<3> {
 
     void add_edge_edge_collisions(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const std::vector<EdgeEdgeCandidate>& candidates,
         const SmoothContactParameters& params,
         const std::function<double(const index_t)>& vert_dhat,
@@ -63,7 +63,7 @@ template <> class SmoothCollisionsBuilder<3> {
 
     void add_face_vertex_collisions(
         const CollisionMesh& mesh,
-        const Eigen::MatrixXd& vertices,
+        Eigen::ConstRef<Eigen::MatrixXd> vertices,
         const std::vector<FaceVertexCandidate>& candidates,
         const SmoothContactParameters& params,
         const std::function<double(const index_t)>& vert_dhat,
diff --git a/src/ipc/utils/eigen_ext.hpp b/src/ipc/utils/eigen_ext.hpp
index e923e4093..3a74f8597 100644
--- a/src/ipc/utils/eigen_ext.hpp
+++ b/src/ipc/utils/eigen_ext.hpp
@@ -32,15 +32,17 @@ using MatrixXb = Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>;
 /// @tparam T The type of the vector elements.
 /// @tparam dim The size of the vector.
 /// @tparam max_dim The maximum size of the vector.
-template <typename T, int dim, int max_dim = dim>
-using Vector = Eigen::Matrix<T, dim, 1, Eigen::ColMajor, max_dim, 1>;
+template <typename T, int max_dim>
+using VectorMax =
+    Eigen::Matrix<T, Eigen::Dynamic, 1, Eigen::ColMajor, max_dim, 1>;
 
 /// @brief A dynamic size row vector with a fixed maximum size.
 /// @tparam T The type of the vector elements.
 /// @tparam dim The size of the vector.
 /// @tparam max_dim The maximum size of the vector.
-template <typename T, int dim, int max_dim = dim>
-using RowVector = Eigen::Matrix<T, 1, dim, Eigen::RowMajor, 1, max_dim>;
+template <typename T, int max_dim>
+using RowVectorMax =
+    Eigen::Matrix<T, 1, Eigen::Dynamic, Eigen::RowMajor, 1, max_dim>;
 
 /// @brief A static size matrix of size of 1×1
 using Vector1d = Eigen::Vector<double, 1>;
@@ -63,17 +65,17 @@ using Matrix12d = Eigen::Matrix<double, 12, 12>;
 using Matrix15d = Eigen::Matrix<double, 15, 15>;
 
 /// @brief A dynamic size matrix with a fixed maximum size of 3×1
-template <typename T> using VectorMax2 = Vector<T, Eigen::Dynamic, 2>;
+template <typename T> using VectorMax2 = VectorMax<T, 2>;
 /// @brief A dynamic size matrix with a fixed maximum size of 3×1
-template <typename T> using VectorMax3 = Vector<T, Eigen::Dynamic, 3>;
+template <typename T> using VectorMax3 = VectorMax<T, 3>;
 /// @brief A dynamic size matrix with a fixed maximum size of 4×1
-template <typename T> using VectorMax4 = Vector<T, Eigen::Dynamic, 4>;
+template <typename T> using VectorMax4 = VectorMax<T, 4>;
 /// @brief A dynamic size matrix with a fixed maximum size of 6×1
-template <typename T> using VectorMax6 = Vector<T, Eigen::Dynamic, 6>;
+template <typename T> using VectorMax6 = VectorMax<T, 6>;
 /// @brief A dynamic size matrix with a fixed maximum size of 9×1
-template <typename T> using VectorMax9 = Vector<T, Eigen::Dynamic, 9>;
+template <typename T> using VectorMax9 = VectorMax<T, 9>;
 /// @brief A dynamic size matrix with a fixed maximum size of 12×1
-template <typename T> using VectorMax12 = Vector<T, Eigen::Dynamic, 12>;
+template <typename T> using VectorMax12 = VectorMax<T, 12>;
 
 /// @brief A dynamic size matrix with a fixed maximum size of 2×1
 using VectorMax2d = VectorMax2<double>;
@@ -95,20 +97,20 @@ using VectorMax9d = VectorMax9<double>;
 using VectorMax12d = VectorMax12<double>;
 
 /// @brief A dynamic size matrix with a fixed maximum size of 1×2
-template <typename T> using RowVectorMax2 = RowVector<T, Eigen::Dynamic, 2>;
+template <typename T> using RowVectorMax2 = RowVectorMax<T, 2>;
 /// @brief A dynamic size matrix with a fixed maximum size of 1×3
-template <typename T> using RowVectorMax3 = RowVector<T, Eigen::Dynamic, 3>;
+template <typename T> using RowVectorMax3 = RowVectorMax<T, 3>;
 
 /// @brief A dynamic size matrix with a fixed maximum size of 1×2
 using RowVectorMax2d = RowVectorMax2<double>;
 /// @brief A dynamic size matrix with a fixed maximum size of 1×3
 using RowVectorMax3d = RowVectorMax3<double>;
 /// @brief A dynamic size matrix with a fixed maximum size of 6×1
-using RowVectorMax6d = RowVector<double, Eigen::Dynamic, 6>;
+using RowVectorMax6d = RowVectorMax<double, 6>;
 /// @brief A dynamic size matrix with a fixed maximum size of 9×1
-using RowVectorMax9d = RowVector<double, Eigen::Dynamic, 9>;
+using RowVectorMax9d = RowVectorMax<double, 9>;
 /// @brief A dynamic size matrix with a fixed maximum size of 12×1
-using RowVectorMax12d = RowVector<double, Eigen::Dynamic, 12>;
+using RowVectorMax12d = RowVectorMax<double, 12>;
 
 template <typename T, int max_rows, int max_cols>
 using MatrixMax = Eigen::Matrix<
@@ -167,10 +169,11 @@ using ArrayMax4d = ArrayMax4<double>;
 /// @brief A dynamic size array with a fixed maximum size of 4×1
 using ArrayMax4i = ArrayMax4<int>;
 
-template <int dim> using GradType = std::tuple<double, Vector<double, dim>>;
 template <int dim>
-using HessianType =
-    std::tuple<double, Vector<double, dim>, Eigen::Matrix<double, dim, dim>>;
+using GradientType = std::tuple<double, Eigen::Vector<double, dim>>;
+template <int dim>
+using HessianType = std::
+    tuple<double, Eigen::Vector<double, dim>, Eigen::Matrix<double, dim, dim>>;
 
 /**@}*/
 
diff --git a/tests/src/tests/barrier/test_barrier.cpp b/tests/src/tests/barrier/test_barrier.cpp
index 1416f647c..75ec106d8 100644
--- a/tests/src/tests/barrier/test_barrier.cpp
+++ b/tests/src/tests/barrier/test_barrier.cpp
@@ -139,8 +139,8 @@ TEST_CASE("Normalize vector derivatives", "[deriv]")
     using T = ipc::ADHessian<3>;
     for (int i = 1; i <= n_samples; i++) {
         Eigen::Vector3d x = Eigen::Vector3d::Random();
-        ipc::Vector<T, 3> x_ad = ipc::slice_positions<T, 3, 1>(x);
-        ipc::Vector<T, 3> y_ad = x_ad / x_ad.norm();
+        Eigen::Vector<T, 3> x_ad = ipc::slice_positions<T, 3, 1>(x);
+        Eigen::Vector<T, 3> y_ad = x_ad / x_ad.norm();
 
         const auto [y, grad, hess] =
             ipc::normalization_and_jacobian_and_hessian(x);
@@ -165,11 +165,11 @@ TEST_CASE("line-line closest direction derivatives", "[deriv]")
     using T = ipc::ADHessian<12>;
     for (int i = 1; i <= n_samples; i++) {
         ipc::Vector6d ea = ipc::Vector6d::Random();
-        ipc::Vector<T, 6> eaT = ipc::slice_positions<T, 6, 1>(ea);
+        Eigen::Vector<T, 6> eaT = ipc::slice_positions<T, 6, 1>(ea);
         ipc::Vector6d eb = ipc::Vector6d::Random();
-        ipc::Vector<T, 6> ebT = ipc::slice_positions<T, 6, 1>(eb, 6);
+        Eigen::Vector<T, 6> ebT = ipc::slice_positions<T, 6, 1>(eb, 6);
 
-        ipc::Vector<T, 3> dT = ipc::line_line_closest_point_direction<T>(
+        Eigen::Vector<T, 3> dT = ipc::line_line_closest_point_direction<T>(
             eaT.head<3>(), eaT.tail<3>(), ebT.head<3>(), ebT.tail<3>());
 
         const auto [d1, grad1] =
@@ -201,7 +201,7 @@ TEST_CASE("opposite_direction_penalty derivatives", "[deriv]")
     const double beta = 1;
     for (int i = 1; i <= n_samples; i++) {
         ipc::Vector6d x = ipc::Vector6d::Random();
-        ipc::Vector<T, 6> x_ad = ipc::slice_positions<T, 6, 1>(x);
+        Eigen::Vector<T, 6> x_ad = ipc::slice_positions<T, 6, 1>(x);
         T y_ad = ipc::Math<T>::smooth_heaviside(
             x_ad.tail(3).dot(x_ad.head(3)) / x_ad.head(3).norm(), alpha, beta);
 
@@ -222,8 +222,8 @@ TEST_CASE("negative_orientation_penalty derivatives", "[deriv]")
         ScalarBase::setVariableCount(9);
         using T = ipc::ADHessian<9>;
         ipc::Vector9d x = ipc::Vector9d::Random();
-        ipc::Vector<T, 9> x_ad = ipc::slice_positions<T, 9, 1>(x);
-        ipc::Vector<T, 3> t = x_ad.head<3>().cross(x_ad.segment<3>(3));
+        Eigen::Vector<T, 9> x_ad = ipc::slice_positions<T, 9, 1>(x);
+        Eigen::Vector<T, 3> t = x_ad.head<3>().cross(x_ad.segment<3>(3));
         T y_ad = ipc::Math<T>::smooth_heaviside(
             x_ad.tail(3).dot(t) / t.norm(), alpha, beta);
 
diff --git a/tests/src/tests/distance/test_edge_edge.cpp b/tests/src/tests/distance/test_edge_edge.cpp
index 5a9a35d0e..9c2d5c7a6 100644
--- a/tests/src/tests/distance/test_edge_edge.cpp
+++ b/tests/src/tests/distance/test_edge_edge.cpp
@@ -294,7 +294,7 @@ TEST_CASE("Edge-edge distance gradient", "[distance][edge-edge][gradient]")
     CAPTURE(e0x, e0y, e0z, edge_edge_distance_type(e00, e01, e10, e11));
     CHECK(fd::compare_gradient(grad, fgrad));
 
-    Vector<ADGrad<12>, 12> X1 = slice_positions<ADGrad<12>, 12, 1>(
+    Eigen::Vector<ADGrad<12>, 12> X1 = slice_positions<ADGrad<12>, 12, 1>(
         (Vector12d() << e00, e01, e10, e11).finished());
     ADGrad<12> dist1 =
         PrimitiveDistanceTemplate<Edge3, Edge3, ADGrad<12>>::compute_distance(
@@ -397,8 +397,8 @@ TEST_CASE(
                 x, dtype);
         CHECK(distance == Catch::Approx(s * s).margin(1e-15));
 
-        Vector<ADGrad<12>, 12> X = slice_positions<ADGrad<12>, 12, 1>(x);
-        Vector<ADGrad<12>, 3> direc = PrimitiveDistanceTemplate<
+        Eigen::Vector<ADGrad<12>, 12> X = slice_positions<ADGrad<12>, 12, 1>(x);
+        Eigen::Vector<ADGrad<12>, 3> direc = PrimitiveDistanceTemplate<
             Edge3, Edge3, ADGrad<12>>::compute_closest_direction(X, dtype);
         CHECK(direc.squaredNorm().val == Catch::Approx(s * s).margin(1e-15));
 
diff --git a/tests/src/tests/distance/test_point_edge.cpp b/tests/src/tests/distance/test_point_edge.cpp
index 5adc739b8..bec5f7b95 100644
--- a/tests/src/tests/distance/test_point_edge.cpp
+++ b/tests/src/tests/distance/test_point_edge.cpp
@@ -253,9 +253,9 @@ TEMPLATE_TEST_CASE_SIG(
         fd::finite_jacobian(
             x,
             [](const Eigen::VectorXd& x) {
-                Vector<double, dim> p = x.segment<dim>(0);
-                Vector<double, dim> e0 = x.segment<dim>(dim);
-                Vector<double, dim> e1 = x.segment<dim>(2 * dim);
+                Eigen::Vector<double, dim> p = x.segment<dim>(0);
+                Eigen::Vector<double, dim> e0 = x.segment<dim>(dim);
+                Eigen::Vector<double, dim> e1 = x.segment<dim>(2 * dim);
                 return PointEdgeDistance<
                     double, dim>::point_edge_closest_point_direction(p, e0, e1);
             },
@@ -269,7 +269,7 @@ TEMPLATE_TEST_CASE_SIG(
         VectorMax9d x(3 * dim);
         x << e0, e1, p;
 
-        Vector<ADGrad<3 * dim>, 3 * dim> X =
+        Eigen::Vector<ADGrad<3 * dim>, 3 * dim> X =
             slice_positions<ADGrad<3 * dim>, 3 * dim, 1>(x);
         ADGrad<3 * dim> dist;
         if constexpr (dim == 2) {
@@ -285,9 +285,9 @@ TEMPLATE_TEST_CASE_SIG(
         fd::finite_gradient(
             x,
             [](const Eigen::VectorXd& x) {
-                Vector<double, dim> e0 = x.segment<dim>(0);
-                Vector<double, dim> e1 = x.segment<dim>(dim);
-                Vector<double, dim> p = x.segment<dim>(2 * dim);
+                Eigen::Vector<double, dim> e0 = x.segment<dim>(0);
+                Eigen::Vector<double, dim> e1 = x.segment<dim>(dim);
+                Eigen::Vector<double, dim> p = x.segment<dim>(2 * dim);
                 return PointEdgeDistance<double, dim>::point_edge_sqr_distance(
                     p, e0, e1);
             },
@@ -308,9 +308,9 @@ TEMPLATE_TEST_CASE_SIG(
             fd::finite_hessian(
                 x,
                 [i, dtype](const Eigen::VectorXd& x) {
-                    Vector<double, dim> p = x.segment<dim>(0);
-                    Vector<double, dim> e0 = x.segment<dim>(dim);
-                    Vector<double, dim> e1 = x.segment<dim>(2 * dim);
+                    Eigen::Vector<double, dim> p = x.segment<dim>(0);
+                    Eigen::Vector<double, dim> e0 = x.segment<dim>(dim);
+                    Eigen::Vector<double, dim> e1 = x.segment<dim>(2 * dim);
                     return PointEdgeDistance<double, dim>::
                         point_edge_closest_point_direction(p, e0, e1, dtype)(i);
                 },
diff --git a/tests/src/tests/distance/test_point_point.cpp b/tests/src/tests/distance/test_point_point.cpp
index 1a0dc11bd..9c4cc2f19 100644
--- a/tests/src/tests/distance/test_point_point.cpp
+++ b/tests/src/tests/distance/test_point_point.cpp
@@ -66,7 +66,7 @@ TEMPLATE_TEST_CASE_SIG(
     CHECK(fd::compare_gradient(grad, fgrad));
 
     constexpr int n_dofs = 2 * dim;
-    Vector<ADGrad<n_dofs>, n_dofs> X =
+    Eigen::Vector<ADGrad<n_dofs>, n_dofs> X =
         slice_positions<ADGrad<n_dofs>, n_dofs, 1>(x);
     ADGrad<n_dofs> dist;
     if constexpr (dim == 2) {
diff --git a/tests/src/tests/distance/test_point_triangle.cpp b/tests/src/tests/distance/test_point_triangle.cpp
index 3ed910814..7d16f11ad 100644
--- a/tests/src/tests/distance/test_point_triangle.cpp
+++ b/tests/src/tests/distance/test_point_triangle.cpp
@@ -191,10 +191,10 @@ TEST_CASE(
 
     CHECK(fd::compare_gradient(grad, fgrad));
 
-    Vector<ADGrad<12>, 12> X = slice_positions<ADGrad<12>, 12, 1>(
+    Eigen::Vector<ADGrad<12>, 12> X = slice_positions<ADGrad<12>, 12, 1>(
         (Vector12d() << p, t0, t1, t2).finished());
     {
-        Vector<ADGrad<12>, 3> tmp = X.head<3>();
+        Eigen::Vector<ADGrad<12>, 3> tmp = X.head<3>();
         X.head<9>() = X.tail<9>().eval();
         X.tail<3>() = tmp;
     }
@@ -288,9 +288,10 @@ TEST_CASE("Point-triangle distance hessian", "[distance][point-triangle][hess]")
     CAPTURE(dtype);
     CHECK(fd::compare_hessian(hess, fhess, 1e-2));
 
-    Vector<ADHessian<12>, 12> X = slice_positions<ADHessian<12>, 12, 1>(x);
+    Eigen::Vector<ADHessian<12>, 12> X =
+        slice_positions<ADHessian<12>, 12, 1>(x);
     {
-        Vector<ADHessian<12>, 3> tmp = X.head<3>();
+        Eigen::Vector<ADHessian<12>, 3> tmp = X.head<3>();
         X.head<9>() = X.tail<9>().eval();
         X.tail<3>() = tmp;
     }