From bd4fa95ebfa3ce5e8b1a9ab68c1d332a046d1934 Mon Sep 17 00:00:00 2001 From: zachferguson Date: Tue, 21 Jan 2025 10:29:32 -0500 Subject: [PATCH 01/10] Add cubic barrier functions --- docs/source/cpp-api/barrier.rst | 12 +++ docs/source/python-api/barrier.rst | 17 ++- notebooks/physical_barrier.ipynb | 97 +++++++++++++---- python/src/barrier/barrier.cpp | 47 ++++++++- python/src/candidates/collision_stencil.cpp | 64 ++++++++++- src/ipc/barrier/barrier.cpp | 74 +++++++++++++ src/ipc/barrier/barrier.hpp | 111 +++++++++++++++++++- src/ipc/candidates/collision_stencil.hpp | 39 +++++++ src/ipc/candidates/edge_edge.hpp | 4 + src/ipc/candidates/edge_vertex.hpp | 4 + src/ipc/candidates/face_vertex.hpp | 4 + src/ipc/candidates/vertex_vertex.hpp | 4 + src/ipc/potentials/normal_potential.cpp | 6 +- tests/src/tests/barrier/test_barrier.cpp | 5 + 14 files changed, 458 insertions(+), 30 deletions(-) diff --git a/docs/source/cpp-api/barrier.rst b/docs/source/cpp-api/barrier.rst index 5834f8e2c..3b243d377 100644 --- a/docs/source/cpp-api/barrier.rst +++ b/docs/source/cpp-api/barrier.rst @@ -33,4 +33,16 @@ Normalized Clamped Log Barrier ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. doxygenclass:: ipc::NormalizedClampedLogBarrier + :allow-dot-graphs: + +Clamped Log Squared Barrier +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. doxygenclass:: ipc::ClampedLogSqBarrier + :allow-dot-graphs: + +Cubic Barrier +~~~~~~~~~~~~~ + +.. doxygenclass:: ipc::CubicBarrier :allow-dot-graphs: \ No newline at end of file diff --git a/docs/source/python-api/barrier.rst b/docs/source/python-api/barrier.rst index 2ad4b74eb..fa0cc329b 100644 --- a/docs/source/python-api/barrier.rst +++ b/docs/source/python-api/barrier.rst @@ -25,4 +25,19 @@ Barrier Class Clamped Log Barrier ~~~~~~~~~~~~~~~~~~~ -.. autoclass:: ipctk.ClampedLogBarrier \ No newline at end of file +.. autoclass:: ipctk.ClampedLogBarrier + +Normalized Clamped Log Barrier +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ipctk.NormalizedClampedLogBarrier + +Clamped Log Squared Barrier +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. autoclass:: ipctk.ClampedLogSqBarrier + +Cubic Barrier +~~~~~~~~~~~~~ + +.. autoclass:: ipctk.CubicBarrier \ No newline at end of file diff --git a/notebooks/physical_barrier.ipynb b/notebooks/physical_barrier.ipynb index 5e256f2b0..6138e8c89 100644 --- a/notebooks/physical_barrier.ipynb +++ b/notebooks/physical_barrier.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 13, + "execution_count": 64, "metadata": {}, "outputs": [], "source": [ @@ -12,18 +12,18 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 65, "metadata": {}, "outputs": [], "source": [ - "dhat = sympy.Symbol('\\hat{d}')\n", + "dhat = sympy.Symbol(r'\\hat{d}')\n", "d = sympy.Symbol('d')\n", "dmin = sympy.Symbol('d_{\\\\text{min}}')" ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 66, "metadata": {}, "outputs": [ { @@ -61,24 +61,48 @@ }, "metadata": {}, "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\displaystyle \\left(- \\hat{d} + d\\right)^{2} \\log{\\left(\\frac{d}{\\hat{d}} \\right)}^{2}$" + ], + "text/plain": [ + "(-\\hat{d} + d)**2*log(d/\\hat{d})**2" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ "def barrier(d, dhat):\n", " return -(d - dhat)**2 * sympy.log(d/dhat)\n", "\n", + "\n", "def normalized_barrier(d, dhat):\n", " return -(d/dhat - 1)**2 * sympy.log(d/dhat)\n", "\n", + "\n", "def physical_barrier(d, p_dhat):\n", " return dhat * normalized_barrier(d, p_dhat)\n", "\n", - "display(barrier(d, dhat), normalized_barrier(d, dhat), physical_barrier(d, dhat))" + "\n", + "def log_sq_barrier(d, dhat):\n", + " return (d - dhat)**2 * sympy.log(d/dhat)**2\n", + "\n", + "\n", + "display(\n", + " barrier(d, dhat),\n", + " normalized_barrier(d, dhat),\n", + " physical_barrier(d, dhat),\n", + " log_sq_barrier(d, dhat)\n", + ")" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 67, "metadata": {}, "outputs": [ { @@ -90,7 +114,7 @@ "1/\\hat{d}" ] }, - "execution_count": 16, + "execution_count": 67, "metadata": {}, "output_type": "execute_result" } @@ -101,7 +125,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 68, "metadata": {}, "outputs": [ { @@ -113,7 +137,7 @@ "\\hat{d}**(-3)" ] }, - "execution_count": 17, + "execution_count": 68, "metadata": {}, "output_type": "execute_result" } @@ -124,7 +148,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 69, "metadata": {}, "outputs": [ { @@ -136,7 +160,7 @@ "1/(\\hat{d}*(\\hat{d} + 2*d_{\\text{min}})**2)" ] }, - "execution_count": 18, + "execution_count": 69, "metadata": {}, "output_type": "execute_result" } @@ -147,9 +171,16 @@ "(physical_barrier(d_param, dhat_param) / barrier(d_param, dhat_param)).simplify()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Code Generation" + ] + }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 70, "metadata": {}, "outputs": [ { @@ -157,17 +188,18 @@ "output_type": "stream", "text": [ "const auto t0 = d/dhat;\n", - "r = -std::pow(1 - t0, 2)*std::log(t0);\n" + "return -std::pow(1 - t0, 2)*std::log(t0);\n" ] } ], "source": [ - "print(generate_code(normalized_barrier(d, dhat), \"r\").replace(\"\\hat{d}\", \"dhat\"))" + "print(generate_code(normalized_barrier(\n", + " d, dhat)).replace(r\"\\hat{d}\", \"dhat\"))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 71, "metadata": {}, "outputs": [ { @@ -177,17 +209,18 @@ "const auto t0 = 1.0/dhat;\n", "const auto t1 = d*t0;\n", "const auto t2 = 1 - t1;\n", - "r = t2*(2*t0*std::log(t1) - t2/d);\n" + "return t2*(2*t0*std::log(t1) - t2/d);\n" ] } ], "source": [ - "print(generate_code(normalized_barrier(d, dhat).diff(d), \"r\").replace(\"\\hat{d}\", \"dhat\"))" + "print(generate_code(normalized_barrier(d, dhat).diff(\n", + " d)).replace(r\"\\hat{d}\", \"dhat\"))" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 72, "metadata": {}, "outputs": [ { @@ -197,12 +230,34 @@ "const auto t0 = 1.0/dhat;\n", "const auto t1 = d*t0;\n", "const auto t2 = 1 - t1;\n", - "r = 4*t0*t2/d + std::pow(t2, 2)/std::pow(d, 2) - 2*std::log(t1)/std::pow(dhat, 2);\n" + "return 4*t0*t2/d + std::pow(t2, 2)/std::pow(d, 2) - 2*std::log(t1)/std::pow(dhat, 2);\n" + ] + } + ], + "source": [ + "print(generate_code(normalized_barrier(d, dhat).diff(\n", + " d).diff(d)).replace(r\"\\hat{d}\", \"dhat\"))" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "const auto t0 = std::log(d/dhat);\n", + "const auto t1 = dhat - d;\n", + "const auto t2 = std::pow(t1, 2)/std::pow(d, 2);\n", + "return 2*std::pow(t0, 2) - 2*t0*t2 + 2*t2 - 8*t0*t1/d;\n" ] } ], "source": [ - "print(generate_code(normalized_barrier(d, dhat).diff(d).diff(d), \"r\").replace(\"\\hat{d}\", \"dhat\"))" + "print(generate_code(log_sq_barrier(\n", + " d, dhat).diff(d).diff(d)).replace(r\"\\hat{d}\", \"dhat\"))" ] }, { @@ -229,7 +284,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.2" + "version": "3.12.7" }, "orig_nbformat": 4 }, diff --git a/python/src/barrier/barrier.cpp b/python/src/barrier/barrier.cpp index b47d4be54..6f34817f8 100644 --- a/python/src/barrier/barrier.cpp +++ b/python/src/barrier/barrier.cpp @@ -72,7 +72,52 @@ void define_barrier(py::module_& m) py::class_>( m, "ClampedLogBarrier", - "Smoothly clamped log barrier functions from [Li et al. 2020].") + R"ipc_Qu8mg5v7( + Smoothly clamped log barrier functions from [Li et al. 2020]. + + .. math:: + + b(d) = -(d-\hat{d})^2\ln\left(\frac{d}{\hat{d}}\right) + + )ipc_Qu8mg5v7") + .def(py::init()); + + py::class_< + NormalizedClampedLogBarrier, ClampedLogBarrier, + std::shared_ptr>( + m, "NormalizedClampedLogBarrier", + R"ipc_Qu8mg5v7( + Normalized barrier function from [Li et al. 2023]. + + .. math:: + + b(d) = -\left(\frac{d}{\hat{d}}-1\right)^2\ln\left(\frac{d}{\hat{d}}\right) + + )ipc_Qu8mg5v7") + .def(py::init()); + + py::class_< + ClampedLogSqBarrier, Barrier, std::shared_ptr>( + m, "ClampedLogSqBarrier", + R"ipc_Qu8mg5v7( + Clamped log barrier with a quadratic log term from [Huang et al. 2024]. + + .. math:: + + b(d) = (d-\hat{d})^2\ln^2\left(\frac{d}{\hat{d}}\right) + + )ipc_Qu8mg5v7") + .def(py::init()); + + py::class_>( + m, "CubicBarrier", R"ipc_Qu8mg5v7( + Cubic barrier function from [Ando 2024]. + + .. math:: + + b(d) = -\frac{2}{3\hat{d}} (d - \hat{d})^3 + + )ipc_Qu8mg5v7") .def(py::init()); m.def( diff --git a/python/src/candidates/collision_stencil.cpp b/python/src/candidates/collision_stencil.cpp index 65091560a..3a4aeb8b3 100644 --- a/python/src/candidates/collision_stencil.cpp +++ b/python/src/candidates/collision_stencil.cpp @@ -69,7 +69,63 @@ void define_collision_stencil(py::module_& m) )ipc_Qu8mg5v7", py::arg("X"), py::arg("edges"), py::arg("faces")) .def( - "compute_distance", &CollisionStencil::compute_distance, + "compute_distance", + py::overload_cast< + const Eigen::MatrixXd&, const Eigen::MatrixXi&, + const Eigen::MatrixXi&>( + &CollisionStencil::compute_distance, py::const_), + R"ipc_Qu8mg5v7( + Compute the distance of the stencil. + + Parameters: + vertices: Collision mesh vertices. + edges: Collision mesh edges. + faces: Collision mesh faces. + + Returns: + Distance of the stencil. + )ipc_Qu8mg5v7", + py::arg("vertices"), py::arg("edges"), py::arg("faces")) + .def( + "compute_distance_gradient", + py::overload_cast< + const Eigen::MatrixXd&, const Eigen::MatrixXi&, + const Eigen::MatrixXi&>( + &CollisionStencil::compute_distance_gradient, py::const_), + R"ipc_Qu8mg5v7( + Compute the distance gradient of the stencil w.r.t. the stencil's vertex positions. + + Parameters: + vertices: Collision mesh vertices. + edges: Collision mesh edges. + faces: Collision mesh faces. + + Returns: + Distance gradient of the stencil w.r.t. the stencil's vertex positions. + )ipc_Qu8mg5v7", + py::arg("vertices"), py::arg("edges"), py::arg("faces")) + .def( + "compute_distance_hessian", + py::overload_cast< + const Eigen::MatrixXd&, const Eigen::MatrixXi&, + const Eigen::MatrixXi&>( + &CollisionStencil::compute_distance_hessian, py::const_), + R"ipc_Qu8mg5v7( + Compute the distance Hessian of the stencil w.r.t. the stencil's vertex positions. + + Parameters: + vertices: Collision mesh vertices. + edges: Collision mesh edges. + faces: Collision mesh faces. + + Returns: + Distance Hessian of the stencil w.r.t. the stencil's vertex positions. + )ipc_Qu8mg5v7", + py::arg("vertices"), py::arg("edges"), py::arg("faces")) + .def( + "compute_distance", + py::overload_cast( + &CollisionStencil::compute_distance, py::const_), R"ipc_Qu8mg5v7( Compute the distance of the stencil. @@ -85,7 +141,8 @@ void define_collision_stencil(py::module_& m) py::arg("positions")) .def( "compute_distance_gradient", - &CollisionStencil::compute_distance_gradient, + py::overload_cast( + &CollisionStencil::compute_distance_gradient, py::const_), R"ipc_Qu8mg5v7( Compute the distance gradient of the stencil w.r.t. the stencil's vertex positions. @@ -101,7 +158,8 @@ void define_collision_stencil(py::module_& m) py::arg("positions")) .def( "compute_distance_hessian", - &CollisionStencil::compute_distance_hessian, + py::overload_cast( + &CollisionStencil::compute_distance_hessian, py::const_), R"ipc_Qu8mg5v7( Compute the distance Hessian of the stencil w.r.t. the stencil's vertex positions. diff --git a/src/ipc/barrier/barrier.cpp b/src/ipc/barrier/barrier.cpp index dee4e30dd..73d6450fa 100644 --- a/src/ipc/barrier/barrier.cpp +++ b/src/ipc/barrier/barrier.cpp @@ -42,4 +42,78 @@ double barrier_second_derivative(const double d, const double dhat) return (dhat_d + 2) * dhat_d - 2 * log(d / dhat) - 3; } +// ============================================================================ + +double ClampedLogSqBarrier::operator()(const double d, const double dhat) const +{ + if (d <= 0.0) { + return std::numeric_limits::infinity(); + } + if (d >= dhat) { + return 0; + } + // b(d) = (d-d̂)²ln²(d / d̂) + const double d_minus_dhat = (d - dhat); + const double log_d_dhat = log(d / dhat); + return d_minus_dhat * d_minus_dhat * log_d_dhat * log_d_dhat; +} + +double +ClampedLogSqBarrier::first_derivative(const double d, const double dhat) const +{ + if (d <= 0.0 || d >= dhat) { + return 0.0; + } + // 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] + const double d_minus_dhat = (d - dhat); + const double log_d_dhat = log(d / dhat); + return 2 * d_minus_dhat * log_d_dhat * (log_d_dhat + d_minus_dhat / d); +} + +double +ClampedLogSqBarrier::second_derivative(const double d, const double dhat) const +{ + if (d <= 0.0 || d >= dhat) { + return 0.0; + } + const double t0 = dhat - d; + const double t1 = log(d / dhat); + const double t2 = (t0 * t0) / (d * d); + return 2 * ((t1 * t1) - (t1 - 1) * t2 - 4 * t1 * t0 / d); +} + +// ============================================================================ + +double CubicBarrier::operator()(const double d, const double dhat) const +{ + if (d < dhat) { + // b(d) = (d - d̂)³ + const double d_minus_dhat = (d - dhat); + return -2.0 / 3.0 / dhat * d_minus_dhat * d_minus_dhat * d_minus_dhat; + } else { + return 0; + } +} + +double CubicBarrier::first_derivative(const double d, const double dhat) const +{ + if (d < dhat) { + const double d_minus_dhat = (d - dhat); + return -2 / dhat * d_minus_dhat * d_minus_dhat; + } else { + return 0; + } +} + +double CubicBarrier::second_derivative(const double d, const double dhat) const +{ + if (d < dhat) { + return 4 * (1 - d / dhat); + } else { + return 0; + } +} + } // namespace ipc diff --git a/src/ipc/barrier/barrier.hpp b/src/ipc/barrier/barrier.hpp index 50f40a1a6..014a970c6 100644 --- a/src/ipc/barrier/barrier.hpp +++ b/src/ipc/barrier/barrier.hpp @@ -131,7 +131,7 @@ class ClampedLogBarrier : public Barrier { /// @return The units of the barrier function. double units(const double dhat) const override { - // (d - d̂)² = d̂² (d/d̂ - 1)² + // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² return dhat * dhat; } }; @@ -198,4 +198,113 @@ class NormalizedClampedLogBarrier : public ClampedLogBarrier { double units(const double dhat) const override { return 1.0; } }; +// ============================================================================ +// Quadratic log barrier functions from [Huang et al. 2024] +// ============================================================================ + +/// @brief Clamped log barrier with a quadratic log term from [Huang et al. 2024]. +class ClampedLogSqBarrier : public Barrier { +public: + ClampedLogSqBarrier() = default; + + /// @brief Function that grows to infinity as d approaches 0 from the right. + /// + /// \f\[ + /// b(d) = (d-\hat{d})^2\ln^2\left(\frac{d}{\hat{d}}\right) + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The value of the barrier function at d. + double operator()(const double d, const double dhat) const override; + + /// @brief Derivative of the barrier function. + /// + /// \f\[ + /// b'(d) = 2 (d - \hat{d}) \ln\left(\frac{d}{\hat{d}}\right) + /// \left[\ln\left(\frac{d}{\hat{d}}\right) + \frac{d - + /// \hat{d}}{d}\right] + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The derivative of the barrier wrt d. + double first_derivative(const double d, const double dhat) const override; + + /// @brief Second derivative of the barrier function. + /// + /// \f\[ + /// b''(d) = 2 \left(\ln^2\left(\frac{d}{\hat{d}}\right) - \left( + /// \ln\left(\frac{d}{\hat{d}}\right) - 1\right) \frac{(\hat{d} - + /// d)^2}{d^2} - 4 \ln\left(\frac{d}{\hat{d}}\right) \frac{\hat{d} - + /// d}{d}\right) + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The second derivative of the barrier wrt d. + double second_derivative(const double d, const double dhat) const override; + + /// @brief Get the units of the barrier function. + /// @param dhat The activation distance of the barrier. + /// @return The units of the barrier function. + double units(const double dhat) const override + { + // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² + return dhat * dhat; + } +}; + +// ============================================================================ +// Cubic barrier from [Ando 2024] +// ============================================================================ + +/// @brief Cubic barrier function from [Ando 2024]. +class CubicBarrier : public Barrier { +public: + CubicBarrier() = default; + + /// @brief Weak barrier function. + /// + /// \f\[ + /// b(d) = -\frac{2}{3\hat{d}} (d - \hat{d})^3 + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The value of the barrier function at d. + double operator()(const double d, const double dhat) const override; + + /// @brief Derivative of the barrier function. + /// + /// \f\[ + /// b'(d) = -2 (d - \hat{d})^2 + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The derivative of the barrier wrt d. + double first_derivative(const double d, const double dhat) const override; + + /// @brief Second derivative of the barrier function. + /// + /// \f\[ + /// b''(d) = -4 (d - \hat{d}) + /// \f\] + /// + /// @param d The distance. + /// @param dhat Activation distance of the barrier. + /// @return The second derivative of the barrier wrt d. + double second_derivative(const double d, const double dhat) const override; + + /// @brief Get the units of the barrier function. + /// @param dhat The activation distance of the barrier. + /// @return The units of the barrier function. + double units(const double dhat) const override + { + // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² + return dhat * dhat; + } +}; + } // namespace ipc diff --git a/src/ipc/candidates/collision_stencil.hpp b/src/ipc/candidates/collision_stencil.hpp index be63ce3f1..c49446cf4 100644 --- a/src/ipc/candidates/collision_stencil.hpp +++ b/src/ipc/candidates/collision_stencil.hpp @@ -80,6 +80,45 @@ class CollisionStencil { return x; } + /// @brief Compute the distance of the stencil. + /// @param vertices Collision mesh vertices + /// @param edges Collision mesh edges + /// @param faces Collision mesh faces + /// @return Distance of the stencil. + double compute_distance( + const Eigen::MatrixXd& vertices, + const Eigen::MatrixXi& edges, + const Eigen::MatrixXi& faces) const + { + return compute_distance(dof(vertices, edges, faces)); + } + + /// @brief Compute the distance gradient of the stencil w.r.t. the stencil's vertex positions. + /// @param vertices Collision mesh vertices + /// @param edges Collision mesh edges + /// @param faces Collision mesh faces + /// @return Distance gradient of the stencil w.r.t. the stencil's vertex positions. + VectorMax12d compute_distance_gradient( + const Eigen::MatrixXd& vertices, + const Eigen::MatrixXi& edges, + const Eigen::MatrixXi& faces) const + { + return compute_distance_gradient(dof(vertices, edges, faces)); + } + + /// @brief Compute the distance Hessian of the stencil w.r.t. the stencil's vertex positions. + /// @param vertices Collision mesh vertices + /// @param edges Collision mesh edges + /// @param faces Collision mesh faces + /// @return Distance Hessian of the stencil w.r.t. the stencil's vertex positions. + MatrixMax12d compute_distance_hessian( + const Eigen::MatrixXd& vertices, + const Eigen::MatrixXi& edges, + const Eigen::MatrixXi& faces) const + { + return compute_distance_hessian(dof(vertices, edges, faces)); + } + // ---------------------------------------------------------------------- // NOTE: The following functions take stencil vertices as output by dof() // ---------------------------------------------------------------------- diff --git a/src/ipc/candidates/edge_edge.hpp b/src/ipc/candidates/edge_edge.hpp index 30fea302e..a2e729072 100644 --- a/src/ipc/candidates/edge_edge.hpp +++ b/src/ipc/candidates/edge_edge.hpp @@ -26,6 +26,10 @@ class EdgeEdgeCandidate : public ContinuousCollisionCandidate { edges(edge1_id, 0), edges(edge1_id, 1) } }; } + using CollisionStencil::compute_distance; + using CollisionStencil::compute_distance_gradient; + using CollisionStencil::compute_distance_hessian; + double compute_distance(const VectorMax12d& positions) const override; VectorMax12d diff --git a/src/ipc/candidates/edge_vertex.hpp b/src/ipc/candidates/edge_vertex.hpp index 0115ab034..464904fca 100644 --- a/src/ipc/candidates/edge_vertex.hpp +++ b/src/ipc/candidates/edge_vertex.hpp @@ -25,6 +25,10 @@ class EdgeVertexCandidate : public ContinuousCollisionCandidate { return { { vertex_id, edges(edge_id, 0), edges(edge_id, 1), -1 } }; } + using CollisionStencil::compute_distance; + using CollisionStencil::compute_distance_gradient; + using CollisionStencil::compute_distance_hessian; + double compute_distance(const VectorMax12d& positions) const override; VectorMax12d diff --git a/src/ipc/candidates/face_vertex.hpp b/src/ipc/candidates/face_vertex.hpp index 1b4dbaab8..89a098b8f 100644 --- a/src/ipc/candidates/face_vertex.hpp +++ b/src/ipc/candidates/face_vertex.hpp @@ -26,6 +26,10 @@ class FaceVertexCandidate : public ContinuousCollisionCandidate { faces(face_id, 2) } }; } + using CollisionStencil::compute_distance; + using CollisionStencil::compute_distance_gradient; + using CollisionStencil::compute_distance_hessian; + double compute_distance(const VectorMax12d& positions) const override; VectorMax12d diff --git a/src/ipc/candidates/vertex_vertex.hpp b/src/ipc/candidates/vertex_vertex.hpp index 767baa94e..e4062a038 100644 --- a/src/ipc/candidates/vertex_vertex.hpp +++ b/src/ipc/candidates/vertex_vertex.hpp @@ -29,6 +29,10 @@ class VertexVertexCandidate : public ContinuousCollisionCandidate { return { { vertex0_id, vertex1_id, -1, -1 } }; } + using CollisionStencil::compute_distance; + using CollisionStencil::compute_distance_gradient; + using CollisionStencil::compute_distance_hessian; + double compute_distance(const VectorMax12d& positions) const override; VectorMax12d diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index d4dcc260b..b468b171c 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -161,9 +161,9 @@ void NormalPotential::shape_derivative( } if (collision.weight_gradient.nonZeros()) { - VectorMax12d grad_b = gradient(collision, positions); + VectorMax12d grad_f = gradient(collision, positions); assert(collision.weight != 0); - grad_b.array() /= collision.weight; // remove weight + grad_f.array() /= collision.weight; // remove weight for (int i = 0; i < collision.num_vertices(); i++) { for (int d = 0; d < dim; d++) { @@ -171,7 +171,7 @@ void NormalPotential::shape_derivative( for (Itr j(collision.weight_gradient); j; ++j) { out.emplace_back( vertex_ids[i] * dim + d, j.index(), - grad_b[dim * i + d] * j.value()); + grad_f[dim * i + d] * j.value()); } } } diff --git a/tests/src/tests/barrier/test_barrier.cpp b/tests/src/tests/barrier/test_barrier.cpp index 5c304c1c3..e85ab5eea 100644 --- a/tests/src/tests/barrier/test_barrier.cpp +++ b/tests/src/tests/barrier/test_barrier.cpp @@ -67,6 +67,11 @@ TEST_CASE("Barrier derivatives", "[barrier]") { barrier = std::make_unique(use_dist_sqr); } + SECTION("ClampedLogSq") + { + barrier = std::make_unique(); + } + SECTION("Cubic") { barrier = std::make_unique(); } if (use_dist_sqr) { d_vec *= d; From 219f10b0a9a014384c96f48c04b0860798c8bb7b Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 24 Jan 2025 16:28:16 -0500 Subject: [PATCH 02/10] Add compute_coefficients to CollisionStencil --- python/src/utils/interval.cpp | 2 + src/ipc/barrier/barrier.hpp | 13 +-- src/ipc/candidates/collision_stencil.hpp | 19 +++++ src/ipc/candidates/edge_edge.cpp | 58 +++++++++++++ src/ipc/candidates/edge_edge.hpp | 4 + src/ipc/candidates/edge_vertex.cpp | 35 ++++++++ src/ipc/candidates/edge_vertex.hpp | 4 + src/ipc/candidates/face_vertex.cpp | 51 +++++++++++ src/ipc/candidates/face_vertex.hpp | 4 + src/ipc/candidates/vertex_vertex.cpp | 8 ++ src/ipc/candidates/vertex_vertex.hpp | 4 + src/ipc/collisions/normal/plane_vertex.cpp | 8 ++ src/ipc/collisions/normal/plane_vertex.hpp | 6 ++ src/ipc/utils/eigen_ext.hpp | 31 +++++-- tests/src/tests/candidates/CMakeLists.txt | 1 + .../tests/candidates/test_coefficients.cpp | 84 +++++++++++++++++++ tests/src/tests/ccd/test_nonlinear_ccd.cpp | 4 + 17 files changed, 324 insertions(+), 12 deletions(-) create mode 100644 tests/src/tests/candidates/test_coefficients.cpp diff --git a/python/src/utils/interval.cpp b/python/src/utils/interval.cpp index 11dfeb8f0..9301cd9d0 100644 --- a/python/src/utils/interval.cpp +++ b/python/src/utils/interval.cpp @@ -4,7 +4,9 @@ #include namespace py = pybind11; +#ifdef IPC_TOOLKIT_WITH_FILIB using namespace filib; +#endif // ============================================================================ diff --git a/src/ipc/barrier/barrier.hpp b/src/ipc/barrier/barrier.hpp index 014a970c6..9c4c50fbc 100644 --- a/src/ipc/barrier/barrier.hpp +++ b/src/ipc/barrier/barrier.hpp @@ -141,9 +141,9 @@ class ClampedLogBarrier : public Barrier { // ============================================================================ /// @brief Normalized barrier function from [Li et al. 2023]. -class NormalizedClampedLogBarrier : public ClampedLogBarrier { +template class NormalizedBarrier : public BarrierT { public: - NormalizedClampedLogBarrier() = default; + NormalizedBarrier() = default; /// @brief Function that grows to infinity as d approaches 0 from the right. /// @@ -157,7 +157,7 @@ class NormalizedClampedLogBarrier : public ClampedLogBarrier { /// @return The value of the barrier function at d. double operator()(const double d, const double dhat) const override { - return ClampedLogBarrier::operator()(d / dhat, 1.0); + return BarrierT::operator()(d / dhat, 1.0); } /// @brief Derivative of the barrier function. @@ -173,7 +173,7 @@ class NormalizedClampedLogBarrier : public ClampedLogBarrier { /// @return The derivative of the barrier wrt d. double first_derivative(const double d, const double dhat) const override { - return ClampedLogBarrier::first_derivative(d / dhat, 1.0) / dhat; + return BarrierT::first_derivative(d / dhat, 1.0) / dhat; } /// @brief Second derivative of the barrier function. @@ -188,8 +188,7 @@ class NormalizedClampedLogBarrier : public ClampedLogBarrier { /// @return The second derivative of the barrier wrt d. double second_derivative(const double d, const double dhat) const override { - return ClampedLogBarrier::second_derivative(d / dhat, 1.0) - / (dhat * dhat); + return BarrierT::second_derivative(d / dhat, 1.0) / (dhat * dhat); } /// @brief Get the units of the barrier function. @@ -198,6 +197,8 @@ class NormalizedClampedLogBarrier : public ClampedLogBarrier { double units(const double dhat) const override { return 1.0; } }; +using NormalizedClampedLogBarrier = NormalizedBarrier; + // ============================================================================ // Quadratic log barrier functions from [Huang et al. 2024] // ============================================================================ diff --git a/src/ipc/candidates/collision_stencil.hpp b/src/ipc/candidates/collision_stencil.hpp index c49446cf4..67eabd454 100644 --- a/src/ipc/candidates/collision_stencil.hpp +++ b/src/ipc/candidates/collision_stencil.hpp @@ -119,6 +119,19 @@ class CollisionStencil { return compute_distance_hessian(dof(vertices, edges, faces)); } + /// @brief Compute the coefficients of the stencil s.t. d(x) = ‖∑ cᵢ xᵢ‖². + /// @param vertices Collision mesh vertices + /// @param edges Collision mesh edges + /// @param faces Collision mesh faces + /// @return Coefficients of the stencil. + VectorMax4d compute_coefficients( + const Eigen::MatrixXd& vertices, + const Eigen::MatrixXi& edges, + const Eigen::MatrixXi& faces) const + { + return compute_coefficients(dof(vertices, edges, faces)); + } + // ---------------------------------------------------------------------- // NOTE: The following functions take stencil vertices as output by dof() // ---------------------------------------------------------------------- @@ -142,6 +155,12 @@ class CollisionStencil { /// @return Distance Hessian of the stencil w.r.t. the stencil's vertex positions. virtual MatrixMax12d compute_distance_hessian(const VectorMax12d& positions) const = 0; + + /// @brief Compute the coefficients of the stencil s.t. d(x) = ‖∑ cᵢ xᵢ‖². + /// @param positions Stencil's vertex positions. + /// @return Coefficients of the stencil. + virtual VectorMax4d + compute_coefficients(const VectorMax12d& positions) const = 0; }; } // namespace ipc \ No newline at end of file diff --git a/src/ipc/candidates/edge_edge.cpp b/src/ipc/candidates/edge_edge.cpp index 768ab3e94..7036d736c 100644 --- a/src/ipc/candidates/edge_edge.cpp +++ b/src/ipc/candidates/edge_edge.cpp @@ -1,6 +1,7 @@ #include "edge_edge.hpp" #include +#include namespace ipc { @@ -36,6 +37,63 @@ EdgeEdgeCandidate::compute_distance_hessian(const VectorMax12d& positions) const positions.tail<3>(), known_dtype()); } +VectorMax4d +EdgeEdgeCandidate::compute_coefficients(const VectorMax12d& positions) const +{ + assert(positions.size() == 12); + const Eigen::Ref ea0 = positions.head<3>(); + const Eigen::Ref ea1 = positions.segment<3>(3); + const Eigen::Ref eb0 = positions.segment<3>(6); + const Eigen::Ref eb1 = positions.tail<3>(); + + // Project the point inside the triangle + auto dtype = known_dtype(); + if (dtype == EdgeEdgeDistanceType::AUTO) { + dtype = edge_edge_distance_type(ea0, ea1, eb0, eb1); + } + + VectorMax4d coeffs(4); + switch (dtype) { + case EdgeEdgeDistanceType::EA0_EB0: + coeffs << 1.0, 0.0, -1.0, 0.0; + break; + case EdgeEdgeDistanceType::EA0_EB1: + coeffs << 1.0, 0.0, 0.0, -1.0; + break; + case EdgeEdgeDistanceType::EA1_EB0: + coeffs << 0.0, 1.0, -1.0, 0.0; + break; + case EdgeEdgeDistanceType::EA1_EB1: + coeffs << 0.0, 1.0, 0.0, -1.0; + break; + case EdgeEdgeDistanceType::EA_EB0: { + const double alpha = point_edge_closest_point(eb0, ea0, ea1); + coeffs << 1.0 - alpha, alpha, -1.0, 0; + } break; + case EdgeEdgeDistanceType::EA_EB1: { + const double alpha = point_edge_closest_point(eb1, ea0, ea1); + coeffs << 1.0 - alpha, alpha, 0, -1.0; + } break; + case EdgeEdgeDistanceType::EA0_EB: { + const double alpha = point_edge_closest_point(ea0, eb0, eb1); + coeffs << 1.0, 0, -1.0 + alpha, -alpha; + } break; + case EdgeEdgeDistanceType::EA1_EB: { + const double alpha = point_edge_closest_point(ea1, eb0, eb1); + coeffs << 0, 1.0, -1.0 + alpha, -alpha; + } break; + case EdgeEdgeDistanceType::EA_EB: { + const Eigen::Vector2d ab = edge_edge_closest_point(ea0, ea1, eb0, eb1); + coeffs << 1.0 - ab[0], ab[0], -1.0 + ab[1], -ab[1]; + } break; + default: + assert(false); + break; + } + + return coeffs; +} + bool EdgeEdgeCandidate::ccd( const VectorMax12d& vertices_t0, const VectorMax12d& vertices_t1, diff --git a/src/ipc/candidates/edge_edge.hpp b/src/ipc/candidates/edge_edge.hpp index a2e729072..650a38485 100644 --- a/src/ipc/candidates/edge_edge.hpp +++ b/src/ipc/candidates/edge_edge.hpp @@ -26,6 +26,7 @@ class EdgeEdgeCandidate : public ContinuousCollisionCandidate { edges(edge1_id, 0), edges(edge1_id, 1) } }; } + using CollisionStencil::compute_coefficients; using CollisionStencil::compute_distance; using CollisionStencil::compute_distance_gradient; using CollisionStencil::compute_distance_hessian; @@ -38,6 +39,9 @@ class EdgeEdgeCandidate : public ContinuousCollisionCandidate { MatrixMax12d compute_distance_hessian(const VectorMax12d& positions) const override; + VectorMax4d + compute_coefficients(const VectorMax12d& positions) const override; + // ------------------------------------------------------------------------ // ContinuousCollisionCandidate diff --git a/src/ipc/candidates/edge_vertex.cpp b/src/ipc/candidates/edge_vertex.cpp index bc54add49..15e6d1ac2 100644 --- a/src/ipc/candidates/edge_vertex.cpp +++ b/src/ipc/candidates/edge_vertex.cpp @@ -1,6 +1,7 @@ #include "edge_vertex.hpp" #include +#include #include @@ -42,6 +43,40 @@ MatrixMax12d EdgeVertexCandidate::compute_distance_hessian( known_dtype()); } +VectorMax4d +EdgeVertexCandidate::compute_coefficients(const VectorMax12d& positions) const +{ + assert(positions.size() == 6 || positions.size() == 9); + const int dim = this->dim(positions.size()); + const Eigen::Ref p = positions.head(dim); + const Eigen::Ref t0 = positions.segment(dim, dim); + const Eigen::Ref t1 = positions.tail(dim); + + auto dtype = known_dtype(); + if (dtype == PointEdgeDistanceType::AUTO) { + dtype = point_edge_distance_type(p, t0, t1); + } + + VectorMax4d coeffs(3); + switch (dtype) { + case PointEdgeDistanceType::P_E0: + coeffs << 1.0, -1.0, 0.0; + break; + case PointEdgeDistanceType::P_E1: + coeffs << 1.0, 0.0, -1.0; + break; + case PointEdgeDistanceType::P_E: { + const double alpha = point_edge_closest_point(p, t0, t1); + coeffs << 1.0, -1.0 + alpha, -alpha; + } break; + default: + assert(false); + break; + } + + return coeffs; +} + bool EdgeVertexCandidate::ccd( const VectorMax12d& vertices_t0, const VectorMax12d& vertices_t1, diff --git a/src/ipc/candidates/edge_vertex.hpp b/src/ipc/candidates/edge_vertex.hpp index 464904fca..78ba306aa 100644 --- a/src/ipc/candidates/edge_vertex.hpp +++ b/src/ipc/candidates/edge_vertex.hpp @@ -25,6 +25,7 @@ class EdgeVertexCandidate : public ContinuousCollisionCandidate { return { { vertex_id, edges(edge_id, 0), edges(edge_id, 1), -1 } }; } + using CollisionStencil::compute_coefficients; using CollisionStencil::compute_distance; using CollisionStencil::compute_distance_gradient; using CollisionStencil::compute_distance_hessian; @@ -37,6 +38,9 @@ class EdgeVertexCandidate : public ContinuousCollisionCandidate { MatrixMax12d compute_distance_hessian(const VectorMax12d& positions) const override; + VectorMax4d + compute_coefficients(const VectorMax12d& positions) const override; + // ------------------------------------------------------------------------ // ContinuousCollisionCandidate diff --git a/src/ipc/candidates/face_vertex.cpp b/src/ipc/candidates/face_vertex.cpp index b18fa6fb6..81349c1af 100644 --- a/src/ipc/candidates/face_vertex.cpp +++ b/src/ipc/candidates/face_vertex.cpp @@ -1,6 +1,7 @@ #include "face_vertex.hpp" #include +#include #include @@ -39,6 +40,56 @@ MatrixMax12d FaceVertexCandidate::compute_distance_hessian( positions.tail<3>(), known_dtype()); } +VectorMax4d +FaceVertexCandidate::compute_coefficients(const VectorMax12d& positions) const +{ + assert(positions.size() == 12); + const Eigen::Ref p = positions.head<3>(); + const Eigen::Ref t0 = positions.segment<3>(3); + const Eigen::Ref t1 = positions.segment<3>(6); + const Eigen::Ref t2 = positions.tail<3>(); + + // Project the point inside the triangle + auto dtype = known_dtype(); + if (dtype == PointTriangleDistanceType::AUTO) { + dtype = point_triangle_distance_type(p, t0, t1, t2); + } + + VectorMax4d coeffs(4); + switch (dtype) { + case PointTriangleDistanceType::P_T0: + coeffs << 1.0, -1.0, 0.0, 0.0; + break; + case PointTriangleDistanceType::P_T1: + coeffs << 1.0, 0.0, -1.0, 0.0; + break; + case PointTriangleDistanceType::P_T2: + coeffs << 1.0, 0.0, 0.0, -1.0; + break; + case PointTriangleDistanceType::P_E0: { + const double alpha = point_edge_closest_point(p, t0, t1); + coeffs << 1.0, -alpha, -1 + alpha, 0.0; + } break; + case PointTriangleDistanceType::P_E1: { + const double alpha = point_edge_closest_point(p, t1, t2); + coeffs << 1.0, 0.0, -alpha, -1 + alpha; + } break; + case PointTriangleDistanceType::P_E2: { + const double alpha = point_edge_closest_point(p, t2, t0); + coeffs << 1.0, -1 + alpha, 0.0, -alpha; + } break; + case PointTriangleDistanceType::P_T: { + const Eigen::Vector2d uv = point_triangle_closest_point(p, t0, t1, t2); + coeffs << 1.0, -1.0 + uv[0] + uv[1], -uv[0], -uv[1]; + } break; + default: + assert(false); + break; + } + + return coeffs; +} + bool FaceVertexCandidate::ccd( const VectorMax12d& vertices_t0, const VectorMax12d& vertices_t1, diff --git a/src/ipc/candidates/face_vertex.hpp b/src/ipc/candidates/face_vertex.hpp index 89a098b8f..9b508e012 100644 --- a/src/ipc/candidates/face_vertex.hpp +++ b/src/ipc/candidates/face_vertex.hpp @@ -26,6 +26,7 @@ class FaceVertexCandidate : public ContinuousCollisionCandidate { faces(face_id, 2) } }; } + using CollisionStencil::compute_coefficients; using CollisionStencil::compute_distance; using CollisionStencil::compute_distance_gradient; using CollisionStencil::compute_distance_hessian; @@ -38,6 +39,9 @@ class FaceVertexCandidate : public ContinuousCollisionCandidate { MatrixMax12d compute_distance_hessian(const VectorMax12d& positions) const override; + VectorMax4d + compute_coefficients(const VectorMax12d& positions) const override; + // ------------------------------------------------------------------------ // ContinuousCollisionCandidate diff --git a/src/ipc/candidates/vertex_vertex.cpp b/src/ipc/candidates/vertex_vertex.cpp index 3f4a039d8..14b26ccc7 100644 --- a/src/ipc/candidates/vertex_vertex.cpp +++ b/src/ipc/candidates/vertex_vertex.cpp @@ -38,6 +38,14 @@ MatrixMax12d VertexVertexCandidate::compute_distance_hessian( positions.head(dim), positions.tail(dim)); } +VectorMax4d +VertexVertexCandidate::compute_coefficients(const VectorMax12d& positions) const +{ + VectorMax4d coeffs(2); + coeffs << 1.0, -1.0; + return coeffs; +} + bool VertexVertexCandidate::ccd( const VectorMax12d& vertices_t0, const VectorMax12d& vertices_t1, diff --git a/src/ipc/candidates/vertex_vertex.hpp b/src/ipc/candidates/vertex_vertex.hpp index e4062a038..db23a90e7 100644 --- a/src/ipc/candidates/vertex_vertex.hpp +++ b/src/ipc/candidates/vertex_vertex.hpp @@ -29,6 +29,7 @@ class VertexVertexCandidate : public ContinuousCollisionCandidate { return { { vertex0_id, vertex1_id, -1, -1 } }; } + using CollisionStencil::compute_coefficients; using CollisionStencil::compute_distance; using CollisionStencil::compute_distance_gradient; using CollisionStencil::compute_distance_hessian; @@ -41,6 +42,9 @@ class VertexVertexCandidate : public ContinuousCollisionCandidate { MatrixMax12d compute_distance_hessian(const VectorMax12d& positions) const override; + VectorMax4d + compute_coefficients(const VectorMax12d& positions) const override; + // ------------------------------------------------------------------------ // ContinuousCollisionCandidate diff --git a/src/ipc/collisions/normal/plane_vertex.cpp b/src/ipc/collisions/normal/plane_vertex.cpp index 63ce2882a..27b3a53c8 100644 --- a/src/ipc/collisions/normal/plane_vertex.cpp +++ b/src/ipc/collisions/normal/plane_vertex.cpp @@ -35,4 +35,12 @@ MatrixMax12d PlaneVertexNormalCollision::compute_distance_hessian( return point_plane_distance_hessian(point, plane_origin, plane_normal); } +VectorMax4d PlaneVertexNormalCollision::compute_coefficients( + const VectorMax12d& positions) const +{ + VectorMax4d coeffs(1); + coeffs << 1.0; + return coeffs; +} + } // namespace ipc diff --git a/src/ipc/collisions/normal/plane_vertex.hpp b/src/ipc/collisions/normal/plane_vertex.hpp index 258928e15..1d9f46c94 100644 --- a/src/ipc/collisions/normal/plane_vertex.hpp +++ b/src/ipc/collisions/normal/plane_vertex.hpp @@ -38,6 +38,12 @@ class PlaneVertexNormalCollision : public NormalCollision { MatrixMax12d compute_distance_hessian(const VectorMax12d& point) const override; + /// @brief Compute the coefficients of the stencil. + /// @param positions Vertex positions. + /// @return Coefficients of the stencil. + VectorMax4d + compute_coefficients(const VectorMax12d& positions) const override; + /// @brief The plane's origin. VectorMax3d plane_origin; diff --git a/src/ipc/utils/eigen_ext.hpp b/src/ipc/utils/eigen_ext.hpp index 16d1784d3..87b9acd8e 100644 --- a/src/ipc/utils/eigen_ext.hpp +++ b/src/ipc/utils/eigen_ext.hpp @@ -30,31 +30,39 @@ using Matrix12d = Eigen::Matrix; template using VectorMax2 = Vector; /// @brief A dynamic size matrix with a fixed maximum size of 3×1 template using VectorMax3 = Vector; +/// @brief A dynamic size matrix with a fixed maximum size of 4×1 +template using VectorMax4 = Vector; +/// @brief A dynamic size matrix with a fixed maximum size of 6×1 +template using VectorMax6 = Vector; +/// @brief A dynamic size matrix with a fixed maximum size of 9×1 +template using VectorMax9 = Vector; +/// @brief A dynamic size matrix with a fixed maximum size of 12×1 +template using VectorMax12 = Vector; + /// @brief A dynamic size matrix with a fixed maximum size of 2×1 using VectorMax2d = VectorMax2; /// @brief A dynamic size matrix with a fixed maximum size of 3×1 using VectorMax3d = VectorMax3; /// @brief A dynamic size matrix with a fixed maximum size of 3×1 using VectorMax3i = VectorMax3; -/// @brief A dynamic size matrix with a fixed maximum size of 6×1 -template using VectorMax6 = Vector; +/// @brief A dynamic size matrix with a fixed maximum size of 4×1 +using VectorMax4d = VectorMax4; +/// @brief A dynamic size matrix with a fixed maximum size of 4×1 +using VectorMax4i = VectorMax4; /// @brief A dynamic size matrix with a fixed maximum size of 6×1 using VectorMax6d = VectorMax6; /// @brief A dynamic size matrix with a fixed maximum size of 6×1 using VectorMax6b = VectorMax6; /// @brief A dynamic size matrix with a fixed maximum size of 9×1 -template using VectorMax9 = Vector; -/// @brief A dynamic size matrix with a fixed maximum size of 9×1 using VectorMax9d = VectorMax9; /// @brief A dynamic size matrix with a fixed maximum size of 12×1 -template using VectorMax12 = Vector; -/// @brief A dynamic size matrix with a fixed maximum size of 12×1 using VectorMax12d = VectorMax12; /// @brief A dynamic size matrix with a fixed maximum size of 1×2 template using RowVectorMax2 = RowVector; /// @brief A dynamic size matrix with a fixed maximum size of 1×3 template using RowVectorMax3 = RowVector; + /// @brief A dynamic size matrix with a fixed maximum size of 1×2 using RowVectorMax2d = RowVectorMax2; /// @brief A dynamic size matrix with a fixed maximum size of 1×3 @@ -74,6 +82,7 @@ using MatrixMax = Eigen::Matrix< Eigen::ColMajor, max_rows, max_cols>; + /// @brief A dynamic size matrix with a fixed maximum size of 3×3 template using MatrixMax2 = MatrixMax; /// @brief A dynamic size matrix with a fixed maximum size of 3×3 @@ -84,6 +93,7 @@ template using MatrixMax6 = MatrixMax; template using MatrixMax9 = MatrixMax; /// @brief A dynamic size matrix with a fixed maximum size of 12×12 template using MatrixMax12 = MatrixMax; + /// @brief A dynamic size matrix with a fixed maximum size of 3×3 using MatrixMax2d = MatrixMax2; /// @brief A dynamic size matrix with a fixed maximum size of 3×3 @@ -94,6 +104,7 @@ using MatrixMax6d = MatrixMax6; using MatrixMax9d = MatrixMax9; /// @brief A dynamic size matrix with a fixed maximum size of 12×12 using MatrixMax12d = MatrixMax12; + /// @brief A dynamic size diagonal matrix using DiagonalMatrixXd = Eigen::DiagonalMatrix; /// @brief A dynamic size diagonal matrix with a fixed maximum size of 6×6 @@ -105,12 +116,20 @@ using ArrayMax2 = Eigen::Array; /// @brief A dynamic size array with a fixed maximum size of 2×1 template using ArrayMax3 = Eigen::Array; +/// @brief A dynamic size array with a fixed maximum size of 4×1 +template +using ArrayMax4 = Eigen::Array; + /// @brief A dynamic size array with a fixed maximum size of 2×1 using ArrayMax2d = ArrayMax2; /// @brief A dynamic size array with a fixed maximum size of 3×1 using ArrayMax3d = ArrayMax3; /// @brief A dynamic size array with a fixed maximum size of 3×1 using ArrayMax3i = ArrayMax3; +/// @brief A dynamic size array with a fixed maximum size of 4×1 +using ArrayMax4d = ArrayMax4; +/// @brief A dynamic size array with a fixed maximum size of 4×1 +using ArrayMax4i = ArrayMax4; /// @brief Matrix projection onto positive definite cone /// @param A Symmetric matrix to project diff --git a/tests/src/tests/candidates/CMakeLists.txt b/tests/src/tests/candidates/CMakeLists.txt index 523cc4225..8c3f5a90b 100644 --- a/tests/src/tests/candidates/CMakeLists.txt +++ b/tests/src/tests/candidates/CMakeLists.txt @@ -1,6 +1,7 @@ set(SOURCES # Tests test_candidates.cpp + test_coefficients.cpp # Benchmarks diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp new file mode 100644 index 000000000..4ae0d7cac --- /dev/null +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -0,0 +1,84 @@ +#include +#include + +#include + +#include + +using namespace ipc; + +TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") +{ + Eigen::MatrixXd V(4, 3); + V << 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1; + + SECTION("default") { } + SECTION("random") { V.row(3).setRandom(); } + SECTION("off triangle 0") { V.row(3) << 2, 0.5, 1; } + SECTION("off triangle 1") { V.row(3) << 2, 2, 1; } + SECTION("off triangle 2") { V.row(3) << 0.5, 2, 1; } + SECTION("off triangle 3") { V.row(3) << -1, 2, 1; } + SECTION("off triangle 4") { V.row(3) << -1, 0.5, 1; } + SECTION("off triangle 5") { V.row(3) << -1, -1, 1; } + SECTION("off triangle 6") { V.row(3) << 0.5, -1, 1; } + SECTION("off triangle 7") { V.row(3) << 2, -0.5, 1; } + + Eigen::MatrixXi F(1, 3); + F << 0, 1, 2; + + Eigen::MatrixXi E; + igl::edges(F, E); + + FaceVertexCandidate fv(0, 3); + + VectorMax4d coeffs = fv.compute_coefficients(V, E, F); + CAPTURE(V.row(3), coeffs.transpose()); + + CHECK(-coeffs[1] == Catch::Approx(1 + coeffs[2] + coeffs[3])); + + std::array vertices = fv.vertices(V, E, F); + Eigen::Vector3d n = Eigen::Vector3d::Zero(); + for (int i = 0; i < fv.num_vertices(); i++) { + n += coeffs[i] * vertices[i]; + } + + CHECK(n.squaredNorm() == Catch::Approx(fv.compute_distance(V, E, F))); +} + +TEST_CASE("Edge-edge collision stencil coeffs.", "[ee][stencil][coeffs]") +{ + Eigen::MatrixXd V(4, 3); + V << -1, 0, 0, /**/ 1, 0, 0, /**/ 0, -1, 0, /**/ 0, 1, 0; + + // clang-format off + SECTION("default") { } + SECTION("random") { V.bottomRows<2>().setRandom(); } + SECTION("quad 0") { V.block<2, 1>(2, 0).array() += 2; } + SECTION("quad 1") { V.block<2, 2>(2, 0).array() += 2; } + SECTION("quad 2") { V.block<2, 1>(2, 1).array() += 2; } + SECTION("quad 3") { V.bottomRows<2>().rowwise() += Eigen::RowVector3d(-2, 2, 0); } + SECTION("quad 4") { V.block<2, 1>(2, 0).array() -= 2; } + SECTION("quad 5") { V.block<2, 2>(2, 0).array() -= 2; } + SECTION("quad 6") { V.block<2, 1>(2, 1).array() -= 2; } + SECTION("quad 7") { V.bottomRows<2>().rowwise() += Eigen::RowVector3d(2, -2, 0); } + // clang-format on + + Eigen::MatrixXi E(2, 2), F; + E << 0, 1, /**/ 2, 3; + + EdgeEdgeCandidate ee(0, 1); + + VectorMax4d coeffs = ee.compute_coefficients(V, E, F); + CAPTURE( + V, coeffs.transpose(), + edge_edge_distance_type(V.row(0), V.row(1), V.row(2), V.row(3))); + + std::array vertices = ee.vertices(V, E, F); + Eigen::Vector3d n = Eigen::Vector3d::Zero(); + for (int i = 0; i < ee.num_vertices(); i++) { + n += coeffs[i] * vertices[i]; + } + + CHECK(n.squaredNorm() == Catch::Approx(ee.compute_distance(V, E, F))); + logger().critical("{}", coeffs.transpose()); +} \ No newline at end of file diff --git a/tests/src/tests/ccd/test_nonlinear_ccd.cpp b/tests/src/tests/ccd/test_nonlinear_ccd.cpp index fc15a85f2..4d6bccf30 100644 --- a/tests/src/tests/ccd/test_nonlinear_ccd.cpp +++ b/tests/src/tests/ccd/test_nonlinear_ccd.cpp @@ -66,6 +66,7 @@ class RotationalTrajectory : virtual public NonlinearTrajectory { } }; +#ifdef IPC_TOOLKIT_WITH_FILIB class IntervalRotationalTrajectory : public RotationalTrajectory, public IntervalNonlinearTrajectory { public: @@ -90,6 +91,7 @@ class IntervalRotationalTrajectory : public RotationalTrajectory, return IntervalNonlinearTrajectory::max_distance_from_linear(t0, t1); } }; +#endif class StaticTrajectory : public NonlinearTrajectory { public: @@ -165,11 +167,13 @@ TEST_CASE("Nonlinear Point-Point CCD", "[ccd][nonlinear][point-point]") p1 = std::make_unique( Eigen::Vector2d(1, 0), Eigen::Vector2d::Zero(), igl::PI); } +#ifdef IPC_TOOLKIT_WITH_FILIB SECTION("Interval max") { p1 = std::make_unique( Eigen::Vector2d(1, 0), Eigen::Vector2d::Zero(), igl::PI); } +#endif double toi; bool collision = point_point_nonlinear_ccd( From 4be3d07d80642a766b97c801d6772ef3199fb46f Mon Sep 17 00:00:00 2001 From: zachferguson Date: Mon, 3 Feb 2025 16:23:58 -0500 Subject: [PATCH 03/10] Add semi_implicit_stiffness --- src/ipc/barrier/adaptive_stiffness.cpp | 113 ++++++++++++++++++ src/ipc/barrier/adaptive_stiffness.hpp | 35 ++++++ .../tests/candidates/test_coefficients.cpp | 1 - 3 files changed, 148 insertions(+), 1 deletion(-) diff --git a/src/ipc/barrier/adaptive_stiffness.cpp b/src/ipc/barrier/adaptive_stiffness.cpp index e9fc00baa..38f53223f 100644 --- a/src/ipc/barrier/adaptive_stiffness.cpp +++ b/src/ipc/barrier/adaptive_stiffness.cpp @@ -3,6 +3,8 @@ #include "adaptive_stiffness.hpp" #include +#include +#include #include // std::min/max #include @@ -72,4 +74,115 @@ double update_barrier_stiffness( return barrier_stiffness; } +// ----------------------------------------------------------------------------- +// +// Based on `compute_stiffness()` in `src/cpp/barrier/barrier.cu` of +// (ppf-contact-solver)[https://github.com/st-tech/ppf-contact-solver] +// +// Original license: +// File: barrier.cu +// Author: Ryoichi Ando (ryoichi.ando@zozo.com) +// License: Apache v2.0 +// + +double semi_implicit_stiffness( + const CollisionStencil& stencil, + const std::array& vertex_ids, + const VectorMax12d& vertices, + const VectorMax4d& mass, + const Eigen::SparseMatrix& hess, + const double dmin) +{ + unsigned N = stencil.num_vertices(); + assert(vertices.size() % N == 0); + unsigned dim = vertices.size() / N; + + const VectorMax4d value = stencil.compute_coefficients(vertices); + + // Compute the contact normal (i.e., the vector from the ) + VectorMax3d normal = VectorMax3d::Zero(dim); + for (unsigned i = 0; i < N; ++i) { + normal += value[i] * vertices.segment(dim * i, dim); + } + + MatrixMax12d local_hess = MatrixMax12d::Zero(dim * N, dim * N); + + // H + for (unsigned i = 0; i < N; ++i) { + for (unsigned j = 0; j < N; ++j) { + for (unsigned k = 0; k < dim; ++k) { + for (unsigned l = 0; l < dim; ++l) { + local_hess(dim * i + k, dim * j + l) = hess.coeff( + dim * vertex_ids[i] + k, dim * vertex_ids[j] + l); + } + } + } + } + + // mᵢ / d² + const double distance = normal.norm() - dmin; + const double distance_sqr = distance * distance; + for (unsigned i = 0; i < N; ++i) { + local_hess.diagonal().segment(dim * i, dim).array() += + mass[i] / distance_sqr; + } + + VectorMax12d w = VectorMax12d::Zero(dim * N); + for (unsigned i = 0; i < N; ++i) { + w.segment(dim * i, dim) = value[i] * normal; + } + w.normalize(); + + return (local_hess * w).dot(w); +} + +template +Eigen::VectorXd semi_implicit_stiffness( + const CollisionMesh& mesh, + const Eigen::MatrixXd& vertices, + const StencilsT& collisions, + const Eigen::VectorXd& vertex_masses, + const Eigen::SparseMatrix& hess, + const double dmin) +{ + Eigen::VectorXd stiffnesses(collisions.size()); + + for (size_t i = 0; i < collisions.size(); i++) { + const CollisionStencil& collision = collisions[i]; + + const std::array vertex_ids = + collision.vertex_ids(mesh.edges(), mesh.faces()); + const VectorMax12d positions = + collision.dof(vertices, mesh.edges(), mesh.faces()); + + VectorMax4d mass(collision.num_vertices()); + for (unsigned j = 0; j < collision.num_vertices(); j++) { + mass[j] = vertex_masses[vertex_ids[j]]; + } + + stiffnesses[i] = semi_implicit_stiffness( + collision, vertex_ids, positions, mass, hess, dmin); + } + + return stiffnesses; +} + +template Eigen::VectorXd semi_implicit_stiffness( + const CollisionMesh& mesh, + const Eigen::MatrixXd& vertices, + const NormalCollisions& collisions, + const Eigen::VectorXd& vertex_masses, + const Eigen::SparseMatrix& hess, + const double dmin); + +template Eigen::VectorXd semi_implicit_stiffness( + const CollisionMesh& mesh, + const Eigen::MatrixXd& vertices, + const Candidates& collisions, + const Eigen::VectorXd& vertex_masses, + const Eigen::SparseMatrix& hess, + const double dmin); + +// ----------------------------------------------------------------------------- + } // namespace ipc diff --git a/src/ipc/barrier/adaptive_stiffness.hpp b/src/ipc/barrier/adaptive_stiffness.hpp index 1ffb49318..dfed422d4 100644 --- a/src/ipc/barrier/adaptive_stiffness.hpp +++ b/src/ipc/barrier/adaptive_stiffness.hpp @@ -2,7 +2,9 @@ #pragma once +#include #include +#include #include @@ -48,4 +50,37 @@ double update_barrier_stiffness( const double dhat_epsilon_scale = 1e-9, const double dmin = 0); +/// @brief Compute the semi-implicit stiffness for a single collision. +/// @param stencil Collision stencil. +/// @param vertex_ids Vertex indices of the collision. +/// @param vertices Vertex positions. +/// @param mass Vertex masses. +/// @param hess Hessian of the energy function. +/// @param dmin Minimum distance between elements. +/// @return The semi-implicit stiffness. +double semi_implicit_stiffness( + const CollisionStencil& stencil, + const std::array& vertex_ids, + const VectorMax12d& vertices, + const VectorMax4d& mass, + const Eigen::SparseMatrix& hess, + double dmin); + +/// @brief Compute the semi-implicit stiffnesses for all collisions. +/// @param mesh Collision mesh. +/// @param vertices Vertex positions. +/// @param collisions Normal collisions. +/// @param vertex_masses Lumped vertex masses. +/// @param hess Hessian of the elasticity energy function. +/// @param dmin Minimum distance between elements. +/// @return The semi-implicit stiffnesses. +template +Eigen::VectorXd semi_implicit_stiffness( + const CollisionMesh& mesh, + const Eigen::MatrixXd& vertices, + const StencilsT& collisions, + const Eigen::VectorXd& vertex_masses, + const Eigen::SparseMatrix& hess, + double dmin); + } // namespace ipc diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp index 4ae0d7cac..65be5984b 100644 --- a/tests/src/tests/candidates/test_coefficients.cpp +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -80,5 +80,4 @@ TEST_CASE("Edge-edge collision stencil coeffs.", "[ee][stencil][coeffs]") } CHECK(n.squaredNorm() == Catch::Approx(ee.compute_distance(V, E, F))); - logger().critical("{}", coeffs.transpose()); } \ No newline at end of file From 3527ad62e36218b57d7e923bbbe7285313a4f15e Mon Sep 17 00:00:00 2001 From: zachferguson Date: Mon, 3 Feb 2025 16:28:30 -0500 Subject: [PATCH 04/10] Revert math type --- src/ipc/barrier/barrier.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ipc/barrier/barrier.hpp b/src/ipc/barrier/barrier.hpp index 9c4c50fbc..a83f48fef 100644 --- a/src/ipc/barrier/barrier.hpp +++ b/src/ipc/barrier/barrier.hpp @@ -131,7 +131,7 @@ class ClampedLogBarrier : public Barrier { /// @return The units of the barrier function. double units(const double dhat) const override { - // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² + // (d - d̂)² = d̂² (d/d̂ - 1)² return dhat * dhat; } }; @@ -251,7 +251,7 @@ class ClampedLogSqBarrier : public Barrier { /// @return The units of the barrier function. double units(const double dhat) const override { - // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² + // (d - d̂)² = d̂² (d/d̂ - 1)² return dhat * dhat; } }; @@ -303,7 +303,7 @@ class CubicBarrier : public Barrier { /// @return The units of the barrier function. double units(const double dhat) const override { - // (d - \hat{d})² = d̂² (d/\hat{d} - 1)² + // (d - d̂)² = d̂² (d/d̂ - 1)² return dhat * dhat; } }; From 0367a1522c5a1d8fc735a088d46ff44fef9d6473 Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 13:47:15 -0500 Subject: [PATCH 05/10] Fix FaceVertexCandidate::compute_coefficients --- src/ipc/candidates/face_vertex.cpp | 6 +++--- tests/src/tests/candidates/test_coefficients.cpp | 5 +++++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/ipc/candidates/face_vertex.cpp b/src/ipc/candidates/face_vertex.cpp index 81349c1af..9a4c8863e 100644 --- a/src/ipc/candidates/face_vertex.cpp +++ b/src/ipc/candidates/face_vertex.cpp @@ -68,15 +68,15 @@ FaceVertexCandidate::compute_coefficients(const VectorMax12d& positions) const break; case PointTriangleDistanceType::P_E0: { const double alpha = point_edge_closest_point(p, t0, t1); - coeffs << 1.0, -alpha, -1 + alpha, 0.0; + coeffs << 1.0, -1 + alpha, -alpha, 0.0; } break; case PointTriangleDistanceType::P_E1: { const double alpha = point_edge_closest_point(p, t1, t2); - coeffs << 1.0, 0.0, -alpha, -1 + alpha; + coeffs << 1.0, 0.0, -1 + alpha, -alpha; } break; case PointTriangleDistanceType::P_E2: { const double alpha = point_edge_closest_point(p, t2, t0); - coeffs << 1.0, -1 + alpha, 0.0, -alpha; + coeffs << 1.0, -alpha, 0.0, -1 + alpha; } break; case PointTriangleDistanceType::P_T: { const Eigen::Vector2d uv = point_triangle_closest_point(p, t0, t1, t2); diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp index 65be5984b..b134aded9 100644 --- a/tests/src/tests/candidates/test_coefficients.cpp +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -22,6 +22,8 @@ TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") SECTION("off triangle 5") { V.row(3) << -1, -1, 1; } SECTION("off triangle 6") { V.row(3) << 0.5, -1, 1; } SECTION("off triangle 7") { V.row(3) << 2, -0.5, 1; } + SECTION("failure case 1") { V.row(3) << 0.680375, -0.211234, 0.566198; } + SECTION("failure case 1") { V.row(3) << -0.997497, 0.127171, -0.613392; } Eigen::MatrixXi F(1, 3); F << 0, 1, 2; @@ -32,6 +34,8 @@ TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") FaceVertexCandidate fv(0, 3); VectorMax4d coeffs = fv.compute_coefficients(V, E, F); + CAPTURE( + point_triangle_distance_type(V.row(3), V.row(0), V.row(1), V.row(2))); CAPTURE(V.row(3), coeffs.transpose()); CHECK(-coeffs[1] == Catch::Approx(1 + coeffs[2] + coeffs[3])); @@ -41,6 +45,7 @@ TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") for (int i = 0; i < fv.num_vertices(); i++) { n += coeffs[i] * vertices[i]; } + CAPTURE(n.transpose()); CHECK(n.squaredNorm() == Catch::Approx(fv.compute_distance(V, E, F))); } From 370cf080aa26e7321fc894ff1dfeba6efca3e15d Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 16:21:00 -0500 Subject: [PATCH 06/10] Add semi-implicit stiffness to python --- docs/requirements.txt | 10 +- docs/source/Doxyfile | 4 +- docs/source/conf.py | 2 +- docs/source/cpp-api/barrier.rst | 5 + docs/source/python-api/barrier.rst | 6 ++ docs/source/references.bib | 110 ++++++++++++++++++++++ docs/source/refs.bib | 88 ----------------- python/src/barrier/adaptive_stiffness.cpp | 78 +++++++++++++++ src/ipc/barrier/adaptive_stiffness.cpp | 78 ++++++++------- src/ipc/barrier/adaptive_stiffness.hpp | 16 ++-- 10 files changed, 264 insertions(+), 133 deletions(-) create mode 100644 docs/source/references.bib delete mode 100644 docs/source/refs.bib diff --git a/docs/requirements.txt b/docs/requirements.txt index 06339c84c..7977ef7dd 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,13 +1,17 @@ sphinx sphinx-sitemap -sphinx-immaterial==0.11.14 +sphinx-immaterial autoclasstoc # sphinx-autodoc-toolbox breathe nbsphinx pandoc -myst_parser==1.0.0 +myst_parser sphinxcontrib-bibtex sphinxemoji sphinx-last-updated-by-git -sphinx-autobuild \ No newline at end of file +sphinx-autobuild +nbconvert==5.6.1 +ipython_genutils +jinja2==3.0.3 +setuptools \ No newline at end of file diff --git a/docs/source/Doxyfile b/docs/source/Doxyfile index e11032475..29108378d 100644 --- a/docs/source/Doxyfile +++ b/docs/source/Doxyfile @@ -806,7 +806,7 @@ LAYOUT_FILE = # LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the # search path. See also \cite for info how to create references. -CITE_BIB_FILES = +CITE_BIB_FILES = references.bib #--------------------------------------------------------------------------- # Configuration options related to warning and progress messages @@ -2054,7 +2054,7 @@ LATEX_HIDE_INDICES = NO # The default value is: plain. # This tag requires that the tag GENERATE_LATEX is set to YES. -LATEX_BIB_STYLE = plain +LATEX_BIB_STYLE = plainnat # If the LATEX_TIMESTAMP tag is set to YES then the footer of each generated # page will contain the date and time when the page was generated. Setting this diff --git a/docs/source/conf.py b/docs/source/conf.py index 66aab18e0..2bc86127b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -64,7 +64,7 @@ "sphinx_last_updated_by_git", ] -bibtex_bibfiles = ['refs.bib'] +bibtex_bibfiles = ['references.bib'] bibtex_reference_style = 'author_year' bibtex_default_style = 'plain' diff --git a/docs/source/cpp-api/barrier.rst b/docs/source/cpp-api/barrier.rst index 3b243d377..ebebc2851 100644 --- a/docs/source/cpp-api/barrier.rst +++ b/docs/source/cpp-api/barrier.rst @@ -17,6 +17,11 @@ Adaptive Barrier Stiffness .. doxygenfunction:: ipc::initial_barrier_stiffness .. doxygenfunction:: ipc::update_barrier_stiffness +Semi-Implicit Stiffness +~~~~~~~~~~~~~~~~~~~~~~~ + +.. doxygenfunction:: ipc::semi_implicit_stiffness(const CollisionMesh&, const Eigen::MatrixXd&, const StencilsT&, const Eigen::VectorXd&, const Eigen::SparseMatrix&, const double) + Barrier Class ------------- diff --git a/docs/source/python-api/barrier.rst b/docs/source/python-api/barrier.rst index fa0cc329b..8296a9c3d 100644 --- a/docs/source/python-api/barrier.rst +++ b/docs/source/python-api/barrier.rst @@ -17,6 +17,12 @@ Adaptive Barrier Stiffness .. autofunction:: ipctk.initial_barrier_stiffness .. autofunction:: ipctk.update_barrier_stiffness + +Semi-Implicit Stiffness +~~~~~~~~~~~~~~~~~~~~~~~ + +.. autofunction:: ipctk.semi_implicit_stiffness + Barrier Class ------------- diff --git a/docs/source/references.bib b/docs/source/references.bib new file mode 100644 index 000000000..fcaf9a07d --- /dev/null +++ b/docs/source/references.bib @@ -0,0 +1,110 @@ +@article{Li2020IPC, + title = {Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics}, + author = {Minchen Li and Zachary Ferguson and Teseo Schneider and Timothy Langlois and Denis Zorin and Daniele Panozzo and Chenfanfu Jiang and Danny M. Kaufman}, + year = 2020, + journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, + volume = 39, + number = 4, + note = {\url{https://ipc-sim.github.io}}, + articleno = 49 +} +@article{Li2021CIPC, + title = {Codimensional Incremental Potential Contact}, + author = {Minchen Li and Danny M. Kaufman and Chenfanfu Jiang}, + year = 2021, + journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, + volume = 40, + number = 4, + note = {\url{https://ipc-sim.github.io/C-IPC}}, + articleno = 170 +} +@misc{Li2023Convergent, + title = {Convergent Incremental Potential Contact}, + author = {Minchen Li and Zachary Ferguson and Teseo Schneider and Timothy Langlois and Denis Zorin and Daniele Panozzo and Chenfanfu Jiang and Danny M. Kaufman}, + year = 2023, + eprint = {2307.15908}, + archiveprefix = {arXiv}, + primaryclass = {math.NA} +} +@article{Wang2021TightInclusion, + title = {A Large Scale Benchmark and an Inclusion-Based Algorithm for Continuous Collision Detection}, + author = {Bolun Wang and Zachary Ferguson and Teseo Schneider and Xin Jiang and Marco Attene and Daniele Panozzo}, + year = 2021, + month = oct, + journal = {{ACM} Transactions on Graphics}, + volume = 40, + number = 5, + note = {\url{https://continuous-collision-detection.github.io/tight_inclusion/}}, + articleno = 188, + numpages = 16 +} +@misc{Belgrod2023Time, + title = {Time of Impact Dataset for Continuous Collision Detection and a Scalable Conservative Algorithm}, + author = {David Belgrod and Bolun Wang and Zachary Ferguson and Xin Zhao and Marco Attene and Daniele Panozzo and Teseo Schneider}, + year = 2023, + eprint = {2112.06300}, + archiveprefix = {arXiv}, + primaryclass = {cs.GR} +} +@inproceedings{Ferguson2023HighOrderIPC, + title = {High-Order Incremental Potential Contact for Elastodynamic Simulation on Curved Meshes}, + author = {Zachary Ferguson and Pranav Jain and Denis Zorin and Teseo Schneider and Daniele Panozzo}, + year = 2023, + booktitle = {{ACM} {SIGGRAPH} 2023 Conference Proceedings}, + location = {Los Angeles, CA, USA}, + publisher = {Association for Computing Machinery}, + address = {New York, NY, USA}, + series = {SIGGRAPH '23}, + note = {\url{https://zferg.us/research/high-order-ipc/}}, + numpages = 11 +} +@article{Ferguson2021RigidIPC, + title = {Intersection-free Rigid Body Dynamics}, + author = {Zachary Ferguson and Minchen Li and Teseo Schneider and Francisca Gil-Ureta and Timothy Langlois and Chenfanfu Jiang and Denis Zorin and Danny M. Kaufman and Daniele Panozzo}, + year = 2021, + journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, + volume = 40, + number = 4, + note = {\url{https://ipc-sim.github.io/rigid-ipc}}, + articleno = 183 +} +@inproceedings{Chen2024Stabler, + title = {Stabler Neo-Hookean Simulation: Absolute Eigenvalue Filtering for Projected Newton}, + author = {Honglin Chen and Hsueh-Ti Derek Liu and David I.W. Levin and Changxi Zheng and Alec Jacobson}, + year = 2024, + booktitle = {{ACM} {SIGGRAPH} 2024 Conference Proceedings}, + note = {\url{https://www.cs.columbia.edu/cg/abs-psd/}} +} +@article{Fang2023AugmentedStickyInteractions, + title = {Augmented Incremental Potential Contact for Sticky Interactions}, + author = {Fang, Yu and Li, Minchen and Cao, Yadi and Li, Xuan and Wolper, Joshuah and Yang, Yin and Jiang, Chenfanfu}, + year = 2024, + journal = {IEEE Transactions on Visualization and Computer Graphics}, + volume = 30, + number = 8, + pages = {5596--5608}, + doi = {10.1109/TVCG.2023.3295656}, + keywords = {Adhesives;Friction;Computational modeling;Deformation;Solids;Force;Deformable models;Physically based animation;optimization time integration} +} +@article{Ando2024Cubic, + title = {A Cubic Barrier with Elasticity-Inclusive Dynamic Stiffness}, + author = {Ando, Ryoichi}, + year = 2024, + month = dec, + journal = {{ACM} Transactions on Graphics ({SIGGRAPH} Asia)}, + volume = 43, + number = 6, + pages = {1--13}, + note = {\url{https://github.com/st-tech/ppf-contact-solver}} +} +@article{Huang2024GIPC, + title = {{GIPC}: Fast and Stable Gauss-Newton Optimization of IPC Barrier Energy}, + author = {Huang, Kemeng and Chitalu, Floyd M. and Lin, Huancheng and Komura, Taku}, + year = 2024, + month = apr, + journal = {ACM Transactions on Graphics}, + volume = 43, + number = 2, + pages = {1--18}, + note = {\url{https://dl.acm.org/doi/10.1145/3643028}} +} diff --git a/docs/source/refs.bib b/docs/source/refs.bib deleted file mode 100644 index accb20bb5..000000000 --- a/docs/source/refs.bib +++ /dev/null @@ -1,88 +0,0 @@ -@article{Li2020IPC, - title = {Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics}, - author = {Minchen Li and Zachary Ferguson and Teseo Schneider and Timothy Langlois and Denis Zorin and Daniele Panozzo and Chenfanfu Jiang and Danny M. Kaufman}, - year = 2020, - journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, - volume = 39, - number = 4, - articleno = 49, - note = {\url{https://ipc-sim.github.io}} -} -@article{Li2021CIPC, - title = {Codimensional Incremental Potential Contact}, - author = {Minchen Li and Danny M. Kaufman and Chenfanfu Jiang}, - year = 2021, - journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, - volume = 40, - number = 4, - articleno = 170, - note = {\url{https://ipc-sim.github.io/C-IPC}} -} -@misc{Li2023Convergent, - title = {Convergent Incremental Potential Contact}, - author = {Minchen Li and Zachary Ferguson and Teseo Schneider and Timothy Langlois and Denis Zorin and Daniele Panozzo and Chenfanfu Jiang and Danny M. Kaufman}, - year = 2023, - eprint = {2307.15908}, - archiveprefix = {arXiv}, - primaryclass = {math.NA} -} -@article{Wang2021TightInclusion, - title = {A Large Scale Benchmark and an Inclusion-Based Algorithm for Continuous Collision Detection}, - author = {Bolun Wang and Zachary Ferguson and Teseo Schneider and Xin Jiang and Marco Attene and Daniele Panozzo}, - year = 2021, - month = oct, - journal = {{ACM} Transactions on Graphics}, - volume = 40, - number = 5, - articleno = 188, - numpages = 16, - note = {\url{https://continuous-collision-detection.github.io/tight_inclusion/}} -} -@misc{Belgrod2023Time, - title = {Time of Impact Dataset for Continuous Collision Detection and a Scalable Conservative Algorithm}, - author = {David Belgrod and Bolun Wang and Zachary Ferguson and Xin Zhao and Marco Attene and Daniele Panozzo and Teseo Schneider}, - year = 2023, - eprint = {2112.06300}, - archiveprefix = {arXiv}, - primaryclass = {cs.GR} -} -@inproceedings{Ferguson2023HighOrderIPC, - title = {High-Order Incremental Potential Contact for Elastodynamic Simulation on Curved Meshes}, - author = {Zachary Ferguson and Pranav Jain and Denis Zorin and Teseo Schneider and Daniele Panozzo}, - year = 2023, - booktitle = {{ACM} {SIGGRAPH} 2023 Conference Proceedings}, - location = {Los Angeles, CA, USA}, - publisher = {Association for Computing Machinery}, - address = {New York, NY, USA}, - series = {SIGGRAPH '23}, - numpages = 11, - note = {\url{https://zferg.us/research/high-order-ipc/}} -} -@article{Ferguson2021RigidIPC, - title = {Intersection-free Rigid Body Dynamics}, - author = {Zachary Ferguson and Minchen Li and Teseo Schneider and Francisca Gil-Ureta and Timothy Langlois and Chenfanfu Jiang and Denis Zorin and Danny M. Kaufman and Daniele Panozzo}, - year = 2021, - journal = {{ACM} Transactions on Graphics ({SIGGRAPH})}, - volume = 40, - number = 4, - articleno = 183, - note = {\url{https://ipc-sim.github.io/rigid-ipc}} -} -@inproceedings{Chen2024Stabler, - title = {Stabler Neo-Hookean Simulation: Absolute Eigenvalue Filtering for Projected Newton}, - author = {Honglin Chen and Hsueh-Ti Derek Liu and David I.W. Levin and Changxi Zheng and Alec Jacobson}, - year = 2024, - booktitle = {{ACM} {SIGGRAPH} 2024 Conference Proceedings}, - note = {\url{https://www.cs.columbia.edu/cg/abs-psd/}} -} -@article{Fang2023AugmentedStickyInteractions, - title = {Augmented Incremental Potential Contact for Sticky Interactions}, - author = {Fang, Yu and Li, Minchen and Cao, Yadi and Li, Xuan and Wolper, Joshuah and Yang, Yin and Jiang, Chenfanfu}, - journal = {IEEE Transactions on Visualization and Computer Graphics}, - year = {2024}, - volume = {30}, - number = {8}, - pages = {5596-5608}, - keywords = {Adhesives;Friction;Computational modeling;Deformation;Solids;Force;Deformable models;Physically based animation;optimization time integration}, - doi = {10.1109/TVCG.2023.3295656} -} diff --git a/python/src/barrier/adaptive_stiffness.cpp b/python/src/barrier/adaptive_stiffness.cpp index 37e1acec9..4d7319528 100644 --- a/python/src/barrier/adaptive_stiffness.cpp +++ b/python/src/barrier/adaptive_stiffness.cpp @@ -1,6 +1,8 @@ #include #include +#include +#include namespace py = pybind11; using namespace ipc; @@ -66,4 +68,80 @@ void define_adaptive_stiffness(py::module_& m) py::arg("max_barrier_stiffness"), py::arg("barrier_stiffness"), py::arg("bbox_diagonal"), py::arg("dhat_epsilon_scale") = 1e-9, py::arg("dmin") = 0); + + m.def( + "semi_implicit_stiffness", + static_cast&, + const VectorMax12d&, const VectorMax4d&, const MatrixMax12d&, + const double)>(&semi_implicit_stiffness), + R"ipc_Qu8mg5v7( + Compute the semi-implicit stiffness for a single collision. + + See [Ando 2024] for details. + + Parameters: + stencil: Collision stencil. + vertex_ids: Vertex indices of the collision. + vertices: Vertex positions. + mass: Vertex masses. + local_hess: Local hessian of the elasticity energy function. + dmin: Minimum distance between elements. + + Returns: + The semi-implicit stiffness. + )ipc_Qu8mg5v7", + py::arg("stencil"), py::arg("vertex_ids"), py::arg("vertices"), + py::arg("mass"), py::arg("local_hess"), py::arg("dmin")); + + m.def( + "semi_implicit_stiffness", + static_cast&, const double)>( + &semi_implicit_stiffness), + R"ipc_Qu8mg5v7( + Compute the semi-implicit stiffness's for all collisions. + + See [Ando 2024] for details. + + Parameters: + mesh: Collision mesh. + vertices: Vertex positions. + collisions: Normal collisions. + vertex_masses: Lumped vertex masses. + hess: Hessian of the elasticity energy function. + dmin: Minimum distance between elements. + + Returns: + The semi-implicit stiffness's. + )ipc_Qu8mg5v7", + py::arg("mesh"), py::arg("vertices"), py::arg("collisions"), + py::arg("vertex_masses"), py::arg("hess"), py::arg("dmin")); + + m.def( + "semi_implicit_stiffness", + static_cast&, + const double)>(&semi_implicit_stiffness), + R"ipc_Qu8mg5v7( + Compute the semi-implicit stiffness's for all collisions. + + See [Ando 2024] for details. + + Parameters: + mesh: Collision mesh. + vertices: Vertex positions. + collisions: Collisions candidates. + vertex_masses: Lumped vertex masses. + hess: Hessian of the elasticity energy function. + dmin: Minimum distance between elements. + + Returns: + The semi-implicit stiffness's. + )ipc_Qu8mg5v7", + py::arg("mesh"), py::arg("vertices"), py::arg("collisions"), + py::arg("vertex_masses"), py::arg("hess"), py::arg("dmin")); } diff --git a/src/ipc/barrier/adaptive_stiffness.cpp b/src/ipc/barrier/adaptive_stiffness.cpp index 38f53223f..60c7cb0de 100644 --- a/src/ipc/barrier/adaptive_stiffness.cpp +++ b/src/ipc/barrier/adaptive_stiffness.cpp @@ -90,12 +90,12 @@ double semi_implicit_stiffness( const std::array& vertex_ids, const VectorMax12d& vertices, const VectorMax4d& mass, - const Eigen::SparseMatrix& hess, + const MatrixMax12d& local_hess, const double dmin) { - unsigned N = stencil.num_vertices(); + const unsigned N = stencil.num_vertices(); assert(vertices.size() % N == 0); - unsigned dim = vertices.size() / N; + const unsigned dim = stencil.dim(vertices.size()); const VectorMax4d value = stencil.compute_coefficients(vertices); @@ -105,27 +105,13 @@ double semi_implicit_stiffness( normal += value[i] * vertices.segment(dim * i, dim); } - MatrixMax12d local_hess = MatrixMax12d::Zero(dim * N, dim * N); - - // H - for (unsigned i = 0; i < N; ++i) { - for (unsigned j = 0; j < N; ++j) { - for (unsigned k = 0; k < dim; ++k) { - for (unsigned l = 0; l < dim; ++l) { - local_hess(dim * i + k, dim * j + l) = hess.coeff( - dim * vertex_ids[i] + k, dim * vertex_ids[j] + l); - } - } - } - } - - // mᵢ / d² + // d² const double distance = normal.norm() - dmin; const double distance_sqr = distance * distance; - for (unsigned i = 0; i < N; ++i) { - local_hess.diagonal().segment(dim * i, dim).array() += - mass[i] / distance_sqr; - } + + // average mass: mᵢ = cᵀMc / ‖c‖² + const double avg_mass = + value.dot(mass.asDiagonal() * value) / value.squaredNorm(); VectorMax12d w = VectorMax12d::Zero(dim * N); for (unsigned i = 0; i < N; ++i) { @@ -133,7 +119,7 @@ double semi_implicit_stiffness( } w.normalize(); - return (local_hess * w).dot(w); + return avg_mass / distance_sqr + w.dot(local_hess * w); } template @@ -145,23 +131,51 @@ Eigen::VectorXd semi_implicit_stiffness( const Eigen::SparseMatrix& hess, const double dmin) { + const int dim = mesh.dim(); + assert(vertices.cols() == dim); // Vertex positions must be 3D + assert(hess.rows() == hess.cols()); // Hessian must be square + // Hess and vertex_masses must have the same number of rows + assert(hess.rows() == vertex_masses.size() * dim); + // Hess can be either for the reduced or full mesh + assert(hess.rows() == mesh.ndof() || hess.rows() == mesh.full_ndof()); + Eigen::VectorXd stiffnesses(collisions.size()); - for (size_t i = 0; i < collisions.size(); i++) { - const CollisionStencil& collision = collisions[i]; + for (size_t ci = 0; ci < collisions.size(); ci++) { + const CollisionStencil& collision = collisions[ci]; + const unsigned N = collision.num_vertices(); - const std::array vertex_ids = - collision.vertex_ids(mesh.edges(), mesh.faces()); const VectorMax12d positions = collision.dof(vertices, mesh.edges(), mesh.faces()); - VectorMax4d mass(collision.num_vertices()); - for (unsigned j = 0; j < collision.num_vertices(); j++) { - mass[j] = vertex_masses[vertex_ids[j]]; + std::array vertex_ids = + collision.vertex_ids(mesh.edges(), mesh.faces()); + if (hess.rows() == mesh.full_ndof()) { + for (int i = 0; i < N; i++) { + vertex_ids[i] = mesh.to_full_vertex_id(vertex_ids[i]); + } + } + + VectorMax4d local_mass(collision.num_vertices()); + for (unsigned i = 0; i < collision.num_vertices(); i++) { + local_mass[i] = vertex_masses[vertex_ids[i]]; + } + + MatrixMax12d local_hess = MatrixMax12d::Zero(dim * N, dim * N); + for (unsigned i = 0; i < N; ++i) { + for (unsigned j = 0; j < N; ++j) { + for (unsigned k = 0; k < dim; ++k) { + for (unsigned l = 0; l < dim; ++l) { + // NOTE: Assumes DOF are flattened in row-major order + local_hess(dim * i + k, dim * j + l) = hess.coeff( + dim * vertex_ids[i] + k, dim * vertex_ids[j] + l); + } + } + } } - stiffnesses[i] = semi_implicit_stiffness( - collision, vertex_ids, positions, mass, hess, dmin); + stiffnesses[ci] = semi_implicit_stiffness( + collision, vertex_ids, positions, local_mass, local_hess, dmin); } return stiffnesses; diff --git a/src/ipc/barrier/adaptive_stiffness.hpp b/src/ipc/barrier/adaptive_stiffness.hpp index dfed422d4..51a155903 100644 --- a/src/ipc/barrier/adaptive_stiffness.hpp +++ b/src/ipc/barrier/adaptive_stiffness.hpp @@ -51,11 +51,12 @@ double update_barrier_stiffness( const double dmin = 0); /// @brief Compute the semi-implicit stiffness for a single collision. +/// See [Ando 2024] for details. /// @param stencil Collision stencil. /// @param vertex_ids Vertex indices of the collision. /// @param vertices Vertex positions. /// @param mass Vertex masses. -/// @param hess Hessian of the energy function. +/// @param local_hess Local hessian of the elasticity energy function. /// @param dmin Minimum distance between elements. /// @return The semi-implicit stiffness. double semi_implicit_stiffness( @@ -63,17 +64,18 @@ double semi_implicit_stiffness( const std::array& vertex_ids, const VectorMax12d& vertices, const VectorMax4d& mass, - const Eigen::SparseMatrix& hess, - double dmin); + const MatrixMax12d& local_hess, + const double dmin); -/// @brief Compute the semi-implicit stiffnesses for all collisions. +/// @brief Compute the semi-implicit stiffness's for all collisions. +/// See [Ando 2024] for details. /// @param mesh Collision mesh. /// @param vertices Vertex positions. -/// @param collisions Normal collisions. +/// @param collisions Normal collisions or collision candidates. /// @param vertex_masses Lumped vertex masses. /// @param hess Hessian of the elasticity energy function. /// @param dmin Minimum distance between elements. -/// @return The semi-implicit stiffnesses. +/// @return The semi-implicit stiffness's. template Eigen::VectorXd semi_implicit_stiffness( const CollisionMesh& mesh, @@ -81,6 +83,6 @@ Eigen::VectorXd semi_implicit_stiffness( const StencilsT& collisions, const Eigen::VectorXd& vertex_masses, const Eigen::SparseMatrix& hess, - double dmin); + const double dmin); } // namespace ipc From ce54d5e8f3a198bd35ee61a0ab359b034f9456d2 Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 21:41:56 -0500 Subject: [PATCH 07/10] Test all coeffs cases --- .../tests/candidates/test_coefficients.cpp | 108 +++++++++++++----- 1 file changed, 78 insertions(+), 30 deletions(-) diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp index b134aded9..972836713 100644 --- a/tests/src/tests/candidates/test_coefficients.cpp +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -7,47 +7,52 @@ using namespace ipc; -TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") +TEST_CASE("Vertex-vertex collision stencil coeffs.", "[vv][stencil][coeffs]") { - Eigen::MatrixXd V(4, 3); - V << 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1; + Eigen::MatrixXd V(2, 3); + V << -1, 0, 0, /**/ 1, 0, 0; - SECTION("default") { } - SECTION("random") { V.row(3).setRandom(); } - SECTION("off triangle 0") { V.row(3) << 2, 0.5, 1; } - SECTION("off triangle 1") { V.row(3) << 2, 2, 1; } - SECTION("off triangle 2") { V.row(3) << 0.5, 2, 1; } - SECTION("off triangle 3") { V.row(3) << -1, 2, 1; } - SECTION("off triangle 4") { V.row(3) << -1, 0.5, 1; } - SECTION("off triangle 5") { V.row(3) << -1, -1, 1; } - SECTION("off triangle 6") { V.row(3) << 0.5, -1, 1; } - SECTION("off triangle 7") { V.row(3) << 2, -0.5, 1; } - SECTION("failure case 1") { V.row(3) << 0.680375, -0.211234, 0.566198; } - SECTION("failure case 1") { V.row(3) << -0.997497, 0.127171, -0.613392; } + Eigen::MatrixXi E, F; - Eigen::MatrixXi F(1, 3); - F << 0, 1, 2; + VertexVertexCandidate vv(0, 1); - Eigen::MatrixXi E; - igl::edges(F, E); + VectorMax4d coeffs = vv.compute_coefficients(V, E, F); - FaceVertexCandidate fv(0, 3); + std::array vertices = vv.vertices(V, E, F); + Eigen::Vector3d n = Eigen::Vector3d::Zero(); + for (int i = 0; i < vv.num_vertices(); i++) { + n += coeffs[i] * vertices[i]; + } - VectorMax4d coeffs = fv.compute_coefficients(V, E, F); - CAPTURE( - point_triangle_distance_type(V.row(3), V.row(0), V.row(1), V.row(2))); - CAPTURE(V.row(3), coeffs.transpose()); + CHECK(n.squaredNorm() == Catch::Approx(vv.compute_distance(V, E, F))); + CHECK(coeffs[0] == 1); + CHECK(coeffs[1] == -1); +} - CHECK(-coeffs[1] == Catch::Approx(1 + coeffs[2] + coeffs[3])); +TEST_CASE("Edge-vertex collision stencil coeffs.", "[ev][stencil][coeffs]") +{ + Eigen::MatrixXd V(3, 3); + V << -1, 0, 0, /**/ 1, 0, 0, /**/ 0, 1, 0; - std::array vertices = fv.vertices(V, E, F); + SECTION("default") { } + SECTION("random") { V.row(2).setRandom(); } + SECTION("e0") { V(0, 2) = -2; } + SECTION("e1") { V(0, 2) = 2; } + + Eigen::MatrixXi E(1, 2), F; + E << 0, 1; + + EdgeVertexCandidate ev(0, 2); + + VectorMax4d coeffs = ev.compute_coefficients(V, E, F); + + std::array vertices = ev.vertices(V, E, F); Eigen::Vector3d n = Eigen::Vector3d::Zero(); - for (int i = 0; i < fv.num_vertices(); i++) { + for (int i = 0; i < ev.num_vertices(); i++) { n += coeffs[i] * vertices[i]; } - CAPTURE(n.transpose()); - CHECK(n.squaredNorm() == Catch::Approx(fv.compute_distance(V, E, F))); + CHECK(n.squaredNorm() == Catch::Approx(ev.compute_distance(V, E, F))); } TEST_CASE("Edge-edge collision stencil coeffs.", "[ee][stencil][coeffs]") @@ -85,4 +90,47 @@ TEST_CASE("Edge-edge collision stencil coeffs.", "[ee][stencil][coeffs]") } CHECK(n.squaredNorm() == Catch::Approx(ee.compute_distance(V, E, F))); -} \ No newline at end of file +} + +TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") +{ + Eigen::MatrixXd V(4, 3); + V << 0, 0, 0, /**/ 1, 0, 0, /**/ 0, 1, 0, /**/ 0.333, 0.333, 1; + + SECTION("default") { } + SECTION("random") { V.row(3).setRandom(); } + SECTION("off triangle 0") { V.row(3) << 2, 0.5, 1; } + SECTION("off triangle 1") { V.row(3) << 2, 2, 1; } + SECTION("off triangle 2") { V.row(3) << 0.5, 2, 1; } + SECTION("off triangle 3") { V.row(3) << -1, 2, 1; } + SECTION("off triangle 4") { V.row(3) << -1, 0.5, 1; } + SECTION("off triangle 5") { V.row(3) << -1, -1, 1; } + SECTION("off triangle 6") { V.row(3) << 0.5, -1, 1; } + SECTION("off triangle 7") { V.row(3) << 2, -0.5, 1; } + SECTION("failure case 1") { V.row(3) << 0.680375, -0.211234, 0.566198; } + SECTION("failure case 1") { V.row(3) << -0.997497, 0.127171, -0.613392; } + + Eigen::MatrixXi F(1, 3); + F << 0, 1, 2; + + Eigen::MatrixXi E; + igl::edges(F, E); + + FaceVertexCandidate fv(0, 3); + + VectorMax4d coeffs = fv.compute_coefficients(V, E, F); + CAPTURE( + point_triangle_distance_type(V.row(3), V.row(0), V.row(1), V.row(2))); + CAPTURE(V.row(3), coeffs.transpose()); + + CHECK(-coeffs[1] == Catch::Approx(1 + coeffs[2] + coeffs[3])); + + std::array vertices = fv.vertices(V, E, F); + Eigen::Vector3d n = Eigen::Vector3d::Zero(); + for (int i = 0; i < fv.num_vertices(); i++) { + n += coeffs[i] * vertices[i]; + } + CAPTURE(n.transpose()); + + CHECK(n.squaredNorm() == Catch::Approx(fv.compute_distance(V, E, F))); +} From 57b20909796aac1853b00f91ae93d739a969acd3 Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 21:48:08 -0500 Subject: [PATCH 08/10] Test plane_vertex coeffs --- src/ipc/collisions/normal/plane_vertex.hpp | 5 +++++ .../tests/candidates/test_coefficients.cpp | 20 ++++++++++++++++++- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/src/ipc/collisions/normal/plane_vertex.hpp b/src/ipc/collisions/normal/plane_vertex.hpp index 1d9f46c94..713689bd5 100644 --- a/src/ipc/collisions/normal/plane_vertex.hpp +++ b/src/ipc/collisions/normal/plane_vertex.hpp @@ -21,6 +21,11 @@ class PlaneVertexNormalCollision : public NormalCollision { return { { vertex_id, -1, -1, -1 } }; } + using NormalCollision::compute_coefficients; + using NormalCollision::compute_distance; + using NormalCollision::compute_distance_gradient; + using NormalCollision::compute_distance_hessian; + /// @brief Compute the distance between the point and plane. /// @param point Point's position. /// @return Distance of the stencil. diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp index 972836713..b6a53eb6a 100644 --- a/tests/src/tests/candidates/test_coefficients.cpp +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -2,6 +2,7 @@ #include #include +#include #include @@ -92,7 +93,7 @@ TEST_CASE("Edge-edge collision stencil coeffs.", "[ee][stencil][coeffs]") CHECK(n.squaredNorm() == Catch::Approx(ee.compute_distance(V, E, F))); } -TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") +TEST_CASE("Face-vertex collision stencil coeffs.", "[fv][stencil][coeffs]") { Eigen::MatrixXd V(4, 3); V << 0, 0, 0, /**/ 1, 0, 0, /**/ 0, 1, 0, /**/ 0.333, 0.333, 1; @@ -134,3 +135,20 @@ TEST_CASE("Face-Vertex collision stencil coeffs.", "[fv][stencil][coeffs]") CHECK(n.squaredNorm() == Catch::Approx(fv.compute_distance(V, E, F))); } + +TEST_CASE("Plane-vertex collision stencil coeffs.", "[pv][stencil][coeffs]") +{ + Eigen::MatrixXd V(1, 3); + V << 0, 1, 0; + + Eigen::Vector3d n(0, 1, 0), o(0, 0, 0); + + Eigen::MatrixXi E, F; + + PlaneVertexNormalCollision pv(o, n, 0); + + VectorMax4d coeffs = pv.compute_coefficients(V, E, F); + + CHECK(coeffs.size() == 1); + CHECK(coeffs[0] == 1); +} \ No newline at end of file From df3db63ce6e46795a04c02505f4a90c74b8f93ed Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 22:02:49 -0500 Subject: [PATCH 09/10] Test semi-implicit stiffness --- .../tests/barrier/test_adaptive_stiffness.cpp | 38 ++++++++++++++++++- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/tests/src/tests/barrier/test_adaptive_stiffness.cpp b/tests/src/tests/barrier/test_adaptive_stiffness.cpp index d0e4f8109..e6c115c30 100644 --- a/tests/src/tests/barrier/test_adaptive_stiffness.cpp +++ b/tests/src/tests/barrier/test_adaptive_stiffness.cpp @@ -1,11 +1,15 @@ #include #include +#include #include +#include + +#include using namespace ipc; -TEST_CASE("Initial barrier stiffness", "[adaptive_stiffness]") +TEST_CASE("Initial barrier stiffness", "[stiffness][adaptive_stiffness]") { const double bbox_diagonal = 1.0; const ClampedLogBarrier barrier; @@ -26,7 +30,7 @@ TEST_CASE("Initial barrier stiffness", "[adaptive_stiffness]") CHECK(max_barrier_stiffness == Catch::Approx(100 * expected_kappa)); } -TEST_CASE("Update barrier stiffness", "[adaptive_stiffness]") +TEST_CASE("Update barrier stiffness", "[stiffness][adaptive_stiffness]") { double prev_min_distance = 1e-10 * 1e-10; double min_distance = 1e-11 * 1e-11; @@ -39,4 +43,34 @@ TEST_CASE("Update barrier stiffness", "[adaptive_stiffness]") barrier_stiffness, bbox_diagonal); REQUIRE(kappa == Catch::Approx(2.0)); +} + +TEST_CASE("Semi-implicit stiffness", "[stiffness]") +{ + Eigen::MatrixXd vertices(4, 3); + Eigen::MatrixXi edges, faces(1, 3); + + const double d = GENERATE(range(0.1, 1.0, 0.1)); + vertices << 0, 0, 0, /**/ 1, 0, 0, /**/ 0, 1, 0, /**/ 0.333, 0.333, d; + faces << 0, 1, 2; + igl::edges(faces, edges); + + const CollisionMesh mesh(vertices, edges, faces); + + NormalCollisions collisions; + collisions.fv_collisions.emplace_back(0, 3); + + const double m = GENERATE(range(1.0, 10.0, 1.0)); + const Eigen::VectorXd vertex_masses = + Eigen::VectorXd::Constant(vertices.rows(), m); + + const double dmin = 0; + + Eigen::SparseMatrix hess(vertices.size(), vertices.size()); // zero + + const Eigen::VectorXd kappa = semi_implicit_stiffness( + mesh, vertices, collisions, vertex_masses, hess, dmin); + + REQUIRE(kappa.size() == 1); + REQUIRE(kappa(0) == Catch::Approx(m / (d * d))); } \ No newline at end of file From f54a6275cf628e32ad3de5f10135e6f97553b3d0 Mon Sep 17 00:00:00 2001 From: zachferguson Date: Fri, 7 Feb 2025 22:05:45 -0500 Subject: [PATCH 10/10] Fix index assignment in edge-vertex collision stencil coefficient tests --- tests/src/tests/candidates/test_coefficients.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/src/tests/candidates/test_coefficients.cpp b/tests/src/tests/candidates/test_coefficients.cpp index b6a53eb6a..4e8dd9125 100644 --- a/tests/src/tests/candidates/test_coefficients.cpp +++ b/tests/src/tests/candidates/test_coefficients.cpp @@ -37,8 +37,8 @@ TEST_CASE("Edge-vertex collision stencil coeffs.", "[ev][stencil][coeffs]") SECTION("default") { } SECTION("random") { V.row(2).setRandom(); } - SECTION("e0") { V(0, 2) = -2; } - SECTION("e1") { V(0, 2) = 2; } + SECTION("e0") { V(2, 0) = -2; } + SECTION("e1") { V(2, 0) = 2; } Eigen::MatrixXi E(1, 2), F; E << 0, 1;