From 511cb01c29ebe32662f69094c85bd9065d921c78 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 8 Sep 2026 00:29:15 -0400 Subject: [PATCH 1/4] Make the low-level math callable from device code Extends the mechanism from the previous commit across distance/, geometry/, tangent/, barrier/, and math/, and adds GPU unit tests that check each function against its host result. Annotating the declarations is most of the diff, but a few functions had to stop using host-only constructs to compile for the device: - distance_type.cpp is split, moving the logger and the exceptions into distance_type_host.cpp. nvcc cannot parse spdlog, so a shared device source must not include it; the header picks between the host helpers and a device trap. - barrier.cpp keeps only the three free functions. The virtual Barrier hierarchy moves to barrier_classes.cpp, which stays host-only -- virtual dispatch has nothing to offer device code. - morton.hpp gets a device-safe clamp. MSVC's debug std::clamp checks bounds with _STL_VERIFY, which expands to __debugbreak(), and nvcc's NVVM backend then emits invalid IR. The comparison order matches std::clamp exactly. - angle.cpp scatters the normal-Jacobian blocks with middleCols instead of dn1_dx(Eigen::all, idx), since Eigen's index slicing is unavailable in device code. It also returns 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, so a pair built at run time inside a kernel comes back zero-filled. The xsimd instantiation blocks are guarded on IPC_TOOLKIT_WITH_SIMD alone. The parent commit retracts that macro for nvcc, so the blocks no longer need an IPC_TOOLKIT_INSTANTIATE_HOST_SCALARS term; only the autodiff instantiations in edge_edge.cpp and point_triangle.cpp still use it. Co-Authored-By: Claude Opus 5 --- src/ipc/barrier/CMakeLists.txt | 12 +- src/ipc/barrier/barrier.cpp | 151 +---- src/ipc/barrier/barrier.hpp | 7 +- src/ipc/barrier/barrier_classes.cpp | 159 +++++ src/ipc/distance/CMakeLists.txt | 17 +- src/ipc/distance/distance_type.cpp | 33 +- src/ipc/distance/distance_type.hpp | 74 ++- src/ipc/distance/distance_type_host.cpp | 41 ++ src/ipc/distance/edge_edge.cpp | 10 +- src/ipc/distance/edge_edge.hpp | 13 +- src/ipc/distance/edge_edge_mollifier.cpp | 16 +- src/ipc/distance/edge_edge_mollifier.hpp | 68 +- src/ipc/distance/line_line.cpp | 6 +- src/ipc/distance/line_line.hpp | 13 +- src/ipc/distance/point_edge.cpp | 5 +- src/ipc/distance/point_edge.hpp | 15 +- src/ipc/distance/point_plane.cpp | 6 +- src/ipc/distance/point_plane.hpp | 33 +- src/ipc/distance/point_point.hpp | 15 +- src/ipc/distance/point_triangle.cpp | 11 +- src/ipc/distance/point_triangle.hpp | 15 +- src/ipc/distance/signed/CMakeLists.txt | 13 +- src/ipc/distance/signed/line_line.cpp | 19 +- src/ipc/distance/signed/line_line.hpp | 15 +- src/ipc/distance/signed/point_line.cpp | 19 +- src/ipc/distance/signed/point_line.hpp | 15 +- src/ipc/distance/signed/point_plane.cpp | 21 +- src/ipc/distance/signed/point_plane.hpp | 15 +- src/ipc/geometry/CMakeLists.txt | 11 +- src/ipc/geometry/angle.cpp | 36 +- src/ipc/geometry/angle.hpp | 13 +- src/ipc/geometry/area.cpp | 4 +- src/ipc/geometry/area.hpp | 20 +- src/ipc/geometry/normal.cpp | 29 +- src/ipc/geometry/normal.hpp | 113 ++-- src/ipc/math/morton.hpp | 39 +- src/ipc/math/scalar_math.hpp | 12 +- src/ipc/tangent/CMakeLists.txt | 9 +- src/ipc/tangent/closest_point.cpp | 18 +- src/ipc/tangent/closest_point.hpp | 60 +- src/ipc/tangent/relative_velocity.hpp | 59 +- src/ipc/tangent/tangent_basis.cpp | 34 +- src/ipc/tangent/tangent_basis.hpp | 53 +- src/ipc/utils/simd.hpp | 22 +- tests/src/tests/barrier/CMakeLists.txt | 6 + tests/src/tests/barrier/test_gpu_barrier.cu | 93 +++ tests/src/tests/distance/CMakeLists.txt | 6 + tests/src/tests/distance/test_gpu_distance.cu | 589 ++++++++++++++++++ tests/src/tests/geometry/CMakeLists.txt | 6 + tests/src/tests/geometry/test_gpu_geometry.cu | 362 +++++++++++ tests/src/tests/tangent/CMakeLists.txt | 6 + tests/src/tests/tangent/test_gpu_tangent.cu | 354 +++++++++++ 52 files changed, 2274 insertions(+), 517 deletions(-) create mode 100644 src/ipc/barrier/barrier_classes.cpp create mode 100644 src/ipc/distance/distance_type_host.cpp create mode 100644 tests/src/tests/barrier/test_gpu_barrier.cu create mode 100644 tests/src/tests/distance/test_gpu_distance.cu create mode 100644 tests/src/tests/geometry/test_gpu_geometry.cu create mode 100644 tests/src/tests/tangent/test_gpu_tangent.cu diff --git a/src/ipc/barrier/CMakeLists.txt b/src/ipc/barrier/CMakeLists.txt index 70f573cca..913951abf 100644 --- a/src/ipc/barrier/CMakeLists.txt +++ b/src/ipc/barrier/CMakeLists.txt @@ -1,10 +1,18 @@ set(SOURCES adaptive_stiffness.cpp adaptive_stiffness.hpp + barrier_classes.cpp 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 (the free functions) is shared between host C++ and CUDA device +# code; under CUDA the helper generates a .cu wrapper (see +# ipc_toolkit_shared_device_sources). The other barrier .cpp files are ordinary +# host-only sources listed in SOURCES above. +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..cdf30e368 100644 --- a/src/ipc/barrier/barrier.cpp +++ b/src/ipc/barrier/barrier.cpp @@ -15,189 +15,58 @@ namespace ipc { // 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 T barrier(const T d, const T dhat) +template IPC_TOOLKIT_HOST_DEVICE T barrier(const T d, const T dhat) { using namespace ipc::numext; // log // b(d) = -(d-d̂)²ln(d / d̂) return select_lazy( d <= T(0), [&] { return infinity(); }, // - d < dhat, [&] { return -sqr(d - dhat) * log(d / dhat); }, + d < dhat, [&] { return -sqr(d - dhat) * ipc::numext::log(d / dhat); }, [&] { return T(0); }); } -template T barrier_first_derivative(const T d, const T dhat) +template +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̂) // b'(d) = -2(d - d̂)ln(d / d̂) - (d-d̂)²(1 / d) // = (d - d̂) * (-2ln(d/d̂) - (d - d̂) / d) // = (d̂ - d) * (2ln(d/d̂) - d̂/d + 1) - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < dhat, - [&] { return (dhat - d) * (2 * log(d / dhat) - dhat / d + 1); }, - [&] { return T(0); }); -} - -template T barrier_second_derivative(const T d, const T dhat) -{ - using namespace ipc::numext; // log - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < dhat, - [&] { - const T dhat_d = dhat / d; - return (dhat_d + 2) * dhat_d - 2 * log(d / dhat) - 3; - }, - [&] { return T(0); }); -} - -// ============================================================================ - -template -T ClampedLogSqBarrier::operator()(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - // b(d) = (d-d̂)²ln²(d / d̂) - return select_lazy( - d <= T(0), [&] { return infinity(); }, // - d < dhat, - [&] { - const T log_d_dhat = log(d / dhat); - return sqr(d - dhat) * sqr(log_d_dhat); - }, - [&] { return T(0); }); -} - -template -T ClampedLogSqBarrier::first_derivative(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - // b(d) = (d - d̂)²ln²(d / d̂) - // b'(d) = 2 (d - d̂) ln²(d / d̂) + 2 (d - d̂)² ln(d / d̂) / d - // = 2 (d - d̂) ln(d / d̂) [ln(d / d̂) + (d - d̂) / d] return select_lazy( d <= T(0), [&] { return T(0); }, // d < dhat, [&] { - const T d_minus_dhat = (d - dhat); - const T log_d_dhat = log(d / dhat); - return T(2) * d_minus_dhat * log_d_dhat - * (log_d_dhat + d_minus_dhat / d); + return (dhat - d) * (2 * ipc::numext::log(d / dhat) - dhat / d + 1); }, [&] { return T(0); }); } template -T ClampedLogSqBarrier::second_derivative(const T d, const T dhat) const +IPC_TOOLKIT_HOST_DEVICE T barrier_second_derivative(const T d, const T dhat) { using namespace ipc::numext; // log return select_lazy( d <= T(0), [&] { return T(0); }, // d < dhat, [&] { - const T t0 = dhat - d; - const T t1 = log(d / dhat); - const T t2 = sqr(t0) / sqr(d); - return T(2) * (sqr(t1) - (t1 - T(1)) * t2 - T(4) * t1 * t0 / d); + const T dhat_d = dhat / d; + return (dhat_d + 2) * dhat_d - 2 * ipc::numext::log(d / dhat) - 3; }, [&] { return T(0); }); } - -// ============================================================================ - -template -T CubicBarrier::operator()(const T d, const T dhat) const -{ - // b(d) = (d - d̂)³ - // - // A polynomial: finite at d <= 0, so unlike the log barriers it needs no - // penetration guard. - return select_lazy( - d < dhat, [&] { return -T(2) / T(3) / dhat * cubic(d - dhat); }, - [&] { return T(0); }); -} - -template -T CubicBarrier::first_derivative(const T d, const T dhat) const -{ - return select_lazy( - d < dhat, [&] { return T(-2) / dhat * sqr(d - dhat); }, - [&] { return T(0); }); -} - -template -T CubicBarrier::second_derivative(const T d, const T dhat) const -{ - return select_lazy( - d < dhat, [&] { return T(4) * (T(1) - d / dhat); }, - [&] { return T(0); }); -} - -// ============================================================================ - -template -T TwoStageBarrier::operator()(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - return select_lazy( - d <= T(0), [&] { return infinity(); }, // - d < T(0.5) * dhat, - [&] { return T(-0.25) * sqr(dhat) * (log(T(2) * d / dhat) - T(0.5)); }, - d < dhat, [&] { return T(0.5) * sqr(dhat - d); }, // - [&] { return T(0); }); -} - -template -T TwoStageBarrier::first_derivative(const T d, const T dhat) const -{ - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < T(0.5) * dhat, [&] { return T(-0.25) * sqr(dhat) / d; }, // - d < dhat, [&] { return d - dhat; }, // - [&] { return T(0); }); -} - -template -T TwoStageBarrier::second_derivative(const T d, const T dhat) const -{ - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < T(0.5) * dhat, [&] { return (T(0.25) * sqr(dhat)) / sqr(d); }, // - d < dhat, [&] { return T(1); }, // - [&] { return T(0); }); -} - // ============================================================================ // Explicit template instantiations /// @cond DOXYGEN_SKIP -template class BarrierBase; -template class BarrierBase; -template class ClampedLogBarrier; -template class ClampedLogBarrier; -template class ClampedLogSqBarrier; -template class ClampedLogSqBarrier; -template class CubicBarrier; -template class CubicBarrier; -template class TwoStageBarrier; -template class TwoStageBarrier; +#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 class BarrierBase>; -template class BarrierBase>; -template class ClampedLogBarrier>; -template class ClampedLogBarrier>; -template class ClampedLogSqBarrier>; -template class ClampedLogSqBarrier>; -template class CubicBarrier>; -template class CubicBarrier>; -template class TwoStageBarrier>; -template class TwoStageBarrier>; template SimdBatch barrier(const SimdBatch d, const SimdBatch dhat); template SimdBatch 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 T barrier(const T d, const T dhat); +template +IPC_TOOLKIT_HOST_DEVICE T barrier(const T d, const T dhat); /// @brief Derivative of the barrier function. /// @@ -71,7 +72,7 @@ template T barrier(const T d, const T dhat); /// @param dhat Activation distance of the barrier. /// @return The derivative of the barrier wrt d. template -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 -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 class ClampedLogBarrier : public BarrierBase { diff --git a/src/ipc/barrier/barrier_classes.cpp b/src/ipc/barrier/barrier_classes.cpp new file mode 100644 index 000000000..2be54fdf5 --- /dev/null +++ b/src/ipc/barrier/barrier_classes.cpp @@ -0,0 +1,159 @@ +// Out-of-line definitions for the Barrier class hierarchy. +// +// These live apart from barrier.cpp because that file is compiled as a CUDA +// translation unit when CUDA is enabled, and a virtual hierarchy has nothing to +// offer device code anyway. Keeping the classes here lets barrier.cpp stay +// small and device-callable, holding only the three free barrier functions. +#include "barrier.hpp" + +#include +#include +#include + +namespace ipc { + +// ============================================================================ + +template +T ClampedLogSqBarrier::operator()(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + // b(d) = (d-d̂)²ln²(d / d̂) + return select_lazy( + d <= T(0), [&] { return infinity(); }, // + d < dhat, + [&] { + const T log_d_dhat = log(d / dhat); + return sqr(d - dhat) * sqr(log_d_dhat); + }, + [&] { return T(0); }); +} + +template +T ClampedLogSqBarrier::first_derivative(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + // b(d) = (d - d̂)²ln²(d / d̂) + // b'(d) = 2 (d - d̂) ln²(d / d̂) + 2 (d - d̂)² ln(d / d̂) / d + // = 2 (d - d̂) ln(d / d̂) [ln(d / d̂) + (d - d̂) / d] + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < dhat, + [&] { + const T d_minus_dhat = (d - dhat); + const T log_d_dhat = log(d / dhat); + return T(2) * d_minus_dhat * log_d_dhat + * (log_d_dhat + d_minus_dhat / d); + }, + [&] { return T(0); }); +} + +template +T ClampedLogSqBarrier::second_derivative(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < dhat, + [&] { + const T t0 = dhat - d; + const T t1 = log(d / dhat); + const T t2 = sqr(t0) / sqr(d); + return T(2) * (sqr(t1) - (t1 - T(1)) * t2 - T(4) * t1 * t0 / d); + }, + [&] { return T(0); }); +} + +// ============================================================================ + +template +T CubicBarrier::operator()(const T d, const T dhat) const +{ + // b(d) = (d - d̂)³ + // + // A polynomial: finite at d <= 0, so unlike the log barriers it needs no + // penetration guard. + return select_lazy( + d < dhat, [&] { return -T(2) / T(3) / dhat * cubic(d - dhat); }, + [&] { return T(0); }); +} + +template +T CubicBarrier::first_derivative(const T d, const T dhat) const +{ + return select_lazy( + d < dhat, [&] { return T(-2) / dhat * sqr(d - dhat); }, + [&] { return T(0); }); +} + +template +T CubicBarrier::second_derivative(const T d, const T dhat) const +{ + return select_lazy( + d < dhat, [&] { return T(4) * (T(1) - d / dhat); }, + [&] { return T(0); }); +} + +// ============================================================================ + +template +T TwoStageBarrier::operator()(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + return select_lazy( + d <= T(0), [&] { return infinity(); }, // + d < T(0.5) * dhat, + [&] { return T(-0.25) * sqr(dhat) * (log(T(2) * d / dhat) - T(0.5)); }, + d < dhat, [&] { return T(0.5) * sqr(dhat - d); }, // + [&] { return T(0); }); +} + +template +T TwoStageBarrier::first_derivative(const T d, const T dhat) const +{ + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < T(0.5) * dhat, [&] { return T(-0.25) * sqr(dhat) / d; }, // + d < dhat, [&] { return d - dhat; }, // + [&] { return T(0); }); +} + +template +T TwoStageBarrier::second_derivative(const T d, const T dhat) const +{ + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < T(0.5) * dhat, [&] { return (T(0.25) * sqr(dhat)) / sqr(d); }, // + d < dhat, [&] { return T(1); }, // + [&] { return T(0); }); +} + +// ============================================================================ +// Explicit template instantiations +/// @cond DOXYGEN_SKIP +template class BarrierBase; +template class BarrierBase; +template class ClampedLogBarrier; +template class ClampedLogBarrier; +template class ClampedLogSqBarrier; +template class ClampedLogSqBarrier; +template class CubicBarrier; +template class CubicBarrier; +template class TwoStageBarrier; +template class TwoStageBarrier; +#ifdef IPC_TOOLKIT_WITH_SIMD +template class BarrierBase>; +template class BarrierBase>; +template class ClampedLogBarrier>; +template class ClampedLogBarrier>; +template class ClampedLogSqBarrier>; +template class ClampedLogSqBarrier>; +template class CubicBarrier>; +template class CubicBarrier>; +template class TwoStageBarrier>; +template class TwoStageBarrier>; +#endif +/// @endcond +// ============================================================================ + +} // namespace ipc diff --git a/src/ipc/distance/CMakeLists.txt b/src/ipc/distance/CMakeLists.txt index f327d575f..e24df9796 100644 --- a/src/ipc/distance/CMakeLists.txt +++ b/src/ipc/distance/CMakeLists.txt @@ -1,31 +1,32 @@ set(SOURCES - distance_type.cpp distance_type.hpp - edge_edge.cpp + distance_type_host.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..3c5b7beda 100644 --- a/src/ipc/distance/distance_type.cpp +++ b/src/ipc/distance/distance_type.cpp @@ -1,40 +1,15 @@ #include "distance_type.hpp" #include -#include #include -#include #include -#include namespace ipc::detail { -void warn_degenerate_point_edge() noexcept -{ - logger().warn("Degenerate edge in point_edge_distance_type!"); -} - -void throw_invalid_distance_type(const char* function) -{ - throw std::invalid_argument( - fmt::format("{}: invalid distance type", function)); -} - -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)); -} - template -PointTriangleDistanceType point_triangle_distance_type( +IPC_TOOLKIT_HOST_DEVICE PointTriangleDistanceType point_triangle_distance_type( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -95,7 +70,7 @@ PointTriangleDistanceType point_triangle_distance_type( // A more robust implementation of http://geomalgorithms.com/a07-_distance.html template -EdgeEdgeDistanceType edge_edge_distance_type( +IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_distance_type( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -199,7 +174,7 @@ EdgeEdgeDistanceType edge_edge_distance_type( } template -EdgeEdgeDistanceType edge_edge_parallel_distance_type( +IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_parallel_distance_type( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -237,12 +212,14 @@ EdgeEdgeDistanceType edge_edge_parallel_distance_type( } // clang-format off +#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS template PointTriangleDistanceType point_triangle_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); template PointTriangleDistanceType point_triangle_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); template EdgeEdgeDistanceType edge_edge_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); template EdgeEdgeDistanceType edge_edge_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); template EdgeEdgeDistanceType edge_edge_parallel_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); template EdgeEdgeDistanceType edge_edge_parallel_distance_type(Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef, Eigen::ConstRef); +#endif // clang-format on } // 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 #include #include @@ -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 - inline PointEdgeDistanceType point_edge_distance_type( + IPC_TOOLKIT_HOST_DEVICE inline PointEdgeDistanceType + point_edge_distance_type( const Eigen::Vector& p, const Eigen::Vector& e0, const Eigen::Vector& 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 - PointTriangleDistanceType point_triangle_distance_type( + IPC_TOOLKIT_HOST_DEVICE PointTriangleDistanceType + point_triangle_distance_type( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - EdgeEdgeDistanceType edge_edge_distance_type( + IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType edge_edge_distance_type( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - EdgeEdgeDistanceType edge_edge_parallel_distance_type( + IPC_TOOLKIT_HOST_DEVICE EdgeEdgeDistanceType + edge_edge_parallel_distance_type( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 -inline PointEdgeDistanceType point_edge_distance_type( +IPC_TOOLKIT_HOST_DEVICE inline PointEdgeDistanceType point_edge_distance_type( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& e1) @@ -190,7 +240,8 @@ template < typename DerivedT1, typename DerivedT2, std::enable_if_t, int> = 0> -inline PointTriangleDistanceType point_triangle_distance_type( +IPC_TOOLKIT_HOST_DEVICE inline PointTriangleDistanceType +point_triangle_distance_type( const Eigen::MatrixBase& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& t1, @@ -212,7 +263,7 @@ template < typename DerivedEB0, typename DerivedEB1, std::enable_if_t, int> = 0> -inline EdgeEdgeDistanceType edge_edge_distance_type( +IPC_TOOLKIT_HOST_DEVICE inline EdgeEdgeDistanceType edge_edge_distance_type( const Eigen::MatrixBase& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& eb0, @@ -234,7 +285,8 @@ template < typename DerivedEB0, typename DerivedEB1, std::enable_if_t, int> = 0> -inline EdgeEdgeDistanceType edge_edge_parallel_distance_type( +IPC_TOOLKIT_HOST_DEVICE inline EdgeEdgeDistanceType +edge_edge_parallel_distance_type( const Eigen::MatrixBase& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& eb0, diff --git a/src/ipc/distance/distance_type_host.cpp b/src/ipc/distance/distance_type_host.cpp new file mode 100644 index 000000000..d3f995e82 --- /dev/null +++ b/src/ipc/distance/distance_type_host.cpp @@ -0,0 +1,41 @@ +#include "distance_type.hpp" + +#include + +#include + +#include + +// Host-only definitions for the error reporting helpers declared in +// distance_type.hpp. +// +// These live apart from distance_type.cpp because that file is compiled as a +// CUDA translation unit when CUDA is enabled, and nvcc cannot parse spdlog. +// Keeping the logger and the exceptions here means the shared device sources +// never have to include either one. The inline wrappers in the header pick +// between these and a device-side trap. +namespace ipc::detail { + +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)); +} + +} // namespace ipc::detail 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 -T edge_edge_distance( +IPC_TOOLKIT_HOST_DEVICE T edge_edge_distance( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -59,7 +59,7 @@ T edge_edge_distance( } template -Eigen::Vector edge_edge_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE Eigen::Vector edge_edge_distance_gradient( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -143,7 +143,7 @@ Eigen::Vector edge_edge_distance_gradient( } template -Eigen::Matrix edge_edge_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix edge_edge_distance_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -290,14 +290,18 @@ Eigen::Matrix edge_edge_distance_hessian( Eigen::ConstRef>, \ Eigen::ConstRef>, 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); IPC_INSTANTIATE_EDGE_EDGE(SimdBatch); 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 #include #include @@ -15,7 +16,7 @@ namespace detail { /// @param dtype The point edge distance type to compute. /// @return The distance between the two edges. template - T edge_edge_distance( + IPC_TOOLKIT_HOST_DEVICE T edge_edge_distance( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - Eigen::Vector edge_edge_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector edge_edge_distance_gradient( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - Eigen::Matrix edge_edge_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix edge_edge_distance_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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 - Eigen::Vector edge_edge_mollifier_gradient_wrt_x( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector + edge_edge_mollifier_gradient_wrt_x( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> eb0_rest, @@ -36,7 +37,8 @@ namespace detail { } template - Eigen::Matrix edge_edge_mollifier_gradient_jacobian_wrt_x( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + edge_edge_mollifier_gradient_jacobian_wrt_x( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> eb0_rest, @@ -90,8 +92,10 @@ namespace detail { Eigen::ConstRef>, \ Eigen::ConstRef>) +#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); IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER(SimdBatch); @@ -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 - 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 - 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 - 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[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); IPC_INSTANTIATE_EDGE_EDGE_MOLLIFIER_AUTOGEN(SimdBatch); 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 #include #include #include @@ -10,13 +11,13 @@ namespace ipc { namespace autogen { // clang-format off template - 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 - 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 - 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(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 inline T edge_edge_mollifier(const T x, const T eps_x) +template +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 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 -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 -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 -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 -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 - inline T edge_edge_cross_squarednorm( + IPC_TOOLKIT_HOST_DEVICE inline T edge_edge_cross_squarednorm( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -114,7 +119,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_cross_squarednorm_gradient front end. template - inline Eigen::Vector edge_edge_cross_squarednorm_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + edge_edge_cross_squarednorm_gradient( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -129,7 +135,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_cross_squarednorm_hessian front end. template - inline Eigen::Matrix edge_edge_cross_squarednorm_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_cross_squarednorm_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -144,7 +151,7 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier front end. template - inline T edge_edge_mollifier( + IPC_TOOLKIT_HOST_DEVICE inline T edge_edge_mollifier( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -157,7 +164,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_gradient front end. template - inline Eigen::Vector edge_edge_mollifier_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + edge_edge_mollifier_gradient( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -178,7 +186,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_hessian front end. template - inline Eigen::Matrix edge_edge_mollifier_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_mollifier_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -205,7 +214,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_gradient_wrt_x front end. template - Eigen::Vector edge_edge_mollifier_gradient_wrt_x( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector + edge_edge_mollifier_gradient_wrt_x( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> eb0_rest, @@ -218,7 +228,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_gradient_jacobian_wrt_x front /// end. template - Eigen::Matrix edge_edge_mollifier_gradient_jacobian_wrt_x( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + edge_edge_mollifier_gradient_jacobian_wrt_x( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> eb0_rest, @@ -230,7 +241,7 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_threshold front end. template - T edge_edge_mollifier_threshold( + IPC_TOOLKIT_HOST_DEVICE T edge_edge_mollifier_threshold( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> eb0_rest, @@ -242,7 +253,8 @@ namespace detail { /// @note Prefer the ipc::edge_edge_mollifier_threshold_gradient front end. template - Eigen::Vector edge_edge_mollifier_threshold_gradient( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector + edge_edge_mollifier_threshold_gradient( Eigen::ConstRef> ea0_rest, Eigen::ConstRef> ea1_rest, Eigen::ConstRef> 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0_rest, const Eigen::MatrixBase& ea1_rest, const Eigen::MatrixBase& 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& ea0_rest, const Eigen::MatrixBase& ea1_rest, const Eigen::MatrixBase& 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& ea0_rest, const Eigen::MatrixBase& ea1_rest, const Eigen::MatrixBase& 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& ea0_rest, const Eigen::MatrixBase& ea1_rest, const Eigen::MatrixBase& 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 -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 -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[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); 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 #include #include @@ -10,10 +11,10 @@ namespace ipc { namespace autogen { // clang-format off template - 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 - 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 - inline T line_line_distance( + IPC_TOOLKIT_HOST_DEVICE inline T line_line_distance( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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 - Eigen::Matrix point_edge_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_edge_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1, @@ -68,10 +69,12 @@ namespace detail { Eigen::ConstRef>, \ Eigen::ConstRef>, 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, 2); IPC_INSTANTIATE_POINT_EDGE_DISTANCE_HESSIAN(SimdBatch, 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 #include #include #include @@ -21,7 +22,7 @@ namespace detail { /// @param dtype The point edge distance type to compute. /// @return The distance between the point and edge. template - inline T point_edge_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_edge_distance( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 - inline Eigen::Vector point_edge_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_edge_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 - Eigen::Matrix point_edge_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_edge_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 -inline auto point_edge_distance( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_edge_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance_gradient( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_edge_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_distance_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -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 -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[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); IPC_INSTANTIATE_POINT_PLANE_AUTOGEN(SimdBatch); 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 #include #include @@ -9,10 +10,10 @@ namespace ipc { namespace autogen { // clang-format off template - 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 - 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 - inline T point_plane_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_plane_distance( Eigen::ConstRef> p, Eigen::ConstRef> origin, Eigen::ConstRef> normal) @@ -42,7 +43,7 @@ namespace detail { /// @param t2 The third vertex of the triangle. /// @return The distance between the point and plane. template - inline T point_plane_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_plane_distance( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -61,7 +62,8 @@ namespace detail { /// @param normal The normal of the plane. /// @return The gradient of the distance wrt p. template - inline Eigen::Vector3 point_plane_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 + point_plane_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> origin, Eigen::ConstRef> 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 - inline Eigen::Vector point_plane_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_plane_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -98,7 +101,8 @@ namespace detail { /// @param normal The normal of the plane. /// @return The hessian of the distance wrt p. template - inline Eigen::Matrix3 point_plane_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix3 + point_plane_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> origin, Eigen::ConstRef> 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 - inline Eigen::Matrix point_plane_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_plane_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -135,7 +140,7 @@ namespace detail { /// @param normal The normal of the plane. /// @return The distance between the point and plane. template -inline auto point_plane_distance( +IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance( const Eigen::MatrixBase& p, const Eigen::MatrixBase& origin, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 -inline auto point_plane_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_gradient( const Eigen::MatrixBase& p, const Eigen::MatrixBase& origin, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 -inline auto point_plane_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_plane_distance_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& origin, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 #include #include @@ -15,7 +16,7 @@ namespace detail { /// @param p1 The second point. /// @return The distance between p0 and p1. template - inline T point_point_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_point_distance( Eigen::ConstRef> p0, Eigen::ConstRef> p1) { @@ -31,7 +32,8 @@ namespace detail { /// @param p1 The second point. /// @return The computed gradient. template - inline Eigen::Vector point_point_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_point_distance_gradient( Eigen::ConstRef> p0, Eigen::ConstRef> p1) { @@ -50,7 +52,8 @@ namespace detail { /// @param p1 The second point. /// @return The computed hessian. template - inline Eigen::Matrix point_point_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_point_distance_hessian( Eigen::ConstRef> /*p0*/, Eigen::ConstRef> /*p1*/) { @@ -73,7 +76,7 @@ namespace detail { /// @param p1 The second point. /// @return The distance between p0 and p1. template -inline auto point_point_distance( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance( const Eigen::MatrixBase& p0, const Eigen::MatrixBase& p1) { @@ -98,7 +101,7 @@ inline auto point_point_distance( /// @param p1 The second point. /// @return The computed gradient. template -inline auto point_point_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance_gradient( const Eigen::MatrixBase& p0, const Eigen::MatrixBase& p1) { @@ -124,7 +127,7 @@ inline auto point_point_distance_gradient( /// @param p1 The second point. /// @return The computed hessian. template -inline auto point_point_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_distance_hessian( const Eigen::MatrixBase& p0, const Eigen::MatrixBase& 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 -T point_triangle_distance( +IPC_TOOLKIT_HOST_DEVICE T point_triangle_distance( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -53,7 +53,7 @@ T point_triangle_distance( } template -Eigen::Vector point_triangle_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE Eigen::Vector point_triangle_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -124,7 +124,8 @@ Eigen::Vector point_triangle_distance_gradient( } template -Eigen::Matrix point_triangle_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_triangle_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -253,12 +254,16 @@ Eigen::Matrix point_triangle_distance_hessian( Eigen::ConstRef>, \ Eigen::ConstRef>, 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 #include #include @@ -15,7 +16,7 @@ namespace detail { /// @param dtype The point-triangle distance type to compute. /// @return The distance between the point and triangle. template - T point_triangle_distance( + IPC_TOOLKIT_HOST_DEVICE T point_triangle_distance( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - Eigen::Vector point_triangle_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector + point_triangle_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - Eigen::Matrix point_triangle_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_triangle_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 -Eigen::Matrix line_line_signed_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +line_line_signed_distance_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -28,8 +29,18 @@ Eigen::Matrix 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> hess_n_3_144( + hess_n.data()); + const Eigen::Matrix contracted = + (hess_n_3_144.transpose() * v).eval(); + hess = Eigen::Map>(contracted.data()); + } // --------------------------------------------------------- // 2. Add Jacobian Terms (Product Rule Corrections) @@ -80,8 +91,10 @@ Eigen::Matrix line_line_signed_distance_hessian( Eigen::ConstRef>, \ Eigen::ConstRef>) +#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); IPC_INSTANTIATE_LINE_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch); 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 #include #include @@ -26,7 +27,7 @@ namespace detail { /// represent segment endpoints). Behavior is undefined if ea0 == ea1 or eb0 /// == eb1. template - inline T line_line_signed_distance( + IPC_TOOLKIT_HOST_DEVICE inline T line_line_signed_distance( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -55,7 +56,8 @@ namespace detail { /// /// @see line_line_signed_distance, line_line_normal template - inline Eigen::Vector line_line_signed_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + line_line_signed_distance_gradient( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -90,7 +92,8 @@ namespace detail { /// /// @see line_line_signed_distance, line_line_signed_distance_gradient template - Eigen::Matrix line_line_signed_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + line_line_signed_distance_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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 -Eigen::Matrix point_line_signed_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_line_signed_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -25,8 +26,18 @@ Eigen::Matrix 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> hess_n_2_36( + hess_n.data()); + const Eigen::Matrix contracted = + (hess_n_2_36.transpose() * v).eval(); + hess = Eigen::Map>(contracted.data()); + } // --------------------------------------------------------- // 2. Add Jacobian Terms (Product Rule Corrections) @@ -80,8 +91,10 @@ Eigen::Matrix point_line_signed_distance_hessian( Eigen::ConstRef>, \ Eigen::ConstRef>) +#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); IPC_INSTANTIATE_POINT_LINE_SIGNED_DISTANCE_HESSIAN(SimdBatch); 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 #include #include @@ -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 - inline T point_line_signed_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_line_signed_distance( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 - inline Eigen::Vector point_line_signed_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_line_signed_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 - Eigen::Matrix point_line_signed_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_line_signed_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 -inline auto point_line_signed_distance( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_signed_distance_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance_gradient( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_signed_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_signed_distance_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -Eigen::Matrix point_plane_signed_distance_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_plane_signed_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -42,9 +43,19 @@ Eigen::Matrix 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> hess_n_3_81( + hess_n.data()); + const Eigen::Matrix contracted = + (hess_n_3_81.transpose() * v).eval(); + hess.template block<9, 9>(3, 3) = + Eigen::Map>(contracted.data()); + } // B. Subtract first derivative terms (Product Rule corrections) // Extract 3x3 Jacobian blocks for t0, t1, t2 @@ -79,8 +90,10 @@ Eigen::Matrix point_plane_signed_distance_hessian( Eigen::ConstRef>, \ Eigen::ConstRef>) +#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); IPC_INSTANTIATE_POINT_PLANE_SIGNED_DISTANCE_HESSIAN(SimdBatch); 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 #include #include @@ -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 - inline T point_plane_signed_distance( + IPC_TOOLKIT_HOST_DEVICE inline T point_plane_signed_distance( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - inline Eigen::Vector point_plane_signed_distance_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_plane_signed_distance_gradient( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - Eigen::Matrix point_plane_signed_distance_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_plane_signed_distance_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 #include -#include 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 struct DihedralNormalJacobians { + /// @brief Jacobian of n0 = normal(x0, x1, x2) w.r.t. all 12 DOFs. + Eigen::Matrix dn0_dx; + /// @brief Jacobian of n1 = normal(x1, x0, x3) w.r.t. all 12 DOFs. + Eigen::Matrix 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 - inline std::pair, Eigen::Matrix> + IPC_TOOLKIT_HOST_DEVICE inline DihedralNormalJacobians dihedral_normal_jacobians( Eigen::ConstRef> x0, Eigen::ConstRef> 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 dn1 = triangle_normal_jacobian(x1, x0, x3); Eigen::Matrix dn1_dx; - const std::array 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 -T dihedral_angle( +IPC_TOOLKIT_HOST_DEVICE T dihedral_angle( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> x2, @@ -56,7 +76,7 @@ T dihedral_angle( } template -Eigen::Vector dihedral_angle_gradient( +IPC_TOOLKIT_HOST_DEVICE Eigen::Vector dihedral_angle_gradient( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> x2, @@ -90,7 +110,7 @@ Eigen::Vector dihedral_angle_gradient( } template -Eigen::Matrix dihedral_angle_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix dihedral_angle_hessian( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> x2, @@ -368,8 +388,10 @@ Eigen::Matrix dihedral_angle_hessian( template Eigen::Vector dihedral_angle_gradient(Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>); \ template Eigen::Matrix dihedral_angle_hessian(Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>) +#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS IPC_INSTANTIATE_ANGLE(float); IPC_INSTANTIATE_ANGLE(double); +#endif #ifdef IPC_TOOLKIT_WITH_SIMD IPC_INSTANTIATE_ANGLE(SimdBatch); IPC_INSTANTIATE_ANGLE(SimdBatch); 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 #include namespace ipc { @@ -7,7 +8,7 @@ namespace ipc { namespace detail { /// @note Prefer the ipc::dihedral_angle front end. template - T dihedral_angle( + IPC_TOOLKIT_HOST_DEVICE T dihedral_angle( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> x2, @@ -15,7 +16,7 @@ namespace detail { /// @note Prefer the ipc::dihedral_angle_gradient front end. template - Eigen::Vector dihedral_angle_gradient( + IPC_TOOLKIT_HOST_DEVICE Eigen::Vector dihedral_angle_gradient( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> x2, @@ -23,7 +24,7 @@ namespace detail { /// @note Prefer the ipc::dihedral_angle_hessian front end. template - Eigen::Matrix dihedral_angle_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix dihedral_angle_hessian( Eigen::ConstRef> x0, Eigen::ConstRef> x1, Eigen::ConstRef> 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& x0, const Eigen::MatrixBase& x1, const Eigen::MatrixBase& 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& x0, const Eigen::MatrixBase& x1, const Eigen::MatrixBase& 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& x0, const Eigen::MatrixBase& x1, const Eigen::MatrixBase& 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 -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[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); IPC_INSTANTIATE_AREA_AUTOGEN(SimdBatch); 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 #include #include @@ -15,7 +16,7 @@ namespace autogen { // clang-format off /// dA is (9×1) flattened in column-major order template - 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 - inline T edge_length( + IPC_TOOLKIT_HOST_DEVICE inline T edge_length( Eigen::ConstRef> e0, Eigen::ConstRef> 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 - inline Eigen::Vector edge_length_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + edge_length_gradient( Eigen::ConstRef> e0, Eigen::ConstRef> e1) { @@ -69,7 +71,7 @@ namespace detail { /// @param t2 The third vertex of the triangle. /// @return The area of the triangle. template - inline T triangle_area( + IPC_TOOLKIT_HOST_DEVICE inline T triangle_area( Eigen::ConstRef> t0, Eigen::ConstRef> t1, Eigen::ConstRef> 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 - inline Eigen::Vector triangle_area_gradient( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector triangle_area_gradient( Eigen::ConstRef> t0, Eigen::ConstRef> t1, Eigen::ConstRef> t2) @@ -102,7 +104,7 @@ namespace detail { /// @param e1 The second vertex of the edge. /// @return The length of the edge. template -inline auto edge_length( +IPC_TOOLKIT_HOST_DEVICE inline auto edge_length( const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto edge_length_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto edge_length_gradient( const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto triangle_area( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_area( const Eigen::MatrixBase& t0, const Eigen::MatrixBase& t1, const Eigen::MatrixBase& 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 -inline auto triangle_area_gradient( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_area_gradient( const Eigen::MatrixBase& t0, const Eigen::MatrixBase& t1, const Eigen::MatrixBase& 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 #include @@ -8,7 +8,7 @@ namespace ipc::detail { // --- point-line normal functions ------------------------------------------- template -VectorMax3 point_line_unnormalized_normal( +IPC_TOOLKIT_HOST_DEVICE VectorMax3 point_line_unnormalized_normal( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -31,7 +31,8 @@ VectorMax3 point_line_unnormalized_normal( } template -MatrixMax point_line_unnormalized_normal_jacobian( +IPC_TOOLKIT_HOST_DEVICE MatrixMax +point_line_unnormalized_normal_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -67,7 +68,8 @@ MatrixMax point_line_unnormalized_normal_jacobian( } template -MatrixMax point_line_unnormalized_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE MatrixMax +point_line_unnormalized_normal_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -173,7 +175,7 @@ MatrixMax point_line_unnormalized_normal_hessian( } template -MatrixMax point_line_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE MatrixMax point_line_normal_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -217,7 +219,7 @@ MatrixMax point_line_normal_hessian( namespace { template - void set_cross_product_matrix_jacobian( + IPC_TOOLKIT_HOST_DEVICE void set_cross_product_matrix_jacobian( Eigen::Ref> 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 Eigen::Matrix cross_product_matrix_jacobian() +template +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix cross_product_matrix_jacobian() { Eigen::Matrix J = Eigen::Matrix::Zero(); J(2, 1) = J(3, 2) = J(7, 0) = T(-1); @@ -234,7 +237,8 @@ template Eigen::Matrix cross_product_matrix_jacobian() } template -Eigen::Matrix triangle_unnormalized_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +triangle_unnormalized_normal_hessian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> c) @@ -263,7 +267,7 @@ Eigen::Matrix triangle_unnormalized_normal_hessian( } template -Eigen::Matrix triangle_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix triangle_normal_hessian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> c) @@ -300,7 +304,8 @@ Eigen::Matrix triangle_normal_hessian( // --- line-line normal functions --------------------------------------------- template -Eigen::Matrix line_line_unnormalized_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +line_line_unnormalized_normal_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -340,7 +345,7 @@ Eigen::Matrix line_line_unnormalized_normal_hessian( } template -Eigen::Matrix line_line_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix line_line_normal_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -390,8 +395,10 @@ Eigen::Matrix line_line_normal_hessian( template Eigen::Matrix line_line_unnormalized_normal_hessian(Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>); \ template Eigen::Matrix line_line_normal_hessian(Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>, Eigen::ConstRef>) +#if IPC_TOOLKIT_INSTANTIATE_DEVICE_SCALARS IPC_INSTANTIATE_NORMAL(float); IPC_INSTANTIATE_NORMAL(double); +#endif #ifdef IPC_TOOLKIT_WITH_SIMD IPC_INSTANTIATE_NORMAL(SimdBatch); IPC_INSTANTIATE_NORMAL(SimdBatch); 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 #include #include @@ -21,8 +22,7 @@ struct NormalizationAndJacobian { Eigen::Matrix 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 - inline NormalizationAndJacobian + IPC_TOOLKIT_HOST_DEVICE inline NormalizationAndJacobian normalization_and_jacobian(const Eigen::Vector& 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 - inline NormalizationAndJacobianAndHessian + IPC_TOOLKIT_HOST_DEVICE inline NormalizationAndJacobianAndHessian normalization_and_jacobian_and_hessian(const Eigen::Vector& 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 -inline auto normalization_and_jacobian(const Eigen::MatrixBase& x) +IPC_TOOLKIT_HOST_DEVICE inline auto +normalization_and_jacobian(const Eigen::MatrixBase& x) { using T = typename DerivedX::Scalar; using DynamicResult = NormalizationAndJacobian; @@ -121,7 +122,8 @@ inline auto normalization_and_jacobian(const Eigen::MatrixBase& x) /// @param x The input vector. /// @return The Jacobian of the normalization operation. template -inline auto normalization_jacobian(const Eigen::MatrixBase& x) +IPC_TOOLKIT_HOST_DEVICE inline auto +normalization_jacobian(const Eigen::MatrixBase& x) { using T = typename DerivedX::Scalar; @@ -143,7 +145,7 @@ inline auto normalization_jacobian(const Eigen::MatrixBase& x) /// @param x The input vector. /// @return The normalized vector, its Jacobian, and its Hessian. template -inline auto +IPC_TOOLKIT_HOST_DEVICE inline auto normalization_and_jacobian_and_hessian(const Eigen::MatrixBase& x) { using T = typename DerivedX::Scalar; @@ -180,7 +182,8 @@ normalization_and_jacobian_and_hessian(const Eigen::MatrixBase& x) /// @param x The input vector. /// @return The Hessian of the normalization operation. template -inline auto normalization_hessian(const Eigen::MatrixBase& x) +IPC_TOOLKIT_HOST_DEVICE inline auto +normalization_hessian(const Eigen::MatrixBase& 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 - inline Eigen::Matrix3 + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix3 cross_product_matrix(Eigen::ConstRef> v) { Eigen::Matrix3 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 - Eigen::Matrix cross_product_matrix_jacobian(); + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + cross_product_matrix_jacobian(); // ========================================================================= @@ -224,7 +228,7 @@ namespace detail { /// @param e1 The end position of the line. /// @return The unnormalized normal vector. template - VectorMax3 point_line_unnormalized_normal( + IPC_TOOLKIT_HOST_DEVICE VectorMax3 point_line_unnormalized_normal( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1); @@ -235,7 +239,7 @@ namespace detail { /// @param e1 The end position of the line. /// @return The normal vector. template - inline VectorMax3 point_line_normal( + IPC_TOOLKIT_HOST_DEVICE inline VectorMax3 point_line_normal( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -250,7 +254,8 @@ namespace detail { /// @param e1 The end position of the line. /// @return The Jacobian of the unnormalized normal vector. template - MatrixMax point_line_unnormalized_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE MatrixMax + point_line_unnormalized_normal_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1); @@ -262,7 +267,8 @@ namespace detail { /// @return The Hessian of the unnormalized normal vector of the point-line /// pair. template - MatrixMax point_line_unnormalized_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE MatrixMax + point_line_unnormalized_normal_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1); @@ -273,7 +279,8 @@ namespace detail { /// @param e1 The end position of the line. /// @return The Jacobian of the normal vector. template - inline MatrixMax point_line_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline MatrixMax + point_line_normal_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -289,7 +296,7 @@ namespace detail { /// @param e1 The end position of the line. /// @return The Hessian of the normal vector. template - MatrixMax point_line_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE MatrixMax point_line_normal_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1); @@ -310,7 +317,8 @@ namespace detail { /// @param c The third vertex of the triangle. /// @return The unnormalized normal vector of the triangle. template - inline Eigen::Vector3 triangle_unnormalized_normal( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 + triangle_unnormalized_normal( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> c) @@ -324,7 +332,7 @@ namespace detail { /// @param c The third vertex of the triangle. /// @return The normal vector of the triangle. template - inline Eigen::Vector3 triangle_normal( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 triangle_normal( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> 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 - inline Eigen::Matrix triangle_unnormalized_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + triangle_unnormalized_normal_jacobian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> 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 - Eigen::Matrix triangle_unnormalized_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + triangle_unnormalized_normal_hessian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> 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 - inline Eigen::Matrix triangle_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + triangle_normal_jacobian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> 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 - Eigen::Matrix triangle_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix triangle_normal_hessian( Eigen::ConstRef> a, Eigen::ConstRef> b, Eigen::ConstRef> 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 - inline Eigen::Vector3 line_line_unnormalized_normal( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 + line_line_unnormalized_normal( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - inline Eigen::Vector3 line_line_normal( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 line_line_normal( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - inline Eigen::Matrix line_line_unnormalized_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + line_line_unnormalized_normal_jacobian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - inline Eigen::Matrix line_line_normal_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + line_line_normal_jacobian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - Eigen::Matrix line_line_unnormalized_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + line_line_unnormalized_normal_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - Eigen::Matrix line_line_normal_hessian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix line_line_normal_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 -inline Eigen::Matrix cross_product_matrix_jacobian() +IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix +cross_product_matrix_jacobian() { return detail::cross_product_matrix_jacobian(); } @@ -518,7 +534,8 @@ inline Eigen::Matrix cross_product_matrix_jacobian() /// @param v Vector to create the cross product matrix for. /// @return The cross product matrix of the vector. template -inline auto cross_product_matrix(const Eigen::MatrixBase& v) +IPC_TOOLKIT_HOST_DEVICE inline auto +cross_product_matrix(const Eigen::MatrixBase& v) { using T = typename DerivedX::Scalar; return detail::cross_product_matrix(v); @@ -530,7 +547,7 @@ inline auto cross_product_matrix(const Eigen::MatrixBase& v) /// @param e1 The end position of the line. /// @return The unnormalized normal vector. template -inline auto point_line_unnormalized_normal( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& e1) @@ -545,7 +562,7 @@ inline auto point_line_unnormalized_normal( /// @param e1 The end position of the line. /// @return The normal vector. template -inline auto point_line_normal( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_unnormalized_normal_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal_jacobian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_unnormalized_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_unnormalized_normal_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_normal_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal_jacobian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_line_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_line_normal_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto triangle_unnormalized_normal( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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 -inline auto triangle_normal( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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 -inline auto triangle_unnormalized_normal_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal_jacobian( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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 -inline auto triangle_unnormalized_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_unnormalized_normal_hessian( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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 -inline auto triangle_normal_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal_jacobian( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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 -inline auto triangle_normal_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto triangle_normal_hessian( const Eigen::MatrixBase& a, const Eigen::MatrixBase& b, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& eb0, diff --git a/src/ipc/math/morton.hpp b/src/ipc/math/morton.hpp index ccc201cda..02b9c5619 100644 --- a/src/ipc/math/morton.hpp +++ b/src/ipc/math/morton.hpp @@ -1,14 +1,31 @@ #pragma once -#include // for std::clamp -#include // for uint64_t +#include // for IPC_TOOLKIT_HOST_DEVICE + +#include // for uint64_t namespace ipc { +namespace detail { + /// @brief Clamp a double to [lo, hi] (a device-safe std::clamp). + /// + /// std::clamp cannot be called from device code: MSVC's debug STL checks + /// the bounds with _STL_VERIFY, which expands to __debugbreak(); 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 input passing through unchanged -- is + /// identical. + IPC_TOOLKIT_HOST_DEVICE inline double + clamp(const double v, const double lo, const double hi) + { + return v < lo ? lo : (hi < v ? hi : v); + } +} // namespace detail + /// @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 +38,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 +52,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 = detail::clamp(x * scale, 0.0, scale - 1); + y = detail::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 +67,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 = detail::clamp(x * scale, 0.0, scale - 1); + y = detail::clamp(y * scale, 0.0, scale - 1); + z = detail::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 && ...)>> \ - 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 inline T sqr(const T& x) { return x * x; } +template 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 inline T cubic(const T& x) { return x * x * x; } +template 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 -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 -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 -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 -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 -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 -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 -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 -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[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); IPC_INSTANTIATE_CLOSEST_POINT_AUTOGEN(SimdBatch); 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 #include #include @@ -13,12 +14,12 @@ namespace ipc { namespace autogen { /// hess is (6×6) flattened in column-major order template - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - inline Eigen::Vector2 solve_spd_2x2( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2 solve_spd_2x2( Eigen::ConstRef> A, Eigen::ConstRef> b) { @@ -215,7 +216,7 @@ namespace detail { /// @param e1 Second edge point /// @return Barycentric coordinate of the closest point template - inline T point_edge_closest_point( + IPC_TOOLKIT_HOST_DEVICE inline T point_edge_closest_point( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -233,7 +234,8 @@ namespace detail { /// @param e1 Second edge point /// @return Jacobian of the closest point template - inline Eigen::Vector point_edge_closest_point_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_edge_closest_point_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -261,7 +263,8 @@ namespace detail { /// @param e1 Second edge point /// @return Hessian of the closest point template - inline Eigen::Matrix point_edge_closest_point_hessian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_edge_closest_point_hessian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -291,7 +294,7 @@ namespace detail { /// @param eb1 Second point of the second edge /// @return Barycentric coordinates of the closest points template - inline Eigen::Vector2 edge_edge_closest_point( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2 edge_edge_closest_point( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -321,7 +324,8 @@ namespace detail { /// @param eb1 Second point of the second edge /// @return Jacobian of the closest points template - inline Eigen::Matrix edge_edge_closest_point_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_closest_point_jacobian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -342,7 +346,7 @@ namespace detail { /// @param eb1 Second point of the second edge /// @return Hessian of the closest points (2x12x12 tensor) template - inline std::array, 2> + IPC_TOOLKIT_HOST_DEVICE inline std::array, 2> edge_edge_closest_point_hessian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, @@ -370,7 +374,8 @@ namespace detail { /// @param t2 Triangle's third vertex /// @return Barycentric coordinates of the closest point template - inline Eigen::Vector2 point_triangle_closest_point( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector2 + point_triangle_closest_point( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -392,7 +397,8 @@ namespace detail { /// @param t2 Triangle's third vertex /// @return Jacobian of the closest point template - inline Eigen::Matrix point_triangle_closest_point_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_triangle_closest_point_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -413,7 +419,7 @@ namespace detail { /// @param t2 Triangle's third vertex /// @return Hessian of the closest point (2x12x12 tensor) template - inline std::array, 2> + IPC_TOOLKIT_HOST_DEVICE inline std::array, 2> point_triangle_closest_point_hessian( Eigen::ConstRef> p, Eigen::ConstRef> t0, @@ -440,7 +446,7 @@ namespace detail { /// @param e1 Second edge point /// @return barycentric coordinates of the closest point template -inline auto point_edge_closest_point( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& e1) @@ -465,7 +471,7 @@ inline auto point_edge_closest_point( /// @param e1 Second edge point /// @return Jacobian of the closest point template -inline auto point_edge_closest_point_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point_jacobian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& e1) @@ -492,7 +498,7 @@ inline auto point_edge_closest_point_jacobian( /// @param e1 Second edge point /// @return Hessian of the closest point template -inline auto point_edge_closest_point_hessian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_closest_point_hessian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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 #include #include @@ -17,7 +18,8 @@ namespace detail { /// @param dp1 Velocity of the second point /// @return The relative velocity of the two points template - inline Eigen::Vector point_point_relative_velocity( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_point_relative_velocity( Eigen::ConstRef> dp0, Eigen::ConstRef> dp1) { @@ -31,7 +33,7 @@ namespace detail { /// @tparam dim The dimension (2 or 3). /// @return The relative velocity Jacobian du/dx template - inline Eigen::Matrix + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix 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 - inline Eigen::Vector + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector 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 - inline Eigen::Vector point_edge_relative_velocity( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector + point_edge_relative_velocity( Eigen::ConstRef> dp, Eigen::ConstRef> de0, Eigen::ConstRef> 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 - inline Eigen::Matrix + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix 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 - inline Eigen::Vector + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector 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 - inline Eigen::Vector3 edge_edge_relative_velocity( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 + edge_edge_relative_velocity( Eigen::ConstRef> dea0, Eigen::ConstRef> dea1, Eigen::ConstRef> deb0, @@ -151,7 +155,8 @@ namespace detail { /// @param coords Two parametric coordinates of the closest points /// @return The relative velocity Jacobian du/dx template - inline Eigen::Matrix edge_edge_relative_velocity_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_relative_velocity_jacobian( Eigen::ConstRef> coords) { Eigen::Matrix J = Eigen::Matrix::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 - inline Eigen::Matrix edge_edge_relative_velocity_dx_dbeta( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_relative_velocity_dx_dbeta( Eigen::ConstRef> 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 - inline Eigen::Vector3 point_triangle_relative_velocity( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Vector3 + point_triangle_relative_velocity( Eigen::ConstRef> dp, Eigen::ConstRef> dt0, Eigen::ConstRef> dt1, @@ -213,7 +220,8 @@ namespace detail { /// @param coords Barycentric coordinates of the closest point /// @return The relative velocity Jacobian du/dx template - inline Eigen::Matrix point_triangle_relative_velocity_jacobian( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_triangle_relative_velocity_jacobian( Eigen::ConstRef> coords) { Eigen::Matrix J = Eigen::Matrix::Zero(); @@ -231,7 +239,8 @@ namespace detail { /// @param coords Baricentric coordinates of the closest point /// @return The vectorized tensor of d²u/dxdβ template - inline Eigen::Matrix point_triangle_relative_velocity_dx_dbeta( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_triangle_relative_velocity_dx_dbeta( Eigen::ConstRef> 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 -inline auto point_point_relative_velocity( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_relative_velocity( const Eigen::MatrixBase& dp0, const Eigen::MatrixBase& 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 -inline MatrixMax point_point_relative_velocity_jacobian(const int dim) +IPC_TOOLKIT_HOST_DEVICE inline MatrixMax +point_point_relative_velocity_jacobian(const int dim) { if (dim == 2) { return MatrixMax( @@ -299,7 +309,8 @@ inline MatrixMax point_point_relative_velocity_jacobian(const int dim) /// @param dim Dimension (2 or 3) /// @return The vectorized tensor of d²u/dxdβ template -inline VectorMax point_point_relative_velocity_dx_dbeta(const int dim) +IPC_TOOLKIT_HOST_DEVICE inline VectorMax +point_point_relative_velocity_dx_dbeta(const int dim) { if (dim == 2) { return VectorMax( @@ -321,7 +332,7 @@ inline VectorMax 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 -inline auto point_edge_relative_velocity( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_relative_velocity( const Eigen::MatrixBase& dp, const Eigen::MatrixBase& de0, const Eigen::MatrixBase& 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 -inline MatrixMax +IPC_TOOLKIT_HOST_DEVICE inline MatrixMax 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 -inline VectorMax +IPC_TOOLKIT_HOST_DEVICE inline VectorMax 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& dea0, const Eigen::MatrixBase& dea1, const Eigen::MatrixBase& 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 -inline auto edge_edge_relative_velocity_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_relative_velocity_jacobian( const Eigen::MatrixBase& 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 -inline auto edge_edge_relative_velocity_dx_dbeta( +IPC_TOOLKIT_HOST_DEVICE inline auto edge_edge_relative_velocity_dx_dbeta( const Eigen::MatrixBase& 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& dp, const Eigen::MatrixBase& dt0, const Eigen::MatrixBase& 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 -inline auto point_triangle_relative_velocity_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_relative_velocity_jacobian( const Eigen::MatrixBase& 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 -inline auto point_triangle_relative_velocity_dx_dbeta( +IPC_TOOLKIT_HOST_DEVICE inline auto point_triangle_relative_velocity_dx_dbeta( const Eigen::MatrixBase& 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 -Eigen::Matrix point_point_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_point_tangent_basis_jacobian( Eigen::ConstRef> p0, Eigen::ConstRef> p1) { @@ -32,7 +33,8 @@ Eigen::Matrix point_point_tangent_basis_jacobian( // Point - Edge template -Eigen::Matrix point_edge_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_edge_tangent_basis_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -56,7 +58,8 @@ Eigen::Matrix point_edge_tangent_basis_jacobian( // Edge - Edge template -Eigen::Matrix edge_edge_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +edge_edge_tangent_basis_jacobian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -73,7 +76,8 @@ Eigen::Matrix edge_edge_tangent_basis_jacobian( // Point - Triangle template -Eigen::Matrix point_triangle_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix +point_triangle_tangent_basis_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -107,10 +111,12 @@ Eigen::Matrix point_triangle_tangent_basis_jacobian( Eigen::ConstRef>, \ Eigen::ConstRef>) +#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, 2); IPC_INSTANTIATE_TANGENT_BASIS_ND(SimdBatch, 3); @@ -142,8 +148,10 @@ IPC_INSTANTIATE_TANGENT_BASIS_ND(SimdBatch, 3); Eigen::ConstRef>, \ Eigen::ConstRef>) +#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); IPC_INSTANTIATE_TANGENT_BASIS_3D(SimdBatch); @@ -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 inline T pow_1_5(T x) + template 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 -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 -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 -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 -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 -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 -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[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); IPC_INSTANTIATE_TANGENT_BASIS_AUTOGEN(SimdBatch); 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 #include #include @@ -18,7 +19,8 @@ namespace detail { /// @param p1 Second point /// @return A dim×(dim-1) matrix whose columns are the basis vectors. template - inline Eigen::Matrix point_point_tangent_basis( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_point_tangent_basis( Eigen::ConstRef> p0, Eigen::ConstRef> p1) { @@ -72,7 +74,8 @@ namespace detail { /// @param p1 Second point /// @return A (dim*(dim-1))×(2*dim) matrix. template - Eigen::Matrix point_point_tangent_basis_jacobian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_point_tangent_basis_jacobian( Eigen::ConstRef> p0, Eigen::ConstRef> 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 - inline Eigen::Matrix point_edge_tangent_basis( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_edge_tangent_basis( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> e1) @@ -115,7 +119,8 @@ namespace detail { /// @param e1 Second edge point /// @return A (dim*(dim-1))×(3*dim) matrix. template - Eigen::Matrix point_edge_tangent_basis_jacobian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_edge_tangent_basis_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> e0, Eigen::ConstRef> 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 - inline Eigen::Matrix edge_edge_tangent_basis( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + edge_edge_tangent_basis( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> 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 - Eigen::Matrix edge_edge_tangent_basis_jacobian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + edge_edge_tangent_basis_jacobian( Eigen::ConstRef> ea0, Eigen::ConstRef> ea1, Eigen::ConstRef> eb0, @@ -178,7 +185,8 @@ namespace detail { /// @param t2 Triangle's third vertex /// @return A 3x2 matrix whose columns are the basis vectors. template - inline Eigen::Matrix point_triangle_tangent_basis( + IPC_TOOLKIT_HOST_DEVICE inline Eigen::Matrix + point_triangle_tangent_basis( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> 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 - Eigen::Matrix point_triangle_tangent_basis_jacobian( + IPC_TOOLKIT_HOST_DEVICE Eigen::Matrix + point_triangle_tangent_basis_jacobian( Eigen::ConstRef> p, Eigen::ConstRef> t0, Eigen::ConstRef> t1, @@ -224,7 +233,7 @@ namespace detail { /// @param p1 Second point /// @return A 3x2 matrix whose columns are the basis vectors. template -inline auto point_point_tangent_basis( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_tangent_basis( const Eigen::MatrixBase& p0, const Eigen::MatrixBase& 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 -inline auto point_point_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_point_tangent_basis_jacobian( const Eigen::MatrixBase& p0, const Eigen::MatrixBase& 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 -inline auto point_edge_tangent_basis( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_tangent_basis( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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 -inline auto point_edge_tangent_basis_jacobian( +IPC_TOOLKIT_HOST_DEVICE inline auto point_edge_tangent_basis_jacobian( const Eigen::MatrixBase& p, const Eigen::MatrixBase& e0, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& ea0, const Eigen::MatrixBase& ea1, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& 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& p, const Eigen::MatrixBase& t0, const Eigen::MatrixBase& t1, @@ -437,22 +446,22 @@ namespace autogen { // J is (2×4) flattened in column-major order template - 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 - 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 - 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 - 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 - 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 - 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..7c4dc701f 100644 --- a/src/ipc/utils/simd.hpp +++ b/src/ipc/utils/simd.hpp @@ -4,9 +4,6 @@ #include #include -#ifdef IPC_TOOLKIT_WITH_SIMD -#include -#endif #include #include @@ -40,7 +37,7 @@ template using scalar_of_t = typename ScalarOf::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 inline T literal(const double c) +template IPC_TOOLKIT_HOST_DEVICE inline T literal(const double c) { return T(static_cast>(c)); } @@ -54,19 +51,20 @@ template 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 inline T select(const bool mask, const T& a, const T& b) +template +IPC_TOOLKIT_HOST_DEVICE inline T select(const bool mask, const T& a, const T& b) { return mask ? a : b; } /// @brief `+infinity` for any scalar the library templates on. -template inline T infinity() +template IPC_TOOLKIT_HOST_DEVICE inline T infinity() { return T(std::numeric_limits>::infinity()); } @@ -92,13 +90,15 @@ template 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 inline auto select_lazy(F&& else_value) +template +IPC_TOOLKIT_HOST_DEVICE inline auto select_lazy(F&& else_value) { return else_value(); } template -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, bool>) { return mask ? value() : select_lazy(std::forward(rest)...); @@ -111,6 +111,8 @@ inline auto select_lazy(const Mask& mask, F&& value, Rest&&... rest) #ifdef IPC_TOOLKIT_WITH_SIMD +#include + namespace Eigen { /// @brief Let Eigen treat an xsimd batch as a scalar type. @@ -210,7 +212,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 -inline typename Derived::PlainObject +IPC_TOOLKIT_HOST_DEVICE inline typename Derived::PlainObject normalized(const Eigen::MatrixBase& v) { using T = typename Derived::Scalar; 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..62e7f1319 --- /dev/null +++ b/tests/src/tests/barrier/test_gpu_barrier.cu @@ -0,0 +1,93 @@ +// 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 + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include + +#include +#include + +#include + +#include +#include + +namespace { + +#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) + +__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]") +{ + int device_count = 0; + const cudaError_t err = cudaGetDeviceCount(&device_count); + if (err != cudaSuccess || device_count == 0) { + SKIP("No CUDA device available; kernel compiled but not executed."); + } + + const double dhat = 1e-2; + const std::vector in = { dhat, 0.5 * dhat, 2 * dhat, -1.0 }; + + std::vector 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)); + } + + double *d_in = nullptr, *d_out = nullptr; + REQUIRE_CUDA(cudaMalloc(&d_in, in.size() * sizeof(double))); + REQUIRE_CUDA(cudaMalloc(&d_out, expected.size() * sizeof(double))); + REQUIRE_CUDA(cudaMemcpy( + d_in, in.data(), in.size() * sizeof(double), cudaMemcpyHostToDevice)); + + barrier_kernel<<<1, 1>>>(d_in, d_out); + REQUIRE_CUDA(cudaGetLastError()); + REQUIRE_CUDA(cudaDeviceSynchronize()); + + std::vector out(expected.size()); + REQUIRE_CUDA(cudaMemcpy( + out.data(), d_out, out.size() * sizeof(double), + cudaMemcpyDeviceToHost)); + cudaFree(d_in); + cudaFree(d_out); + + 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)); + } + } +} + +#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..30320ed6a --- /dev/null +++ b/tests/src/tests/distance/test_gpu_distance.cu @@ -0,0 +1,589 @@ +// 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 + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include + +namespace { + +#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) + +/// @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 +__device__ int +write_out(const Eigen::MatrixBase& 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 +void write_expected( + const Eigen::MatrixBase& m, std::vector& expected) +{ + for (int r = 0; r < m.rows(); ++r) { + for (int c = 0; c < m.cols(); ++c) { + expected.push_back(m(r, c)); + } + } +} + +/// @brief Launch `kernel` with the given inputs; returns the device outputs. +std::vector run_gpu_kernel( + GpuTestKernel kernel, const std::vector& 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 out(n_out); + REQUIRE_CUDA(cudaMemcpy( + out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); + cudaFree(d_in); + cudaFree(d_out); + return out; +} + +/// @brief SKIP the test if no CUDA device is present (kernels have still been +/// compiled as device code, which is the primary verification). +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."); + } +} + +void check_gpu_matches_host( + GpuTestKernel kernel, + const std::vector& in, + const std::vector& expected) +{ + const std::vector out = run_gpu_kernel(kernel, in, expected.size()); + for (size_t i = 0; i < expected.size(); ++i) { + CAPTURE(i); + CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); + } +} + +// --------------------------------------------------------------------------- + +// Shared fixtures (kept in one place so kernels and host references agree). +// clang-format off +const std::vector PP_IN = { 0, 0, 0, 1, 2, 2 }; +const std::vector PL_IN = { 0.1, 1, 0.2, -1, 0, 0, 1, 0, 0.1 }; +const std::vector LL_IN = { -1, 0, 0, 1, 0, 0, 0, 1, -1, 0.1, 1, 1 }; +const std::vector PLANE_IN = { 0, 1, 0, 0, 0, 0, 0, 1, 0, + -1, 0, 1, 1, 0, 1, 0, 0, -1 }; +const std::vector PE_IN = { 0.5, 1, 0, -1, 0, 0, 1, 0, 0 }; +const std::vector EE_IN = { -1, 0, 0, 1, 0, 0, -0.5, 1, -1, 0.5, 1, 1 }; +const std::vector PT_IN = { 0, 1, 0, -1, 0, 1, 1, 0, 1, 0, 0, -1 }; +const std::vector 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(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(ipc::edge_edge_distance_type(ea0, ea1, eb0, eb1)); + out[offset++] = static_cast( + 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(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 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 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 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 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 expected; + expected.push_back( + static_cast(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 expected; + expected.push_back( + static_cast(ipc::edge_edge_distance_type(ea0, ea1, eb0, eb1))); + expected.push_back( + static_cast( + 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 expected; + expected.push_back( + static_cast(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 in = EE_IN; + in.push_back(eps_x); + + std::vector 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 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 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..054833eac --- /dev/null +++ b/tests/src/tests/geometry/test_gpu_geometry.cu @@ -0,0 +1,362 @@ +// 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 + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include + +#include +#include + +#include + +#include + +namespace { + +#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) + +/// @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 +__device__ int +write_out(const Eigen::MatrixBase& 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 +void write_expected( + const Eigen::MatrixBase& m, std::vector& expected) +{ + for (int r = 0; r < m.rows(); ++r) { + for (int c = 0; c < m.cols(); ++c) { + expected.push_back(m(r, c)); + } + } +} + +/// @brief Launch `kernel` with the given inputs; returns the device outputs. +/// Skips the surrounding test if no CUDA device is available. +std::vector run_gpu_kernel( + GpuTestKernel kernel, const std::vector& 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 out(n_out); + REQUIRE_CUDA(cudaMemcpy( + out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); + cudaFree(d_in); + cudaFree(d_out); + return out; +} + +/// @brief SKIP the test if no CUDA device is present (kernels have still been +/// compiled as device code, which is the primary verification). +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."); + } +} + +void check_gpu_matches_host( + GpuTestKernel kernel, + const std::vector& in, + const std::vector& expected) +{ + const std::vector out = run_gpu_kernel(kernel, in, expected.size()); + for (size_t i = 0; i < expected.size(); ++i) { + CAPTURE(i); + CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); + } +} + +// --------------------------------------------------------------------------- + +__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(), 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 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 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 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(), 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 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 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 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/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..1f5e6b158 --- /dev/null +++ b/tests/src/tests/tangent/test_gpu_tangent.cu @@ -0,0 +1,354 @@ +// 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 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 + +#ifdef IPC_TOOLKIT_WITH_CUDA + +#include +#include +#include + +#include +#include + +#include + +#include + +namespace { + +#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) + +/// @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 +__device__ int +write_out(const Eigen::MatrixBase& 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 +void write_expected( + const Eigen::MatrixBase& m, std::vector& expected) +{ + for (int r = 0; r < m.rows(); ++r) { + for (int c = 0; c < m.cols(); ++c) { + expected.push_back(m(r, c)); + } + } +} + +/// @brief Launch `kernel` with the given inputs; returns the device outputs. +std::vector run_gpu_kernel( + GpuTestKernel kernel, const std::vector& 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 out(n_out); + REQUIRE_CUDA(cudaMemcpy( + out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); + cudaFree(d_in); + cudaFree(d_out); + return out; +} + +/// @brief SKIP the test if no CUDA device is present (kernels have still been +/// compiled as device code, which is the primary verification). +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."); + } +} + +void check_gpu_matches_host( + GpuTestKernel kernel, + const std::vector& in, + const std::vector& expected) +{ + const std::vector out = run_gpu_kernel(kernel, in, expected.size()); + for (size_t i = 0; i < expected.size(); ++i) { + CAPTURE(i); + CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); + } +} + +// --------------------------------------------------------------------------- + +// Shared fixtures (kept in one place so kernels and host references agree). +// clang-format off +const std::vector PE_IN = { 0.5, 1, 0.2, -1, 0, 0, 1, 0, 0.1 }; +const std::vector EE_IN = { -1, 0, 0, 1, 0, 0, -0.5, 1, -1, 0.5, 1, 1 }; +const std::vector 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 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 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 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 tangent_fixture() +{ + std::vector 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 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 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 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 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 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 From 3825f3c0a83b83a84f60033d7c89e16ef44c64ae Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 8 Sep 2026 16:02:14 -0400 Subject: [PATCH 2/4] Fold the host-only sources back into their shared device siblings barrier_classes.cpp and distance_type_host.cpp existed only to keep code away from nvcc, which left barrier.hpp and distance_type.hpp each implemented across two .cpp files. The dual-compilation mechanism does not require a separate file, only that the host-only code be invisible to the device pass, so both are merged back behind `#ifndef __CUDACC__`. - barrier.cpp regains the virtual Barrier hierarchy. Virtual dispatch cannot cross the host/device boundary -- a host-built vtable 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 -- so the guard also keeps every class symbol, float and double included, in the host object. - distance_type.cpp regains the error reporting helpers. Here the blocker is that nvcc cannot parse spdlog, so the guard covers the includes as well as the definitions; an #include in a false branch is never processed. Every use of the logger and of fmt must stay inside it. The free functions in both files are unchanged and still emit relocatable device code: cuobjdump confirms all six barrier functions in the fatbin, and the device object holds zero class symbols and zero spdlog symbols. Co-Authored-By: Claude Opus 5 (1M context) --- src/ipc/barrier/CMakeLists.txt | 7 +- src/ipc/barrier/barrier.cpp | 168 +++++++++++++++++++++++- src/ipc/barrier/barrier_classes.cpp | 159 ---------------------- src/ipc/distance/CMakeLists.txt | 1 - src/ipc/distance/distance_type.cpp | 45 +++++++ src/ipc/distance/distance_type_host.cpp | 41 ------ 6 files changed, 210 insertions(+), 211 deletions(-) delete mode 100644 src/ipc/barrier/barrier_classes.cpp delete mode 100644 src/ipc/distance/distance_type_host.cpp diff --git a/src/ipc/barrier/CMakeLists.txt b/src/ipc/barrier/CMakeLists.txt index 913951abf..6cfa28ae1 100644 --- a/src/ipc/barrier/CMakeLists.txt +++ b/src/ipc/barrier/CMakeLists.txt @@ -1,7 +1,6 @@ set(SOURCES adaptive_stiffness.cpp adaptive_stiffness.hpp - barrier_classes.cpp barrier_force_magnitude.cpp barrier_force_magnitude.hpp barrier.hpp @@ -9,10 +8,8 @@ set(SOURCES target_sources(ipc_toolkit PRIVATE ${SOURCES}) -# barrier.cpp (the free functions) is shared between host C++ and CUDA device -# code; under CUDA the helper generates a .cu wrapper (see -# ipc_toolkit_shared_device_sources). The other barrier .cpp files are ordinary -# host-only sources listed in SOURCES above. +# 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 cdf30e368..5538430f2 100644 --- a/src/ipc/barrier/barrier.cpp +++ b/src/ipc/barrier/barrier.cpp @@ -10,6 +10,10 @@ 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 @@ -21,7 +25,7 @@ template IPC_TOOLKIT_HOST_DEVICE T barrier(const T d, const T dhat) // b(d) = -(d-d̂)²ln(d / d̂) return select_lazy( d <= T(0), [&] { return infinity(); }, // - d < dhat, [&] { return -sqr(d - dhat) * ipc::numext::log(d / dhat); }, + d < dhat, [&] { return -sqr(d - dhat) * log(d / dhat); }, [&] { return T(0); }); } @@ -36,9 +40,7 @@ IPC_TOOLKIT_HOST_DEVICE T barrier_first_derivative(const T d, const T dhat) return select_lazy( d <= T(0), [&] { return T(0); }, // d < dhat, - [&] { - return (dhat - d) * (2 * ipc::numext::log(d / dhat) - dhat / d + 1); - }, + [&] { return (dhat - d) * (2 * log(d / dhat) - dhat / d + 1); }, [&] { return T(0); }); } @@ -51,7 +53,7 @@ IPC_TOOLKIT_HOST_DEVICE T barrier_second_derivative(const T d, const T dhat) d < dhat, [&] { const T dhat_d = dhat / d; - return (dhat_d + 2) * dhat_d - 2 * ipc::numext::log(d / dhat) - 3; + return (dhat_d + 2) * dhat_d - 2 * log(d / dhat) - 3; }, [&] { return T(0); }); } @@ -81,6 +83,162 @@ template SimdBatch barrier_second_derivative( const SimdBatch d, const SimdBatch 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 +T ClampedLogSqBarrier::operator()(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + // b(d) = (d-d̂)²ln²(d / d̂) + return select_lazy( + d <= T(0), [&] { return infinity(); }, // + d < dhat, + [&] { + const T log_d_dhat = log(d / dhat); + return sqr(d - dhat) * sqr(log_d_dhat); + }, + [&] { return T(0); }); +} + +template +T ClampedLogSqBarrier::first_derivative(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + // b(d) = (d - d̂)²ln²(d / d̂) + // b'(d) = 2 (d - d̂) ln²(d / d̂) + 2 (d - d̂)² ln(d / d̂) / d + // = 2 (d - d̂) ln(d / d̂) [ln(d / d̂) + (d - d̂) / d] + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < dhat, + [&] { + const T d_minus_dhat = (d - dhat); + const T log_d_dhat = log(d / dhat); + return T(2) * d_minus_dhat * log_d_dhat + * (log_d_dhat + d_minus_dhat / d); + }, + [&] { return T(0); }); +} + +template +T ClampedLogSqBarrier::second_derivative(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < dhat, + [&] { + const T t0 = dhat - d; + const T t1 = log(d / dhat); + const T t2 = sqr(t0) / sqr(d); + return T(2) * (sqr(t1) - (t1 - T(1)) * t2 - T(4) * t1 * t0 / d); + }, + [&] { return T(0); }); +} + +// ============================================================================ + +template +T CubicBarrier::operator()(const T d, const T dhat) const +{ + // b(d) = (d - d̂)³ + // + // A polynomial: finite at d <= 0, so unlike the log barriers it needs no + // penetration guard. + return select_lazy( + d < dhat, [&] { return -T(2) / T(3) / dhat * cubic(d - dhat); }, + [&] { return T(0); }); +} + +template +T CubicBarrier::first_derivative(const T d, const T dhat) const +{ + return select_lazy( + d < dhat, [&] { return T(-2) / dhat * sqr(d - dhat); }, + [&] { return T(0); }); +} + +template +T CubicBarrier::second_derivative(const T d, const T dhat) const +{ + return select_lazy( + d < dhat, [&] { return T(4) * (T(1) - d / dhat); }, + [&] { return T(0); }); +} + +// ============================================================================ + +template +T TwoStageBarrier::operator()(const T d, const T dhat) const +{ + using namespace ipc::numext; // log + return select_lazy( + d <= T(0), [&] { return infinity(); }, // + d < T(0.5) * dhat, + [&] { return T(-0.25) * sqr(dhat) * (log(T(2) * d / dhat) - T(0.5)); }, + d < dhat, [&] { return T(0.5) * sqr(dhat - d); }, // + [&] { return T(0); }); +} + +template +T TwoStageBarrier::first_derivative(const T d, const T dhat) const +{ + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < T(0.5) * dhat, [&] { return T(-0.25) * sqr(dhat) / d; }, // + d < dhat, [&] { return d - dhat; }, // + [&] { return T(0); }); +} + +template +T TwoStageBarrier::second_derivative(const T d, const T dhat) const +{ + return select_lazy( + d <= T(0), [&] { return T(0); }, // + d < T(0.5) * dhat, [&] { return (T(0.25) * sqr(dhat)) / sqr(d); }, // + d < dhat, [&] { return T(1); }, // + [&] { return T(0); }); +} + +// ============================================================================ +// Explicit template instantiations +/// @cond DOXYGEN_SKIP +template class BarrierBase; +template class BarrierBase; +template class ClampedLogBarrier; +template class ClampedLogBarrier; +template class ClampedLogSqBarrier; +template class ClampedLogSqBarrier; +template class CubicBarrier; +template class CubicBarrier; +template class TwoStageBarrier; +template class TwoStageBarrier; +#ifdef IPC_TOOLKIT_WITH_SIMD +template class BarrierBase>; +template class BarrierBase>; +template class ClampedLogBarrier>; +template class ClampedLogBarrier>; +template class ClampedLogSqBarrier>; +template class ClampedLogSqBarrier>; +template class CubicBarrier>; +template class CubicBarrier>; +template class TwoStageBarrier>; +template class TwoStageBarrier>; +#endif +/// @endcond // ============================================================================ +#endif // !__CUDACC__ + } // namespace ipc diff --git a/src/ipc/barrier/barrier_classes.cpp b/src/ipc/barrier/barrier_classes.cpp deleted file mode 100644 index 2be54fdf5..000000000 --- a/src/ipc/barrier/barrier_classes.cpp +++ /dev/null @@ -1,159 +0,0 @@ -// Out-of-line definitions for the Barrier class hierarchy. -// -// These live apart from barrier.cpp because that file is compiled as a CUDA -// translation unit when CUDA is enabled, and a virtual hierarchy has nothing to -// offer device code anyway. Keeping the classes here lets barrier.cpp stay -// small and device-callable, holding only the three free barrier functions. -#include "barrier.hpp" - -#include -#include -#include - -namespace ipc { - -// ============================================================================ - -template -T ClampedLogSqBarrier::operator()(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - // b(d) = (d-d̂)²ln²(d / d̂) - return select_lazy( - d <= T(0), [&] { return infinity(); }, // - d < dhat, - [&] { - const T log_d_dhat = log(d / dhat); - return sqr(d - dhat) * sqr(log_d_dhat); - }, - [&] { return T(0); }); -} - -template -T ClampedLogSqBarrier::first_derivative(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - // b(d) = (d - d̂)²ln²(d / d̂) - // b'(d) = 2 (d - d̂) ln²(d / d̂) + 2 (d - d̂)² ln(d / d̂) / d - // = 2 (d - d̂) ln(d / d̂) [ln(d / d̂) + (d - d̂) / d] - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < dhat, - [&] { - const T d_minus_dhat = (d - dhat); - const T log_d_dhat = log(d / dhat); - return T(2) * d_minus_dhat * log_d_dhat - * (log_d_dhat + d_minus_dhat / d); - }, - [&] { return T(0); }); -} - -template -T ClampedLogSqBarrier::second_derivative(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < dhat, - [&] { - const T t0 = dhat - d; - const T t1 = log(d / dhat); - const T t2 = sqr(t0) / sqr(d); - return T(2) * (sqr(t1) - (t1 - T(1)) * t2 - T(4) * t1 * t0 / d); - }, - [&] { return T(0); }); -} - -// ============================================================================ - -template -T CubicBarrier::operator()(const T d, const T dhat) const -{ - // b(d) = (d - d̂)³ - // - // A polynomial: finite at d <= 0, so unlike the log barriers it needs no - // penetration guard. - return select_lazy( - d < dhat, [&] { return -T(2) / T(3) / dhat * cubic(d - dhat); }, - [&] { return T(0); }); -} - -template -T CubicBarrier::first_derivative(const T d, const T dhat) const -{ - return select_lazy( - d < dhat, [&] { return T(-2) / dhat * sqr(d - dhat); }, - [&] { return T(0); }); -} - -template -T CubicBarrier::second_derivative(const T d, const T dhat) const -{ - return select_lazy( - d < dhat, [&] { return T(4) * (T(1) - d / dhat); }, - [&] { return T(0); }); -} - -// ============================================================================ - -template -T TwoStageBarrier::operator()(const T d, const T dhat) const -{ - using namespace ipc::numext; // log - return select_lazy( - d <= T(0), [&] { return infinity(); }, // - d < T(0.5) * dhat, - [&] { return T(-0.25) * sqr(dhat) * (log(T(2) * d / dhat) - T(0.5)); }, - d < dhat, [&] { return T(0.5) * sqr(dhat - d); }, // - [&] { return T(0); }); -} - -template -T TwoStageBarrier::first_derivative(const T d, const T dhat) const -{ - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < T(0.5) * dhat, [&] { return T(-0.25) * sqr(dhat) / d; }, // - d < dhat, [&] { return d - dhat; }, // - [&] { return T(0); }); -} - -template -T TwoStageBarrier::second_derivative(const T d, const T dhat) const -{ - return select_lazy( - d <= T(0), [&] { return T(0); }, // - d < T(0.5) * dhat, [&] { return (T(0.25) * sqr(dhat)) / sqr(d); }, // - d < dhat, [&] { return T(1); }, // - [&] { return T(0); }); -} - -// ============================================================================ -// Explicit template instantiations -/// @cond DOXYGEN_SKIP -template class BarrierBase; -template class BarrierBase; -template class ClampedLogBarrier; -template class ClampedLogBarrier; -template class ClampedLogSqBarrier; -template class ClampedLogSqBarrier; -template class CubicBarrier; -template class CubicBarrier; -template class TwoStageBarrier; -template class TwoStageBarrier; -#ifdef IPC_TOOLKIT_WITH_SIMD -template class BarrierBase>; -template class BarrierBase>; -template class ClampedLogBarrier>; -template class ClampedLogBarrier>; -template class ClampedLogSqBarrier>; -template class ClampedLogSqBarrier>; -template class CubicBarrier>; -template class CubicBarrier>; -template class TwoStageBarrier>; -template class TwoStageBarrier>; -#endif -/// @endcond -// ============================================================================ - -} // namespace ipc diff --git a/src/ipc/distance/CMakeLists.txt b/src/ipc/distance/CMakeLists.txt index e24df9796..9ede9e1e7 100644 --- a/src/ipc/distance/CMakeLists.txt +++ b/src/ipc/distance/CMakeLists.txt @@ -1,6 +1,5 @@ set(SOURCES distance_type.hpp - distance_type_host.cpp edge_edge.hpp edge_edge_mollifier.hpp line_line.hpp diff --git a/src/ipc/distance/distance_type.cpp b/src/ipc/distance/distance_type.cpp index 3c5b7beda..f5bf59d84 100644 --- a/src/ipc/distance/distance_type.cpp +++ b/src/ipc/distance/distance_type.cpp @@ -6,6 +6,19 @@ #include +// 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 + +#include + +#include +#endif + namespace ipc::detail { template @@ -222,4 +235,36 @@ template EdgeEdgeDistanceType edge_edge_parallel_distance_type(Eigen::Co #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_host.cpp b/src/ipc/distance/distance_type_host.cpp deleted file mode 100644 index d3f995e82..000000000 --- a/src/ipc/distance/distance_type_host.cpp +++ /dev/null @@ -1,41 +0,0 @@ -#include "distance_type.hpp" - -#include - -#include - -#include - -// Host-only definitions for the error reporting helpers declared in -// distance_type.hpp. -// -// These live apart from distance_type.cpp because that file is compiled as a -// CUDA translation unit when CUDA is enabled, and nvcc cannot parse spdlog. -// Keeping the logger and the exceptions here means the shared device sources -// never have to include either one. The inline wrappers in the header pick -// between these and a device-side trap. -namespace ipc::detail { - -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)); -} - -} // namespace ipc::detail From cb59cded50443dd5ab7270415152ce6a053db24c Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 8 Sep 2026 16:02:14 -0400 Subject: [PATCH 3/4] Move clamp into simd.hpp and add a batch overload morton.hpp carried a private detail::clamp, a device-safe std::clamp usable only from that header. It now sits in simd.hpp beside select, all_of, and infinity, which pair a scalar and a batch overload the same way, and covers float, double, and an xsimd batch. - The scalar overload is IPC_TOOLKIT_HOST_DEVICE, so it works on the GPU. std::clamp does not: MSVC's debug STL checks the bounds with _STL_VERIFY, which expands to __debugbreak(), and nvcc's NVVM backend then emits invalid IR. The comparison order matches std::clamp exactly, so a NaN v still passes through unchanged. - The batch overload blends the same first-match-wins order with two xsimd::select calls, since a batch comparison answers per-lane rather than with one bool. It is the more specialized template, so it wins for a batch argument. Host only, as config.hpp retracts IPC_TOOLKIT_WITH_SIMD for a CUDA translation unit. Co-Authored-By: Claude Opus 5 (1M context) --- src/ipc/math/morton.hpp | 29 +++++++---------------------- src/ipc/utils/simd.hpp | 33 +++++++++++++++++++++++++++++++-- 2 files changed, 38 insertions(+), 24 deletions(-) diff --git a/src/ipc/math/morton.hpp b/src/ipc/math/morton.hpp index 02b9c5619..d43c06bac 100644 --- a/src/ipc/math/morton.hpp +++ b/src/ipc/math/morton.hpp @@ -1,27 +1,12 @@ #pragma once -#include // for IPC_TOOLKIT_HOST_DEVICE +#include // for IPC_TOOLKIT_HOST_DEVICE +#include // for clamp #include // for uint64_t namespace ipc { -namespace detail { - /// @brief Clamp a double to [lo, hi] (a device-safe std::clamp). - /// - /// std::clamp cannot be called from device code: MSVC's debug STL checks - /// the bounds with _STL_VERIFY, which expands to __debugbreak(); 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 input passing through unchanged -- is - /// identical. - IPC_TOOLKIT_HOST_DEVICE inline double - clamp(const double v, const double lo, const double hi) - { - return v < lo ? lo : (hi < v ? hi : v); - } -} // namespace detail - /// @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. @@ -55,8 +40,8 @@ IPC_TOOLKIT_HOST_DEVICE inline uint64_t expand_bits_2(uint64_t v) IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_2D(double x, double y) { constexpr double scale = 1ULL << 32; - x = detail::clamp(x * scale, 0.0, scale - 1); - y = detail::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; @@ -70,9 +55,9 @@ IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_2D(double x, double y) IPC_TOOLKIT_HOST_DEVICE inline uint64_t morton_3D(double x, double y, double z) { constexpr double scale = 1ULL << 21; - x = detail::clamp(x * scale, 0.0, scale - 1); - y = detail::clamp(y * scale, 0.0, scale - 1); - z = detail::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/utils/simd.hpp b/src/ipc/utils/simd.hpp index 7c4dc701f..72c099635 100644 --- a/src/ipc/utils/simd.hpp +++ b/src/ipc/utils/simd.hpp @@ -4,6 +4,9 @@ #include #include +#ifdef IPC_TOOLKIT_WITH_SIMD +#include +#endif #include #include @@ -63,6 +66,19 @@ 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 +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 IPC_TOOLKIT_HOST_DEVICE inline T infinity() { @@ -111,8 +127,6 @@ select_lazy(const Mask& mask, F&& value, Rest&&... rest) #ifdef IPC_TOOLKIT_WITH_SIMD -#include - namespace Eigen { /// @brief Let Eigen treat an xsimd batch as a scalar type. @@ -200,6 +214,21 @@ inline bool all_of(const xsimd::batch_bool& 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 +inline xsimd::batch clamp( + const xsimd::batch& v, + const xsimd::batch& lo, + const xsimd::batch& hi) +{ + return xsimd::select(v < lo, lo, xsimd::select(hi < v, hi, v)); +} + } // namespace ipc #endif From 785174beb0204792f20ea143c2e35ee9c0627f69 Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Tue, 8 Sep 2026 16:16:05 -0400 Subject: [PATCH 4/4] Share the GPU parity test scaffolding across the test_gpu_*.cu files The distance, geometry, and tangent GPU tests each carried a byte-identical copy of the same seven helpers -- REQUIRE_CUDA, GpuTestKernel, write_out, write_expected, run_gpu_kernel, skip_if_no_cuda_device, and check_gpu_matches_host -- and the barrier test hand-rolled a fourth version of the same launch, transfer, and comparison. They now live in tests/src/tests/gpu_utils.hpp, beside simd_utils.hpp and following its include and namespace conventions, leaving only the kernels and their fixtures in the individual files. - check_gpu_matches_host compares an infinite reference by sign rather than by margin. That was the only thing keeping the barrier test from sharing the comparison, and it is a no-op for the finite values the other three produce. - The header defines __device__ code, so it #errors when included outside a .cu rather than failing on an unrecognized __device__. - geometry's copy documented run_gpu_kernel as skipping the test when no device is present, which it never did; the comment now sits on skip_if_no_cuda_device. Net -152 lines. Verified with the full suite (345/345) and with the [gpu] filter, whose 21 cases execute rather than skip on this machine. --- tests/src/tests/CMakeLists.txt | 5 + tests/src/tests/barrier/test_gpu_barrier.cu | 41 +------ tests/src/tests/distance/test_gpu_distance.cu | 81 +----------- tests/src/tests/geometry/test_gpu_geometry.cu | 82 +------------ tests/src/tests/gpu_utils.hpp | 116 ++++++++++++++++++ tests/src/tests/tangent/test_gpu_tangent.cu | 81 +----------- 6 files changed, 138 insertions(+), 268 deletions(-) create mode 100644 tests/src/tests/gpu_utils.hpp 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/test_gpu_barrier.cu b/tests/src/tests/barrier/test_gpu_barrier.cu index 62e7f1319..7f2aede1d 100644 --- a/tests/src/tests/barrier/test_gpu_barrier.cu +++ b/tests/src/tests/barrier/test_gpu_barrier.cu @@ -14,17 +14,15 @@ #include -#include -#include +#include #include -#include #include -namespace { +using namespace ipc::tests; -#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) +namespace { __global__ void barrier_kernel(const double* in, double* out) { @@ -45,11 +43,7 @@ __global__ void barrier_kernel(const double* in, double* out) TEST_CASE("GPU barrier", "[barrier][gpu]") { - int device_count = 0; - const cudaError_t err = cudaGetDeviceCount(&device_count); - if (err != cudaSuccess || device_count == 0) { - SKIP("No CUDA device available; kernel compiled but not executed."); - } + skip_if_no_cuda_device(); const double dhat = 1e-2; const std::vector in = { dhat, 0.5 * dhat, 2 * dhat, -1.0 }; @@ -62,32 +56,7 @@ TEST_CASE("GPU barrier", "[barrier][gpu]") expected.push_back(ipc::barrier_second_derivative(d, dhat)); } - double *d_in = nullptr, *d_out = nullptr; - REQUIRE_CUDA(cudaMalloc(&d_in, in.size() * sizeof(double))); - REQUIRE_CUDA(cudaMalloc(&d_out, expected.size() * sizeof(double))); - REQUIRE_CUDA(cudaMemcpy( - d_in, in.data(), in.size() * sizeof(double), cudaMemcpyHostToDevice)); - - barrier_kernel<<<1, 1>>>(d_in, d_out); - REQUIRE_CUDA(cudaGetLastError()); - REQUIRE_CUDA(cudaDeviceSynchronize()); - - std::vector out(expected.size()); - REQUIRE_CUDA(cudaMemcpy( - out.data(), d_out, out.size() * sizeof(double), - cudaMemcpyDeviceToHost)); - cudaFree(d_in); - cudaFree(d_out); - - 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)); - } - } + check_gpu_matches_host(barrier_kernel, in, expected); } #endif // IPC_TOOLKIT_WITH_CUDA diff --git a/tests/src/tests/distance/test_gpu_distance.cu b/tests/src/tests/distance/test_gpu_distance.cu index 30320ed6a..3f2289b95 100644 --- a/tests/src/tests/distance/test_gpu_distance.cu +++ b/tests/src/tests/distance/test_gpu_distance.cu @@ -28,88 +28,15 @@ #include #include +#include + #include #include -namespace { - -#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) - -/// @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 -__device__ int -write_out(const Eigen::MatrixBase& 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 -void write_expected( - const Eigen::MatrixBase& m, std::vector& expected) -{ - for (int r = 0; r < m.rows(); ++r) { - for (int c = 0; c < m.cols(); ++c) { - expected.push_back(m(r, c)); - } - } -} +using namespace ipc::tests; -/// @brief Launch `kernel` with the given inputs; returns the device outputs. -std::vector run_gpu_kernel( - GpuTestKernel kernel, const std::vector& 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 out(n_out); - REQUIRE_CUDA(cudaMemcpy( - out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); - cudaFree(d_in); - cudaFree(d_out); - return out; -} - -/// @brief SKIP the test if no CUDA device is present (kernels have still been -/// compiled as device code, which is the primary verification). -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."); - } -} - -void check_gpu_matches_host( - GpuTestKernel kernel, - const std::vector& in, - const std::vector& expected) -{ - const std::vector out = run_gpu_kernel(kernel, in, expected.size()); - for (size_t i = 0; i < expected.size(); ++i) { - CAPTURE(i); - CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); - } -} +namespace { // --------------------------------------------------------------------------- diff --git a/tests/src/tests/geometry/test_gpu_geometry.cu b/tests/src/tests/geometry/test_gpu_geometry.cu index 054833eac..5838ff30a 100644 --- a/tests/src/tests/geometry/test_gpu_geometry.cu +++ b/tests/src/tests/geometry/test_gpu_geometry.cu @@ -18,89 +18,15 @@ #include #include +#include + #include #include -namespace { - -#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) - -/// @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 -__device__ int -write_out(const Eigen::MatrixBase& 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 -void write_expected( - const Eigen::MatrixBase& m, std::vector& expected) -{ - for (int r = 0; r < m.rows(); ++r) { - for (int c = 0; c < m.cols(); ++c) { - expected.push_back(m(r, c)); - } - } -} +using namespace ipc::tests; -/// @brief Launch `kernel` with the given inputs; returns the device outputs. -/// Skips the surrounding test if no CUDA device is available. -std::vector run_gpu_kernel( - GpuTestKernel kernel, const std::vector& 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 out(n_out); - REQUIRE_CUDA(cudaMemcpy( - out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); - cudaFree(d_in); - cudaFree(d_out); - return out; -} - -/// @brief SKIP the test if no CUDA device is present (kernels have still been -/// compiled as device code, which is the primary verification). -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."); - } -} - -void check_gpu_matches_host( - GpuTestKernel kernel, - const std::vector& in, - const std::vector& expected) -{ - const std::vector out = run_gpu_kernel(kernel, in, expected.size()); - for (size_t i = 0; i < expected.size(); ++i) { - CAPTURE(i); - CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); - } -} +namespace { // --------------------------------------------------------------------------- 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 +#include + +#include + +#include + +#include +#include + +#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 +__device__ int +write_out(const Eigen::MatrixBase& 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 +inline void write_expected( + const Eigen::MatrixBase& m, std::vector& 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 run_gpu_kernel( + GpuTestKernel kernel, const std::vector& 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 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& in, + const std::vector& expected) +{ + const std::vector 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/test_gpu_tangent.cu b/tests/src/tests/tangent/test_gpu_tangent.cu index 1f5e6b158..6b682a502 100644 --- a/tests/src/tests/tangent/test_gpu_tangent.cu +++ b/tests/src/tests/tangent/test_gpu_tangent.cu @@ -19,88 +19,15 @@ #include #include +#include + #include #include -namespace { - -#define REQUIRE_CUDA(expr) REQUIRE((expr) == cudaSuccess) - -/// @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 -__device__ int -write_out(const Eigen::MatrixBase& 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 -void write_expected( - const Eigen::MatrixBase& m, std::vector& expected) -{ - for (int r = 0; r < m.rows(); ++r) { - for (int c = 0; c < m.cols(); ++c) { - expected.push_back(m(r, c)); - } - } -} +using namespace ipc::tests; -/// @brief Launch `kernel` with the given inputs; returns the device outputs. -std::vector run_gpu_kernel( - GpuTestKernel kernel, const std::vector& 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 out(n_out); - REQUIRE_CUDA(cudaMemcpy( - out.data(), d_out, n_out * sizeof(double), cudaMemcpyDeviceToHost)); - cudaFree(d_in); - cudaFree(d_out); - return out; -} - -/// @brief SKIP the test if no CUDA device is present (kernels have still been -/// compiled as device code, which is the primary verification). -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."); - } -} - -void check_gpu_matches_host( - GpuTestKernel kernel, - const std::vector& in, - const std::vector& expected) -{ - const std::vector out = run_gpu_kernel(kernel, in, expected.size()); - for (size_t i = 0; i < expected.size(); ++i) { - CAPTURE(i); - CHECK(out[i] == Catch::Approx(expected[i]).margin(1e-12)); - } -} +namespace { // ---------------------------------------------------------------------------