Update Nonlinear CCD Docs

Author: zfergus

docs/source/cpp-api/interval.rst

@@ -0,0 +1,31 @@
+Interval Arithmetic
+===================
+
+For interval arithmetic, we use the `Filib <https://github.com/zfergus/filib>`_ library.
+
+.. doxygenclass:: filib::Interval
+
+Helper Functions
+----------------
+
+.. doxygenfunction:: ipc::norm
+.. doxygenfunction:: ipc::squared_norm
+
+Type Definitions
+----------------
+
+.. doxygentypedef:: ipc::Vector2I
+.. doxygentypedef:: ipc::Vector3I
+.. doxygentypedef:: ipc::VectorXI
+
+.. doxygentypedef:: ipc::RowVector2I
+.. doxygentypedef:: ipc::RowVector3I
+.. doxygentypedef:: ipc::RowVectorXI
+
+.. doxygentypedef:: ipc::Matrix2I
+.. doxygentypedef:: ipc::Matrix3I
+.. doxygentypedef:: ipc::MatrixXI
+
+.. doxygentypedef:: ipc::VectorMax3I
+.. doxygentypedef:: ipc::RowVectorMax3I
+.. doxygentypedef:: ipc::MatrixMax3I

docs/source/cpp.rst

@@ -31,6 +31,8 @@ C++
     :start-after: <!--- BEGIN C++ README 2 --->
     :end-before: <!--- FILIB DEPENDENCY NOTE --->
 
+.. _filib_dependency_note:
+
 .. warning::
     ``filib`` is licensed under `LGPL-2.1 <https://github.com/zfergus/filib/blob/main/LICENSE>`_ and as such it is required to be dynamically linked. Doing so automatically is a challenge, so by default we use static linkage. Enabling dynaic linkage requires copying the ``.so``/``.dylib``/``.dll`` file to the binary directory or system path. To enable this, set the CMake option :cmake:`FILIB_BUILD_SHARED_LIBS` to :cmake:`ON` and add this CMake code to copy the shared libaray object to the binary directory:
 

docs/source/index.rst

@@ -37,6 +37,7 @@
     cpp-api/barrier.rst
     cpp-api/adhesion.rst
     cpp-api/intersections.rst
+    cpp-api/interval.rst
     cpp-api/utils.rst
 
 .. toctree::
@@ -57,6 +58,7 @@
     python-api/barrier.rst
     python-api/adhesion.rst
     python-api/intersections.rst
+    python-api/interval.rst
     python-api/utils.rst
 
 .. toctree::

docs/source/python-api/interval.rst

@@ -0,0 +1,8 @@
+Interval Arithmetic
+===================
+
+.. automodule:: ipctk.filib
+    :members:
+    :undoc-members:
+    :show-inheritance:
+    :inherited-members:

docs/source/tutorial/adhesion.rst

@@ -3,7 +3,7 @@ Adhesion
 
 Adhesion simulates the sticking interaction between surfaces, such as glue bonding or material contact.
 
-We provide functionality to compute normal adhesion (perpendicular forces) and tangential adhesion (parallel forces) using the model of :cite:`Fang2023AugmentedStickyInteractions`.
+We provide functionality to compute normal adhesion (perpendicular forces) and tangential adhesion (parallel forces) using the model of :cite:t:`Fang2023AugmentedStickyInteractions`.
 
 Normal Adhesion
 ---------------

docs/source/tutorial/nonlinear_ccd.rst

@@ -1,7 +1,9 @@
 Nonlinear CCD
 =============
 
-We can also perform CCD of nonlinear trajectories (of linear geometry) using the method of :cite:t:`Ferguson2021RigidIPC`. While :cite:t:`Ferguson2021RigidIPC` introduce their method in the context of rigid bodies, it can be applied to any nonlinear trajectory of linear geometry. The method works by transforming the nonlinear trajectories into (adaptive) piecewise linear trajectories with an envelope/minimum separation around each piece, enclosing the nonlinear trajectory. The method then performs CCD on the piecewise linear trajectories to find the earliest time of impact.
+We also implement CCD of nonlinear trajectories (of linear geometry) using the method of :cite:t:`Ferguson2021RigidIPC`. While :cite:t:`Ferguson2021RigidIPC` introduce their method in the context of rigid bodies, it can be applied to any nonlinear trajectory of linear geometry.
+
+The method works by transforming the nonlinear trajectories into (adaptive) piecewise linear trajectories with an envelope/minimum separation around each piece, enclosing the nonlinear trajectory. The method then performs CCD on the piecewise linear trajectories to find the earliest time of impact.
 
 We provide the following functions to perform nonlinear CCD:
 
@@ -14,7 +16,7 @@ We provide the following functions to perform nonlinear CCD:
         * :cpp:func:`ipc::edge_edge_nonlinear_ccd`, and
         * :cpp:func:`ipc::point_triangle_nonlinear_ccd`.
 
-        Each of these functions take as input a :cpp:func:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry.
+        Each of these functions take as input a :cpp:class:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry.
 
     .. md-tab-item:: Python
 
@@ -23,7 +25,7 @@ We provide the following functions to perform nonlinear CCD:
         * :py:func:`ipctk.edge_edge_nonlinear_ccd`, and
         * :py:func:`ipctk.point_triangle_nonlinear_ccd`.
 
-        Each of these functions take as input a :py:func:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry.
+        Each of these functions take as input a :py:class:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry.
 
 For example, the following code defines a rigid trajectory in 2D in order to perform nonlinear CCD between a point and edge:
 
@@ -47,7 +49,7 @@ For example, the following code defines a rigid trajectory in 2D in order to per
 Defining the Trajectory
 -----------------------
 
-Let's dive deeper by breaking down the implementation of Rigid2DTrajectory. The first function we need to implement is the call operator:
+Let's dive deeper by breaking down the implementation of ``Rigid2DTrajectory``. The first function we need to implement is the call operator:
 
 .. md-tab-set::
 
@@ -63,7 +65,7 @@ Let's dive deeper by breaking down the implementation of Rigid2DTrajectory. The
 
         .. literalinclude:: ../../../python/tests/test_ccd.py
             :language: python
-            :lines: 91-95
+            :lines: 90-94
             :dedent: 8
 
 This function computes the position of the point at a time :math:`t \in [0, 1]`. This defines the trajectory of the point. In this case, we have a rigid body with a center of mass (COM) at the origin. The trajectory of the point is given by:
@@ -92,7 +94,7 @@ The second function we need to implement is ``max_distance_from_linear``.
 
         .. literalinclude:: ../../../python/tests/test_ccd.py
             :language: python
-            :lines: 97-103
+            :lines: 96-102
             :dedent: 8
 
 This function computes the maximum distance over a time interval :math:`[t_0, t_1]` between the nonlinear trajectory and a line segment from :math:`x(t_0)` to :math:`x(t_1)`. Mathematically this function computes
@@ -135,4 +137,96 @@ Last, we use the ``Rigid2DTrajectory`` to perform nonlinear CCD between a point
             :dedent: 4
 
 .. note::
-    We adjust the ``conservative_rescaling`` parameter to get a more accurate time of impact (TOI), but in practice, this is not needed as a more conservative estimate of the TOI is sufficient to avoid penetrations.
+    We adjust the ``conservative_rescaling`` parameter to get a more accurate time of impact (TOI), but in practice, this is not needed as a more conservative estimate of the TOI is sufficient to avoid penetrations.
+
+Interval-Based Nonlinear CCD
+----------------------------
+
+.. warning::
+    The following section requires enabling the ``filib`` interval arithmetic library. ``filib`` is licensed under L-GPL 2.1, so special care must be taken when using it. See the `filib dependency note <../../cpp.html#filib-dependency-note>`_ for more information.
+
+If an analytic expression for the ``max_distance_from_linear`` function is not available, we can use interval arithmetic to compute a conservative envelope.
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        The :cpp:class:`ipc::IntervalNonlinearTrajectory` class does this for us and all we need to provide is a function to compute the point's position over a time interval.
+
+    .. md-tab-item:: Python
+
+        The :py:class:`ipctk.IntervalNonlinearTrajectory` class does this for us and all we need to provide is a function to compute the point's position over a time interval.
+
+Conservative Envelope
+~~~~~~~~~~~~~~~~~~~~~
+
+Our implementation of the ``max_distance_from_linear`` function is as follows.
+
+Start by defining a linear interpolation function:
+
+.. math::
+    \operatorname{lerp}(a, b, t) := (b - a) t + a,
+
+which interpolates between two points :math:`a` and :math:`b` at time :math:`t`.
+
+The exact envelope from above is bounded by a interval approximation:
+
+.. math::
+    \begin{align}
+    &\max_{t \in [0, 1]} \| p(\operatorname{lerp}(t_0, t_1, t)) - \operatorname{lerp}(p(t_0), p(t_1), t) \|_2\\
+    &\leq \sup(\| p_{\Box}([t_0, t_1]) - ((p(t_1) - p(t_0)) \cdot [0, 1] + p(t_0)) \|_2).
+    \end{align}
+
+Therefore, we can compute the a conservative estimate of the envelope by implementing :math:`p_{\Box}([t_0, t_1])` and :math:`p(t)`.
+
+.. note::
+    In practice we subdivide the interval into smaller intervals and compute the envelope over each subinterval. This is done to create a more accurate estimate.
+
+Interval Arithmetic
+~~~~~~~~~~~~~~~~~~~
+
+`Interval arithmetic <https://en.wikipedia.org/wiki/Interval_arithmetic>`_ is a method to compute bounds on the range of a function over an interval. We can construct a vector of intervals to represent the position of the point over a time interval.
+
+The following code snippet shows an example of how to use interval arithmetic to compute the position of a point over a time interval:
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: cpp
+
+            #include <ipc/utils/interval.hpp>
+
+            using namespace ipc;
+
+            Vector2I position(
+                const VectorMax3d& center,
+                const VectorMax3d& point,
+                const double omega,
+                const Interval& t)
+            {
+                // 2×2 matrix of intervals representing the rotation matrix
+                Matrix2I R;
+                R << cos(omega * t), -sin(omega * t),
+                     sin(omega * t),  cos(omega * t);
+                return R * (point - center) + center;
+            }
+
+        The full documentation for the ``Interval`` class can be found `here <../../cpp-api/interval.html>`_.
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            import numpy as np
+
+            from ipctk.filib import Interval, sin, cos
+
+            def position(center, point, omega, t : Interval):
+                R = np.array([
+                    [cos(omega * t), -sin(omega * t)],
+                    [sin(omega * t),  cos(omega * t)]
+                ])
+                return R @ (point - center) + center
+
+        The full documentation for the filib python bindings can be found `here <../python-api/interval.html>`_.

python/src/utils/interval.cpp

@@ -36,7 +36,9 @@ using namespace filib;
 void define_interval(py::module_& m)
 {
 #ifdef IPC_TOOLKIT_WITH_FILIB
-    auto m_filib = m.def_submodule("filib", "Fast Interval Library");
+    auto m_filib = m.def_submodule(
+        "filib",
+        "Python bindings for the `Fast Interval Library (Filib) <https://github.com/zfergus/filib>`_ by Werner Hofschuster and Walter Kraemer.");
 
     py::class_<Interval>(m_filib, "Interval")
         .def(py::init())

src/ipc/candidates/collision_stencil.hpp

@@ -39,7 +39,7 @@ class CollisionStencil {
     /// @return The vertex positions of the collision stencil. Size is always 4, but elements i > num_vertices() are NaN.
     template <typename T>
     std::array<VectorMax3<T>, 4> vertices(
-        const MatrixX<T>& vertices,
+        const Eigen::MatrixX<T>& vertices,
         const Eigen::MatrixXi& edges,
         const Eigen::MatrixXi& faces) const
     {
@@ -67,7 +67,7 @@ class CollisionStencil {
     /// @return This stencil's DOF.
     template <typename T>
     VectorMax12<T>
-    dof(const MatrixX<T>& X,
+    dof(const Eigen::MatrixX<T>& X,
         const Eigen::MatrixXi& edges,
         const Eigen::MatrixXi& faces) const
     {

src/ipc/collisions/normal/normal_collisions.cpp

@@ -136,6 +136,8 @@ namespace {
 
         return ev_candidates;
     }
+
+    inline double sqr(double x) { return x * x; }
 } // namespace
 
 void NormalCollisions::build(
@@ -147,7 +149,7 @@ void NormalCollisions::build(
 {
     assert(vertices.rows() == mesh.num_vertices());
 
-    double inflation_radius = (dhat + dmin) / 2;
+    const double inflation_radius = 0.5 * (dhat + dmin);
 
     Candidates candidates;
     candidates.build(mesh, vertices, inflation_radius, broad_phase_method);
@@ -168,8 +170,7 @@ void NormalCollisions::build(
 
     // Cull the candidates by measuring the distance and dropping those that are
     // greater than dhat.
-    const double offset_sqr = (dmin + dhat) * (dmin + dhat);
-    auto is_active = [&](double distance_sqr) {
+    auto is_active = [offset_sqr = sqr(dmin + dhat)](double distance_sqr) {
         return distance_sqr < offset_sqr;
     };
 

src/ipc/tangent/closest_point.cpp

@@ -40,7 +40,7 @@ VectorMax9d point_edge_closest_point_jacobian(
 // ============================================================================
 // Edge - Edge
 
-Vector2<double> edge_edge_closest_point(
+Eigen::Vector2d edge_edge_closest_point(
     const Eigen::Ref<const Eigen::Vector3d>& ea0,
     const Eigen::Ref<const Eigen::Vector3d>& ea1,
     const Eigen::Ref<const Eigen::Vector3d>& eb0,
@@ -55,7 +55,7 @@ Vector2<double> edge_edge_closest_point(
     coefMtr(0, 1) = coefMtr(1, 0) = -eb.dot(ea);
     coefMtr(1, 1) = eb.squaredNorm();
 
-    Vector2<double> rhs;
+    Eigen::Vector2d rhs;
     rhs[0] = -eb_to_ea.dot(ea);
     rhs[1] = eb_to_ea.dot(eb);
 
@@ -82,15 +82,15 @@ Eigen::Matrix<double, 2, 12> edge_edge_closest_point_jacobian(
 // ============================================================================
 // Point - Triangle
 
-Vector2<double> point_triangle_closest_point(
+Eigen::Vector2d point_triangle_closest_point(
     const Eigen::Ref<const Eigen::Vector3d>& p,
     const Eigen::Ref<const Eigen::Vector3d>& t0,
     const Eigen::Ref<const Eigen::Vector3d>& t1,
     const Eigen::Ref<const Eigen::Vector3d>& t2)
 {
     Eigen::Matrix<double, 2, 3> basis;
-    basis.row(0) = RowVector3<double>(t1 - t0); // edge 0
-    basis.row(1) = RowVector3<double>(t2 - t0); // edge 1
+    basis.row(0) = Eigen::RowVector3d(t1 - t0); // edge 0
+    basis.row(1) = Eigen::RowVector3d(t2 - t0); // edge 1
     const Eigen::Matrix2d A = basis * basis.transpose();
     const Eigen::Vector2d b = basis * (p - t0);
     const Eigen::Vector2d x = A.ldlt().solve(b);

src/ipc/utils/eigen_ext.hpp

@@ -8,31 +8,20 @@
 namespace ipc {
 
 // Boolean scalar
-typedef Eigen::Array<bool, Eigen::Dynamic, 1> ArrayXb;
-typedef Eigen::Matrix<bool, Eigen::Dynamic, 1> VectorXb;
-typedef Eigen::Matrix<bool, 3, 1> Vector3b;
-typedef Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic> MatrixXb;
+using ArrayXb = Eigen::Array<bool, Eigen::Dynamic, 1>;
+using VectorXb = Eigen::Matrix<bool, Eigen::Dynamic, 1>;
+using Vector3b = Eigen::Matrix<bool, 3, 1>;
+using MatrixXb = Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>;
 
 template <typename T, int dim, int max_dim = dim>
 using Vector = Eigen::Matrix<T, dim, 1, Eigen::ColMajor, max_dim, 1>;
 template <typename T, int dim, int max_dim = dim>
 using RowVector = Eigen::Matrix<T, 1, dim, Eigen::RowMajor, 1, max_dim>;
-template <typename T> using Vector1 = Vector<T, 1>;
-template <typename T> using Vector2 = Vector<T, 2>;
-template <typename T> using Vector3 = Vector<T, 3>;
-template <typename T> using VectorX = Vector<T, Eigen::Dynamic>;
-template <typename T> using RowVector2 = RowVector<T, 2>;
-template <typename T> using RowVector3 = RowVector<T, 3>;
-template <typename T> using RowVectorX = RowVector<T, Eigen::Dynamic>;
-template <typename T> using Matrix2 = Eigen::Matrix<T, 2, 2>;
-template <typename T> using Matrix3 = Eigen::Matrix<T, 3, 3>;
-template <typename T>
-using MatrixX = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
 
-using Vector1d = Vector1<double>;
-using Vector6d = Vector<double, 6>;
-using Vector9d = Vector<double, 9>;
-using Vector12d = Vector<double, 12>;
+using Vector1d = Eigen::Vector<double, 1>;
+using Vector6d = Eigen::Vector<double, 6>;
+using Vector9d = Eigen::Vector<double, 9>;
+using Vector12d = Eigen::Vector<double, 12>;
 using Matrix6d = Eigen::Matrix<double, 6, 6>;
 using Matrix9d = Eigen::Matrix<double, 9, 9>;
 using Matrix12d = Eigen::Matrix<double, 12, 12>;