Make the low-level math callable from device code

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})
+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
+)

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

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

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
 )
 
 ################################################################################

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