Move NormalizedClampedLogBarrier into src (#143)

Author: zfergus

docs/source/cpp-api/barrier.rst

@@ -27,4 +27,10 @@ Clamped Log Barrier
 ~~~~~~~~~~~~~~~~~~~
 
 .. doxygenclass:: ipc::ClampedLogBarrier
+    :allow-dot-graphs:
+
+Normalized Clamped Log Barrier
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. doxygenclass:: ipc::NormalizedClampedLogBarrier
     :allow-dot-graphs:

docs/source/tutorial/getting_started.rst

@@ -186,7 +186,7 @@ you will get the gradient of size :math:`|V|d \times 1` with the order
     \frac{\partial B}{\partial x_n} &
     \frac{\partial B}{\partial y_n} &
     \frac{\partial B}{\partial z_n}
-    \end{bmatrix}^T,
+    \end{bmatrix}^\top,
 
 and the Hessian of size :math:`|V|d \times |V|d` with the order
 
@@ -417,7 +417,7 @@ This returns a scalar value ``friction_potential`` which is the sum of the indiv
 Mathematically this is defined as
 
 .. math::
-   D(x) = \sum_{k \in C} \mu\lambda_k^nf_0\left(\|T_k^Tv\|, \epsilon_v\right),
+   D(v) = \sum_{k \in C} \mu\lambda_k^nf_0\left(\|T_k^\top v\|, \epsilon_v\right),
 
 where :math:`C` is the lagged collisions, :math:`\lambda_k^n` is the normal force magnitude for the :math:`k`-th collision, :math:`T_k` is the tangential basis for the :math:`k`-th collision, and :math:`f_0` is the smooth friction function used to approximate the non-smooth transition from dynamic to static friction.
 

src/ipc/barrier/barrier.hpp

@@ -136,4 +136,66 @@ class ClampedLogBarrier : public Barrier {
     }
 };
 
+// ============================================================================
+// Normalized Barrier functions from [Li et al. 2023]
+// ============================================================================
+
+/// @brief Normalized barrier function from [Li et al. 2023].
+class NormalizedClampedLogBarrier : public ClampedLogBarrier {
+public:
+    NormalizedClampedLogBarrier() = default;
+
+    /// @brief Function that grows to infinity as d approaches 0 from the right.
+    ///
+    /// \f\[
+    ///     b(d) =
+    ///     -\left(\frac{d}{\hat{d}}-1\right)^2\ln\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
+    {
+        return ClampedLogBarrier::operator()(d / dhat, 1.0);
+    }
+
+    /// @brief Derivative of the barrier function.
+    ///
+    /// \f\[
+    ///     b'(d) =
+    ///     2\frac{1}{\hat{d}}\left(1-\frac{d}{\hat{d}}\right)\ln\left(\frac{d}{\hat{d}}\right)
+    ///             + \left(1-\frac{d}{\hat{d}}\right)^2 \frac{1}{d}
+    /// \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
+    {
+        return ClampedLogBarrier::first_derivative(d / dhat, 1.0) / dhat;
+    }
+
+    /// @brief Second derivative of the barrier function.
+    ///
+    /// \f\[
+    ///     b''(d) = \frac{\hat{d}^2-2 d^2 \ln \left(\frac{d}{\hat{d}}\right)+2
+    ///     \hat{d} d-3 d^2}{\hat{d}^2 d^2}
+    /// \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
+    {
+        return ClampedLogBarrier::second_derivative(d / dhat, 1.0)
+            / (dhat * dhat);
+    }
+
+    /// @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 { return 1.0; }
+};
+
 } // namespace ipc

tests/src/tests/barrier/test_barrier.cpp

@@ -12,71 +12,27 @@
 
 namespace ipc {
 
-class NormalizedClampedLogBarrier : public ipc::Barrier {
-public:
-    double operator()(const double d, const double dhat) const override
-    {
-        if (d <= 0.0) {
-            return std::numeric_limits<double>::infinity();
-        }
-        if (d >= dhat) {
-            return 0;
-        }
-
-        // b(d) = -(d/d̂-1)²ln(d / d̂)
-        const auto d_dhat = d / dhat;
-        const auto d_dhat_minus_1 = d_dhat - 1;
-        return -d_dhat_minus_1 * d_dhat_minus_1 * std::log(d_dhat);
-    }
-
-    double first_derivative(const double d, const double dhat) const override
-    {
-        if (d <= 0.0 || d >= dhat) {
-            return 0.0;
-        }
-        const double t0 = 1.0 / dhat;
-        const double t1 = d * t0;
-        const double t2 = 1 - t1;
-        return t2 * (2 * t0 * std::log(t1) - t2 / d);
-    }
-
-    double second_derivative(const double d, const double dhat) const override
-    {
-        if (d <= 0.0 || d >= dhat) {
-            return 0.0;
-        }
-
-        const double t0 = 1.0 / dhat;
-        const double t1 = d * t0;
-        const double t2 = 1 - t1;
-        return 4 * t0 * t2 / d + (t2 * t2) / (d * d)
-            - 2 * std::log(t1) / (dhat * dhat);
-    }
-
-    double units(const double dhat) const override { return 1; }
-};
-
 /// @warning This implementation will not work with dmin > 0
-class PhysicalBarrier : public ipc::NormalizedClampedLogBarrier {
+class PhysicalBarrier : public NormalizedClampedLogBarrier {
 public:
     PhysicalBarrier(const bool _use_dist_sqr) : use_dist_sqr(_use_dist_sqr) { }
 
     double operator()(const double d, const double dhat) const override
     {
         return (use_dist_sqr ? sqrt(dhat) : dhat)
-            * ipc::NormalizedClampedLogBarrier::operator()(d, dhat);
+            * NormalizedClampedLogBarrier::operator()(d, dhat);
     }
 
     double first_derivative(const double d, const double dhat) const override
     {
         return (use_dist_sqr ? sqrt(dhat) : dhat)
-            * ipc::NormalizedClampedLogBarrier::first_derivative(d, dhat);
+            * NormalizedClampedLogBarrier::first_derivative(d, dhat);
     }
 
     double second_derivative(const double d, const double dhat) const override
     {
         return (use_dist_sqr ? sqrt(dhat) : dhat)
-            * ipc::NormalizedClampedLogBarrier::second_derivative(d, dhat);
+            * NormalizedClampedLogBarrier::second_derivative(d, dhat);
     }
 
 private: