diff --git a/src/ipc/barrier/CMakeLists.txt b/src/ipc/barrier/CMakeLists.txt
index 70f573cca..6cfa28ae1 100644
--- a/src/ipc/barrier/CMakeLists.txt
+++ b/src/ipc/barrier/CMakeLists.txt
@@ -3,8 +3,13 @@ set(SOURCES
   adaptive_stiffness.hpp
   barrier_force_magnitude.cpp
   barrier_force_magnitude.hpp
-  barrier.cpp
   barrier.hpp
 )
 
-target_sources(ipc_toolkit PRIVATE ${SOURCES})
\ No newline at end of file
+target_sources(ipc_toolkit PRIVATE ${SOURCES})
+
+# barrier.cpp is shared between host C++ and CUDA device code; under CUDA the
+# helper generates a .cu wrapper (see ipc_toolkit_shared_device_sources).
+ipc_toolkit_target_shared_device_sources(ipc_toolkit
+  barrier.cpp
+)
diff --git a/src/ipc/barrier/barrier.cpp b/src/ipc/barrier/barrier.cpp
index 38221d975..5538430f2 100644
--- a/src/ipc/barrier/barrier.cpp
+++ b/src/ipc/barrier/barrier.cpp
@@ -10,12 +10,16 @@
 
 namespace ipc {
 
+// ============================================================================
+// Free barrier functions -- shared between host C++ and CUDA device code.
+// ============================================================================
+//
 // Each barrier is one select_lazy cascade, ordered by increasing d so it
 // reads like the piecewise definition in the header. A scalar evaluates only
 // the case it lands in -- so the log below is never reached for d <= 0 -- while
 // a batch evaluates every case and blends per-lane, earlier cases winning.
 
-template <typename T> T barrier(const T d, const T dhat)
+template <typename T> IPC_TOOLKIT_HOST_DEVICE T barrier(const T d, const T dhat)
 {
     using namespace ipc::numext; // log
     // b(d) = -(d-d̂)²ln(d / d̂)
@@ -25,7 +29,8 @@ template <typename T> T barrier(const T d, const T dhat)
         [&] { return T(0); });
 }
 
-template <typename T> T barrier_first_derivative(const T d, const T dhat)
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE T barrier_first_derivative(const T d, const T dhat)
 {
     using namespace ipc::numext; // log
     // b(d) = -(d - d̂)²ln(d / d̂)
@@ -39,7 +44,8 @@ template <typename T> T barrier_first_derivative(const T d, const T dhat)
         [&] { return T(0); });
 }
 
-template <typename T> T barrier_second_derivative(const T d, const T dhat)
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE T barrier_second_derivative(const T d, const T dhat)
 {
     using namespace ipc::numext; // log
     return select_lazy(
@@ -51,8 +57,45 @@ template <typename T> T barrier_second_derivative(const T d, const T dhat)
         },
         [&] { return T(0); });
 }
+// ============================================================================
+// Explicit template instantiations
+/// @cond DOXYGEN_SKIP
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
+template float barrier(const float d, const float dhat);
+template double barrier(const double d, const double dhat);
+template float barrier_first_derivative(const float d, const float dhat);
+template double barrier_first_derivative(const double d, const double dhat);
+template float barrier_second_derivative(const float d, const float dhat);
+template double barrier_second_derivative(const double d, const double dhat);
+#endif
+#ifdef IPC_TOOLKIT_WITH_SIMD
+template SimdBatch<float>
+barrier(const SimdBatch<float> d, const SimdBatch<float> dhat);
+template SimdBatch<double>
+barrier(const SimdBatch<double> d, const SimdBatch<double> dhat);
+template SimdBatch<float>
+barrier_first_derivative(const SimdBatch<float> d, const SimdBatch<float> dhat);
+template SimdBatch<double> barrier_first_derivative(
+    const SimdBatch<double> d, const SimdBatch<double> dhat);
+template SimdBatch<float> barrier_second_derivative(
+    const SimdBatch<float> d, const SimdBatch<float> dhat);
+template SimdBatch<double> barrier_second_derivative(
+    const SimdBatch<double> d, const SimdBatch<double> dhat);
+#endif
+/// @endcond
 
 // ============================================================================
+// Barrier class hierarchy -- host only.
+// ============================================================================
+//
+// The classes are a virtual dispatch layer over the free functions above, and
+// virtual dispatch cannot cross the host/device boundary: a vtable built on the
+// host holds host code addresses, CUDA forbids passing an object of a class
+// with virtual functions to a __global__ function, and BarrierPotential owns
+// its barrier through a host-only std::shared_ptr. Skipping the hierarchy in
+// the device pass also keeps every class symbol -- including the float and
+// double ones -- in the host object, so each is emitted exactly once.
+#ifndef __CUDACC__
 
 template <typename T>
 T ClampedLogSqBarrier<T>::operator()(const T d, const T dhat) const
@@ -181,12 +224,6 @@ template class CubicBarrier<float>;
 template class CubicBarrier<double>;
 template class TwoStageBarrier<float>;
 template class TwoStageBarrier<double>;
-template float barrier(const float d, const float dhat);
-template double barrier(const double d, const double dhat);
-template float barrier_first_derivative(const float d, const float dhat);
-template double barrier_first_derivative(const double d, const double dhat);
-template float barrier_second_derivative(const float d, const float dhat);
-template double barrier_second_derivative(const double d, const double dhat);
 #ifdef IPC_TOOLKIT_WITH_SIMD
 template class BarrierBase<SimdBatch<float>>;
 template class BarrierBase<SimdBatch<double>>;
@@ -198,20 +235,10 @@ template class CubicBarrier<SimdBatch<float>>;
 template class CubicBarrier<SimdBatch<double>>;
 template class TwoStageBarrier<SimdBatch<float>>;
 template class TwoStageBarrier<SimdBatch<double>>;
-template SimdBatch<float>
-barrier(const SimdBatch<float> d, const SimdBatch<float> dhat);
-template SimdBatch<double>
-barrier(const SimdBatch<double> d, const SimdBatch<double> dhat);
-template SimdBatch<float>
-barrier_first_derivative(const SimdBatch<float> d, const SimdBatch<float> dhat);
-template SimdBatch<double> barrier_first_derivative(
-    const SimdBatch<double> d, const SimdBatch<double> dhat);
-template SimdBatch<float> barrier_second_derivative(
-    const SimdBatch<float> d, const SimdBatch<float> dhat);
-template SimdBatch<double> barrier_second_derivative(
-    const SimdBatch<double> d, const SimdBatch<double> dhat);
 #endif
 /// @endcond
 // ============================================================================
 
+#endif // !__CUDACC__
+
 } // namespace ipc
diff --git a/src/ipc/barrier/barrier.hpp b/src/ipc/barrier/barrier.hpp
index 58c869eff..de7456e32 100644
--- a/src/ipc/barrier/barrier.hpp
+++ b/src/ipc/barrier/barrier.hpp
@@ -58,7 +58,8 @@ using Barrier = BarrierBase<>;
 /// @param d The distance.
 /// @param dhat Activation distance of the barrier.
 /// @return The value of the barrier function at d.
-template <typename T = double> T barrier(const T d, const T dhat);
+template <typename T = double>
+IPC_TOOLKIT_HOST_DEVICE T barrier(const T d, const T dhat);
 
 /// @brief Derivative of the barrier function.
 ///
@@ -71,7 +72,7 @@ template <typename T = double> T barrier(const T d, const T dhat);
 /// @param dhat Activation distance of the barrier.
 /// @return The derivative of the barrier wrt d.
 template <typename T = double>
-T barrier_first_derivative(const T d, const T dhat);
+IPC_TOOLKIT_HOST_DEVICE T barrier_first_derivative(const T d, const T dhat);
 
 /// @brief Second derivative of the barrier function.
 ///
@@ -84,7 +85,7 @@ T barrier_first_derivative(const T d, const T dhat);
 /// @param dhat Activation distance of the barrier.
 /// @return The second derivative of the barrier wrt d.
 template <typename T = double>
-T barrier_second_derivative(const T d, const T dhat);
+IPC_TOOLKIT_HOST_DEVICE T barrier_second_derivative(const T d, const T dhat);
 
 /// @brief Smoothly clamped log barrier functions from [Li et al. 2020].
 template <typename T = double> class ClampedLogBarrier : public BarrierBase<T> {
diff --git a/src/ipc/distance/CMakeLists.txt b/src/ipc/distance/CMakeLists.txt
index f327d575f..9ede9e1e7 100644
--- a/src/ipc/distance/CMakeLists.txt
+++ b/src/ipc/distance/CMakeLists.txt
@@ -1,31 +1,31 @@
 set(SOURCES
-  distance_type.cpp
   distance_type.hpp
-  edge_edge.cpp
   edge_edge.hpp
-  edge_edge_mollifier.cpp
   edge_edge_mollifier.hpp
-  line_line.cpp
   line_line.hpp
-  point_edge.cpp
   point_edge.hpp
   point_line.hpp
-  point_plane.cpp
   point_plane.hpp
   point_point.hpp
-  point_triangle.cpp
   point_triangle.hpp
 )
 
 target_sources(ipc_toolkit PRIVATE ${SOURCES})
 
 # These definitions are shared between host C++ and CUDA device code. The
-# header only declares them (IPC_TOOLKIT_HOST_DEVICE), so editing a definition
+# headers only declare them (IPC_TOOLKIT_HOST_DEVICE), so editing a definition
 # rebuilds one TU plus a link step rather than every TU that includes the
 # header. Under CUDA the helper generates a .cu wrapper per file so nvcc emits
 # relocatable device code linkable from other TUs.
 ipc_toolkit_target_shared_device_sources(ipc_toolkit
+  distance_type.cpp
+  edge_edge.cpp
+  edge_edge_mollifier.cpp
+  line_line.cpp
+  point_edge.cpp
   point_line.cpp
+  point_plane.cpp
+  point_triangle.cpp
 )
 
 ################################################################################
diff --git a/src/ipc/distance/distance_type.cpp b/src/ipc/distance/distance_type.cpp
index f9b281e42..f5bf59d84 100644
--- a/src/ipc/distance/distance_type.cpp
+++ b/src/ipc/distance/distance_type.cpp
@@ -1,40 +1,28 @@
 #include "distance_type.hpp"
 
 #include <ipc/utils/eigen_ext.hpp>
-#include <ipc/utils/logger.hpp>
 
 #include <Eigen/Geometry>
-#include <spdlog/spdlog.h>
 
 #include <limits>
-#include <stdexcept>
 
-namespace ipc::detail {
+// The error reporting helpers at the bottom of this file are host only: nvcc
+// cannot parse spdlog, so their definitions and the headers they need both sit
+// behind `#ifndef __CUDACC__` and are compiled by the host pass alone. Keep
+// every use of the logger and of fmt inside those blocks -- an unguarded one
+// fails as a parse error deep inside fmt rather than as a missing declaration.
+#ifndef __CUDACC__
+#include <ipc/utils/logger.hpp>
 
-void warn_degenerate_point_edge() noexcept
-{
-    logger().warn("Degenerate edge in point_edge_distance_type!");
-}
+#include <spdlog/spdlog.h>
 
-void throw_invalid_distance_type(const char* function)
-{
-    throw std::invalid_argument(
-        fmt::format("{}: invalid distance type", function));
-}
+#include <stdexcept>
+#endif
 
-void throw_auto_requires_explicit_dtype(const char* function)
-{
-    throw std::invalid_argument(
-        fmt::format(
-            "{}: an explicit distance type is required for non-floating-point "
-            "scalars; resolving AUTO means comparing single ordered values, "
-            "which an autodiff, SIMD batch, or interval scalar does not "
-            "provide",
-            function));
-}
+namespace ipc::detail {
 
 template <typename T>
-PointTriangleDistanceType point_triangle_distance_type(
+IPC_TOOLKIT_HOST_DEVICE PointTriangleDistanceType point_triangle_distance_type(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -95,7 +83,7 @@ PointTriangleDistanceType point_triangle_distance_type(
 
 // A more robust implementation of http://geomalgorithms.com/a07-_distance.html
 template <typename T>
-EdgeEdgeDistanceType edge_edge_distance_type(
+IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_distance_type(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -199,7 +187,7 @@ EdgeEdgeDistanceType edge_edge_distance_type(
 }
 
 template <typename T>
-EdgeEdgeDistanceType edge_edge_parallel_distance_type(
+IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_parallel_distance_type(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -237,12 +225,46 @@ EdgeEdgeDistanceType edge_edge_parallel_distance_type(
 }
 
 // clang-format off
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 template PointTriangleDistanceType point_triangle_distance_type<float>(Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>);
 template PointTriangleDistanceType point_triangle_distance_type<double>(Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>);
 template EdgeEdgeDistanceType edge_edge_distance_type<float>(Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>);
 template EdgeEdgeDistanceType edge_edge_distance_type<double>(Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>);
 template EdgeEdgeDistanceType edge_edge_parallel_distance_type<float>(Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>, Eigen::ConstRef<Eigen::Vector3f>);
 template EdgeEdgeDistanceType edge_edge_parallel_distance_type<double>(Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>, Eigen::ConstRef<Eigen::Vector3d>);
+#endif
 // clang-format on
 
+// ============================================================================
+// Error reporting -- host only.
+// ============================================================================
+//
+// The inline wrappers in distance_type.hpp call these on the host and trap on
+// the device, so the device pass needs neither the definitions nor spdlog.
+#ifndef __CUDACC__
+
+void warn_degenerate_point_edge_host() noexcept
+{
+    logger().warn("Degenerate edge in point_edge_distance_type!");
+}
+
+void throw_invalid_distance_type_host(const char* function)
+{
+    throw std::invalid_argument(
+        fmt::format("{}: invalid distance type", function));
+}
+
+void throw_auto_requires_explicit_dtype_host(const char* function)
+{
+    throw std::invalid_argument(
+        fmt::format(
+            "{}: an explicit distance type is required for non-floating-point "
+            "scalars; resolving AUTO means comparing single ordered values, "
+            "which an autodiff, SIMD batch, or interval scalar does not "
+            "provide",
+            function));
+}
+
+#endif // !__CUDACC__
+
 } // namespace ipc::detail
diff --git a/src/ipc/distance/distance_type.hpp b/src/ipc/distance/distance_type.hpp
index 786d38ec2..5e4b23440 100644
--- a/src/ipc/distance/distance_type.hpp
+++ b/src/ipc/distance/distance_type.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 #include <cassert>
@@ -60,14 +61,60 @@ enum class EdgeEdgeDistanceType : uint8_t {
 namespace detail {
     /// @brief Warn about a degenerate edge.
     /// @note Out of line to keep the logger (and spdlog) out of this header.
-    void warn_degenerate_point_edge() noexcept;
+    void warn_degenerate_point_edge_host() noexcept;
 
     /// @brief Throw for an invalid distance type.
     /// @note Out of line and [[noreturn]] so that constructing the exception does not consume the caller's inlining budget on the hot path.
-    [[noreturn]] void throw_invalid_distance_type(const char* function);
+    [[noreturn]] void throw_invalid_distance_type_host(const char* function);
 
     /// @brief Throw when AUTO is requested for a scalar that cannot resolve it.
-    [[noreturn]] void throw_auto_requires_explicit_dtype(const char* function);
+    [[noreturn]] void
+    throw_auto_requires_explicit_dtype_host(const char* function);
+
+    // The distance functions below are shared between host C++ and CUDA device
+    // code, but the device has neither the logger nor exceptions. So we funnel
+    // every error report through the wrappers below: on the host they forward
+    // to the out-of-line helpers above, and on the device they assert and trap.
+    // The tradeoff is that a bad distance type takes down the kernel instead of
+    // unwinding to a handler, but that is the best we can do without
+    // exceptions, and it still fails loudly instead of returning garbage.
+
+    /// @brief Warn about a degenerate edge on the host; a no-op on the device.
+    IPC_TOOLKIT_HOST_DEVICE inline void warn_degenerate_point_edge() noexcept
+    {
+#ifdef __CUDA_ARCH__
+        // A degenerate edge is recoverable (we fall back to an arbitrary
+        // end-point), so on the device we stay silent rather than trapping.
+#else
+        warn_degenerate_point_edge_host();
+#endif
+    }
+
+    /// @brief Report an invalid distance type: throws on the host, traps on the device.
+    [[noreturn]] IPC_TOOLKIT_HOST_DEVICE inline void
+    throw_invalid_distance_type(const char* function)
+    {
+#ifdef __CUDA_ARCH__
+        (void)function;
+        assert(false && "Invalid distance type!");
+        __trap();
+#else
+        throw_invalid_distance_type_host(function);
+#endif
+    }
+
+    /// @brief Report an unresolvable AUTO: throws on the host, traps on the device.
+    [[noreturn]] IPC_TOOLKIT_HOST_DEVICE inline void
+    throw_auto_requires_explicit_dtype(const char* function)
+    {
+#ifdef __CUDA_ARCH__
+        (void)function;
+        assert(false && "An explicit distance type is required!");
+        __trap();
+#else
+        throw_auto_requires_explicit_dtype_host(function);
+#endif
+    }
 
     /// @brief Determine the closest pair between a point and edge.
     /// @note Prefer the ipc::point_edge_distance_type front end below, which deduces both the scalar type and the dimension.
@@ -78,7 +125,8 @@ namespace detail {
     /// @param e1 The second vertex of the edge.
     /// @return The distance type of the point-edge pair.
     template <typename T, int dim>
-    inline PointEdgeDistanceType point_edge_distance_type(
+    IPC_TOOLKIT_HOST_DEVICE inline PointEdgeDistanceType
+    point_edge_distance_type(
         const Eigen::Vector<T, dim>& p,
         const Eigen::Vector<T, dim>& e0,
         const Eigen::Vector<T, dim>& e1)
@@ -110,7 +158,8 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The distance type of the point-triangle pair.
     template <typename T>
-    PointTriangleDistanceType point_triangle_distance_type(
+    IPC_TOOLKIT_HOST_DEVICE PointTriangleDistanceType
+    point_triangle_distance_type(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -123,7 +172,7 @@ namespace detail {
     /// @param eb1 The second vertex of the second edge.
     /// @return The distance type of the edge-edge pair.
     template <typename T>
-    EdgeEdgeDistanceType edge_edge_distance_type(
+    IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_distance_type(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -136,7 +185,8 @@ namespace detail {
     /// @param eb1 The second vertex of the second edge.
     /// @return The distance type of the edge-edge pair.
     template <typename T>
-    EdgeEdgeDistanceType edge_edge_parallel_distance_type(
+    IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType
+    edge_edge_parallel_distance_type(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -159,7 +209,7 @@ namespace detail {
 /// @param e1 The second vertex of the edge.
 /// @return The distance type of the point-edge pair.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline PointEdgeDistanceType point_edge_distance_type(
+IPC_TOOLKIT_HOST_DEVICE inline PointEdgeDistanceType point_edge_distance_type(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -190,7 +240,8 @@ template <
     typename DerivedT1,
     typename DerivedT2,
     std::enable_if_t<std::is_class_v<DerivedP>, int> = 0>
-inline PointTriangleDistanceType point_triangle_distance_type(
+IPC_TOOLKIT_HOST_DEVICE inline PointTriangleDistanceType
+point_triangle_distance_type(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -212,7 +263,7 @@ template <
     typename DerivedEB0,
     typename DerivedEB1,
     std::enable_if_t<std::is_class_v<DerivedEA0>, int> = 0>
-inline EdgeEdgeDistanceType edge_edge_distance_type(
+IPC_TOOLKIT_HOST_DEVICE inline EdgeEdgeDistanceType edge_edge_distance_type(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -234,7 +285,8 @@ template <
     typename DerivedEB0,
     typename DerivedEB1,
     std::enable_if_t<std::is_class_v<DerivedEA0>, int> = 0>
-inline EdgeEdgeDistanceType edge_edge_parallel_distance_type(
+IPC_TOOLKIT_HOST_DEVICE inline EdgeEdgeDistanceType
+edge_edge_parallel_distance_type(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
diff --git a/src/ipc/distance/edge_edge.cpp b/src/ipc/distance/edge_edge.cpp
index 1d308ca76..71b3487f3 100644
--- a/src/ipc/distance/edge_edge.cpp
+++ b/src/ipc/distance/edge_edge.cpp
@@ -10,7 +10,7 @@
 namespace ipc::detail {
 
 template <typename T>
-T edge_edge_distance(
+IPC_TOOLKIT_HOST_DEVICE T edge_edge_distance(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -59,7 +59,7 @@ T edge_edge_distance(
 }
 
 template <typename T>
-Eigen::Vector<T, 12> edge_edge_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12> edge_edge_distance_gradient(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -143,7 +143,7 @@ Eigen::Vector<T, 12> edge_edge_distance_gradient(
 }
 
 template <typename T>
-Eigen::Matrix<T, 12, 12> edge_edge_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12> edge_edge_distance_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -290,14 +290,18 @@ Eigen::Matrix<T, 12, 12> edge_edge_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>, EdgeEdgeDistanceType)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_EDGE_EDGE(float);
 IPC_INSTANTIATE_EDGE_EDGE(double);
+#endif
+#if IPC_TOOLKIT_INSTANTIATE_HOST_SCALARS
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADGrad<9>);
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADHessian<9>);
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADGrad<12>);
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADHessian<12>);
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADGrad<13>);
 IPC_INSTANTIATE_EDGE_EDGE_VALUE(ADHessian<13>);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_EDGE_EDGE(SimdBatch<float>);
 IPC_INSTANTIATE_EDGE_EDGE(SimdBatch<double>);
diff --git a/src/ipc/distance/edge_edge.hpp b/src/ipc/distance/edge_edge.hpp
index d07e3b464..df715ab5d 100644
--- a/src/ipc/distance/edge_edge.hpp
+++ b/src/ipc/distance/edge_edge.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/distance/distance_type.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -15,7 +16,7 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The distance between the two edges.
     template <typename T>
-    T edge_edge_distance(
+    IPC_TOOLKIT_HOST_DEVICE T edge_edge_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -31,7 +32,7 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The gradient of the distance wrt ea0, ea1, eb0, and eb1.
     template <typename T>
-    Eigen::Vector<T, 12> edge_edge_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12> edge_edge_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -47,7 +48,7 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The hessian of the distance wrt ea0, ea1, eb0, and eb1.
     template <typename T>
-    Eigen::Matrix<T, 12, 12> edge_edge_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12> edge_edge_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -70,7 +71,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_distance(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -94,7 +95,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_distance_gradient(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -118,7 +119,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_distance_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
diff --git a/src/ipc/distance/edge_edge_mollifier.cpp b/src/ipc/distance/edge_edge_mollifier.cpp
index 44bf82e45..d4229636b 100644
--- a/src/ipc/distance/edge_edge_mollifier.cpp
+++ b/src/ipc/distance/edge_edge_mollifier.cpp
@@ -4,7 +4,8 @@ namespace ipc {
 namespace detail {
 
     template <typename T>
-    Eigen::Vector<T, 12> edge_edge_mollifier_gradient_wrt_x(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12>
+    edge_edge_mollifier_gradient_wrt_x(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -36,7 +37,8 @@ namespace detail {
     }
 
     template <typename T>
-    Eigen::Matrix<T, 12, 12> edge_edge_mollifier_gradient_jacobian_wrt_x(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+    edge_edge_mollifier_gradient_jacobian_wrt_x(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -90,8 +92,10 @@ namespace detail {
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER(float);
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER(SimdBatch<float>);
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER(SimdBatch<double>);
@@ -104,7 +108,7 @@ namespace autogen {
     // This function was generated by the Symbolic Math Toolbox version 8.3.
     // 01-Nov-2019 16:54:23
     template <typename T>
-    void edge_edge_cross_squarednorm_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_cross_squarednorm_gradient(
         T v01,
         T v02,
         T v03,
@@ -160,7 +164,7 @@ namespace autogen {
     // This function was generated by the Symbolic Math Toolbox version 8.3.
     // 01-Nov-2019 16:54:23
     template <typename T>
-    void edge_edge_cross_squarednorm_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_cross_squarednorm_hessian(
         T v01,
         T v02,
         T v03,
@@ -396,7 +400,7 @@ namespace autogen {
     }
 
     template <typename T>
-    void edge_edge_mollifier_threshold_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_mollifier_threshold_gradient(
         T ea0x,
         T ea0y,
         T ea0z,
@@ -449,8 +453,10 @@ namespace autogen {
     template void edge_edge_mollifier_threshold_gradient<T>(                   \
         T, T, T, T, T, T, T, T, T, T, T, T, T[12], T)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER_AUTOGEN(float);
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER_AUTOGEN(SimdBatch<float>);
     IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER_AUTOGEN(SimdBatch<double>);
diff --git a/src/ipc/distance/edge_edge_mollifier.hpp b/src/ipc/distance/edge_edge_mollifier.hpp
index 670a7b18d..5e3b8ea0b 100644
--- a/src/ipc/distance/edge_edge_mollifier.hpp
+++ b/src/ipc/distance/edge_edge_mollifier.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/math/scalar_math.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 #include <ipc/utils/simd.hpp>
@@ -10,13 +11,13 @@ namespace ipc {
 namespace autogen {
     // clang-format off
     template <typename T>
-    void edge_edge_cross_squarednorm_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_cross_squarednorm_gradient(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T g[12]);
     template <typename T>
-    void edge_edge_cross_squarednorm_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_cross_squarednorm_hessian(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T H[144]);
     template <typename T>
-    void edge_edge_mollifier_threshold_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_mollifier_threshold_gradient(
         T ea0x, T ea0y, T ea0z, T ea1x, T ea1y, T ea1z, T eb0x, T eb0y, T eb0z, T eb1x, T eb1y, T eb1z, T grad[12], T scale = literal<T>(1e-3));
     // clang-format on
 } // namespace autogen
@@ -27,7 +28,8 @@ namespace autogen {
 /// @param x Squared norm of the edge-edge cross product.
 /// @param eps_x Mollifier activation threshold.
 /// @return The mollifier coefficient to premultiply the edge-edge distance.
-template <typename T> inline T edge_edge_mollifier(const T x, const T eps_x)
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE inline T edge_edge_mollifier(const T x, const T eps_x)
 {
     return select_lazy(
         x < eps_x,
@@ -43,7 +45,8 @@ template <typename T> inline T edge_edge_mollifier(const T x, const T eps_x)
 /// @param eps_x Mollifier activation threshold.
 /// @return The gradient of the mollifier function for edge-edge distance wrt x.
 template <typename T>
-inline T edge_edge_mollifier_gradient(const T x, const T eps_x)
+IPC_TOOLKIT_HOST_DEVICE inline T
+edge_edge_mollifier_gradient(const T x, const T eps_x)
 {
     using namespace ipc::numext; // fma
     return select_lazy(
@@ -62,7 +65,8 @@ inline T edge_edge_mollifier_gradient(const T x, const T eps_x)
 /// @return The derivative of the mollifier function for edge-edge distance wrt
 ///     eps_x.
 template <typename T>
-inline T edge_edge_mollifier_derivative_wrt_eps_x(const T x, const T eps_x)
+IPC_TOOLKIT_HOST_DEVICE inline T
+edge_edge_mollifier_derivative_wrt_eps_x(const T x, const T eps_x)
 {
     return select_lazy(
         x < eps_x,
@@ -75,7 +79,8 @@ inline T edge_edge_mollifier_derivative_wrt_eps_x(const T x, const T eps_x)
 /// @param eps_x Mollifier activation threshold.
 /// @return The hessian of the mollifier function for edge-edge distance wrt x.
 template <typename T>
-inline T edge_edge_mollifier_hessian(const T x, const T eps_x)
+IPC_TOOLKIT_HOST_DEVICE inline T
+edge_edge_mollifier_hessian(const T x, const T eps_x)
 {
     return select_lazy(
         x < eps_x, [&] { return T(-2) / (eps_x * eps_x); },
@@ -89,7 +94,7 @@ inline T edge_edge_mollifier_hessian(const T x, const T eps_x)
 /// @return The derivative of the gradient of the mollifier function for
 ///     edge-edge distance wrt eps_x.
 template <typename T>
-inline T
+IPC_TOOLKIT_HOST_DEVICE inline T
 edge_edge_mollifier_gradient_derivative_wrt_eps_x(const T x, const T eps_x)
 {
     return select_lazy(
@@ -103,7 +108,7 @@ edge_edge_mollifier_gradient_derivative_wrt_eps_x(const T x, const T eps_x)
 namespace detail {
     /// @note Prefer the ipc::edge_edge_cross_squarednorm front end.
     template <typename T>
-    inline T edge_edge_cross_squarednorm(
+    IPC_TOOLKIT_HOST_DEVICE inline T edge_edge_cross_squarednorm(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -114,7 +119,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_cross_squarednorm_gradient front end.
     template <typename T>
-    inline Eigen::Vector<T, 12> edge_edge_cross_squarednorm_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 12>
+    edge_edge_cross_squarednorm_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -129,7 +135,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_cross_squarednorm_hessian front end.
     template <typename T>
-    inline Eigen::Matrix<T, 12, 12> edge_edge_cross_squarednorm_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 12, 12>
+    edge_edge_cross_squarednorm_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -144,7 +151,7 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier front end.
     template <typename T>
-    inline T edge_edge_mollifier(
+    IPC_TOOLKIT_HOST_DEVICE inline T edge_edge_mollifier(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -157,7 +164,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier_gradient front end.
     template <typename T>
-    inline Eigen::Vector<T, 12> edge_edge_mollifier_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 12>
+    edge_edge_mollifier_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -178,7 +186,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier_hessian front end.
     template <typename T>
-    inline Eigen::Matrix<T, 12, 12> edge_edge_mollifier_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 12, 12>
+    edge_edge_mollifier_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -205,7 +214,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier_gradient_wrt_x front end.
     template <typename T>
-    Eigen::Vector<T, 12> edge_edge_mollifier_gradient_wrt_x(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12>
+    edge_edge_mollifier_gradient_wrt_x(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -218,7 +228,8 @@ namespace detail {
     /// @note Prefer the ipc::edge_edge_mollifier_gradient_jacobian_wrt_x front
     ///     end.
     template <typename T>
-    Eigen::Matrix<T, 12, 12> edge_edge_mollifier_gradient_jacobian_wrt_x(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+    edge_edge_mollifier_gradient_jacobian_wrt_x(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -230,7 +241,7 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier_threshold front end.
     template <typename T>
-    T edge_edge_mollifier_threshold(
+    IPC_TOOLKIT_HOST_DEVICE T edge_edge_mollifier_threshold(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -242,7 +253,8 @@ namespace detail {
 
     /// @note Prefer the ipc::edge_edge_mollifier_threshold_gradient front end.
     template <typename T>
-    Eigen::Vector<T, 12> edge_edge_mollifier_threshold_gradient(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12>
+    edge_edge_mollifier_threshold_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1_rest,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0_rest,
@@ -270,7 +282,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_cross_squarednorm(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_cross_squarednorm(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -294,7 +306,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_cross_squarednorm_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_cross_squarednorm_gradient(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -316,7 +328,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_cross_squarednorm_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_cross_squarednorm_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -341,7 +353,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_mollifier(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -364,7 +376,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_mollifier_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_gradient(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -387,7 +399,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_mollifier_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -418,7 +430,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_mollifier_gradient_wrt_x(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_gradient_wrt_x(
     const Eigen::MatrixBase<DerivedEA0Rest>& ea0_rest,
     const Eigen::MatrixBase<DerivedEA1Rest>& ea1_rest,
     const Eigen::MatrixBase<DerivedEB0Rest>& eb0_rest,
@@ -455,7 +467,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_mollifier_gradient_jacobian_wrt_x(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_gradient_jacobian_wrt_x(
     const Eigen::MatrixBase<DerivedEA0Rest>& ea0_rest,
     const Eigen::MatrixBase<DerivedEA1Rest>& ea1_rest,
     const Eigen::MatrixBase<DerivedEB0Rest>& eb0_rest,
@@ -484,7 +496,7 @@ template <
     typename DerivedEA1Rest,
     typename DerivedEB0Rest,
     typename DerivedEB1Rest>
-inline auto edge_edge_mollifier_threshold(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_threshold(
     const Eigen::MatrixBase<DerivedEA0Rest>& ea0_rest,
     const Eigen::MatrixBase<DerivedEA1Rest>& ea1_rest,
     const Eigen::MatrixBase<DerivedEB0Rest>& eb0_rest,
@@ -510,7 +522,7 @@ template <
     typename DerivedEA1Rest,
     typename DerivedEB0Rest,
     typename DerivedEB1Rest>
-inline auto edge_edge_mollifier_threshold_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_mollifier_threshold_gradient(
     const Eigen::MatrixBase<DerivedEA0Rest>& ea0_rest,
     const Eigen::MatrixBase<DerivedEA1Rest>& ea1_rest,
     const Eigen::MatrixBase<DerivedEB0Rest>& eb0_rest,
diff --git a/src/ipc/distance/line_line.cpp b/src/ipc/distance/line_line.cpp
index 12d7301a3..ca506ec73 100644
--- a/src/ipc/distance/line_line.cpp
+++ b/src/ipc/distance/line_line.cpp
@@ -8,7 +8,7 @@ namespace ipc::autogen {
 // This function was generated by the Symbolic Math Toolbox version 8.3.
 // 14-Jun-2019 13:58:25
 template <typename T>
-void line_line_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE void line_line_distance_gradient(
     T v01,
     T v02,
     T v03,
@@ -75,7 +75,7 @@ void line_line_distance_gradient(
 // This function was generated by the Symbolic Math Toolbox version 8.3.
 // 14-Jun-2019 13:58:38
 template <typename T>
-void line_line_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE void line_line_distance_hessian(
     T v01,
     T v02,
     T v03,
@@ -565,8 +565,10 @@ void line_line_distance_hessian(
     template void line_line_distance_hessian<T>(                               \
         T, T, T, T, T, T, T, T, T, T, T, T, T[144])
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_LINE_LINE_AUTOGEN(float);
 IPC_INSTANTIATE_LINE_LINE_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 // SIMD batches, so that the batched distance gradients/Hessians link.
 IPC_INSTANTIATE_LINE_LINE_AUTOGEN(SimdBatch<float>);
diff --git a/src/ipc/distance/line_line.hpp b/src/ipc/distance/line_line.hpp
index cb0b35f34..d51b7017c 100644
--- a/src/ipc/distance/line_line.hpp
+++ b/src/ipc/distance/line_line.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 #include <Eigen/Geometry>
@@ -10,10 +11,10 @@ namespace ipc {
 namespace autogen {
     // clang-format off
     template <typename T>
-    void line_line_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void line_line_distance_gradient(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T g[12]);
     template <typename T>
-    void line_line_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void line_line_distance_hessian(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T H[144]);
     // clang-format on
 } // namespace autogen
@@ -28,7 +29,7 @@ namespace detail {
     /// @param eb1 The second vertex of the edge defining the second line.
     /// @return The distance between the two lines.
     template <typename T>
-    inline T line_line_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T line_line_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -53,7 +54,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_distance(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -70,7 +71,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_distance_gradient(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -88,7 +89,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_distance_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
diff --git a/src/ipc/distance/point_edge.cpp b/src/ipc/distance/point_edge.cpp
index 31fefdbcf..e93ba6215 100644
--- a/src/ipc/distance/point_edge.cpp
+++ b/src/ipc/distance/point_edge.cpp
@@ -12,7 +12,8 @@ namespace ipc {
 namespace detail {
 
     template <typename T, int dim>
-    Eigen::Matrix<T, 3 * dim, 3 * dim> point_edge_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 3 * dim, 3 * dim>
+    point_edge_distance_hessian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1,
@@ -68,10 +69,12 @@ namespace detail {
         Eigen::ConstRef<Eigen::Vector<T, dim>>,                                \
         Eigen::ConstRef<Eigen::Vector<T, dim>>, PointEdgeDistanceType)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(float, 2);
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(float, 3);
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(double, 2);
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(double, 3);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(SimdBatch<float>, 2);
     IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(SimdBatch<float>, 3);
diff --git a/src/ipc/distance/point_edge.hpp b/src/ipc/distance/point_edge.hpp
index 80875cec6..584795d5a 100644
--- a/src/ipc/distance/point_edge.hpp
+++ b/src/ipc/distance/point_edge.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/distance/distance_type.hpp>
 #include <ipc/distance/point_line.hpp>
 #include <ipc/distance/point_point.hpp>
@@ -21,7 +22,7 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The distance between the point and edge.
     template <typename T, int dim>
-    inline T point_edge_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_edge_distance(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1,
@@ -62,7 +63,8 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The gradient of the distance wrt p, e0, and e1.
     template <typename T, int dim>
-    inline Eigen::Vector<T, 3 * dim> point_edge_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 3 * dim>
+    point_edge_distance_gradient(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1,
@@ -114,7 +116,8 @@ namespace detail {
     /// @param dtype The point edge distance type to compute.
     /// @return The hessian of the distance wrt p, e0, and e1.
     template <typename T, int dim>
-    Eigen::Matrix<T, 3 * dim, 3 * dim> point_edge_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 3 * dim, 3 * dim>
+    point_edge_distance_hessian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1,
@@ -129,7 +132,7 @@ namespace detail {
 /// @param dtype The point edge distance type to compute.
 /// @return The distance between the point and edge.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1,
@@ -157,7 +160,7 @@ inline auto point_edge_distance(
 /// @param dtype The point edge distance type to compute.
 /// @return The gradient of the distance wrt p, e0, and e1.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1,
@@ -186,7 +189,7 @@ inline auto point_edge_distance_gradient(
 /// @param dtype The point edge distance type to compute.
 /// @return The hessian of the distance wrt p, e0, and e1.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1,
diff --git a/src/ipc/distance/point_plane.cpp b/src/ipc/distance/point_plane.cpp
index 70374eaa6..6b168af00 100644
--- a/src/ipc/distance/point_plane.cpp
+++ b/src/ipc/distance/point_plane.cpp
@@ -7,7 +7,7 @@ namespace ipc::autogen {
 // This function was generated by the Symbolic Math Toolbox version 8.3.
 // 10-Jun-2019 17:42:16
 template <typename T>
-void point_plane_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE void point_plane_distance_gradient(
     T v01,
     T v02,
     T v03,
@@ -72,7 +72,7 @@ void point_plane_distance_gradient(
 // This function was generated by the Symbolic Math Toolbox version 8.3.
 // 10-Jun-2019 17:42:25
 template <typename T>
-void point_plane_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE void point_plane_distance_hessian(
     T v01,
     T v02,
     T v03,
@@ -566,8 +566,10 @@ void point_plane_distance_hessian(
     template void point_plane_distance_hessian<T>(                             \
         T, T, T, T, T, T, T, T, T, T, T, T, T[144])
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_POINT_PLANE_AUTOGEN(float);
 IPC_INSTANTIATE_POINT_PLANE_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_POINT_PLANE_AUTOGEN(SimdBatch<float>);
 IPC_INSTANTIATE_POINT_PLANE_AUTOGEN(SimdBatch<double>);
diff --git a/src/ipc/distance/point_plane.hpp b/src/ipc/distance/point_plane.hpp
index 9a2416acb..d200e8547 100644
--- a/src/ipc/distance/point_plane.hpp
+++ b/src/ipc/distance/point_plane.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/geometry/normal.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -9,10 +10,10 @@ namespace ipc {
 namespace autogen {
     // clang-format off
     template <typename T>
-    void point_plane_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void point_plane_distance_gradient(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T g[12]);
     template <typename T>
-    void point_plane_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void point_plane_distance_hessian(
         T v01, T v02, T v03, T v11, T v12, T v13, T v21, T v22, T v23, T v31, T v32, T v33, T H[144]);
     // clang-format on
 } // namespace autogen
@@ -25,7 +26,7 @@ namespace detail {
     /// @param normal The normal of the plane.
     /// @return The distance between the point and plane.
     template <typename T>
-    inline T point_plane_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_plane_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> origin,
         Eigen::ConstRef<Eigen::Vector3<T>> normal)
@@ -42,7 +43,7 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The distance between the point and plane.
     template <typename T>
-    inline T point_plane_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_plane_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -61,7 +62,8 @@ namespace detail {
     /// @param normal The normal of the plane.
     /// @return The gradient of the distance wrt p.
     template <typename T>
-    inline Eigen::Vector3<T> point_plane_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T>
+    point_plane_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> origin,
         Eigen::ConstRef<Eigen::Vector3<T>> normal)
@@ -78,7 +80,8 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The gradient of the distance wrt p, t0, t1, and t2.
     template <typename T>
-    inline Eigen::Vector<T, 12> point_plane_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 12>
+    point_plane_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -98,7 +101,8 @@ namespace detail {
     /// @param normal The normal of the plane.
     /// @return The hessian of the distance wrt p.
     template <typename T>
-    inline Eigen::Matrix3<T> point_plane_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix3<T>
+    point_plane_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> origin,
         Eigen::ConstRef<Eigen::Vector3<T>> normal)
@@ -114,7 +118,8 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The hessian of the distance wrt p, t0, t1, and t2.
     template <typename T>
-    inline Eigen::Matrix<T, 12, 12> point_plane_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 12, 12>
+    point_plane_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -135,7 +140,7 @@ namespace detail {
 /// @param normal The normal of the plane.
 /// @return The distance between the point and plane.
 template <typename DerivedP, typename DerivedOrigin, typename DerivedNormal>
-inline auto point_plane_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedOrigin>& origin,
     const Eigen::MatrixBase<DerivedNormal>& normal)
@@ -158,7 +163,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -175,7 +180,7 @@ inline auto point_plane_distance(
 /// @param normal The normal of the plane.
 /// @return The gradient of the distance wrt p.
 template <typename DerivedP, typename DerivedOrigin, typename DerivedNormal>
-inline auto point_plane_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedOrigin>& origin,
     const Eigen::MatrixBase<DerivedNormal>& normal)
@@ -196,7 +201,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -213,7 +218,7 @@ inline auto point_plane_distance_gradient(
 /// @param normal The normal of the plane.
 /// @return The hessian of the distance wrt p.
 template <typename DerivedP, typename DerivedOrigin, typename DerivedNormal>
-inline auto point_plane_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedOrigin>& origin,
     const Eigen::MatrixBase<DerivedNormal>& normal)
@@ -234,7 +239,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
diff --git a/src/ipc/distance/point_point.hpp b/src/ipc/distance/point_point.hpp
index af6676327..b4e2006fc 100644
--- a/src/ipc/distance/point_point.hpp
+++ b/src/ipc/distance/point_point.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 #include <cassert>
@@ -15,7 +16,7 @@ namespace detail {
     /// @param p1 The second point.
     /// @return The distance between p0 and p1.
     template <typename T, int dim>
-    inline T point_point_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_point_distance(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> p1)
     {
@@ -31,7 +32,8 @@ namespace detail {
     /// @param p1 The second point.
     /// @return The computed gradient.
     template <typename T, int dim>
-    inline Eigen::Vector<T, 2 * dim> point_point_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 2 * dim>
+    point_point_distance_gradient(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> p1)
     {
@@ -50,7 +52,8 @@ namespace detail {
     /// @param p1 The second point.
     /// @return The computed hessian.
     template <typename T, int dim>
-    inline Eigen::Matrix<T, 2 * dim, 2 * dim> point_point_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 2 * dim, 2 * dim>
+    point_point_distance_hessian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> /*p0*/,
         Eigen::ConstRef<Eigen::Vector<T, dim>> /*p1*/)
     {
@@ -73,7 +76,7 @@ namespace detail {
 /// @param p1 The second point.
 /// @return The distance between p0 and p1.
 template <typename DerivedP0, typename DerivedP1>
-inline auto point_point_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance(
     const Eigen::MatrixBase<DerivedP0>& p0,
     const Eigen::MatrixBase<DerivedP1>& p1)
 {
@@ -98,7 +101,7 @@ inline auto point_point_distance(
 /// @param p1 The second point.
 /// @return The computed gradient.
 template <typename DerivedP0, typename DerivedP1>
-inline auto point_point_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance_gradient(
     const Eigen::MatrixBase<DerivedP0>& p0,
     const Eigen::MatrixBase<DerivedP1>& p1)
 {
@@ -124,7 +127,7 @@ inline auto point_point_distance_gradient(
 /// @param p1 The second point.
 /// @return The computed hessian.
 template <typename DerivedP0, typename DerivedP1>
-inline auto point_point_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance_hessian(
     const Eigen::MatrixBase<DerivedP0>& p0,
     const Eigen::MatrixBase<DerivedP1>& p1)
 {
diff --git a/src/ipc/distance/point_triangle.cpp b/src/ipc/distance/point_triangle.cpp
index 484d8ca2f..7803fd5a7 100644
--- a/src/ipc/distance/point_triangle.cpp
+++ b/src/ipc/distance/point_triangle.cpp
@@ -10,7 +10,7 @@
 namespace ipc::detail {
 
 template <typename T>
-T point_triangle_distance(
+IPC_TOOLKIT_HOST_DEVICE T point_triangle_distance(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -53,7 +53,7 @@ T point_triangle_distance(
 }
 
 template <typename T>
-Eigen::Vector<T, 12> point_triangle_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12> point_triangle_distance_gradient(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -124,7 +124,8 @@ Eigen::Vector<T, 12> point_triangle_distance_gradient(
 }
 
 template <typename T>
-Eigen::Matrix<T, 12, 12> point_triangle_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+point_triangle_distance_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -253,12 +254,16 @@ Eigen::Matrix<T, 12, 12> point_triangle_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>, PointTriangleDistanceType)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_POINT_TRIANGLE(float);
 IPC_INSTANTIATE_POINT_TRIANGLE(double);
+#endif
+#if IPC_TOOLKIT_INSTANTIATE_HOST_SCALARS
 IPC_INSTANTIATE_POINT_TRIANGLE_VALUE(ADGrad<12>);
 IPC_INSTANTIATE_POINT_TRIANGLE_VALUE(ADHessian<12>);
 IPC_INSTANTIATE_POINT_TRIANGLE_VALUE(ADGrad<13>);
 IPC_INSTANTIATE_POINT_TRIANGLE_VALUE(ADHessian<13>);
+#endif
 
 #ifdef IPC_TOOLKIT_WITH_SIMD
 // SIMD batches. Only an explicit distance type is supported.
diff --git a/src/ipc/distance/point_triangle.hpp b/src/ipc/distance/point_triangle.hpp
index 1ae668e7e..e534232eb 100644
--- a/src/ipc/distance/point_triangle.hpp
+++ b/src/ipc/distance/point_triangle.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/distance/distance_type.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -15,7 +16,7 @@ namespace detail {
     /// @param dtype The point-triangle distance type to compute.
     /// @return The distance between the point and triangle.
     template <typename T>
-    T point_triangle_distance(
+    IPC_TOOLKIT_HOST_DEVICE T point_triangle_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -31,7 +32,8 @@ namespace detail {
     /// @param dtype The point-triangle distance type to compute.
     /// @return The gradient of the distance wrt p, t0, t1, and t2.
     template <typename T>
-    Eigen::Vector<T, 12> point_triangle_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12>
+    point_triangle_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -47,7 +49,8 @@ namespace detail {
     /// @param dtype The point-triangle distance type to compute.
     /// @return The hessian of the distance wrt p, t0, t1, and t2.
     template <typename T>
-    Eigen::Matrix<T, 12, 12> point_triangle_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+    point_triangle_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -70,7 +73,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -94,7 +97,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -118,7 +121,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
diff --git a/src/ipc/distance/signed/CMakeLists.txt b/src/ipc/distance/signed/CMakeLists.txt
index 55e7737ad..77c3325f9 100644
--- a/src/ipc/distance/signed/CMakeLists.txt
+++ b/src/ipc/distance/signed/CMakeLists.txt
@@ -1,10 +1,15 @@
 set(SOURCES
-  line_line.cpp
   line_line.hpp
-  point_line.cpp
   point_line.hpp
-  point_plane.cpp
   point_plane.hpp
 )
 
-target_sources(ipc_toolkit PRIVATE ${SOURCES})
\ No newline at end of file
+target_sources(ipc_toolkit PRIVATE ${SOURCES})
+
+# Definitions shared between host C++ and CUDA device code; under CUDA the
+# helper generates a .cu wrapper per file (see ipc_toolkit_shared_device_sources).
+ipc_toolkit_target_shared_device_sources(ipc_toolkit
+  line_line.cpp
+  point_line.cpp
+  point_plane.cpp
+)
\ No newline at end of file
diff --git a/src/ipc/distance/signed/line_line.cpp b/src/ipc/distance/signed/line_line.cpp
index c82e8a7aa..6ccdb3db1 100644
--- a/src/ipc/distance/signed/line_line.cpp
+++ b/src/ipc/distance/signed/line_line.cpp
@@ -5,7 +5,8 @@
 namespace ipc::detail {
 
 template <typename T>
-Eigen::Matrix<T, 12, 12> line_line_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+line_line_signed_distance_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -28,8 +29,18 @@ Eigen::Matrix<T, 12, 12> line_line_signed_distance_hessian(
     // Contract the normal Hessian (3x12x12) with vector v (3x1).
     // This computes (v ⋅ d²n/dx²).
     // The result is a 1x12x12 vector, which maps to the 12x12 Hessian matrix.
-    hess = (hess_n.reshaped(Eigen::fix<3>, Eigen::fix<144>).transpose() * v)
-               .reshaped(Eigen::fix<12>, Eigen::fix<12>);
+    // We spell the two reshapes out as Maps rather than calling
+    // .reshaped(Eigen::fix<...>). Both are fixed-size views over the same
+    // column-major storage, so the result is identical, but .reshaped()
+    // returns a nested expression template that nvcc does not handle, and
+    // this contraction has to stay device-callable.
+    {
+        const Eigen::Map<const Eigen::Matrix<T, 3, 144>> hess_n_3_144(
+            hess_n.data());
+        const Eigen::Matrix<T, 144, 1> contracted =
+            (hess_n_3_144.transpose() * v).eval();
+        hess = Eigen::Map<const Eigen::Matrix<T, 12, 12>>(contracted.data());
+    }
 
     // ---------------------------------------------------------
     // 2. Add Jacobian Terms (Product Rule Corrections)
@@ -80,8 +91,10 @@ Eigen::Matrix<T, 12, 12> line_line_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_LINE_LINE_SIGNED_DISTANCE_HESSIAN(float);
 IPC_INSTANTIATE_LINE_LINE_SIGNED_DISTANCE_HESSIAN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_LINE_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch<float>);
 IPC_INSTANTIATE_LINE_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch<double>);
diff --git a/src/ipc/distance/signed/line_line.hpp b/src/ipc/distance/signed/line_line.hpp
index 34fa41a86..27a43ec00 100644
--- a/src/ipc/distance/signed/line_line.hpp
+++ b/src/ipc/distance/signed/line_line.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/geometry/normal.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -26,7 +27,7 @@ namespace detail {
     /// represent segment endpoints). Behavior is undefined if ea0 == ea1 or eb0
     /// == eb1.
     template <typename T>
-    inline T line_line_signed_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T line_line_signed_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -55,7 +56,8 @@ namespace detail {
     ///
     /// @see line_line_signed_distance, line_line_normal
     template <typename T>
-    inline Eigen::Vector<T, 12> line_line_signed_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 12>
+    line_line_signed_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -90,7 +92,8 @@ namespace detail {
     ///
     /// @see line_line_signed_distance, line_line_signed_distance_gradient
     template <typename T>
-    Eigen::Matrix<T, 12, 12> line_line_signed_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+    line_line_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -121,7 +124,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_signed_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_signed_distance(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -156,7 +159,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_signed_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_signed_distance_gradient(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -190,7 +193,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_signed_distance_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
diff --git a/src/ipc/distance/signed/point_line.cpp b/src/ipc/distance/signed/point_line.cpp
index 6000aaf20..d13c728b1 100644
--- a/src/ipc/distance/signed/point_line.cpp
+++ b/src/ipc/distance/signed/point_line.cpp
@@ -5,7 +5,8 @@
 namespace ipc::detail {
 
 template <typename T>
-Eigen::Matrix<T, 6, 6> point_line_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 6>
+point_line_signed_distance_hessian(
     Eigen::ConstRef<Eigen::Vector2<T>> p,
     Eigen::ConstRef<Eigen::Vector2<T>> e0,
     Eigen::ConstRef<Eigen::Vector2<T>> e1)
@@ -25,8 +26,18 @@ Eigen::Matrix<T, 6, 6> point_line_signed_distance_hessian(
     // ---------------------------------------------------------
     // Contract the normal Hessian (2x36) with vector v (2x1).
     // Result is 1x36, mapped to 6x6.
-    hess = (hess_n.reshaped(Eigen::fix<2>, Eigen::fix<36>).transpose() * v)
-               .reshaped(Eigen::fix<6>, Eigen::fix<6>);
+    // We spell the two reshapes out as Maps rather than calling
+    // .reshaped(Eigen::fix<...>). Both are fixed-size views over the same
+    // column-major storage, so the result is identical, but .reshaped()
+    // returns a nested expression template that nvcc does not handle, and
+    // this contraction has to stay device-callable.
+    {
+        const Eigen::Map<const Eigen::Matrix<T, 2, 36>> hess_n_2_36(
+            hess_n.data());
+        const Eigen::Matrix<T, 36, 1> contracted =
+            (hess_n_2_36.transpose() * v).eval();
+        hess = Eigen::Map<const Eigen::Matrix<T, 6, 6>>(contracted.data());
+    }
 
     // ---------------------------------------------------------
     // 2. Add Jacobian Terms (Product Rule Corrections)
@@ -80,8 +91,10 @@ Eigen::Matrix<T, 6, 6> point_line_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector2<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector2<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_POINT_LINE_SIGNED_DISTANCE_HESSIAN(float);
 IPC_INSTANTIATE_POINT_LINE_SIGNED_DISTANCE_HESSIAN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_POINT_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch<float>);
 IPC_INSTANTIATE_POINT_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch<double>);
diff --git a/src/ipc/distance/signed/point_line.hpp b/src/ipc/distance/signed/point_line.hpp
index fc5556e7f..b76ae2096 100644
--- a/src/ipc/distance/signed/point_line.hpp
+++ b/src/ipc/distance/signed/point_line.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/geometry/normal.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -20,7 +21,7 @@ namespace detail {
     /// @return    The signed scalar distance from p to the infinite line through e0 and e1.
     /// @note      The edge must be non-degenerate (e0 != e1).
     template <typename T>
-    inline T point_line_signed_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_line_signed_distance(
         Eigen::ConstRef<Eigen::Vector2<T>> p,
         Eigen::ConstRef<Eigen::Vector2<T>> e0,
         Eigen::ConstRef<Eigen::Vector2<T>> e1)
@@ -41,7 +42,8 @@ namespace detail {
     /// @return    A 6-vector containing the gradient of the signed distance.
     /// @note      The edge must be non-degenerate (e0 != e1).
     template <typename T>
-    inline Eigen::Vector<T, 6> point_line_signed_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 6>
+    point_line_signed_distance_gradient(
         Eigen::ConstRef<Eigen::Vector2<T>> p,
         Eigen::ConstRef<Eigen::Vector2<T>> e0,
         Eigen::ConstRef<Eigen::Vector2<T>> e1)
@@ -70,7 +72,8 @@ namespace detail {
     /// @return    A 6x6 Hessian matrix of the signed distance.
     /// @note      The edge must be non-degenerate (e0 != e1).
     template <typename T>
-    Eigen::Matrix<T, 6, 6> point_line_signed_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 6>
+    point_line_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector2<T>> p,
         Eigen::ConstRef<Eigen::Vector2<T>> e0,
         Eigen::ConstRef<Eigen::Vector2<T>> e1);
@@ -90,7 +93,7 @@ namespace detail {
 /// @return    The signed scalar distance from p to the infinite line through e0 and e1.
 /// @note      The edge must be non-degenerate (e0 != e1).
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_signed_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -114,7 +117,7 @@ inline auto point_line_signed_distance(
 /// @return    A 6-vector containing the gradient of the signed distance.
 /// @note      The edge must be non-degenerate (e0 != e1).
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_signed_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -138,7 +141,7 @@ inline auto point_line_signed_distance_gradient(
 /// @return    A 6x6 Hessian matrix of the signed distance.
 /// @note      The edge must be non-degenerate (e0 != e1).
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
diff --git a/src/ipc/distance/signed/point_plane.cpp b/src/ipc/distance/signed/point_plane.cpp
index 3b5b41c3a..288dd3470 100644
--- a/src/ipc/distance/signed/point_plane.cpp
+++ b/src/ipc/distance/signed/point_plane.cpp
@@ -5,7 +5,8 @@
 namespace ipc::detail {
 
 template <typename T>
-Eigen::Matrix<T, 12, 12> point_plane_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+point_plane_signed_distance_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -42,9 +43,19 @@ Eigen::Matrix<T, 12, 12> point_plane_signed_distance_hessian(
 
     // A. Contraction of the normal Hessian tensor with vector v
     // hess_n is 3x81. v is 3x1. Result is 1x81, which maps to 9x9.
-    hess.template block<9, 9>(3, 3) =
-        (hess_n.reshaped(Eigen::fix<3>, Eigen::fix<81>).transpose() * v)
-            .reshaped(Eigen::fix<9>, Eigen::fix<9>);
+    // We spell the two reshapes out as Maps rather than calling
+    // .reshaped(Eigen::fix<...>). Both are fixed-size views over the same
+    // column-major storage, so the result is identical, but .reshaped()
+    // returns a nested expression template that nvcc does not handle, and
+    // this contraction has to stay device-callable.
+    {
+        const Eigen::Map<const Eigen::Matrix<T, 3, 81>> hess_n_3_81(
+            hess_n.data());
+        const Eigen::Matrix<T, 81, 1> contracted =
+            (hess_n_3_81.transpose() * v).eval();
+        hess.template block<9, 9>(3, 3) =
+            Eigen::Map<const Eigen::Matrix<T, 9, 9>>(contracted.data());
+    }
 
     // B. Subtract first derivative terms (Product Rule corrections)
     // Extract 3x3 Jacobian blocks for t0, t1, t2
@@ -79,8 +90,10 @@ Eigen::Matrix<T, 12, 12> point_plane_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_POINT_PLANE_SIGNED_DISTANCE_HESSIAN(float);
 IPC_INSTANTIATE_POINT_PLANE_SIGNED_DISTANCE_HESSIAN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_POINT_PLANE_SIGNED_DISTANCE_HESSIAN(SimdBatch<float>);
 IPC_INSTANTIATE_POINT_PLANE_SIGNED_DISTANCE_HESSIAN(SimdBatch<double>);
diff --git a/src/ipc/distance/signed/point_plane.hpp b/src/ipc/distance/signed/point_plane.hpp
index 645646f18..0272b7700 100644
--- a/src/ipc/distance/signed/point_plane.hpp
+++ b/src/ipc/distance/signed/point_plane.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/geometry/normal.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
@@ -18,7 +19,7 @@ namespace detail {
     /// @param t2  Third vertex of the triangle (3D).
     /// @return    The signed distance from p to the plane of the triangle.
     template <typename T>
-    inline T point_plane_signed_distance(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_plane_signed_distance(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -38,7 +39,8 @@ namespace detail {
     /// @param t2  Third vertex of the triangle (3D).
     /// @return    A 12-vector containing the gradient of the signed distance.
     template <typename T>
-    inline Eigen::Vector<T, 12> point_plane_signed_distance_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 12>
+    point_plane_signed_distance_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -70,7 +72,8 @@ namespace detail {
     /// @param t2  Third vertex of the triangle (3D).
     /// @return    A 12x12 Hessian matrix of the signed distance.
     template <typename T>
-    Eigen::Matrix<T, 12, 12> point_plane_signed_distance_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12>
+    point_plane_signed_distance_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -93,7 +96,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_signed_distance(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_signed_distance(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -120,7 +123,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_signed_distance_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_signed_distance_gradient(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -146,7 +149,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_plane_signed_distance_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_signed_distance_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
diff --git a/src/ipc/geometry/CMakeLists.txt b/src/ipc/geometry/CMakeLists.txt
index a771e9d6e..4e4405f05 100644
--- a/src/ipc/geometry/CMakeLists.txt
+++ b/src/ipc/geometry/CMakeLists.txt
@@ -1,16 +1,21 @@
 set(SOURCES
-  angle.cpp
   angle.hpp
-  area.cpp
   area.hpp
   intersection.cpp
   intersection.hpp
-  normal.cpp
   normal.hpp
 )
 
 target_sources(ipc_toolkit PRIVATE ${SOURCES})
 
+# Definitions shared between host C++ and CUDA device code; under CUDA the
+# helper generates a .cu wrapper per file (see ipc_toolkit_shared_device_sources).
+ipc_toolkit_target_shared_device_sources(ipc_toolkit
+  angle.cpp
+  area.cpp
+  normal.cpp
+)
+
 ################################################################################
 # Subfolders
 ################################################################################
\ No newline at end of file
diff --git a/src/ipc/geometry/angle.cpp b/src/ipc/geometry/angle.cpp
index 0014899aa..74af3f471 100644
--- a/src/ipc/geometry/angle.cpp
+++ b/src/ipc/geometry/angle.cpp
@@ -6,11 +6,24 @@
 #include <ipc/utils/simd.hpp>
 
 #include <array>
-#include <utility>
 
 namespace ipc::detail {
 
 namespace {
+    /// @brief Return type of dihedral_normal_jacobians below.
+    ///
+    /// An aggregate rather than a std::pair: nvcc admits libstdc++'s constexpr
+    /// std::pair into device code only under --expt-relaxed-constexpr, which
+    /// covers compile-time evaluation, and a pair built at run time inside a
+    /// kernel silently comes back zero-filled. See
+    /// ipc::NormalizationAndJacobian.
+    template <typename T> struct DihedralNormalJacobians {
+        /// @brief Jacobian of n0 = normal(x0, x1, x2) w.r.t. all 12 DOFs.
+        Eigen::Matrix<T, 3, 12> dn0_dx;
+        /// @brief Jacobian of n1 = normal(x1, x0, x3) w.r.t. all 12 DOFs.
+        Eigen::Matrix<T, 3, 12> dn1_dx;
+    };
+
     /// @brief Jacobians of the two triangle normals with respect to all 12 DOFs
     /// (x0, x1, x2, x3).
     ///
@@ -18,7 +31,7 @@ namespace {
     /// zero. n1 = normal(x1, x0, x3) comes back in (x1, x0, x3) order, so we
     /// permute its columns into place and zero the x2 block.
     template <typename T>
-    inline std::pair<Eigen::Matrix<T, 3, 12>, Eigen::Matrix<T, 3, 12>>
+    IPC_TOOLKIT_HOST_DEVICE inline DihedralNormalJacobians<T>
     dihedral_normal_jacobians(
         Eigen::ConstRef<Eigen::Vector3<T>> x0,
         Eigen::ConstRef<Eigen::Vector3<T>> x1,
@@ -29,17 +42,24 @@ namespace {
         dn0_dx.template leftCols<9>() = triangle_normal_jacobian(x0, x1, x2);
         dn0_dx.template rightCols<3>().setZero();
 
+        // We scatter the per-vertex blocks explicitly rather than writing
+        // dn1_dx(Eigen::all, idx). Eigen's index slicing is unavailable in
+        // device code, and this function has to stay device-callable. The
+        // permutation is the same one: x0 <- block 1, x1 <- block 0, x3 <-
+        // block 2, with the x2 block zeroed.
+        const Eigen::Matrix<T, 3, 9> dn1 = triangle_normal_jacobian(x1, x0, x3);
         Eigen::Matrix<T, 3, 12> dn1_dx;
-        const std::array<int, 9> idx = { { 3, 4, 5, 0, 1, 2, 9, 10, 11 } };
-        dn1_dx(Eigen::all, idx) = triangle_normal_jacobian(x1, x0, x3);
+        dn1_dx.template middleCols<3>(0) = dn1.template middleCols<3>(3);
+        dn1_dx.template middleCols<3>(3) = dn1.template middleCols<3>(0);
         dn1_dx.template middleCols<3>(6).setZero();
+        dn1_dx.template middleCols<3>(9) = dn1.template middleCols<3>(6);
 
         return { dn0_dx, dn1_dx };
     }
 } // namespace
 
 template <typename T>
-T dihedral_angle(
+IPC_TOOLKIT_HOST_DEVICE T dihedral_angle(
     Eigen::ConstRef<Eigen::Vector3<T>> x0,
     Eigen::ConstRef<Eigen::Vector3<T>> x1,
     Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -56,7 +76,7 @@ T dihedral_angle(
 }
 
 template <typename T>
-Eigen::Vector<T, 12> dihedral_angle_gradient(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12> dihedral_angle_gradient(
     Eigen::ConstRef<Eigen::Vector3<T>> x0,
     Eigen::ConstRef<Eigen::Vector3<T>> x1,
     Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -90,7 +110,7 @@ Eigen::Vector<T, 12> dihedral_angle_gradient(
 }
 
 template <typename T>
-Eigen::Matrix<T, 12, 12> dihedral_angle_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12> dihedral_angle_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> x0,
     Eigen::ConstRef<Eigen::Vector3<T>> x1,
     Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -368,8 +388,10 @@ Eigen::Matrix<T, 12, 12> dihedral_angle_hessian(
     template Eigen::Vector<T, 12> dihedral_angle_gradient<T>(Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>); \
     template Eigen::Matrix<T, 12, 12> dihedral_angle_hessian<T>(Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_ANGLE(float);
 IPC_INSTANTIATE_ANGLE(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_ANGLE(SimdBatch<float>);
 IPC_INSTANTIATE_ANGLE(SimdBatch<double>);
diff --git a/src/ipc/geometry/angle.hpp b/src/ipc/geometry/angle.hpp
index 6ae1ce2a5..eac624435 100644
--- a/src/ipc/geometry/angle.hpp
+++ b/src/ipc/geometry/angle.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 namespace ipc {
@@ -7,7 +8,7 @@ namespace ipc {
 namespace detail {
     /// @note Prefer the ipc::dihedral_angle front end.
     template <typename T>
-    T dihedral_angle(
+    IPC_TOOLKIT_HOST_DEVICE T dihedral_angle(
         Eigen::ConstRef<Eigen::Vector3<T>> x0,
         Eigen::ConstRef<Eigen::Vector3<T>> x1,
         Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -15,7 +16,7 @@ namespace detail {
 
     /// @note Prefer the ipc::dihedral_angle_gradient front end.
     template <typename T>
-    Eigen::Vector<T, 12> dihedral_angle_gradient(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Vector<T, 12> dihedral_angle_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> x0,
         Eigen::ConstRef<Eigen::Vector3<T>> x1,
         Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -23,7 +24,7 @@ namespace detail {
 
     /// @note Prefer the ipc::dihedral_angle_hessian front end.
     template <typename T>
-    Eigen::Matrix<T, 12, 12> dihedral_angle_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 12, 12> dihedral_angle_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> x0,
         Eigen::ConstRef<Eigen::Vector3<T>> x1,
         Eigen::ConstRef<Eigen::Vector3<T>> x2,
@@ -44,7 +45,7 @@ template <
     typename DerivedX1,
     typename DerivedX2,
     typename DerivedX3>
-inline auto dihedral_angle(
+IPC_TOOLKIT_HOST_DEVICE inline auto dihedral_angle(
     const Eigen::MatrixBase<DerivedX0>& x0,
     const Eigen::MatrixBase<DerivedX1>& x1,
     const Eigen::MatrixBase<DerivedX2>& x2,
@@ -68,7 +69,7 @@ template <
     typename DerivedX1,
     typename DerivedX2,
     typename DerivedX3>
-inline auto dihedral_angle_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto dihedral_angle_gradient(
     const Eigen::MatrixBase<DerivedX0>& x0,
     const Eigen::MatrixBase<DerivedX1>& x1,
     const Eigen::MatrixBase<DerivedX2>& x2,
@@ -92,7 +93,7 @@ template <
     typename DerivedX1,
     typename DerivedX2,
     typename DerivedX3>
-inline auto dihedral_angle_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto dihedral_angle_hessian(
     const Eigen::MatrixBase<DerivedX0>& x0,
     const Eigen::MatrixBase<DerivedX1>& x1,
     const Eigen::MatrixBase<DerivedX2>& x2,
diff --git a/src/ipc/geometry/area.cpp b/src/ipc/geometry/area.cpp
index 150ac7245..c33c9f3e1 100644
--- a/src/ipc/geometry/area.cpp
+++ b/src/ipc/geometry/area.cpp
@@ -6,7 +6,7 @@ namespace ipc::autogen {
 
 // dA is (9×1) flattened in column-major order
 template <typename T>
-void triangle_area_gradient(
+IPC_TOOLKIT_HOST_DEVICE void triangle_area_gradient(
     T t0_x,
     T t0_y,
     T t0_z,
@@ -48,8 +48,10 @@ void triangle_area_gradient(
 #define IPC_INSTANTIATE_AREA_AUTOGEN(T)                                        \
     template void triangle_area_gradient<T>(T, T, T, T, T, T, T, T, T, T[9])
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_AREA_AUTOGEN(float);
 IPC_INSTANTIATE_AREA_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_AREA_AUTOGEN(SimdBatch<float>);
 IPC_INSTANTIATE_AREA_AUTOGEN(SimdBatch<double>);
diff --git a/src/ipc/geometry/area.hpp b/src/ipc/geometry/area.hpp
index b57bb8d85..6b65cd303 100644
--- a/src/ipc/geometry/area.hpp
+++ b/src/ipc/geometry/area.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 #include <Eigen/Geometry>
@@ -15,7 +16,7 @@ namespace autogen {
     // clang-format off
     /// dA is (9×1) flattened in column-major order
     template <typename T>
-    void triangle_area_gradient(
+    IPC_TOOLKIT_HOST_DEVICE void triangle_area_gradient(
         T t0_x, T t0_y, T t0_z, T t1_x, T t1_y, T t1_z, T t2_x, T t2_y, T t2_z, T dA[9]);
     // clang-format on
 
@@ -29,7 +30,7 @@ namespace detail {
     /// @param e1 The second vertex of the edge.
     /// @return The length of the edge.
     template <typename T, int dim>
-    inline T edge_length(
+    IPC_TOOLKIT_HOST_DEVICE inline T edge_length(
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
     {
@@ -44,7 +45,8 @@ namespace detail {
     /// @param e1 The second vertex of the edge.
     /// @return The gradient of the edge's length wrt e0, and e1.
     template <typename T, int dim>
-    inline Eigen::Vector<T, 2 * dim> edge_length_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 2 * dim>
+    edge_length_gradient(
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
     {
@@ -69,7 +71,7 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The area of the triangle.
     template <typename T>
-    inline T triangle_area(
+    IPC_TOOLKIT_HOST_DEVICE inline T triangle_area(
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
         Eigen::ConstRef<Eigen::Vector3<T>> t2)
@@ -84,7 +86,7 @@ namespace detail {
     /// @param t2 The third vertex of the triangle.
     /// @return The gradient of the triangle's area t0, t1, and t2.
     template <typename T>
-    inline Eigen::Vector<T, 9> triangle_area_gradient(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 9> triangle_area_gradient(
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
         Eigen::ConstRef<Eigen::Vector3<T>> t2)
@@ -102,7 +104,7 @@ namespace detail {
 /// @param e1 The second vertex of the edge.
 /// @return The length of the edge.
 template <typename DerivedE0, typename DerivedE1>
-inline auto edge_length(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_length(
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
 {
@@ -125,7 +127,7 @@ inline auto edge_length(
 /// @param e1 The second vertex of the edge.
 /// @return The gradient of the edge's length wrt e0, and e1.
 template <typename DerivedE0, typename DerivedE1>
-inline auto edge_length_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_length_gradient(
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
 {
@@ -152,7 +154,7 @@ inline auto edge_length_gradient(
 /// @param t2 The third vertex of the triangle.
 /// @return The area of the triangle.
 template <typename DerivedT0, typename DerivedT1, typename DerivedT2>
-inline auto triangle_area(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_area(
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
     const Eigen::MatrixBase<DerivedT2>& t2)
@@ -172,7 +174,7 @@ inline auto triangle_area(
 /// @param t2 The third vertex of the triangle.
 /// @return The gradient of the triangle's area t0, t1, and t2.
 template <typename DerivedT0, typename DerivedT1, typename DerivedT2>
-inline auto triangle_area_gradient(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_area_gradient(
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
     const Eigen::MatrixBase<DerivedT2>& t2)
diff --git a/src/ipc/geometry/normal.cpp b/src/ipc/geometry/normal.cpp
index ebba14be2..6a8db9e16 100644
--- a/src/ipc/geometry/normal.cpp
+++ b/src/ipc/geometry/normal.cpp
@@ -1,4 +1,4 @@
-#include "ipc/geometry/normal.hpp"
+#include "normal.hpp"
 
 #include <ipc/config.hpp>
 #include <ipc/utils/simd.hpp>
@@ -8,7 +8,7 @@ namespace ipc::detail {
 // --- point-line normal functions -------------------------------------------
 
 template <typename T>
-VectorMax3<T> point_line_unnormalized_normal(
+IPC_TOOLKIT_HOST_DEVICE VectorMax3<T> point_line_unnormalized_normal(
     Eigen::ConstRef<VectorMax3<T>> p,
     Eigen::ConstRef<VectorMax3<T>> e0,
     Eigen::ConstRef<VectorMax3<T>> e1)
@@ -31,7 +31,8 @@ VectorMax3<T> point_line_unnormalized_normal(
 }
 
 template <typename T>
-MatrixMax<T, 3, 9> point_line_unnormalized_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 3, 9>
+point_line_unnormalized_normal_jacobian(
     Eigen::ConstRef<VectorMax3<T>> p,
     Eigen::ConstRef<VectorMax3<T>> e0,
     Eigen::ConstRef<VectorMax3<T>> e1)
@@ -67,7 +68,8 @@ MatrixMax<T, 3, 9> point_line_unnormalized_normal_jacobian(
 }
 
 template <typename T>
-MatrixMax<T, 27, 9> point_line_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 27, 9>
+point_line_unnormalized_normal_hessian(
     Eigen::ConstRef<VectorMax3<T>> p,
     Eigen::ConstRef<VectorMax3<T>> e0,
     Eigen::ConstRef<VectorMax3<T>> e1)
@@ -173,7 +175,7 @@ MatrixMax<T, 27, 9> point_line_unnormalized_normal_hessian(
 }
 
 template <typename T>
-MatrixMax<T, 27, 9> point_line_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 27, 9> point_line_normal_hessian(
     Eigen::ConstRef<VectorMax3<T>> p,
     Eigen::ConstRef<VectorMax3<T>> e0,
     Eigen::ConstRef<VectorMax3<T>> e1)
@@ -217,7 +219,7 @@ MatrixMax<T, 27, 9> point_line_normal_hessian(
 
 namespace {
     template <typename T>
-    void set_cross_product_matrix_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void set_cross_product_matrix_jacobian(
         Eigen::Ref<Eigen::Matrix<T, 9, 3>> Jx, T chain_rule = T(1.0))
     {
         Jx(2, 1) = Jx(3, 2) = Jx(7, 0) = -chain_rule;
@@ -225,7 +227,8 @@ namespace {
     }
 } // namespace
 
-template <typename T> Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian()
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian()
 {
     Eigen::Matrix<T, 9, 3> J = Eigen::Matrix<T, 9, 3>::Zero();
     J(2, 1) = J(3, 2) = J(7, 0) = T(-1);
@@ -234,7 +237,8 @@ template <typename T> Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian()
 }
 
 template <typename T>
-Eigen::Matrix<T, 27, 9> triangle_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 27, 9>
+triangle_unnormalized_normal_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> a,
     Eigen::ConstRef<Eigen::Vector3<T>> b,
     Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -263,7 +267,7 @@ Eigen::Matrix<T, 27, 9> triangle_unnormalized_normal_hessian(
 }
 
 template <typename T>
-Eigen::Matrix<T, 27, 9> triangle_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 27, 9> triangle_normal_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> a,
     Eigen::ConstRef<Eigen::Vector3<T>> b,
     Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -300,7 +304,8 @@ Eigen::Matrix<T, 27, 9> triangle_normal_hessian(
 // --- line-line normal functions ---------------------------------------------
 
 template <typename T>
-Eigen::Matrix<T, 36, 12> line_line_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 36, 12>
+line_line_unnormalized_normal_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -340,7 +345,7 @@ Eigen::Matrix<T, 36, 12> line_line_unnormalized_normal_hessian(
 }
 
 template <typename T>
-Eigen::Matrix<T, 36, 12> line_line_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 36, 12> line_line_normal_hessian(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -390,8 +395,10 @@ Eigen::Matrix<T, 36, 12> line_line_normal_hessian(
     template Eigen::Matrix<T, 36, 12> line_line_unnormalized_normal_hessian<T>(Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>); \
     template Eigen::Matrix<T, 36, 12> line_line_normal_hessian<T>(Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>, Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_NORMAL(float);
 IPC_INSTANTIATE_NORMAL(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_NORMAL(SimdBatch<float>);
 IPC_INSTANTIATE_NORMAL(SimdBatch<double>);
diff --git a/src/ipc/geometry/normal.hpp b/src/ipc/geometry/normal.hpp
index 66df2305c..5c59478b6 100644
--- a/src/ipc/geometry/normal.hpp
+++ b/src/ipc/geometry/normal.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 #include <ipc/utils/simd.hpp>
 
@@ -21,8 +22,7 @@ struct NormalizationAndJacobian {
     Eigen::Matrix<T, dim, dim, Eigen::ColMajor, max_dim, max_dim> jacobian;
 };
 
-/// @brief A normalized vector paired with the first two derivatives of the
-/// normalization.
+/// @brief A normalized vector paired with the first two derivatives of the normalization.
 /// @tparam T The scalar type.
 /// @tparam dim The dimension (2 or 3).
 /// @tparam max_dim The maximum dimension (2 or 3).
@@ -49,7 +49,7 @@ namespace detail {
     /// @param x The input vector.
     /// @return The normalized vector and its Jacobian.
     template <typename T, int dim>
-    inline NormalizationAndJacobian<T, dim>
+    IPC_TOOLKIT_HOST_DEVICE inline NormalizationAndJacobian<T, dim>
     normalization_and_jacobian(const Eigen::Vector<T, dim>& x)
     {
         static_assert(dim == 2 || dim == 3, "normalization is only 2D or 3D");
@@ -67,7 +67,7 @@ namespace detail {
     /// @param x The input vector.
     /// @return The normalized vector, its Jacobian, and its Hessian.
     template <typename T, int dim>
-    inline NormalizationAndJacobianAndHessian<T, dim>
+    IPC_TOOLKIT_HOST_DEVICE inline NormalizationAndJacobianAndHessian<T, dim>
     normalization_and_jacobian_and_hessian(const Eigen::Vector<T, dim>& x)
     {
         static_assert(dim == 2 || dim == 3, "normalization is only 2D or 3D");
@@ -98,7 +98,8 @@ namespace detail {
 /// @param x The input vector.
 /// @return The normalized vector and its Jacobian.
 template <typename DerivedX>
-inline auto normalization_and_jacobian(const Eigen::MatrixBase<DerivedX>& x)
+IPC_TOOLKIT_HOST_DEVICE inline auto
+normalization_and_jacobian(const Eigen::MatrixBase<DerivedX>& x)
 {
     using T = typename DerivedX::Scalar;
     using DynamicResult = NormalizationAndJacobian<T, Eigen::Dynamic, 3>;
@@ -121,7 +122,8 @@ inline auto normalization_and_jacobian(const Eigen::MatrixBase<DerivedX>& x)
 /// @param x The input vector.
 /// @return The Jacobian of the normalization operation.
 template <typename DerivedX>
-inline auto normalization_jacobian(const Eigen::MatrixBase<DerivedX>& x)
+IPC_TOOLKIT_HOST_DEVICE inline auto
+normalization_jacobian(const Eigen::MatrixBase<DerivedX>& x)
 {
     using T = typename DerivedX::Scalar;
 
@@ -143,7 +145,7 @@ inline auto normalization_jacobian(const Eigen::MatrixBase<DerivedX>& x)
 /// @param x The input vector.
 /// @return The normalized vector, its Jacobian, and its Hessian.
 template <typename DerivedX>
-inline auto
+IPC_TOOLKIT_HOST_DEVICE inline auto
 normalization_and_jacobian_and_hessian(const Eigen::MatrixBase<DerivedX>& x)
 {
     using T = typename DerivedX::Scalar;
@@ -180,7 +182,8 @@ normalization_and_jacobian_and_hessian(const Eigen::MatrixBase<DerivedX>& x)
 /// @param x The input vector.
 /// @return The Hessian of the normalization operation.
 template <typename DerivedX>
-inline auto normalization_hessian(const Eigen::MatrixBase<DerivedX>& x)
+IPC_TOOLKIT_HOST_DEVICE inline auto
+normalization_hessian(const Eigen::MatrixBase<DerivedX>& x)
 {
     return normalization_and_jacobian_and_hessian(x).hessian;
 }
@@ -193,7 +196,7 @@ namespace detail {
     /// @param v Vector to create the cross product matrix for.
     /// @return The cross product matrix of the vector.
     template <typename T>
-    inline Eigen::Matrix3<T>
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix3<T>
     cross_product_matrix(Eigen::ConstRef<Eigen::Vector3<T>> v)
     {
         Eigen::Matrix3<T> m;
@@ -208,7 +211,8 @@ namespace detail {
     /// @brief Computes the Jacobian of the cross product matrix.
     /// @return The Jacobian of the cross product matrix.
     template <typename T>
-    Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian();
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 9, 3>
+    cross_product_matrix_jacobian();
 
     // =========================================================================
 
@@ -224,7 +228,7 @@ namespace detail {
     /// @param e1 The end position of the line.
     /// @return The unnormalized normal vector.
     template <typename T>
-    VectorMax3<T> point_line_unnormalized_normal(
+    IPC_TOOLKIT_HOST_DEVICE VectorMax3<T> point_line_unnormalized_normal(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1);
@@ -235,7 +239,7 @@ namespace detail {
     /// @param e1 The end position of the line.
     /// @return The normal vector.
     template <typename T>
-    inline VectorMax3<T> point_line_normal(
+    IPC_TOOLKIT_HOST_DEVICE inline VectorMax3<T> point_line_normal(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1)
@@ -250,7 +254,8 @@ namespace detail {
     /// @param e1 The end position of the line.
     /// @return The Jacobian of the unnormalized normal vector.
     template <typename T>
-    MatrixMax<T, 3, 9> point_line_unnormalized_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 3, 9>
+    point_line_unnormalized_normal_jacobian(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1);
@@ -262,7 +267,8 @@ namespace detail {
     /// @return The Hessian of the unnormalized normal vector of the point-line
     /// pair.
     template <typename T>
-    MatrixMax<T, 27, 9> point_line_unnormalized_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 27, 9>
+    point_line_unnormalized_normal_hessian(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1);
@@ -273,7 +279,8 @@ namespace detail {
     /// @param e1 The end position of the line.
     /// @return The Jacobian of the normal vector.
     template <typename T>
-    inline MatrixMax<T, 3, 9> point_line_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline MatrixMax<T, 3, 9>
+    point_line_normal_jacobian(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1)
@@ -289,7 +296,7 @@ namespace detail {
     /// @param e1 The end position of the line.
     /// @return The Hessian of the normal vector.
     template <typename T>
-    MatrixMax<T, 27, 9> point_line_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE MatrixMax<T, 27, 9> point_line_normal_hessian(
         Eigen::ConstRef<VectorMax3<T>> p,
         Eigen::ConstRef<VectorMax3<T>> e0,
         Eigen::ConstRef<VectorMax3<T>> e1);
@@ -310,7 +317,8 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The unnormalized normal vector of the triangle.
     template <typename T>
-    inline Eigen::Vector3<T> triangle_unnormalized_normal(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T>
+    triangle_unnormalized_normal(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -324,7 +332,7 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The normal vector of the triangle.
     template <typename T>
-    inline Eigen::Vector3<T> triangle_normal(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T> triangle_normal(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -339,7 +347,8 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The Jacobian of the unnormalized normal vector of the triangle.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 9> triangle_unnormalized_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 9>
+    triangle_unnormalized_normal_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -358,7 +367,8 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The Hessian of the unnormalized normal vector of the triangle.
     template <typename T>
-    Eigen::Matrix<T, 27, 9> triangle_unnormalized_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 27, 9>
+    triangle_unnormalized_normal_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c);
@@ -369,7 +379,8 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The Jacobian of the normal vector of the triangle.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 9> triangle_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 9>
+    triangle_normal_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c)
@@ -385,7 +396,7 @@ namespace detail {
     /// @param c The third vertex of the triangle.
     /// @return The Hessian of the normal vector of the triangle.
     template <typename T>
-    Eigen::Matrix<T, 27, 9> triangle_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 27, 9> triangle_normal_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> a,
         Eigen::ConstRef<Eigen::Vector3<T>> b,
         Eigen::ConstRef<Eigen::Vector3<T>> c);
@@ -407,7 +418,8 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The unnormalized normal vector of the two lines.
     template <typename T>
-    inline Eigen::Vector3<T> line_line_unnormalized_normal(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T>
+    line_line_unnormalized_normal(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -423,7 +435,7 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The normal vector of the two lines.
     template <typename T>
-    inline Eigen::Vector3<T> line_line_normal(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T> line_line_normal(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -440,7 +452,8 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The Jacobian of the unnormalized normal vector of the two lines.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 12> line_line_unnormalized_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 12>
+    line_line_unnormalized_normal_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -461,7 +474,8 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The Jacobian of the normal vector of the two lines.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 12> line_line_normal_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 12>
+    line_line_normal_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -481,7 +495,8 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The Hessian of the unnormalized normal vector of the two lines.
     template <typename T>
-    Eigen::Matrix<T, 36, 12> line_line_unnormalized_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 36, 12>
+    line_line_unnormalized_normal_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -494,7 +509,7 @@ namespace detail {
     /// @param eb1 The second vertex of the second line.
     /// @return The Hessian of the normal vector of the two lines.
     template <typename T>
-    Eigen::Matrix<T, 36, 12> line_line_normal_hessian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 36, 12> line_line_normal_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -509,7 +524,8 @@ namespace detail {
 /// @brief Computes the Jacobian of the cross product matrix.
 /// @return The Jacobian of the cross product matrix.
 template <typename T>
-inline Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian()
+IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 9, 3>
+cross_product_matrix_jacobian()
 {
     return detail::cross_product_matrix_jacobian<T>();
 }
@@ -518,7 +534,8 @@ inline Eigen::Matrix<T, 9, 3> cross_product_matrix_jacobian()
 /// @param v Vector to create the cross product matrix for.
 /// @return The cross product matrix of the vector.
 template <typename DerivedX>
-inline auto cross_product_matrix(const Eigen::MatrixBase<DerivedX>& v)
+IPC_TOOLKIT_HOST_DEVICE inline auto
+cross_product_matrix(const Eigen::MatrixBase<DerivedX>& v)
 {
     using T = typename DerivedX::Scalar;
     return detail::cross_product_matrix<T>(v);
@@ -530,7 +547,7 @@ inline auto cross_product_matrix(const Eigen::MatrixBase<DerivedX>& v)
 /// @param e1 The end position of the line.
 /// @return The unnormalized normal vector.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_unnormalized_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -545,7 +562,7 @@ inline auto point_line_unnormalized_normal(
 /// @param e1 The end position of the line.
 /// @return The normal vector.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -561,7 +578,7 @@ inline auto point_line_normal(
 /// @param e1 The end position of the line.
 /// @return The Jacobian of the unnormalized normal vector.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_unnormalized_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -578,7 +595,7 @@ inline auto point_line_unnormalized_normal_jacobian(
 /// @return The Hessian of the unnormalized normal vector of the point-line
 /// pair.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -593,7 +610,7 @@ inline auto point_line_unnormalized_normal_hessian(
 /// @param e1 The end position of the line.
 /// @return The Jacobian of the normal vector.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -608,7 +625,7 @@ inline auto point_line_normal_jacobian(
 /// @param e1 The end position of the line.
 /// @return The Hessian of the normal vector.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_line_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -623,7 +640,7 @@ inline auto point_line_normal_hessian(
 /// @param c The third vertex of the triangle.
 /// @return The unnormalized normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_unnormalized_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -638,7 +655,7 @@ inline auto triangle_unnormalized_normal(
 /// @param c The third vertex of the triangle.
 /// @return The normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -654,7 +671,7 @@ inline auto triangle_normal(
 /// @param c The third vertex of the triangle.
 /// @return The Jacobian of the unnormalized normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_unnormalized_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal_jacobian(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -670,7 +687,7 @@ inline auto triangle_unnormalized_normal_jacobian(
 /// @param c The third vertex of the triangle.
 /// @return The Hessian of the unnormalized normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal_hessian(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -685,7 +702,7 @@ inline auto triangle_unnormalized_normal_hessian(
 /// @param c The third vertex of the triangle.
 /// @return The Jacobian of the normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal_jacobian(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -700,7 +717,7 @@ inline auto triangle_normal_jacobian(
 /// @param c The third vertex of the triangle.
 /// @return The Hessian of the normal vector of the triangle.
 template <typename DerivedA, typename DerivedB, typename DerivedC>
-inline auto triangle_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal_hessian(
     const Eigen::MatrixBase<DerivedA>& a,
     const Eigen::MatrixBase<DerivedB>& b,
     const Eigen::MatrixBase<DerivedC>& c)
@@ -720,7 +737,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_unnormalized_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_unnormalized_normal(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -741,7 +758,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_normal(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_normal(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -763,7 +780,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_unnormalized_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_unnormalized_normal_jacobian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -785,7 +802,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_normal_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_normal_jacobian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -807,7 +824,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_unnormalized_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_unnormalized_normal_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -828,7 +845,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto line_line_normal_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto line_line_normal_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
diff --git a/src/ipc/math/morton.hpp b/src/ipc/math/morton.hpp
index ccc201cda..d43c06bac 100644
--- a/src/ipc/math/morton.hpp
+++ b/src/ipc/math/morton.hpp
@@ -1,14 +1,16 @@
 #pragma once
 
-#include <algorithm> // for std::clamp
-#include <cstdint>   // for uint64_t
+#include <ipc/config.hpp>     // for IPC_TOOLKIT_HOST_DEVICE
+#include <ipc/utils/simd.hpp> // for clamp
+
+#include <cstdint> // for uint64_t
 
 namespace ipc {
 
 /// @brief Expands a 32-bit integer into 64 bits by inserting 1 zero after each bit.
 /// @param v The 32-bit integer to expand.
 /// @return The expanded 64-bit integer.
-inline uint64_t expand_bits_1(uint64_t v)
+IPC_TOOLKIT_HOST_DEVICE inline uint64_t expand_bits_1(uint64_t v)
 {
     v = (v | (v << 16)) & 0x0000FFFF0000FFFF;
     v = (v | (v << 8)) & 0x00FF00FF00FF00FF;
@@ -21,7 +23,7 @@ inline uint64_t expand_bits_1(uint64_t v)
 /// @brief Expands a 21-bit integer into 63 bits by inserting 2 zeros after each bit.
 /// @param v The 21-bit integer to expand.
 /// @return The expanded 63-bit integer.
-inline uint64_t expand_bits_2(uint64_t v)
+IPC_TOOLKIT_HOST_DEVICE inline uint64_t expand_bits_2(uint64_t v)
 {
     v = (v | v << 32) & 0x1F00000000FFFF;
     v = (v | v << 16) & 0x1F0000FF0000FF;
@@ -35,11 +37,11 @@ inline uint64_t expand_bits_2(uint64_t v)
 /// @param x The x-coordinate of the point.
 /// @param y The y-coordinate of the point.
 /// @return The 64-bit Morton code.
-inline uint64_t morton_2D(double x, double y)
+IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_2D(double x, double y)
 {
     constexpr double scale = 1ULL << 32;
-    x = std::clamp(x * scale, 0.0, scale - 1);
-    y = std::clamp(y * scale, 0.0, scale - 1);
+    x = ipc::clamp(x * scale, 0.0, scale - 1);
+    y = ipc::clamp(y * scale, 0.0, scale - 1);
     uint64_t xx = expand_bits_1(uint64_t(x));
     uint64_t yy = expand_bits_1(uint64_t(y));
     return (xx << 1) | yy;
@@ -50,12 +52,12 @@ inline uint64_t morton_2D(double x, double y)
 /// @param y The y-coordinate of the point.
 /// @param z The z-coordinate of the point.
 /// @return The 63-bit Morton code.
-inline uint64_t morton_3D(double x, double y, double z)
+IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_3D(double x, double y, double z)
 {
     constexpr double scale = 1ULL << 21;
-    x = std::clamp(x * scale, 0.0, scale - 1);
-    y = std::clamp(y * scale, 0.0, scale - 1);
-    z = std::clamp(z * scale, 0.0, scale - 1);
+    x = ipc::clamp(x * scale, 0.0, scale - 1);
+    y = ipc::clamp(y * scale, 0.0, scale - 1);
+    z = ipc::clamp(z * scale, 0.0, scale - 1);
     uint64_t xx = expand_bits_2(uint64_t(x));
     uint64_t yy = expand_bits_2(uint64_t(y));
     uint64_t zz = expand_bits_2(uint64_t(z));
diff --git a/src/ipc/math/scalar_math.hpp b/src/ipc/math/scalar_math.hpp
index b8cf484ae..47160358f 100644
--- a/src/ipc/math/scalar_math.hpp
+++ b/src/ipc/math/scalar_math.hpp
@@ -31,7 +31,7 @@ namespace numext {
     template <                                                                 \
         typename T, typename... Ts,                                            \
         typename = std::enable_if_t<(std::is_same_v<T, Ts> && ...)>>           \
-    inline auto FUNC(const T& x, const Ts&... rest)                            \
+    IPC_TOOLKIT_HOST_DEVICE inline auto FUNC(const T& x, const Ts&... rest)    \
     {                                                                          \
         IPC_TOOLKIT_USING_STD(FUNC)                                            \
         return FUNC(x, rest...);                                               \
@@ -52,9 +52,15 @@ constexpr double MOLLIFIER_THRESHOLD_EPS = 1e-2;
 
 /// @brief Square of `x`, for any scalar the library templates on.
 /// @note Faster than `std::pow(x, 2)`.
-template <typename T> inline T sqr(const T& x) { return x * x; }
+template <typename T> IPC_TOOLKIT_HOST_DEVICE inline T sqr(const T& x)
+{
+    return x * x;
+}
 
 /// @brief Cube of `x`, for any scalar the library templates on.
-template <typename T> inline T cubic(const T& x) { return x * x * x; }
+template <typename T> IPC_TOOLKIT_HOST_DEVICE inline T cubic(const T& x)
+{
+    return x * x * x;
+}
 
 } // namespace ipc
diff --git a/src/ipc/tangent/CMakeLists.txt b/src/ipc/tangent/CMakeLists.txt
index a2d198d00..1603a5110 100644
--- a/src/ipc/tangent/CMakeLists.txt
+++ b/src/ipc/tangent/CMakeLists.txt
@@ -1,13 +1,18 @@
 set(SOURCES
-  closest_point.cpp
   closest_point.hpp
   relative_velocity.hpp
-  tangent_basis.cpp
   tangent_basis.hpp
 )
 
 target_sources(ipc_toolkit PRIVATE ${SOURCES})
 
+# Definitions shared between host C++ and CUDA device code; under CUDA the
+# helper generates a .cu wrapper per file (see ipc_toolkit_shared_device_sources).
+ipc_toolkit_target_shared_device_sources(ipc_toolkit
+  closest_point.cpp
+  tangent_basis.cpp
+)
+
 ################################################################################
 # Subfolders
 ################################################################################
\ No newline at end of file
diff --git a/src/ipc/tangent/closest_point.cpp b/src/ipc/tangent/closest_point.cpp
index a7002076a..0b45b38a2 100644
--- a/src/ipc/tangent/closest_point.cpp
+++ b/src/ipc/tangent/closest_point.cpp
@@ -5,7 +5,7 @@
 namespace ipc::autogen {
 // hess is (6×6) flattened in column-major order
 template <typename T>
-void point_edge_closest_point_2D_hessian(
+IPC_TOOLKIT_HOST_DEVICE void point_edge_closest_point_2D_hessian(
     T p_x, T p_y, T e0_x, T e0_y, T e1_x, T e1_y, T hess[36])
 {
     const T t0 = e0_x - e1_x;
@@ -91,7 +91,7 @@ void point_edge_closest_point_2D_hessian(
 
 // hess is (9×9) flattened in column-major order
 template <typename T>
-void point_edge_closest_point_3D_hessian(
+IPC_TOOLKIT_HOST_DEVICE void point_edge_closest_point_3D_hessian(
     T p_x,
     T p_y,
     T p_z,
@@ -268,7 +268,7 @@ void point_edge_closest_point_3D_hessian(
 
 // J is (2×12) flattened in column-major order
 template <typename T>
-void edge_edge_closest_point_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_jacobian(
     T ea0_x,
     T ea0_y,
     T ea0_z,
@@ -524,7 +524,7 @@ void edge_edge_closest_point_jacobian(
 
 // hess is (144×1) flattened in column-major order
 template <typename T>
-void edge_edge_closest_point_hessian_a(
+IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_hessian_a(
     T ea0_x,
     T ea0_y,
     T ea0_z,
@@ -1369,7 +1369,7 @@ void edge_edge_closest_point_hessian_a(
 
 // hess is (144×1) flattened in column-major order
 template <typename T>
-void edge_edge_closest_point_hessian_b(
+IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_hessian_b(
     T ea0_x,
     T ea0_y,
     T ea0_z,
@@ -2259,7 +2259,7 @@ void edge_edge_closest_point_hessian_b(
 
 // J is (2×12) flattened in column-major order
 template <typename T>
-void point_triangle_closest_point_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_jacobian(
     T p_x,
     T p_y,
     T p_z,
@@ -2524,7 +2524,7 @@ void point_triangle_closest_point_jacobian(
 
 // hess is (144×1) flattened in column-major order
 template <typename T>
-void point_triangle_closest_point_hessian_0(
+IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_hessian_0(
     T p_x,
     T p_y,
     T p_z,
@@ -3445,7 +3445,7 @@ void point_triangle_closest_point_hessian_0(
 
 // hess is (144×1) flattened in column-major order
 template <typename T>
-void point_triangle_closest_point_hessian_1(
+IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_hessian_1(
     T p_x,
     T p_y,
     T p_z,
@@ -4387,8 +4387,10 @@ void point_triangle_closest_point_hessian_1(
     template void point_triangle_closest_point_hessian_1<T>(                   \
         T, T, T, T, T, T, T, T, T, T, T, T, T[144])
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_CLOSEST_POINT_AUTOGEN(float);
 IPC_INSTANTIATE_CLOSEST_POINT_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_CLOSEST_POINT_AUTOGEN(SimdBatch<float>);
 IPC_INSTANTIATE_CLOSEST_POINT_AUTOGEN(SimdBatch<double>);
diff --git a/src/ipc/tangent/closest_point.hpp b/src/ipc/tangent/closest_point.hpp
index 9f5c4603d..35edc707b 100644
--- a/src/ipc/tangent/closest_point.hpp
+++ b/src/ipc/tangent/closest_point.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 #include <ipc/utils/simd.hpp>
 
@@ -13,12 +14,12 @@ namespace ipc {
 namespace autogen {
     /// hess is (6×6) flattened in column-major order
     template <typename T>
-    void point_edge_closest_point_2D_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void point_edge_closest_point_2D_hessian(
         T p_x, T p_y, T e0_x, T e0_y, T e1_x, T e1_y, T hess[36]);
 
     /// hess is (9×9) flattened in column-major order
     template <typename T>
-    void point_edge_closest_point_3D_hessian(
+    IPC_TOOLKIT_HOST_DEVICE void point_edge_closest_point_3D_hessian(
         T p_x,
         T p_y,
         T p_z,
@@ -32,7 +33,7 @@ namespace autogen {
 
     /// J is (2×12) flattened in column-major order
     template <typename T>
-    void edge_edge_closest_point_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_jacobian(
         T ea0_x,
         T ea0_y,
         T ea0_z,
@@ -49,7 +50,7 @@ namespace autogen {
 
     /// hess is (144×1) flattened in column-major order
     template <typename T>
-    void edge_edge_closest_point_hessian_a(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_hessian_a(
         T ea0_x,
         T ea0_y,
         T ea0_z,
@@ -66,7 +67,7 @@ namespace autogen {
 
     /// hess is (144×1) flattened in column-major order
     template <typename T>
-    void edge_edge_closest_point_hessian_b(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_closest_point_hessian_b(
         T ea0_x,
         T ea0_y,
         T ea0_z,
@@ -83,7 +84,7 @@ namespace autogen {
 
     /// J is (2×12) flattened in column-major order
     template <typename T>
-    void point_triangle_closest_point_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_jacobian(
         T p_x,
         T p_y,
         T p_z,
@@ -100,7 +101,7 @@ namespace autogen {
 
     /// hess is (144×1) flattened in column-major order
     template <typename T>
-    void point_triangle_closest_point_hessian_0(
+    IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_hessian_0(
         T p_x,
         T p_y,
         T p_z,
@@ -117,7 +118,7 @@ namespace autogen {
 
     /// hess is (144×1) flattened in column-major order
     template <typename T>
-    void point_triangle_closest_point_hessian_1(
+    IPC_TOOLKIT_HOST_DEVICE void point_triangle_closest_point_hessian_1(
         T p_x,
         T p_y,
         T p_z,
@@ -172,7 +173,7 @@ namespace detail {
     /// @param b The right-hand side vector.
     /// @return The solution vector `x`.
     template <typename T>
-    inline Eigen::Vector2<T> solve_spd_2x2(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2<T> solve_spd_2x2(
         Eigen::ConstRef<Eigen::Matrix2<T>> A,
         Eigen::ConstRef<Eigen::Vector2<T>> b)
     {
@@ -215,7 +216,7 @@ namespace detail {
     /// @param e1 Second edge point
     /// @return Barycentric coordinate of the closest point
     template <typename T, int dim>
-    inline T point_edge_closest_point(
+    IPC_TOOLKIT_HOST_DEVICE inline T point_edge_closest_point(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
@@ -233,7 +234,8 @@ namespace detail {
     /// @param e1 Second edge point
     /// @return Jacobian of the closest point
     template <typename T, int dim>
-    inline Eigen::Vector<T, 3 * dim> point_edge_closest_point_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 3 * dim>
+    point_edge_closest_point_jacobian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
@@ -261,7 +263,8 @@ namespace detail {
     /// @param e1 Second edge point
     /// @return Hessian of the closest point
     template <typename T, int dim>
-    inline Eigen::Matrix<T, 3 * dim, 3 * dim> point_edge_closest_point_hessian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3 * dim, 3 * dim>
+    point_edge_closest_point_hessian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
@@ -291,7 +294,7 @@ namespace detail {
     /// @param eb1 Second point of the second edge
     /// @return Barycentric coordinates of the closest points
     template <typename T>
-    inline Eigen::Vector2<T> edge_edge_closest_point(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2<T> edge_edge_closest_point(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -321,7 +324,8 @@ namespace detail {
     /// @param eb1 Second point of the second edge
     /// @return Jacobian of the closest points
     template <typename T>
-    inline Eigen::Matrix<T, 2, 12> edge_edge_closest_point_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 2, 12>
+    edge_edge_closest_point_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -342,7 +346,7 @@ namespace detail {
     /// @param eb1 Second point of the second edge
     /// @return Hessian of the closest points (2x12x12 tensor)
     template <typename T>
-    inline std::array<Eigen::Matrix<T, 12, 12>, 2>
+    IPC_TOOLKIT_HOST_DEVICE inline std::array<Eigen::Matrix<T, 12, 12>, 2>
     edge_edge_closest_point_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
@@ -370,7 +374,8 @@ namespace detail {
     /// @param t2 Triangle's third vertex
     /// @return Barycentric coordinates of the closest point
     template <typename T>
-    inline Eigen::Vector2<T> point_triangle_closest_point(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2<T>
+    point_triangle_closest_point(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -392,7 +397,8 @@ namespace detail {
     /// @param t2 Triangle's third vertex
     /// @return Jacobian of the closest point
     template <typename T>
-    inline Eigen::Matrix<T, 2, 12> point_triangle_closest_point_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 2, 12>
+    point_triangle_closest_point_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -413,7 +419,7 @@ namespace detail {
     /// @param t2 Triangle's third vertex
     /// @return Hessian of the closest point (2x12x12 tensor)
     template <typename T>
-    inline std::array<Eigen::Matrix<T, 12, 12>, 2>
+    IPC_TOOLKIT_HOST_DEVICE inline std::array<Eigen::Matrix<T, 12, 12>, 2>
     point_triangle_closest_point_hessian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
@@ -440,7 +446,7 @@ namespace detail {
 /// @param e1 Second edge point
 /// @return barycentric coordinates of the closest point
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_closest_point(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -465,7 +471,7 @@ inline auto point_edge_closest_point(
 /// @param e1 Second edge point
 /// @return Jacobian of the closest point
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_closest_point_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -492,7 +498,7 @@ inline auto point_edge_closest_point_jacobian(
 /// @param e1 Second edge point
 /// @return Hessian of the closest point
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_closest_point_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -535,7 +541,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_closest_point(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_closest_point(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -559,7 +565,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_closest_point_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_closest_point_jacobian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -583,7 +589,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_closest_point_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_closest_point_hessian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -610,7 +616,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_closest_point(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_closest_point(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -634,7 +640,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_closest_point_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_closest_point_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -658,7 +664,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_closest_point_hessian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_closest_point_hessian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
diff --git a/src/ipc/tangent/relative_velocity.hpp b/src/ipc/tangent/relative_velocity.hpp
index b540ea852..489bab855 100644
--- a/src/ipc/tangent/relative_velocity.hpp
+++ b/src/ipc/tangent/relative_velocity.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
 #include <cassert>
@@ -17,7 +18,8 @@ namespace detail {
     /// @param dp1 Velocity of the second point
     /// @return The relative velocity of the two points
     template <typename T, int dim>
-    inline Eigen::Vector<T, dim> point_point_relative_velocity(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, dim>
+    point_point_relative_velocity(
         Eigen::ConstRef<Eigen::Vector<T, dim>> dp0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> dp1)
     {
@@ -31,7 +33,7 @@ namespace detail {
     /// @tparam dim The dimension (2 or 3).
     /// @return The relative velocity Jacobian du/dx
     template <typename T, int dim>
-    inline Eigen::Matrix<T, dim, 2 * dim>
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, dim, 2 * dim>
     point_point_relative_velocity_jacobian()
     {
         static_assert(
@@ -47,7 +49,7 @@ namespace detail {
     /// @tparam dim The dimension (2 or 3).
     /// @return The vectorized tensor of d²u/dxdβ
     template <typename T, int dim>
-    inline Eigen::Vector<T, 2 * dim * dim>
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 2 * dim * dim>
     point_point_relative_velocity_dx_dbeta()
     {
         static_assert(
@@ -68,7 +70,8 @@ namespace detail {
     /// @param alpha Parametric coordinate of the closest point on the edge
     /// @return The relative velocity of the point and the edge
     template <typename T, int dim>
-    inline Eigen::Vector<T, dim> point_edge_relative_velocity(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, dim>
+    point_edge_relative_velocity(
         Eigen::ConstRef<Eigen::Vector<T, dim>> dp,
         Eigen::ConstRef<Eigen::Vector<T, dim>> de0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> de1,
@@ -86,7 +89,7 @@ namespace detail {
     /// @param alpha Parametric coordinate of the closest point on the edge
     /// @return The relative velocity Jacobian du/dx
     template <typename T, int dim>
-    inline Eigen::Matrix<T, dim, 3 * dim>
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, dim, 3 * dim>
     point_edge_relative_velocity_jacobian(const T alpha)
     {
         static_assert(
@@ -106,7 +109,7 @@ namespace detail {
     /// @param alpha Parametric coordinate of the closest point on the edge
     /// @return The vectorized tensor of d²u/dxdα
     template <typename T, int dim>
-    inline Eigen::Vector<T, 3 * dim * dim>
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector<T, 3 * dim * dim>
     point_edge_relative_velocity_dx_dbeta(const T alpha)
     {
         static_assert(
@@ -134,7 +137,8 @@ namespace detail {
     /// @param coords Two parametric coordinates of the closest points
     /// @return The relative velocity of the edges
     template <typename T>
-    inline Eigen::Vector3<T> edge_edge_relative_velocity(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T>
+    edge_edge_relative_velocity(
         Eigen::ConstRef<Eigen::Vector3<T>> dea0,
         Eigen::ConstRef<Eigen::Vector3<T>> dea1,
         Eigen::ConstRef<Eigen::Vector3<T>> deb0,
@@ -151,7 +155,8 @@ namespace detail {
     /// @param coords Two parametric coordinates of the closest points
     /// @return The relative velocity Jacobian du/dx
     template <typename T>
-    inline Eigen::Matrix<T, 3, 12> edge_edge_relative_velocity_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 12>
+    edge_edge_relative_velocity_jacobian(
         Eigen::ConstRef<Eigen::Vector2<T>> coords)
     {
         Eigen::Matrix<T, 3, 12> J = Eigen::Matrix<T, 3, 12>::Zero();
@@ -167,7 +172,8 @@ namespace detail {
     /// @param coords Two parametric coordinates of the closest points
     /// @return The vectorized tensor of d²u/dxdβ
     template <typename T>
-    inline Eigen::Matrix<T, 36, 2> edge_edge_relative_velocity_dx_dbeta(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 36, 2>
+    edge_edge_relative_velocity_dx_dbeta(
         Eigen::ConstRef<Eigen::Vector2<T>> coords)
     {
         constexpr int dim = 3;
@@ -195,7 +201,8 @@ namespace detail {
     /// @param coords Baricentric coordinates of the closest point
     /// @return The relative velocity of the point to the triangle
     template <typename T>
-    inline Eigen::Vector3<T> point_triangle_relative_velocity(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3<T>
+    point_triangle_relative_velocity(
         Eigen::ConstRef<Eigen::Vector3<T>> dp,
         Eigen::ConstRef<Eigen::Vector3<T>> dt0,
         Eigen::ConstRef<Eigen::Vector3<T>> dt1,
@@ -213,7 +220,8 @@ namespace detail {
     /// @param coords Barycentric coordinates of the closest point
     /// @return The relative velocity Jacobian du/dx
     template <typename T>
-    inline Eigen::Matrix<T, 3, 12> point_triangle_relative_velocity_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 12>
+    point_triangle_relative_velocity_jacobian(
         Eigen::ConstRef<Eigen::Vector2<T>> coords)
     {
         Eigen::Matrix<T, 3, 12> J = Eigen::Matrix<T, 3, 12>::Zero();
@@ -231,7 +239,8 @@ namespace detail {
     /// @param coords Baricentric coordinates of the closest point
     /// @return The vectorized tensor of d²u/dxdβ
     template <typename T>
-    inline Eigen::Matrix<T, 36, 2> point_triangle_relative_velocity_dx_dbeta(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 36, 2>
+    point_triangle_relative_velocity_dx_dbeta(
         Eigen::ConstRef<Eigen::Vector2<T>> coords)
     {
         constexpr int dim = 3;
@@ -256,7 +265,7 @@ namespace detail {
 /// @param dp1 Velocity of the second point
 /// @return The relative velocity of the two points
 template <typename DerivedDP0, typename DerivedDP1>
-inline auto point_point_relative_velocity(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_relative_velocity(
     const Eigen::MatrixBase<DerivedDP0>& dp0,
     const Eigen::MatrixBase<DerivedDP1>& dp1)
 {
@@ -282,7 +291,8 @@ inline auto point_point_relative_velocity(
 /// @param dim Dimension (2 or 3)
 /// @return The relative velocity Jacobian du/dx
 template <typename T = double>
-inline MatrixMax<T, 3, 6> point_point_relative_velocity_jacobian(const int dim)
+IPC_TOOLKIT_HOST_DEVICE inline MatrixMax<T, 3, 6>
+point_point_relative_velocity_jacobian(const int dim)
 {
     if (dim == 2) {
         return MatrixMax<T, 3, 6>(
@@ -299,7 +309,8 @@ inline MatrixMax<T, 3, 6> point_point_relative_velocity_jacobian(const int dim)
 /// @param dim Dimension (2 or 3)
 /// @return The vectorized tensor of d²u/dxdβ
 template <typename T = double>
-inline VectorMax<T, 18> point_point_relative_velocity_dx_dbeta(const int dim)
+IPC_TOOLKIT_HOST_DEVICE inline VectorMax<T, 18>
+point_point_relative_velocity_dx_dbeta(const int dim)
 {
     if (dim == 2) {
         return VectorMax<T, 18>(
@@ -321,7 +332,7 @@ inline VectorMax<T, 18> point_point_relative_velocity_dx_dbeta(const int dim)
 /// @param alpha Parametric coordinate of the closest point on the edge
 /// @return The relative velocity of the point and the edge
 template <typename DerivedDP, typename DerivedDE0, typename DerivedDE1>
-inline auto point_edge_relative_velocity(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_relative_velocity(
     const Eigen::MatrixBase<DerivedDP>& dp,
     const Eigen::MatrixBase<DerivedDE0>& de0,
     const Eigen::MatrixBase<DerivedDE1>& de1,
@@ -350,7 +361,7 @@ inline auto point_edge_relative_velocity(
 /// @param alpha Parametric coordinate of the closest point on the edge
 /// @return The relative velocity Jacobian du/dx
 template <typename T = double>
-inline MatrixMax<T, 3, 9>
+IPC_TOOLKIT_HOST_DEVICE inline MatrixMax<T, 3, 9>
 point_edge_relative_velocity_jacobian(const int dim, const T alpha)
 {
     if (dim == 2) {
@@ -375,7 +386,7 @@ point_edge_relative_velocity_jacobian(const int dim, const T alpha)
 /// @param alpha Parametric coordinate of the closest point on the edge
 /// @return The vectorized tensor of d²u/dxdα
 template <typename T = double>
-inline VectorMax<T, 27>
+IPC_TOOLKIT_HOST_DEVICE inline VectorMax<T, 27>
 point_edge_relative_velocity_dx_dbeta(const int dim, const T alpha)
 {
     if (dim == 2) {
@@ -404,7 +415,7 @@ template <
     typename DerivedDEB0,
     typename DerivedDEB1,
     typename DerivedCoords>
-inline auto edge_edge_relative_velocity(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_relative_velocity(
     const Eigen::MatrixBase<DerivedDEA0>& dea0,
     const Eigen::MatrixBase<DerivedDEA1>& dea1,
     const Eigen::MatrixBase<DerivedDEB0>& deb0,
@@ -422,7 +433,7 @@ inline auto edge_edge_relative_velocity(
 /// @param coords Two parametric coordinates of the closest points on the edges
 /// @return The relative velocity Jacobian du/dx
 template <typename DerivedCoords>
-inline auto edge_edge_relative_velocity_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_relative_velocity_jacobian(
     const Eigen::MatrixBase<DerivedCoords>& coords)
 {
     using T = typename DerivedCoords::Scalar;
@@ -441,7 +452,7 @@ inline auto edge_edge_relative_velocity_jacobian(
 /// @param coords Two parametric coordinates of the closest points on the edges
 /// @return The vectorized tensor of d²u/dxdβ
 template <typename DerivedCoords>
-inline auto edge_edge_relative_velocity_dx_dbeta(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_relative_velocity_dx_dbeta(
     const Eigen::MatrixBase<DerivedCoords>& coords)
 {
     using T = typename DerivedCoords::Scalar;
@@ -464,7 +475,7 @@ template <
     typename DerivedDT1,
     typename DerivedDT2,
     typename DerivedCoords>
-inline auto point_triangle_relative_velocity(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_relative_velocity(
     const Eigen::MatrixBase<DerivedDP>& dp,
     const Eigen::MatrixBase<DerivedDT0>& dt0,
     const Eigen::MatrixBase<DerivedDT1>& dt1,
@@ -480,7 +491,7 @@ inline auto point_triangle_relative_velocity(
 /// @param coords Barycentric coordinates of the closest point on the triangle
 /// @return The relative velocity Jacobian du/dx
 template <typename DerivedCoords>
-inline auto point_triangle_relative_velocity_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_relative_velocity_jacobian(
     const Eigen::MatrixBase<DerivedCoords>& coords)
 {
     using T = typename DerivedCoords::Scalar;
@@ -499,7 +510,7 @@ inline auto point_triangle_relative_velocity_jacobian(
 /// @param coords Baricentric coordinates of the closest point on the triangle
 /// @return The vectorized tensor of d²u/dxdβ
 template <typename DerivedCoords>
-inline auto point_triangle_relative_velocity_dx_dbeta(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_relative_velocity_dx_dbeta(
     const Eigen::MatrixBase<DerivedCoords>& coords)
 {
     using T = typename DerivedCoords::Scalar;
diff --git a/src/ipc/tangent/tangent_basis.cpp b/src/ipc/tangent/tangent_basis.cpp
index 182bba771..9de6ed245 100644
--- a/src/ipc/tangent/tangent_basis.cpp
+++ b/src/ipc/tangent/tangent_basis.cpp
@@ -10,7 +10,8 @@ namespace ipc::detail {
 // Point - Point
 
 template <typename T, int dim>
-Eigen::Matrix<T, dim*(dim - 1), 2 * dim> point_point_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, dim*(dim - 1), 2 * dim>
+point_point_tangent_basis_jacobian(
     Eigen::ConstRef<Eigen::Vector<T, dim>> p0,
     Eigen::ConstRef<Eigen::Vector<T, dim>> p1)
 {
@@ -32,7 +33,8 @@ Eigen::Matrix<T, dim*(dim - 1), 2 * dim> point_point_tangent_basis_jacobian(
 // Point - Edge
 
 template <typename T, int dim>
-Eigen::Matrix<T, dim*(dim - 1), 3 * dim> point_edge_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, dim*(dim - 1), 3 * dim>
+point_edge_tangent_basis_jacobian(
     Eigen::ConstRef<Eigen::Vector<T, dim>> p,
     Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
     Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
@@ -56,7 +58,8 @@ Eigen::Matrix<T, dim*(dim - 1), 3 * dim> point_edge_tangent_basis_jacobian(
 // Edge - Edge
 
 template <typename T>
-Eigen::Matrix<T, 6, 12> edge_edge_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 12>
+edge_edge_tangent_basis_jacobian(
     Eigen::ConstRef<Eigen::Vector3<T>> ea0,
     Eigen::ConstRef<Eigen::Vector3<T>> ea1,
     Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -73,7 +76,8 @@ Eigen::Matrix<T, 6, 12> edge_edge_tangent_basis_jacobian(
 // Point - Triangle
 
 template <typename T>
-Eigen::Matrix<T, 6, 12> point_triangle_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 12>
+point_triangle_tangent_basis_jacobian(
     Eigen::ConstRef<Eigen::Vector3<T>> p,
     Eigen::ConstRef<Eigen::Vector3<T>> t0,
     Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -107,10 +111,12 @@ Eigen::Matrix<T, 6, 12> point_triangle_tangent_basis_jacobian(
         Eigen::ConstRef<Eigen::Vector<T, dim>>,                                \
         Eigen::ConstRef<Eigen::Vector<T, dim>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_TANGENT_BASIS_ND(float, 2);
 IPC_INSTANTIATE_TANGENT_BASIS_ND(float, 3);
 IPC_INSTANTIATE_TANGENT_BASIS_ND(double, 2);
 IPC_INSTANTIATE_TANGENT_BASIS_ND(double, 3);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_TANGENT_BASIS_ND(SimdBatch<float>, 2);
 IPC_INSTANTIATE_TANGENT_BASIS_ND(SimdBatch<float>, 3);
@@ -142,8 +148,10 @@ IPC_INSTANTIATE_TANGENT_BASIS_ND(SimdBatch<double>, 3);
         Eigen::ConstRef<Eigen::Vector3<T>>,                                    \
         Eigen::ConstRef<Eigen::Vector3<T>>)
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_TANGENT_BASIS_3D(float);
 IPC_INSTANTIATE_TANGENT_BASIS_3D(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_TANGENT_BASIS_3D(SimdBatch<float>);
 IPC_INSTANTIATE_TANGENT_BASIS_3D(SimdBatch<double>);
@@ -159,7 +167,7 @@ namespace {
     /// @brief Compute the power of 1.5 of a number.
     /// @param x Number to compute the power of 1.5
     /// @return x^(1.5)
-    template <typename T> inline T pow_1_5(T x)
+    template <typename T> IPC_TOOLKIT_HOST_DEVICE inline T pow_1_5(T x)
     {
         return x * ipc::numext::sqrt(x);
     }
@@ -167,8 +175,8 @@ namespace {
 
 // J is (2×4) flattened in column-major order
 template <typename T>
-void point_point_tangent_basis_2D_jacobian(
-    T p0_x, T p0_y, T p1_x, T p1_y, T J[8])
+IPC_TOOLKIT_HOST_DEVICE void
+point_point_tangent_basis_2D_jacobian(T p0_x, T p0_y, T p1_x, T p1_y, T J[8])
 {
     const T t0 = p0_x - p1_x;
     const T t1 = p0_y - p1_y;
@@ -193,7 +201,7 @@ void point_point_tangent_basis_2D_jacobian(
 
 // J is (6×6) flattened in column-major order
 template <typename T>
-void point_point_tangent_basis_3D_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void point_point_tangent_basis_3D_jacobian(
     T p0_x, T p0_y, T p0_z, T p1_x, T p1_y, T p1_z, T J[36])
 {
     const T t0 = p0_x - p1_x;
@@ -330,7 +338,7 @@ void point_point_tangent_basis_3D_jacobian(
 
 // J is (2×6) flattened in column-major order
 template <typename T>
-void point_edge_tangent_basis_2D_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void point_edge_tangent_basis_2D_jacobian(
     T p_x, T p_y, T e0_x, T e0_y, T e1_x, T e1_y, T J[12])
 {
     const T t0 = e0_x - e1_x;
@@ -360,7 +368,7 @@ void point_edge_tangent_basis_2D_jacobian(
 
 // J is (6×9) flattened in column-major order
 template <typename T>
-void point_edge_tangent_basis_3D_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void point_edge_tangent_basis_3D_jacobian(
     T p_x,
     T p_y,
     T p_z,
@@ -507,7 +515,7 @@ void point_edge_tangent_basis_3D_jacobian(
 
 // J is (6×12) flattened in column-major order
 template <typename T>
-void edge_edge_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void edge_edge_tangent_basis_jacobian(
     T ea0_x,
     T ea0_y,
     T ea0_z,
@@ -734,7 +742,7 @@ void edge_edge_tangent_basis_jacobian(
 
 // J is (6×12) flattened in column-major order
 template <typename T>
-void point_triangle_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE void point_triangle_tangent_basis_jacobian(
     T p_x,
     T p_y,
     T p_z,
@@ -958,8 +966,10 @@ void point_triangle_tangent_basis_jacobian(
     template void point_triangle_tangent_basis_jacobian<T>(                    \
         T, T, T, T, T, T, T, T, T, T, T, T, T[72])
 
+#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS
 IPC_INSTANTIATE_TANGENT_BASIS_AUTOGEN(float);
 IPC_INSTANTIATE_TANGENT_BASIS_AUTOGEN(double);
+#endif
 #ifdef IPC_TOOLKIT_WITH_SIMD
 IPC_INSTANTIATE_TANGENT_BASIS_AUTOGEN(SimdBatch<float>);
 IPC_INSTANTIATE_TANGENT_BASIS_AUTOGEN(SimdBatch<double>);
diff --git a/src/ipc/tangent/tangent_basis.hpp b/src/ipc/tangent/tangent_basis.hpp
index 14ba361ac..f7ae3541a 100644
--- a/src/ipc/tangent/tangent_basis.hpp
+++ b/src/ipc/tangent/tangent_basis.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <ipc/config.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 #include <ipc/utils/simd.hpp>
 
@@ -18,7 +19,8 @@ namespace detail {
     /// @param p1 Second point
     /// @return A dim×(dim-1) matrix whose columns are the basis vectors.
     template <typename T, int dim>
-    inline Eigen::Matrix<T, dim, dim - 1> point_point_tangent_basis(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, dim, dim - 1>
+    point_point_tangent_basis(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> p1)
     {
@@ -72,7 +74,8 @@ namespace detail {
     /// @param p1 Second point
     /// @return A (dim*(dim-1))×(2*dim) matrix.
     template <typename T, int dim>
-    Eigen::Matrix<T, dim*(dim - 1), 2 * dim> point_point_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, dim*(dim - 1), 2 * dim>
+    point_point_tangent_basis_jacobian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> p1);
 
@@ -87,7 +90,8 @@ namespace detail {
     /// @param e1 Second edge point
     /// @return A dim×(dim-1) matrix whose columns are the basis vectors.
     template <typename T, int dim>
-    inline Eigen::Matrix<T, dim, dim - 1> point_edge_tangent_basis(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, dim, dim - 1>
+    point_edge_tangent_basis(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1)
@@ -115,7 +119,8 @@ namespace detail {
     /// @param e1 Second edge point
     /// @return A (dim*(dim-1))×(3*dim) matrix.
     template <typename T, int dim>
-    Eigen::Matrix<T, dim*(dim - 1), 3 * dim> point_edge_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, dim*(dim - 1), 3 * dim>
+    point_edge_tangent_basis_jacobian(
         Eigen::ConstRef<Eigen::Vector<T, dim>> p,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e0,
         Eigen::ConstRef<Eigen::Vector<T, dim>> e1);
@@ -131,7 +136,8 @@ namespace detail {
     /// @param eb1 Second point of the second edge
     /// @return A 3x2 matrix whose columns are the basis vectors.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 2> edge_edge_tangent_basis(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 2>
+    edge_edge_tangent_basis(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -161,7 +167,8 @@ namespace detail {
     /// @param eb1 Second point of the second edge
     /// @return A (3*2)x12 matrix whose columns are the basis vectors.
     template <typename T>
-    Eigen::Matrix<T, 6, 12> edge_edge_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 12>
+    edge_edge_tangent_basis_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> ea0,
         Eigen::ConstRef<Eigen::Vector3<T>> ea1,
         Eigen::ConstRef<Eigen::Vector3<T>> eb0,
@@ -178,7 +185,8 @@ namespace detail {
     /// @param t2 Triangle's third vertex
     /// @return A 3x2 matrix whose columns are the basis vectors.
     template <typename T>
-    inline Eigen::Matrix<T, 3, 2> point_triangle_tangent_basis(
+    IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix<T, 3, 2>
+    point_triangle_tangent_basis(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -209,7 +217,8 @@ namespace detail {
     /// @param t2 Triangle's third vertex
     /// @return A (3*2)x12 matrix whose columns are the basis vectors.
     template <typename T>
-    Eigen::Matrix<T, 6, 12> point_triangle_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix<T, 6, 12>
+    point_triangle_tangent_basis_jacobian(
         Eigen::ConstRef<Eigen::Vector3<T>> p,
         Eigen::ConstRef<Eigen::Vector3<T>> t0,
         Eigen::ConstRef<Eigen::Vector3<T>> t1,
@@ -224,7 +233,7 @@ namespace detail {
 /// @param p1 Second point
 /// @return A 3x2 matrix whose columns are the basis vectors.
 template <typename DerivedP0, typename DerivedP1>
-inline auto point_point_tangent_basis(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_tangent_basis(
     const Eigen::MatrixBase<DerivedP0>& p0,
     const Eigen::MatrixBase<DerivedP1>& p1)
 {
@@ -250,7 +259,7 @@ inline auto point_point_tangent_basis(
 /// @param p1 Second point
 /// @return A (3*2)x6 matrix whose columns are the basis vectors.
 template <typename DerivedP0, typename DerivedP1>
-inline auto point_point_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_point_tangent_basis_jacobian(
     const Eigen::MatrixBase<DerivedP0>& p0,
     const Eigen::MatrixBase<DerivedP1>& p1)
 {
@@ -280,7 +289,7 @@ inline auto point_point_tangent_basis_jacobian(
 /// @param e1 Second edge point
 /// @return A 3x2 matrix whose columns are the basis vectors.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_tangent_basis(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_tangent_basis(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -308,7 +317,7 @@ inline auto point_edge_tangent_basis(
 /// @param e1 Second edge point
 /// @return A (3*2)x9 matrix whose columns are the basis vectors.
 template <typename DerivedP, typename DerivedE0, typename DerivedE1>
-inline auto point_edge_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_tangent_basis_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedE0>& e0,
     const Eigen::MatrixBase<DerivedE1>& e1)
@@ -344,7 +353,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_tangent_basis(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_tangent_basis(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -367,7 +376,7 @@ template <
     typename DerivedEA1,
     typename DerivedEB0,
     typename DerivedEB1>
-inline auto edge_edge_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_tangent_basis_jacobian(
     const Eigen::MatrixBase<DerivedEA0>& ea0,
     const Eigen::MatrixBase<DerivedEA1>& ea1,
     const Eigen::MatrixBase<DerivedEB0>& eb0,
@@ -399,7 +408,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_tangent_basis(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_tangent_basis(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -420,7 +429,7 @@ template <
     typename DerivedT0,
     typename DerivedT1,
     typename DerivedT2>
-inline auto point_triangle_tangent_basis_jacobian(
+IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_tangent_basis_jacobian(
     const Eigen::MatrixBase<DerivedP>& p,
     const Eigen::MatrixBase<DerivedT0>& t0,
     const Eigen::MatrixBase<DerivedT1>& t1,
@@ -437,22 +446,22 @@ namespace autogen {
 
     // J is (2×4) flattened in column-major order
     template <typename T>
-    void point_point_tangent_basis_2D_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_point_tangent_basis_2D_jacobian(
         T p0_x, T p0_y, T p1_x, T p1_y, T J[8]);
 
     // J is (6×6) flattened in column-major order
     template <typename T>
-    void point_point_tangent_basis_3D_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_point_tangent_basis_3D_jacobian(
         T p0_x, T p0_y, T p0_z, T p1_x, T p1_y, T p1_z, T J[36]);
 
     // J is (2×6) flattened in column-major order
     template <typename T>
-    void point_edge_tangent_basis_2D_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_edge_tangent_basis_2D_jacobian(
         T p_x, T p_y, T e0_x, T e0_y, T e1_x, T e1_y, T J[12]);
 
     // J is (6×9) flattened in column-major order
     template <typename T>
-    void point_edge_tangent_basis_3D_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_edge_tangent_basis_3D_jacobian(
         T p_x,
         T p_y,
         T p_z,
@@ -466,7 +475,7 @@ namespace autogen {
 
     // J is (6×12) flattened in column-major order
     template <typename T>
-    void edge_edge_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void edge_edge_tangent_basis_jacobian(
         T ea0_x,
         T ea0_y,
         T ea0_z,
@@ -483,7 +492,7 @@ namespace autogen {
 
     // J is (6×12) flattened in column-major order
     template <typename T>
-    void point_triangle_tangent_basis_jacobian(
+    IPC_TOOLKIT_HOST_DEVICE void point_triangle_tangent_basis_jacobian(
         T p_x,
         T p_y,
         T p_z,
diff --git a/src/ipc/utils/simd.hpp b/src/ipc/utils/simd.hpp
index e4ad700f5..72c099635 100644
--- a/src/ipc/utils/simd.hpp
+++ b/src/ipc/utils/simd.hpp
@@ -40,7 +40,7 @@ template <typename T> using scalar_of_t = typename ScalarOf<T>::type;
 /// constructor call receiving a `double`, an implicit narrowing that
 /// `-Wfloat-conversion` reports at every instantiation. For a plain scalar this
 /// is the explicit cast the code would have written anyway.
-template <typename T> inline T literal(const double c)
+template <typename T> IPC_TOOLKIT_HOST_DEVICE inline T literal(const double c)
 {
     return T(static_cast<scalar_of_t<T>>(c));
 }
@@ -54,19 +54,33 @@ template <typename T> inline T literal(const double c)
 /// We overload on the two argument types rather than relying on ADL to find
 /// `xsimd::all`. The tradeoff is a little duplication in exchange for keeping
 /// a name this generic from matching arbitrary types elsewhere in `ipc`.
-inline bool all_of(const bool mask) { return mask; }
+IPC_TOOLKIT_HOST_DEVICE inline bool all_of(const bool mask) { return mask; }
 
 /// @brief Pick between `a` and `b`.
 ///
 /// The scalar counterpart of `xsimd::select`, which ADL finds for a batch
 /// `mask`, so one `select(cond, a, b)` compiles for both.
-template <typename T> inline T select(const bool mask, const T& a, const T& b)
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE inline T select(const bool mask, const T& a, const T& b)
 {
     return mask ? a : b;
 }
 
+/// @brief Clamp `v` to `[lo, hi]`.
+///
+/// Not `std::clamp`, which cannot be called from device code: MSVC's debug STL
+/// checks the bounds with `_STL_VERIFY`, which expands to `__debugbreak()`, and
+/// nvcc's NVVM backend then emits invalid IR ("Terminator found in the middle
+/// of a basic block"). The comparison order matches `std::clamp` exactly, so
+/// the result -- including a NaN `v` passing through unchanged -- is identical.
+template <typename T>
+IPC_TOOLKIT_HOST_DEVICE inline T clamp(const T& v, const T& lo, const T& hi)
+{
+    return v < lo ? lo : (hi < v ? hi : v);
+}
+
 /// @brief `+infinity` for any scalar the library templates on.
-template <typename T> inline T infinity()
+template <typename T> IPC_TOOLKIT_HOST_DEVICE inline T infinity()
 {
     return T(std::numeric_limits<scalar_of_t<T>>::infinity());
 }
@@ -92,13 +106,15 @@ template <typename T> inline T infinity()
 /// comes from the order the blend is folded. Masks may overlap, and only that
 /// order decides the winner, so a test must cover lanes that fall in
 /// overlapping cases.
-template <typename F> inline auto select_lazy(F&& else_value)
+template <typename F>
+IPC_TOOLKIT_HOST_DEVICE inline auto select_lazy(F&& else_value)
 {
     return else_value();
 }
 
 template <typename Mask, typename F, typename... Rest>
-inline auto select_lazy(const Mask& mask, F&& value, Rest&&... rest)
+IPC_TOOLKIT_HOST_DEVICE inline auto
+select_lazy(const Mask& mask, F&& value, Rest&&... rest)
 {
     if constexpr (std::is_same_v<std::decay_t<Mask>, bool>) {
         return mask ? value() : select_lazy(std::forward<Rest>(rest)...);
@@ -198,6 +214,21 @@ inline bool all_of(const xsimd::batch_bool<T, A>& mask)
     return xsimd::all(mask);
 }
 
+/// @brief Clamp each lane of `v` to `[lo, hi]`.
+///
+/// The batch counterpart of the scalar `clamp` above: a batch comparison
+/// answers per-lane rather than with one `bool`, so the `?:` cascade becomes
+/// two `xsimd::select` blends in the same order. A NaN lane matches neither
+/// comparison and passes through, as it does for a scalar.
+template <typename T, typename A>
+inline xsimd::batch<T, A> clamp(
+    const xsimd::batch<T, A>& v,
+    const xsimd::batch<T, A>& lo,
+    const xsimd::batch<T, A>& hi)
+{
+    return xsimd::select(v < lo, lo, xsimd::select(hi < v, hi, v));
+}
+
 } // namespace ipc
 
 #endif
@@ -210,7 +241,7 @@ namespace ipc {
 /// `if (squaredNorm() > 0)`, which a batch cannot answer with one bool. This
 /// applies that same rule per-lane and otherwise defers to Eigen.
 template <typename Derived>
-inline typename Derived::PlainObject
+IPC_TOOLKIT_HOST_DEVICE inline typename Derived::PlainObject
 normalized(const Eigen::MatrixBase<Derived>& v)
 {
     using T = typename Derived::Scalar;
diff --git a/tests/src/tests/CMakeLists.txt b/tests/src/tests/CMakeLists.txt
index cca3ed913..c9c1ea62d 100644
--- a/tests/src/tests/CMakeLists.txt
+++ b/tests/src/tests/CMakeLists.txt
@@ -18,6 +18,11 @@ set(SOURCES
   utils.hpp
 )
 
+# Shared scaffolding for the test_gpu_*.cu parity tests; defines device code.
+if(IPC_TOOLKIT_WITH_CUDA)
+  list(APPEND SOURCES gpu_utils.hpp)
+endif()
+
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
 
 ################################################################################
diff --git a/tests/src/tests/barrier/CMakeLists.txt b/tests/src/tests/barrier/CMakeLists.txt
index 0a9268f2a..e8f70839f 100644
--- a/tests/src/tests/barrier/CMakeLists.txt
+++ b/tests/src/tests/barrier/CMakeLists.txt
@@ -9,6 +9,12 @@ set(SOURCES
   # Utilities
 )
 
+if(IPC_TOOLKIT_WITH_CUDA)
+  list(APPEND SOURCES
+    test_gpu_barrier.cu
+  )
+endif()
+
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
 
 ################################################################################
diff --git a/tests/src/tests/barrier/test_gpu_barrier.cu b/tests/src/tests/barrier/test_gpu_barrier.cu
new file mode 100644
index 000000000..7f2aede1d
--- /dev/null
+++ b/tests/src/tests/barrier/test_gpu_barrier.cu
@@ -0,0 +1,62 @@
+// Exercises the __host__ __device__ barrier free functions from a CUDA
+// kernel.
+//
+// Compilation proves the framework works: nvcc must instantiate the barrier
+// free-function trio as device code. On a machine with a CUDA device the
+// kernel also runs and its outputs are checked against the host
+// implementations — host/device parity is the contract. Without a device
+// (e.g. a GPU-less container) the runtime checks are skipped but the device
+// code has still been compiled.
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_CUDA
+
+#include <ipc/barrier/barrier.hpp>
+
+#include <tests/gpu_utils.hpp>
+
+#include <cuda_runtime.h>
+
+#include <vector>
+
+using namespace ipc::tests;
+
+namespace {
+
+__global__ void barrier_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const double dhat = in[0];
+    // in[1..3]: interior, inactive (d >= dhat), and invalid (d <= 0) samples
+    for (int i = 0; i < 3; ++i) {
+        const double d = in[1 + i];
+        out[3 * i + 0] = ipc::barrier(d, dhat);
+        out[3 * i + 1] = ipc::barrier_first_derivative(d, dhat);
+        out[3 * i + 2] = ipc::barrier_second_derivative(d, dhat);
+    }
+}
+
+} // namespace
+
+TEST_CASE("GPU barrier", "[barrier][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const double dhat = 1e-2;
+    const std::vector<double> in = { dhat, 0.5 * dhat, 2 * dhat, -1.0 };
+
+    std::vector<double> expected;
+    for (int i = 0; i < 3; ++i) {
+        const double d = in[1 + i];
+        expected.push_back(ipc::barrier(d, dhat));
+        expected.push_back(ipc::barrier_first_derivative(d, dhat));
+        expected.push_back(ipc::barrier_second_derivative(d, dhat));
+    }
+
+    check_gpu_matches_host(barrier_kernel, in, expected);
+}
+
+#endif // IPC_TOOLKIT_WITH_CUDA
diff --git a/tests/src/tests/distance/CMakeLists.txt b/tests/src/tests/distance/CMakeLists.txt
index 829389de9..3bd2738b3 100644
--- a/tests/src/tests/distance/CMakeLists.txt
+++ b/tests/src/tests/distance/CMakeLists.txt
@@ -19,6 +19,12 @@ set(SOURCES
   # Utilities
 )
 
+if(IPC_TOOLKIT_WITH_CUDA)
+  list(APPEND SOURCES
+    test_gpu_distance.cu
+  )
+endif()
+
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
 
 ################################################################################
diff --git a/tests/src/tests/distance/test_gpu_distance.cu b/tests/src/tests/distance/test_gpu_distance.cu
new file mode 100644
index 000000000..3f2289b95
--- /dev/null
+++ b/tests/src/tests/distance/test_gpu_distance.cu
@@ -0,0 +1,516 @@
+// Exercises the __host__ __device__ distance functions from CUDA kernels.
+//
+// Compilation proves the framework works: nvcc must instantiate every
+// converted distance function (including the distance-type classification and
+// the signed variants) as device code. On a machine with a CUDA device the
+// kernels also run and their outputs are checked against the host
+// implementations — host/device parity is the contract. Without a device
+// (e.g. a GPU-less container) the runtime checks are skipped but the device
+// code has still been compiled.
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_CUDA
+
+#include <ipc/distance/distance_type.hpp>
+#include <ipc/distance/edge_edge.hpp>
+#include <ipc/distance/edge_edge_mollifier.hpp>
+#include <ipc/distance/line_line.hpp>
+#include <ipc/distance/point_edge.hpp>
+#include <ipc/distance/point_line.hpp>
+#include <ipc/distance/point_plane.hpp>
+#include <ipc/distance/point_point.hpp>
+#include <ipc/distance/point_triangle.hpp>
+#include <ipc/distance/signed/line_line.hpp>
+#include <ipc/distance/signed/point_line.hpp>
+#include <ipc/distance/signed/point_plane.hpp>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <tests/gpu_utils.hpp>
+
+#include <cuda_runtime.h>
+
+#include <vector>
+
+using namespace ipc::tests;
+
+namespace {
+
+// ---------------------------------------------------------------------------
+
+// Shared fixtures (kept in one place so kernels and host references agree).
+// clang-format off
+const std::vector<double> PP_IN  = { 0, 0, 0,  1, 2, 2 };
+const std::vector<double> PL_IN  = { 0.1, 1, 0.2,  -1, 0, 0,  1, 0, 0.1 };
+const std::vector<double> LL_IN  = { -1, 0, 0,  1, 0, 0,  0, 1, -1,  0.1, 1, 1 };
+const std::vector<double> PLANE_IN = { 0, 1, 0,  0, 0, 0,  0, 1, 0,
+                                       -1, 0, 1,  1, 0, 1,  0, 0, -1 };
+const std::vector<double> PE_IN  = { 0.5, 1, 0,  -1, 0, 0,  1, 0, 0 };
+const std::vector<double> EE_IN  = { -1, 0, 0,  1, 0, 0,  -0.5, 1, -1,  0.5, 1, 1 };
+const std::vector<double> PT_IN  = { 0, 1, 0,  -1, 0, 1,  1, 0, 1,  0, 0, -1 };
+const std::vector<double> SGN_IN = { 0.1, 1,  -1, 0,  1, 0 }; // 2D point-line
+// clang-format on
+
+__global__ void point_point_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p0 = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d p1 = Eigen::Vector3d(in[3], in[4], in[5]);
+
+    int offset = 0;
+    out[offset++] = ipc::point_point_distance(p0, p1);
+    offset = write_out(ipc::point_point_distance_gradient(p0, p1), out, offset);
+    offset = write_out(ipc::point_point_distance_hessian(p0, p1), out, offset);
+}
+
+__global__ void point_line_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(in[3], in[4], in[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(in[6], in[7], in[8]);
+
+    int offset = 0;
+    out[offset++] = ipc::point_line_distance(p, e0, e1);
+    offset =
+        write_out(ipc::point_line_distance_gradient(p, e0, e1), out, offset);
+    offset =
+        write_out(ipc::point_line_distance_hessian(p, e0, e1), out, offset);
+}
+
+__global__ void line_line_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d ea0(in[0], in[1], in[2]);
+    const Eigen::Vector3d ea1(in[3], in[4], in[5]);
+    const Eigen::Vector3d eb0(in[6], in[7], in[8]);
+    const Eigen::Vector3d eb1(in[9], in[10], in[11]);
+
+    int offset = 0;
+    out[offset++] = ipc::line_line_distance(ea0, ea1, eb0, eb1);
+    offset = write_out(
+        ipc::line_line_distance_gradient(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::line_line_distance_hessian(ea0, ea1, eb0, eb1), out, offset);
+}
+
+__global__ void point_plane_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d p(in[0], in[1], in[2]);
+    const Eigen::Vector3d origin(in[3], in[4], in[5]);
+    const Eigen::Vector3d normal(in[6], in[7], in[8]);
+    const Eigen::Vector3d t0(in[9], in[10], in[11]);
+    const Eigen::Vector3d t1(in[12], in[13], in[14]);
+    const Eigen::Vector3d t2(in[15], in[16], in[17]);
+
+    int offset = 0;
+    out[offset++] = ipc::point_plane_distance(p, origin, normal);
+    offset = write_out(
+        ipc::point_plane_distance_gradient(p, origin, normal), out, offset);
+    offset = write_out(
+        ipc::point_plane_distance_hessian(p, origin, normal), out, offset);
+    out[offset++] = ipc::point_plane_distance(p, t0, t1, t2);
+    offset = write_out(
+        ipc::point_plane_distance_gradient(p, t0, t1, t2), out, offset);
+    offset = write_out(
+        ipc::point_plane_distance_hessian(p, t0, t1, t2), out, offset);
+}
+
+__global__ void point_edge_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(in[3], in[4], in[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(in[6], in[7], in[8]);
+
+    int offset = 0;
+    out[offset++] =
+        static_cast<double>(ipc::point_edge_distance_type(p, e0, e1));
+    out[offset++] = ipc::point_edge_distance(p, e0, e1);
+    offset =
+        write_out(ipc::point_edge_distance_gradient(p, e0, e1), out, offset);
+    offset =
+        write_out(ipc::point_edge_distance_hessian(p, e0, e1), out, offset);
+}
+
+__global__ void edge_edge_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d ea0(in[0], in[1], in[2]);
+    const Eigen::Vector3d ea1(in[3], in[4], in[5]);
+    const Eigen::Vector3d eb0(in[6], in[7], in[8]);
+    const Eigen::Vector3d eb1(in[9], in[10], in[11]);
+
+    int offset = 0;
+    out[offset++] =
+        static_cast<double>(ipc::edge_edge_distance_type(ea0, ea1, eb0, eb1));
+    out[offset++] = static_cast<double>(
+        ipc::edge_edge_parallel_distance_type(ea0, ea1, eb0, eb1));
+    out[offset++] = ipc::edge_edge_distance(ea0, ea1, eb0, eb1);
+    offset = write_out(
+        ipc::edge_edge_distance_gradient(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::edge_edge_distance_hessian(ea0, ea1, eb0, eb1), out, offset);
+}
+
+__global__ void point_triangle_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d p(in[0], in[1], in[2]);
+    const Eigen::Vector3d t0(in[3], in[4], in[5]);
+    const Eigen::Vector3d t1(in[6], in[7], in[8]);
+    const Eigen::Vector3d t2(in[9], in[10], in[11]);
+
+    int offset = 0;
+    out[offset++] =
+        static_cast<double>(ipc::point_triangle_distance_type(p, t0, t1, t2));
+    out[offset++] = ipc::point_triangle_distance(p, t0, t1, t2);
+    offset = write_out(
+        ipc::point_triangle_distance_gradient(p, t0, t1, t2), out, offset);
+    offset = write_out(
+        ipc::point_triangle_distance_hessian(p, t0, t1, t2), out, offset);
+}
+
+__global__ void mollifier_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d ea0(in[0], in[1], in[2]);
+    const Eigen::Vector3d ea1(in[3], in[4], in[5]);
+    const Eigen::Vector3d eb0(in[6], in[7], in[8]);
+    const Eigen::Vector3d eb1(in[9], in[10], in[11]);
+    const double eps_x = in[12];
+
+    int offset = 0;
+    const double c = ipc::edge_edge_cross_squarednorm(ea0, ea1, eb0, eb1);
+    out[offset++] = c;
+    offset = write_out(
+        ipc::edge_edge_cross_squarednorm_gradient(ea0, ea1, eb0, eb1), out,
+        offset);
+    offset = write_out(
+        ipc::edge_edge_cross_squarednorm_hessian(ea0, ea1, eb0, eb1), out,
+        offset);
+
+    const double x = 0.5 * eps_x; // scalar mollifier in the active range
+    out[offset++] = ipc::edge_edge_mollifier(x, eps_x);
+    out[offset++] = ipc::edge_edge_mollifier_gradient(x, eps_x);
+    out[offset++] = ipc::edge_edge_mollifier_hessian(x, eps_x);
+    out[offset++] = ipc::edge_edge_mollifier_derivative_wrt_eps_x(x, eps_x);
+    out[offset++] =
+        ipc::edge_edge_mollifier_gradient_derivative_wrt_eps_x(x, eps_x);
+
+    out[offset++] = ipc::edge_edge_mollifier(ea0, ea1, eb0, eb1, eps_x);
+    offset = write_out(
+        ipc::edge_edge_mollifier_gradient(ea0, ea1, eb0, eb1, eps_x), out,
+        offset);
+    offset = write_out(
+        ipc::edge_edge_mollifier_hessian(ea0, ea1, eb0, eb1, eps_x), out,
+        offset);
+    offset = write_out(
+        ipc::edge_edge_mollifier_gradient_wrt_x(
+            ea0, ea1, eb0, eb1, ea0, ea1, eb0, eb1),
+        out, offset);
+    offset = write_out(
+        ipc::edge_edge_mollifier_gradient_jacobian_wrt_x(
+            ea0, ea1, eb0, eb1, ea0, ea1, eb0, eb1),
+        out, offset);
+
+    out[offset++] = ipc::edge_edge_mollifier_threshold(ea0, ea1, eb0, eb1);
+    offset = write_out(
+        ipc::edge_edge_mollifier_threshold_gradient(ea0, ea1, eb0, eb1), out,
+        offset);
+}
+
+__global__ void signed_distance_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    // 2D point-line (inputs 0-5)
+    const Eigen::Vector2d p2(in[0], in[1]);
+    const Eigen::Vector2d f0(in[2], in[3]);
+    const Eigen::Vector2d f1(in[4], in[5]);
+    // 3D line-line + point-plane (inputs 6-17 and 18-29)
+    const Eigen::Vector3d ea0(in[6], in[7], in[8]);
+    const Eigen::Vector3d ea1(in[9], in[10], in[11]);
+    const Eigen::Vector3d eb0(in[12], in[13], in[14]);
+    const Eigen::Vector3d eb1(in[15], in[16], in[17]);
+    const Eigen::Vector3d p(in[18], in[19], in[20]);
+    const Eigen::Vector3d t0(in[21], in[22], in[23]);
+    const Eigen::Vector3d t1(in[24], in[25], in[26]);
+    const Eigen::Vector3d t2(in[27], in[28], in[29]);
+
+    int offset = 0;
+    out[offset++] = ipc::point_line_signed_distance(p2, f0, f1);
+    offset = write_out(
+        ipc::point_line_signed_distance_gradient(p2, f0, f1), out, offset);
+    offset = write_out(
+        ipc::point_line_signed_distance_hessian(p2, f0, f1), out, offset);
+
+    out[offset++] = ipc::line_line_signed_distance(ea0, ea1, eb0, eb1);
+    offset = write_out(
+        ipc::line_line_signed_distance_gradient(ea0, ea1, eb0, eb1), out,
+        offset);
+    offset = write_out(
+        ipc::line_line_signed_distance_hessian(ea0, ea1, eb0, eb1), out,
+        offset);
+
+    out[offset++] = ipc::point_plane_signed_distance(p, t0, t1, t2);
+    offset = write_out(
+        ipc::point_plane_signed_distance_gradient(p, t0, t1, t2), out, offset);
+    offset = write_out(
+        ipc::point_plane_signed_distance_hessian(p, t0, t1, t2), out, offset);
+}
+
+} // namespace
+
+TEST_CASE("GPU point-point distance", "[distance][point_point][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d p0 = Eigen::Vector3d(0, 0, 0);
+    const ipc::VectorMax3d p1 = Eigen::Vector3d(1, 2, 2);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::point_point_distance(p0, p1));
+    write_expected(ipc::point_point_distance_gradient(p0, p1), expected);
+    write_expected(ipc::point_point_distance_hessian(p0, p1), expected);
+    REQUIRE(expected[0] == Catch::Approx(9.0)); // sanity-check the fixture
+
+    check_gpu_matches_host(point_point_kernel, PP_IN, expected);
+}
+
+TEST_CASE("GPU point-line distance", "[distance][point_line][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d p = Eigen::Vector3d(PL_IN[0], PL_IN[1], PL_IN[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(PL_IN[3], PL_IN[4], PL_IN[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(PL_IN[6], PL_IN[7], PL_IN[8]);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::point_line_distance(p, e0, e1));
+    write_expected(ipc::point_line_distance_gradient(p, e0, e1), expected);
+    write_expected(ipc::point_line_distance_hessian(p, e0, e1), expected);
+
+    check_gpu_matches_host(point_line_kernel, PL_IN, expected);
+}
+
+TEST_CASE("GPU line-line distance", "[distance][line_line][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d ea0(LL_IN[0], LL_IN[1], LL_IN[2]);
+    const Eigen::Vector3d ea1(LL_IN[3], LL_IN[4], LL_IN[5]);
+    const Eigen::Vector3d eb0(LL_IN[6], LL_IN[7], LL_IN[8]);
+    const Eigen::Vector3d eb1(LL_IN[9], LL_IN[10], LL_IN[11]);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::line_line_distance(ea0, ea1, eb0, eb1));
+    write_expected(
+        ipc::line_line_distance_gradient(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::line_line_distance_hessian(ea0, ea1, eb0, eb1), expected);
+
+    check_gpu_matches_host(line_line_kernel, LL_IN, expected);
+}
+
+TEST_CASE("GPU point-plane distance", "[distance][point_plane][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d p(PLANE_IN[0], PLANE_IN[1], PLANE_IN[2]);
+    const Eigen::Vector3d origin(PLANE_IN[3], PLANE_IN[4], PLANE_IN[5]);
+    const Eigen::Vector3d normal(PLANE_IN[6], PLANE_IN[7], PLANE_IN[8]);
+    const Eigen::Vector3d t0(PLANE_IN[9], PLANE_IN[10], PLANE_IN[11]);
+    const Eigen::Vector3d t1(PLANE_IN[12], PLANE_IN[13], PLANE_IN[14]);
+    const Eigen::Vector3d t2(PLANE_IN[15], PLANE_IN[16], PLANE_IN[17]);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::point_plane_distance(p, origin, normal));
+    write_expected(
+        ipc::point_plane_distance_gradient(p, origin, normal), expected);
+    write_expected(
+        ipc::point_plane_distance_hessian(p, origin, normal), expected);
+    expected.push_back(ipc::point_plane_distance(p, t0, t1, t2));
+    write_expected(ipc::point_plane_distance_gradient(p, t0, t1, t2), expected);
+    write_expected(ipc::point_plane_distance_hessian(p, t0, t1, t2), expected);
+
+    check_gpu_matches_host(point_plane_kernel, PLANE_IN, expected);
+}
+
+TEST_CASE("GPU point-edge distance", "[distance][point_edge][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d p = Eigen::Vector3d(PE_IN[0], PE_IN[1], PE_IN[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(PE_IN[3], PE_IN[4], PE_IN[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(PE_IN[6], PE_IN[7], PE_IN[8]);
+
+    std::vector<double> expected;
+    expected.push_back(
+        static_cast<double>(ipc::point_edge_distance_type(p, e0, e1)));
+    expected.push_back(ipc::point_edge_distance(p, e0, e1));
+    write_expected(ipc::point_edge_distance_gradient(p, e0, e1), expected);
+    write_expected(ipc::point_edge_distance_hessian(p, e0, e1), expected);
+
+    check_gpu_matches_host(point_edge_kernel, PE_IN, expected);
+}
+
+TEST_CASE("GPU edge-edge distance", "[distance][edge_edge][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d ea0(EE_IN[0], EE_IN[1], EE_IN[2]);
+    const Eigen::Vector3d ea1(EE_IN[3], EE_IN[4], EE_IN[5]);
+    const Eigen::Vector3d eb0(EE_IN[6], EE_IN[7], EE_IN[8]);
+    const Eigen::Vector3d eb1(EE_IN[9], EE_IN[10], EE_IN[11]);
+
+    std::vector<double> expected;
+    expected.push_back(
+        static_cast<double>(ipc::edge_edge_distance_type(ea0, ea1, eb0, eb1)));
+    expected.push_back(
+        static_cast<double>(
+            ipc::edge_edge_parallel_distance_type(ea0, ea1, eb0, eb1)));
+    expected.push_back(ipc::edge_edge_distance(ea0, ea1, eb0, eb1));
+    write_expected(
+        ipc::edge_edge_distance_gradient(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::edge_edge_distance_hessian(ea0, ea1, eb0, eb1), expected);
+
+    check_gpu_matches_host(edge_edge_kernel, EE_IN, expected);
+}
+
+TEST_CASE("GPU point-triangle distance", "[distance][point_triangle][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d p(PT_IN[0], PT_IN[1], PT_IN[2]);
+    const Eigen::Vector3d t0(PT_IN[3], PT_IN[4], PT_IN[5]);
+    const Eigen::Vector3d t1(PT_IN[6], PT_IN[7], PT_IN[8]);
+    const Eigen::Vector3d t2(PT_IN[9], PT_IN[10], PT_IN[11]);
+
+    std::vector<double> expected;
+    expected.push_back(
+        static_cast<double>(ipc::point_triangle_distance_type(p, t0, t1, t2)));
+    expected.push_back(ipc::point_triangle_distance(p, t0, t1, t2));
+    write_expected(
+        ipc::point_triangle_distance_gradient(p, t0, t1, t2), expected);
+    write_expected(
+        ipc::point_triangle_distance_hessian(p, t0, t1, t2), expected);
+
+    check_gpu_matches_host(point_triangle_kernel, PT_IN, expected);
+}
+
+TEST_CASE("GPU edge-edge mollifier", "[distance][edge_edge_mollifier][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d ea0(EE_IN[0], EE_IN[1], EE_IN[2]);
+    const Eigen::Vector3d ea1(EE_IN[3], EE_IN[4], EE_IN[5]);
+    const Eigen::Vector3d eb0(EE_IN[6], EE_IN[7], EE_IN[8]);
+    const Eigen::Vector3d eb1(EE_IN[9], EE_IN[10], EE_IN[11]);
+    const double eps_x = 1e-3;
+
+    std::vector<double> in = EE_IN;
+    in.push_back(eps_x);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::edge_edge_cross_squarednorm(ea0, ea1, eb0, eb1));
+    write_expected(
+        ipc::edge_edge_cross_squarednorm_gradient(ea0, ea1, eb0, eb1),
+        expected);
+    write_expected(
+        ipc::edge_edge_cross_squarednorm_hessian(ea0, ea1, eb0, eb1), expected);
+
+    const double x = 0.5 * eps_x;
+    expected.push_back(ipc::edge_edge_mollifier(x, eps_x));
+    expected.push_back(ipc::edge_edge_mollifier_gradient(x, eps_x));
+    expected.push_back(ipc::edge_edge_mollifier_hessian(x, eps_x));
+    expected.push_back(ipc::edge_edge_mollifier_derivative_wrt_eps_x(x, eps_x));
+    expected.push_back(
+        ipc::edge_edge_mollifier_gradient_derivative_wrt_eps_x(x, eps_x));
+
+    expected.push_back(ipc::edge_edge_mollifier(ea0, ea1, eb0, eb1, eps_x));
+    write_expected(
+        ipc::edge_edge_mollifier_gradient(ea0, ea1, eb0, eb1, eps_x), expected);
+    write_expected(
+        ipc::edge_edge_mollifier_hessian(ea0, ea1, eb0, eb1, eps_x), expected);
+    write_expected(
+        ipc::edge_edge_mollifier_gradient_wrt_x(
+            ea0, ea1, eb0, eb1, ea0, ea1, eb0, eb1),
+        expected);
+    write_expected(
+        ipc::edge_edge_mollifier_gradient_jacobian_wrt_x(
+            ea0, ea1, eb0, eb1, ea0, ea1, eb0, eb1),
+        expected);
+
+    expected.push_back(ipc::edge_edge_mollifier_threshold(ea0, ea1, eb0, eb1));
+    write_expected(
+        ipc::edge_edge_mollifier_threshold_gradient(ea0, ea1, eb0, eb1),
+        expected);
+
+    check_gpu_matches_host(mollifier_kernel, in, expected);
+}
+
+TEST_CASE("GPU signed distances", "[distance][signed][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    std::vector<double> in = SGN_IN;                 // 2D point-line
+    in.insert(in.end(), LL_IN.begin(), LL_IN.end()); // 3D line-line
+    in.insert(in.end(), PT_IN.begin(), PT_IN.end()); // 3D point-plane
+
+    const Eigen::Vector2d p2(SGN_IN[0], SGN_IN[1]);
+    const Eigen::Vector2d f0(SGN_IN[2], SGN_IN[3]);
+    const Eigen::Vector2d f1(SGN_IN[4], SGN_IN[5]);
+    const Eigen::Vector3d ea0(LL_IN[0], LL_IN[1], LL_IN[2]);
+    const Eigen::Vector3d ea1(LL_IN[3], LL_IN[4], LL_IN[5]);
+    const Eigen::Vector3d eb0(LL_IN[6], LL_IN[7], LL_IN[8]);
+    const Eigen::Vector3d eb1(LL_IN[9], LL_IN[10], LL_IN[11]);
+    const Eigen::Vector3d p(PT_IN[0], PT_IN[1], PT_IN[2]);
+    const Eigen::Vector3d t0(PT_IN[3], PT_IN[4], PT_IN[5]);
+    const Eigen::Vector3d t1(PT_IN[6], PT_IN[7], PT_IN[8]);
+    const Eigen::Vector3d t2(PT_IN[9], PT_IN[10], PT_IN[11]);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::point_line_signed_distance(p2, f0, f1));
+    write_expected(
+        ipc::point_line_signed_distance_gradient(p2, f0, f1), expected);
+    write_expected(
+        ipc::point_line_signed_distance_hessian(p2, f0, f1), expected);
+
+    expected.push_back(ipc::line_line_signed_distance(ea0, ea1, eb0, eb1));
+    write_expected(
+        ipc::line_line_signed_distance_gradient(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::line_line_signed_distance_hessian(ea0, ea1, eb0, eb1), expected);
+
+    expected.push_back(ipc::point_plane_signed_distance(p, t0, t1, t2));
+    write_expected(
+        ipc::point_plane_signed_distance_gradient(p, t0, t1, t2), expected);
+    write_expected(
+        ipc::point_plane_signed_distance_hessian(p, t0, t1, t2), expected);
+
+    check_gpu_matches_host(signed_distance_kernel, in, expected);
+}
+
+#endif // IPC_TOOLKIT_WITH_CUDA
diff --git a/tests/src/tests/geometry/CMakeLists.txt b/tests/src/tests/geometry/CMakeLists.txt
index 8e8c0b4dc..1714058ce 100644
--- a/tests/src/tests/geometry/CMakeLists.txt
+++ b/tests/src/tests/geometry/CMakeLists.txt
@@ -4,6 +4,12 @@ set(SOURCES
   test_simd_geometry.cpp
 )
 
+if(IPC_TOOLKIT_WITH_CUDA)
+  list(APPEND SOURCES
+    test_gpu_geometry.cu
+  )
+endif()
+
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
 
 ################################################################################
diff --git a/tests/src/tests/geometry/test_gpu_geometry.cu b/tests/src/tests/geometry/test_gpu_geometry.cu
new file mode 100644
index 000000000..5838ff30a
--- /dev/null
+++ b/tests/src/tests/geometry/test_gpu_geometry.cu
@@ -0,0 +1,288 @@
+// Exercises the __host__ __device__ geometry functions from CUDA kernels.
+//
+// Compilation proves the framework works: nvcc must instantiate every
+// converted geometry function (and its Eigen/std::tuple/std::array machinery)
+// as device code. On a machine with a CUDA device the kernels also run and
+// their outputs are checked against the host implementations — host/device
+// parity is the contract. Without a device (e.g. a GPU-less container) the
+// runtime checks are skipped but the device code has still been compiled.
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_CUDA
+
+#include <ipc/geometry/angle.hpp>
+#include <ipc/geometry/area.hpp>
+#include <ipc/geometry/normal.hpp>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <tests/gpu_utils.hpp>
+
+#include <cuda_runtime.h>
+
+#include <vector>
+
+using namespace ipc::tests;
+
+namespace {
+
+// ---------------------------------------------------------------------------
+
+__global__ void area_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d e0(in[0], in[1], in[2]);
+    const Eigen::Vector3d e1(in[3], in[4], in[5]);
+    const Eigen::Vector3d t0(in[6], in[7], in[8]);
+    const Eigen::Vector3d t1(in[9], in[10], in[11]);
+    const Eigen::Vector3d t2(in[12], in[13], in[14]);
+
+    int offset = 0;
+    out[offset++] = ipc::edge_length(e0, e1);
+    offset = write_out(ipc::edge_length_gradient(e0, e1), out, offset);
+    out[offset++] = ipc::triangle_area(t0, t1, t2);
+    offset = write_out(ipc::triangle_area_gradient(t0, t1, t2), out, offset);
+}
+
+__global__ void angle_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d x0(in[0], in[1], in[2]);
+    const Eigen::Vector3d x1(in[3], in[4], in[5]);
+    const Eigen::Vector3d x2(in[6], in[7], in[8]);
+    const Eigen::Vector3d x3(in[9], in[10], in[11]);
+
+    int offset = 0;
+    out[offset++] = ipc::dihedral_angle(x0, x1, x2, x3);
+    offset =
+        write_out(ipc::dihedral_angle_gradient(x0, x1, x2, x3), out, offset);
+    offset =
+        write_out(ipc::dihedral_angle_hessian(x0, x1, x2, x3), out, offset);
+}
+
+__global__ void normalization_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d x = Eigen::Vector3d(in[0], in[1], in[2]);
+
+    int offset = 0;
+    const auto [xhat, J] = ipc::normalization_and_jacobian(x);
+    offset = write_out(xhat, out, offset);
+    offset = write_out(J, out, offset);
+
+    const auto [xhat2, J2, H] = ipc::normalization_and_jacobian_and_hessian(x);
+    offset = write_out(xhat2, out, offset);
+    offset = write_out(J2, out, offset);
+    for (const auto& Hi : H) {
+        offset = write_out(Hi, out, offset);
+    }
+
+    offset = write_out(ipc::cross_product_matrix(x), out, offset);
+    offset =
+        write_out(ipc::cross_product_matrix_jacobian<double>(), out, offset);
+}
+
+__global__ void point_line_normal_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(in[3], in[4], in[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(in[6], in[7], in[8]);
+
+    int offset = 0;
+    offset =
+        write_out(ipc::point_line_unnormalized_normal(p, e0, e1), out, offset);
+    offset = write_out(ipc::point_line_normal(p, e0, e1), out, offset);
+    offset = write_out(
+        ipc::point_line_unnormalized_normal_jacobian(p, e0, e1), out, offset);
+    offset = write_out(ipc::point_line_normal_jacobian(p, e0, e1), out, offset);
+    offset = write_out(
+        ipc::point_line_unnormalized_normal_hessian(p, e0, e1), out, offset);
+    offset = write_out(ipc::point_line_normal_hessian(p, e0, e1), out, offset);
+}
+
+__global__ void triangle_normal_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d a(in[0], in[1], in[2]);
+    const Eigen::Vector3d b(in[3], in[4], in[5]);
+    const Eigen::Vector3d c(in[6], in[7], in[8]);
+
+    int offset = 0;
+    offset = write_out(ipc::triangle_unnormalized_normal(a, b, c), out, offset);
+    offset = write_out(ipc::triangle_normal(a, b, c), out, offset);
+    offset = write_out(
+        ipc::triangle_unnormalized_normal_jacobian(a, b, c), out, offset);
+    offset = write_out(ipc::triangle_normal_jacobian(a, b, c), out, offset);
+    offset = write_out(
+        ipc::triangle_unnormalized_normal_hessian(a, b, c), out, offset);
+    offset = write_out(ipc::triangle_normal_hessian(a, b, c), out, offset);
+}
+
+__global__ void line_line_normal_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d ea0(in[0], in[1], in[2]);
+    const Eigen::Vector3d ea1(in[3], in[4], in[5]);
+    const Eigen::Vector3d eb0(in[6], in[7], in[8]);
+    const Eigen::Vector3d eb1(in[9], in[10], in[11]);
+
+    int offset = 0;
+    offset = write_out(
+        ipc::line_line_unnormalized_normal(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(ipc::line_line_normal(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::line_line_unnormalized_normal_jacobian(ea0, ea1, eb0, eb1), out,
+        offset);
+    offset = write_out(
+        ipc::line_line_normal_jacobian(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::line_line_unnormalized_normal_hessian(ea0, ea1, eb0, eb1), out,
+        offset);
+    offset = write_out(
+        ipc::line_line_normal_hessian(ea0, ea1, eb0, eb1), out, offset);
+}
+
+} // namespace
+
+TEST_CASE("GPU edge length and triangle area", "[geometry][area][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d e0(0, 0, 0), e1(1, 2, 2);
+    const Eigen::Vector3d t0(-1, 0, 1), t1(1, 0, 1), t2(0, 0, -1);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::edge_length(e0, e1));
+    write_expected(ipc::edge_length_gradient(e0, e1), expected);
+    expected.push_back(ipc::triangle_area(t0, t1, t2));
+    write_expected(ipc::triangle_area_gradient(t0, t1, t2), expected);
+
+    check_gpu_matches_host(
+        area_kernel, { 0, 0, 0, 1, 2, 2, -1, 0, 1, 1, 0, 1, 0, 0, -1 },
+        expected);
+}
+
+TEST_CASE("GPU dihedral angle", "[geometry][angle][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    // Two triangles sharing the edge (x0, x1), folded at 90°.
+    const Eigen::Vector3d x0(0, -1, 0), x1(0, 1, 0);
+    const Eigen::Vector3d x2(0.5, 0, 0.5), x3(-0.5, 0, 0.5);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::dihedral_angle(x0, x1, x2, x3));
+    write_expected(ipc::dihedral_angle_gradient(x0, x1, x2, x3), expected);
+    write_expected(ipc::dihedral_angle_hessian(x0, x1, x2, x3), expected);
+
+    check_gpu_matches_host(
+        angle_kernel, { 0, -1, 0, 0, 1, 0, 0.5, 0, 0.5, -0.5, 0, 0.5 },
+        expected);
+}
+
+TEST_CASE("GPU normalization and cross product", "[geometry][normal][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d x = Eigen::Vector3d(1, 2, 2);
+
+    std::vector<double> expected;
+    const auto [xhat, J] = ipc::normalization_and_jacobian(x);
+    write_expected(xhat, expected);
+    write_expected(J, expected);
+    const auto [xhat2, J2, H] = ipc::normalization_and_jacobian_and_hessian(x);
+    write_expected(xhat2, expected);
+    write_expected(J2, expected);
+    for (const auto& Hi : H) {
+        write_expected(Hi, expected);
+    }
+    write_expected(
+        ipc::cross_product_matrix(Eigen::Vector3d(1, 2, 2)), expected);
+    write_expected(ipc::cross_product_matrix_jacobian<double>(), expected);
+
+    check_gpu_matches_host(normalization_kernel, { 1, 2, 2 }, expected);
+}
+
+TEST_CASE("GPU point-line normal", "[geometry][normal][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d p(0.1, 1, 0.2), e0(-1, 0, 0), e1(1, 0, 0.1);
+
+    std::vector<double> expected;
+    write_expected(ipc::point_line_unnormalized_normal(p, e0, e1), expected);
+    write_expected(ipc::point_line_normal(p, e0, e1), expected);
+    write_expected(
+        ipc::point_line_unnormalized_normal_jacobian(p, e0, e1), expected);
+    write_expected(ipc::point_line_normal_jacobian(p, e0, e1), expected);
+    write_expected(
+        ipc::point_line_unnormalized_normal_hessian(p, e0, e1), expected);
+    write_expected(ipc::point_line_normal_hessian(p, e0, e1), expected);
+
+    check_gpu_matches_host(
+        point_line_normal_kernel, { 0.1, 1, 0.2, -1, 0, 0, 1, 0, 0.1 },
+        expected);
+}
+
+TEST_CASE("GPU triangle normal", "[geometry][normal][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d a(-1, 0, 1), b(1, 0, 1), c(0, 0.1, -1);
+
+    std::vector<double> expected;
+    write_expected(ipc::triangle_unnormalized_normal(a, b, c), expected);
+    write_expected(ipc::triangle_normal(a, b, c), expected);
+    write_expected(
+        ipc::triangle_unnormalized_normal_jacobian(a, b, c), expected);
+    write_expected(ipc::triangle_normal_jacobian(a, b, c), expected);
+    write_expected(
+        ipc::triangle_unnormalized_normal_hessian(a, b, c), expected);
+    write_expected(ipc::triangle_normal_hessian(a, b, c), expected);
+
+    check_gpu_matches_host(
+        triangle_normal_kernel, { -1, 0, 1, 1, 0, 1, 0, 0.1, -1 }, expected);
+}
+
+TEST_CASE("GPU line-line normal", "[geometry][normal][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d ea0(-1, 0, 0), ea1(1, 0, 0);
+    const Eigen::Vector3d eb0(0, 1, -1), eb1(0.1, 1, 1);
+
+    std::vector<double> expected;
+    write_expected(
+        ipc::line_line_unnormalized_normal(ea0, ea1, eb0, eb1), expected);
+    write_expected(ipc::line_line_normal(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::line_line_unnormalized_normal_jacobian(ea0, ea1, eb0, eb1),
+        expected);
+    write_expected(
+        ipc::line_line_normal_jacobian(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::line_line_unnormalized_normal_hessian(ea0, ea1, eb0, eb1),
+        expected);
+    write_expected(ipc::line_line_normal_hessian(ea0, ea1, eb0, eb1), expected);
+
+    check_gpu_matches_host(
+        line_line_normal_kernel, { -1, 0, 0, 1, 0, 0, 0, 1, -1, 0.1, 1, 1 },
+        expected);
+}
+
+#endif // IPC_TOOLKIT_WITH_CUDA
diff --git a/tests/src/tests/gpu_utils.hpp b/tests/src/tests/gpu_utils.hpp
new file mode 100644
index 000000000..0f841b665
--- /dev/null
+++ b/tests/src/tests/gpu_utils.hpp
@@ -0,0 +1,116 @@
+#pragma once
+
+// Shared scaffolding for the GPU parity tests
+// (tests/src/tests/**/test_gpu_*.cu).
+//
+// Each of those files launches one-thread kernels over a flat double buffer and
+// checks the results against the host implementations. The launch, transfer,
+// and comparison are identical everywhere, so they live here; only the kernels
+// and their fixtures belong in the individual files.
+//
+// This header defines __device__ code, so include it from a .cu, inside the
+// file's `#ifdef IPC_TOOLKIT_WITH_CUDA` guard.
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <Eigen/Core>
+
+#include <cuda_runtime.h>
+
+#include <cmath>
+#include <vector>
+
+#ifndef __CUDACC__
+#error "tests/gpu_utils.hpp defines device code; include it from a .cu file."
+#endif
+
+#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess)
+
+namespace ipc::tests {
+
+/// @brief Uniform kernel signature: unpack inputs from `in`, write results
+/// flattened (row-major) to `out`.
+using GpuTestKernel = void (*)(const double* in, double* out);
+
+/// @brief Copy a fixed-size Eigen object into `out` row-major; returns the
+/// number of doubles written.
+template <typename Derived>
+__device__ int
+write_out(const Eigen::MatrixBase<Derived>& m, double* out, int offset)
+{
+    for (int r = 0; r < m.rows(); ++r) {
+        for (int c = 0; c < m.cols(); ++c) {
+            out[offset++] = m(r, c);
+        }
+    }
+    return offset;
+}
+
+/// @brief Append a fixed-size Eigen object to a host vector row-major.
+template <typename Derived>
+inline void write_expected(
+    const Eigen::MatrixBase<Derived>& m, std::vector<double>& expected)
+{
+    for (int r = 0; r < m.rows(); ++r) {
+        for (int c = 0; c < m.cols(); ++c) {
+            expected.push_back(m(r, c));
+        }
+    }
+}
+
+/// @brief SKIP the test if no CUDA device is present (kernels have still been
+/// compiled as device code, which is the primary verification).
+inline void skip_if_no_cuda_device()
+{
+    int device_count = 0;
+    const cudaError_t err = cudaGetDeviceCount(&device_count);
+    if (err != cudaSuccess || device_count == 0) {
+        SKIP("No CUDA device available; kernels compiled but not executed.");
+    }
+}
+
+/// @brief Launch `kernel` with the given inputs; returns the device outputs.
+inline std::vector<double> run_gpu_kernel(
+    GpuTestKernel kernel, const std::vector<double>& in, const size_t n_out)
+{
+    double *d_in = nullptr, *d_out = nullptr;
+    REQUIRE_CUDA(cudaMalloc(&d_in, in.size() * sizeof(double)));
+    REQUIRE_CUDA(cudaMalloc(&d_out, n_out * sizeof(double)));
+    REQUIRE_CUDA(cudaMemcpy(
+        d_in, in.data(), in.size() * sizeof(double), cudaMemcpyHostToDevice));
+
+    kernel<<<1, 1>>>(d_in, d_out);
+    REQUIRE_CUDA(cudaGetLastError());
+    REQUIRE_CUDA(cudaDeviceSynchronize());
+
+    std::vector<double> out(n_out);
+    REQUIRE_CUDA(cudaMemcpy(
+        out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost));
+    cudaFree(d_in);
+    cudaFree(d_out);
+    return out;
+}
+
+/// @brief Run `kernel` and check every output against the host reference.
+///
+/// An infinite reference is compared by sign rather than by margin, so a
+/// barrier evaluated at d <= 0 can be checked alongside finite values.
+inline void check_gpu_matches_host(
+    GpuTestKernel kernel,
+    const std::vector<double>& in,
+    const std::vector<double>& expected)
+{
+    const std::vector<double> out = run_gpu_kernel(kernel, in, expected.size());
+    for (size_t i = 0; i < expected.size(); ++i) {
+        CAPTURE(i);
+        if (std::isinf(expected[i])) {
+            CHECK(std::isinf(out[i]));
+            CHECK((out[i] > 0) == (expected[i] > 0));
+        } else {
+            CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12));
+        }
+    }
+}
+
+} // namespace ipc::tests
diff --git a/tests/src/tests/tangent/CMakeLists.txt b/tests/src/tests/tangent/CMakeLists.txt
index 6a265b934..5a9326648 100644
--- a/tests/src/tests/tangent/CMakeLists.txt
+++ b/tests/src/tests/tangent/CMakeLists.txt
@@ -12,6 +12,12 @@ set(SOURCES
   # Utilities
 )
 
+if(IPC_TOOLKIT_WITH_CUDA)
+  list(APPEND SOURCES
+    test_gpu_tangent.cu
+  )
+endif()
+
 target_sources(ipc_toolkit_tests PRIVATE ${SOURCES})
 
 ################################################################################
diff --git a/tests/src/tests/tangent/test_gpu_tangent.cu b/tests/src/tests/tangent/test_gpu_tangent.cu
new file mode 100644
index 000000000..6b682a502
--- /dev/null
+++ b/tests/src/tests/tangent/test_gpu_tangent.cu
@@ -0,0 +1,281 @@
+// Exercises the __host__ __device__ tangent functions from CUDA kernels.
+//
+// Compilation proves the framework works: nvcc must instantiate every
+// converted tangent function (tangent bases, closest points — including the
+// std::array<Matrix12d, 2> hessians — and relative velocities) as device
+// code. On a machine with a CUDA device the kernels also run and their
+// outputs are checked against the host implementations — host/device parity
+// is the contract. Without a device (e.g. a GPU-less container) the runtime
+// checks are skipped but the device code has still been compiled.
+
+#include <ipc/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_CUDA
+
+#include <ipc/tangent/closest_point.hpp>
+#include <ipc/tangent/relative_velocity.hpp>
+#include <ipc/tangent/tangent_basis.hpp>
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <tests/gpu_utils.hpp>
+
+#include <cuda_runtime.h>
+
+#include <vector>
+
+using namespace ipc::tests;
+
+namespace {
+
+// ---------------------------------------------------------------------------
+
+// Shared fixtures (kept in one place so kernels and host references agree).
+// clang-format off
+const std::vector<double> PE_IN = { 0.5, 1, 0.2,  -1, 0, 0,  1, 0, 0.1 };
+const std::vector<double> EE_IN = { -1, 0, 0,  1, 0, 0,  -0.5, 1, -1,  0.5, 1, 1 };
+const std::vector<double> PT_IN = { 0, 1, 0.1,  -1, 0, 1,  1, 0, 1,  0, 0, -1 };
+// Velocities for the relative-velocity kernels (4 x 3D) + coords (2).
+const std::vector<double> RV_IN = { 0.1, -0.2, 0.3,  0.4, 0.5, -0.6,
+                                    -0.7, 0.8, 0.9,  1.0, -1.1, 1.2,
+                                    0.25, 0.5 };
+// clang-format on
+
+__global__ void tangent_basis_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(in[3], in[4], in[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(in[6], in[7], in[8]);
+    const Eigen::Vector3d ea0(in[9], in[10], in[11]);
+    const Eigen::Vector3d ea1(in[12], in[13], in[14]);
+    const Eigen::Vector3d eb0(in[15], in[16], in[17]);
+    const Eigen::Vector3d eb1(in[18], in[19], in[20]);
+    const Eigen::Vector3d q(in[21], in[22], in[23]);
+    const Eigen::Vector3d t0(in[24], in[25], in[26]);
+    const Eigen::Vector3d t1(in[27], in[28], in[29]);
+    const Eigen::Vector3d t2(in[30], in[31], in[32]);
+
+    int offset = 0;
+    offset = write_out(ipc::point_point_tangent_basis(p, e0), out, offset);
+    offset =
+        write_out(ipc::point_point_tangent_basis_jacobian(p, e0), out, offset);
+    offset = write_out(ipc::point_edge_tangent_basis(p, e0, e1), out, offset);
+    offset = write_out(
+        ipc::point_edge_tangent_basis_jacobian(p, e0, e1), out, offset);
+    offset = write_out(
+        ipc::edge_edge_tangent_basis(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::edge_edge_tangent_basis_jacobian(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::point_triangle_tangent_basis(q, t0, t1, t2), out, offset);
+    offset = write_out(
+        ipc::point_triangle_tangent_basis_jacobian(q, t0, t1, t2), out, offset);
+}
+
+__global__ void closest_point_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const ipc::VectorMax3d p = Eigen::Vector3d(in[0], in[1], in[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(in[3], in[4], in[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(in[6], in[7], in[8]);
+    const Eigen::Vector3d ea0(in[9], in[10], in[11]);
+    const Eigen::Vector3d ea1(in[12], in[13], in[14]);
+    const Eigen::Vector3d eb0(in[15], in[16], in[17]);
+    const Eigen::Vector3d eb1(in[18], in[19], in[20]);
+    const Eigen::Vector3d q(in[21], in[22], in[23]);
+    const Eigen::Vector3d t0(in[24], in[25], in[26]);
+    const Eigen::Vector3d t1(in[27], in[28], in[29]);
+    const Eigen::Vector3d t2(in[30], in[31], in[32]);
+
+    int offset = 0;
+    out[offset++] = ipc::point_edge_closest_point(p, e0, e1);
+    offset = write_out(
+        ipc::point_edge_closest_point_jacobian(p, e0, e1), out, offset);
+    offset = write_out(
+        ipc::point_edge_closest_point_hessian(p, e0, e1), out, offset);
+
+    offset = write_out(
+        ipc::edge_edge_closest_point(ea0, ea1, eb0, eb1), out, offset);
+    offset = write_out(
+        ipc::edge_edge_closest_point_jacobian(ea0, ea1, eb0, eb1), out, offset);
+    const std::array<ipc::Matrix12d, 2> ee_hess =
+        ipc::edge_edge_closest_point_hessian(ea0, ea1, eb0, eb1);
+    offset = write_out(ee_hess[0], out, offset);
+    offset = write_out(ee_hess[1], out, offset);
+
+    offset = write_out(
+        ipc::point_triangle_closest_point(q, t0, t1, t2), out, offset);
+    offset = write_out(
+        ipc::point_triangle_closest_point_jacobian(q, t0, t1, t2), out, offset);
+    const std::array<ipc::Matrix12d, 2> pt_hess =
+        ipc::point_triangle_closest_point_hessian(q, t0, t1, t2);
+    offset = write_out(pt_hess[0], out, offset);
+    offset = write_out(pt_hess[1], out, offset);
+}
+
+__global__ void relative_velocity_kernel(const double* in, double* out)
+{
+    if (blockIdx.x != 0 || threadIdx.x != 0) {
+        return;
+    }
+    const Eigen::Vector3d dp(in[0], in[1], in[2]);
+    const Eigen::Vector3d dv0(in[3], in[4], in[5]);
+    const Eigen::Vector3d dv1(in[6], in[7], in[8]);
+    const Eigen::Vector3d dv2(in[9], in[10], in[11]);
+    const Eigen::Vector2d coords(in[12], in[13]);
+    const double alpha = in[12];
+
+    int offset = 0;
+    offset =
+        write_out(ipc::point_point_relative_velocity(dp, dv0), out, offset);
+    offset =
+        write_out(ipc::point_point_relative_velocity_jacobian(3), out, offset);
+    offset =
+        write_out(ipc::point_point_relative_velocity_dx_dbeta(3), out, offset);
+
+    offset = write_out(
+        ipc::point_edge_relative_velocity(dp, dv0, dv1, alpha), out, offset);
+    offset = write_out(
+        ipc::point_edge_relative_velocity_jacobian(3, alpha), out, offset);
+    offset = write_out(
+        ipc::point_edge_relative_velocity_dx_dbeta(3, alpha), out, offset);
+
+    offset = write_out(
+        ipc::edge_edge_relative_velocity(dp, dv0, dv1, dv2, coords), out,
+        offset);
+    offset = write_out(
+        ipc::edge_edge_relative_velocity_jacobian(coords), out, offset);
+    offset = write_out(
+        ipc::edge_edge_relative_velocity_dx_dbeta(coords), out, offset);
+
+    offset = write_out(
+        ipc::point_triangle_relative_velocity(dp, dv0, dv1, dv2, coords), out,
+        offset);
+    offset = write_out(
+        ipc::point_triangle_relative_velocity_jacobian(coords), out, offset);
+    offset = write_out(
+        ipc::point_triangle_relative_velocity_dx_dbeta(coords), out, offset);
+}
+
+std::vector<double> tangent_fixture()
+{
+    std::vector<double> in = PE_IN;
+    in.insert(in.end(), EE_IN.begin(), EE_IN.end());
+    in.insert(in.end(), PT_IN.begin(), PT_IN.end());
+    return in;
+}
+
+} // namespace
+
+TEST_CASE("GPU tangent basis", "[tangent][tangent_basis][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d p = Eigen::Vector3d(PE_IN[0], PE_IN[1], PE_IN[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(PE_IN[3], PE_IN[4], PE_IN[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(PE_IN[6], PE_IN[7], PE_IN[8]);
+    const Eigen::Vector3d ea0(EE_IN[0], EE_IN[1], EE_IN[2]);
+    const Eigen::Vector3d ea1(EE_IN[3], EE_IN[4], EE_IN[5]);
+    const Eigen::Vector3d eb0(EE_IN[6], EE_IN[7], EE_IN[8]);
+    const Eigen::Vector3d eb1(EE_IN[9], EE_IN[10], EE_IN[11]);
+    const Eigen::Vector3d q(PT_IN[0], PT_IN[1], PT_IN[2]);
+    const Eigen::Vector3d t0(PT_IN[3], PT_IN[4], PT_IN[5]);
+    const Eigen::Vector3d t1(PT_IN[6], PT_IN[7], PT_IN[8]);
+    const Eigen::Vector3d t2(PT_IN[9], PT_IN[10], PT_IN[11]);
+
+    std::vector<double> expected;
+    write_expected(ipc::point_point_tangent_basis(p, e0), expected);
+    write_expected(ipc::point_point_tangent_basis_jacobian(p, e0), expected);
+    write_expected(ipc::point_edge_tangent_basis(p, e0, e1), expected);
+    write_expected(ipc::point_edge_tangent_basis_jacobian(p, e0, e1), expected);
+    write_expected(ipc::edge_edge_tangent_basis(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::edge_edge_tangent_basis_jacobian(ea0, ea1, eb0, eb1), expected);
+    write_expected(ipc::point_triangle_tangent_basis(q, t0, t1, t2), expected);
+    write_expected(
+        ipc::point_triangle_tangent_basis_jacobian(q, t0, t1, t2), expected);
+
+    check_gpu_matches_host(tangent_basis_kernel, tangent_fixture(), expected);
+}
+
+TEST_CASE("GPU closest point", "[tangent][closest_point][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const ipc::VectorMax3d p = Eigen::Vector3d(PE_IN[0], PE_IN[1], PE_IN[2]);
+    const ipc::VectorMax3d e0 = Eigen::Vector3d(PE_IN[3], PE_IN[4], PE_IN[5]);
+    const ipc::VectorMax3d e1 = Eigen::Vector3d(PE_IN[6], PE_IN[7], PE_IN[8]);
+    const Eigen::Vector3d ea0(EE_IN[0], EE_IN[1], EE_IN[2]);
+    const Eigen::Vector3d ea1(EE_IN[3], EE_IN[4], EE_IN[5]);
+    const Eigen::Vector3d eb0(EE_IN[6], EE_IN[7], EE_IN[8]);
+    const Eigen::Vector3d eb1(EE_IN[9], EE_IN[10], EE_IN[11]);
+    const Eigen::Vector3d q(PT_IN[0], PT_IN[1], PT_IN[2]);
+    const Eigen::Vector3d t0(PT_IN[3], PT_IN[4], PT_IN[5]);
+    const Eigen::Vector3d t1(PT_IN[6], PT_IN[7], PT_IN[8]);
+    const Eigen::Vector3d t2(PT_IN[9], PT_IN[10], PT_IN[11]);
+
+    std::vector<double> expected;
+    expected.push_back(ipc::point_edge_closest_point(p, e0, e1));
+    write_expected(ipc::point_edge_closest_point_jacobian(p, e0, e1), expected);
+    write_expected(ipc::point_edge_closest_point_hessian(p, e0, e1), expected);
+    write_expected(ipc::edge_edge_closest_point(ea0, ea1, eb0, eb1), expected);
+    write_expected(
+        ipc::edge_edge_closest_point_jacobian(ea0, ea1, eb0, eb1), expected);
+    const std::array<ipc::Matrix12d, 2> ee_hess =
+        ipc::edge_edge_closest_point_hessian(ea0, ea1, eb0, eb1);
+    write_expected(ee_hess[0], expected);
+    write_expected(ee_hess[1], expected);
+    write_expected(ipc::point_triangle_closest_point(q, t0, t1, t2), expected);
+    write_expected(
+        ipc::point_triangle_closest_point_jacobian(q, t0, t1, t2), expected);
+    const std::array<ipc::Matrix12d, 2> pt_hess =
+        ipc::point_triangle_closest_point_hessian(q, t0, t1, t2);
+    write_expected(pt_hess[0], expected);
+    write_expected(pt_hess[1], expected);
+
+    check_gpu_matches_host(closest_point_kernel, tangent_fixture(), expected);
+}
+
+TEST_CASE("GPU relative velocity", "[tangent][relative_velocity][gpu]")
+{
+    skip_if_no_cuda_device();
+
+    const Eigen::Vector3d dp(RV_IN[0], RV_IN[1], RV_IN[2]);
+    const Eigen::Vector3d dv0(RV_IN[3], RV_IN[4], RV_IN[5]);
+    const Eigen::Vector3d dv1(RV_IN[6], RV_IN[7], RV_IN[8]);
+    const Eigen::Vector3d dv2(RV_IN[9], RV_IN[10], RV_IN[11]);
+    const Eigen::Vector2d coords(RV_IN[12], RV_IN[13]);
+    const double alpha = RV_IN[12];
+
+    std::vector<double> expected;
+    write_expected(ipc::point_point_relative_velocity(dp, dv0), expected);
+    write_expected(ipc::point_point_relative_velocity_jacobian(3), expected);
+    write_expected(ipc::point_point_relative_velocity_dx_dbeta(3), expected);
+    write_expected(
+        ipc::point_edge_relative_velocity(dp, dv0, dv1, alpha), expected);
+    write_expected(
+        ipc::point_edge_relative_velocity_jacobian(3, alpha), expected);
+    write_expected(
+        ipc::point_edge_relative_velocity_dx_dbeta(3, alpha), expected);
+    write_expected(
+        ipc::edge_edge_relative_velocity(dp, dv0, dv1, dv2, coords), expected);
+    write_expected(ipc::edge_edge_relative_velocity_jacobian(coords), expected);
+    write_expected(ipc::edge_edge_relative_velocity_dx_dbeta(coords), expected);
+    write_expected(
+        ipc::point_triangle_relative_velocity(dp, dv0, dv1, dv2, coords),
+        expected);
+    write_expected(
+        ipc::point_triangle_relative_velocity_jacobian(coords), expected);
+    write_expected(
+        ipc::point_triangle_relative_velocity_dx_dbeta(coords), expected);
+
+    check_gpu_matches_host(relative_velocity_kernel, RV_IN, expected);
+}
+
+#endif // IPC_TOOLKIT_WITH_CUDA