Replace smooth mu with a piecewise quadratic function (#179)

Author: zfergus

docs/source/_static/img/int_mu_f1_dx.png

@@ -201,4 +201,4 @@ We integrate their product to obtain a smooth tangential adhesion mollifier:
 .. figure:: /_static/img/int_mu_f1a_dx.png
    :align: center
 
-   Integrated mollifier :math:`\int \mu(y) f_1^a(y) \mathrm{d} y` where :math:`\mu_s = 1`, :math:`\mu_k = 0.1`, and :math:`\epsilon_v = 0.001`.
+   Integrated mollifier :math:`\int \mu(y) f_1^a(y) \mathrm{d} y` where :math:`y=\|\mathbf{u}\|`, :math:`\mu_s = 1`, :math:`\mu_k = 0.1`, and :math:`\epsilon_v = 0.001`.

docs/source/_static/img/int_mu_f1a_dx.png

@@ -108,8 +108,10 @@ Smooth :math:`\mu`
 When adding separate coefficients for static and kinetic friction, we need to maintain the :math:`C^1` continuity of the friction force. This lead us to define a smooth coefficient of friction :math:`\mu(y)` that transitions between the static and kinetic coefficients based on the magnitude of the relative velocity :math:`y = \|\mathbf{u}\|`. The smooth coefficient of friction is defined as
 
 .. math::
-    \mu(y) = \begin{cases} \mu_{s} + \left(\mu_{k} - \mu_{s}\right) \left(\frac{3 y^{2}}{\epsilon_{v}^{2}} - \frac{2 y^{3}}{\epsilon_{v}^{3}}\right) & \text{for}\: 0 \leq y \leq \epsilon_{v}\\
-        \mu_{k} & \text{otherwise}.
+    \mu(y) = \begin{cases}
+        2(\mu_{k} - \mu_{s})\frac{y^{2}}{\epsilon_{v}^{2}} + \mu_{s} & \text{for}\: y \leq \frac{\epsilon_{v}}{2} \\
+        -2(\mu_{k} - \mu_{s})\frac{\left(\epsilon_{v} - y\right)^{2}}{\epsilon_{v}^{2}} + \mu_k & \text{for}\: y \leq \epsilon_{v} \\
+        \mu_{k} & \text{for}\: y > \epsilon_{v}
     \end{cases}
 
 We plot the smooth coefficient of friction :math:`\mu(y)` below:
@@ -132,13 +134,9 @@ Replacing the constant coefficient of friction :math:`\mu` with a smooth functio
 However, :math:`\frac{\mathrm{d}}{\mathrm{d}y} \mu(y) f_0(y) \neq \mu(y) f_1(y)`, so we need to adjust the integrated mollifier by integrating the product of the smooth coefficient of friction and the mollifier:
 
 .. math::
-    \int \mu(y) f_1(y) \mathrm{d} y = \begin{cases}
-        \frac{\epsilon_{v}^{6} \left(17 \mu_{k} - 7 \mu_{s}\right) + 30 \epsilon_{v}^{4} \mu_{s} y^{2} - 10 \epsilon_{v}^{3} \mu_{s} y^{3} + 45 \epsilon_{v}^{2} \left(\mu_{k} - \mu_{s}\right) y^{4} + 42 \epsilon_{v} \left(- \mu_{k} + \mu_{s}\right) y^{5} + 10 \left(\mu_{k} - \mu_{s}\right) y^{6}}{30 \epsilon_{v}^{5}} & \text{for}\: \epsilon_{v} > y \\
-        \mu_{k} y & \text{otherwise}
-    \end{cases}.
-
-
-
+    \int \mu(y) f_1(y) \mathrm{d} y = \begin{cases} \frac{\frac{\epsilon_{v}^{5} \left(27 \mu_{k} - 11 \mu_{s}\right)}{48} + \epsilon_{v}^{3} \mu_{s} y^{2} - \frac{\epsilon_{v}^{2} \mu_{s} y^{3}}{3} + \epsilon_{v} y^{4} \left(\mu_{k} - \mu_{s}\right) + \frac{2 y^{5} \left(- \mu_{k} + \mu_{s}\right)}{5}}{\epsilon_{v}^{4}} & \text{for}\: y \leq \frac{\epsilon_{v}}{2} \\
+    \frac{\epsilon_{v}^{5} \left(9 \mu_{k} - 4 \mu_{s}\right) + 15 \epsilon_{v}^{3} y^{2} \left(- \mu_{k} + 2 \mu_{s}\right) + 5 \epsilon_{v}^{2} y^{3} \left(9 \mu_{k} - 10 \mu_{s}\right) - 30 \epsilon_{v} y^{4} \left(\mu_{k} - \mu_{s}\right) + 6 y^{5} \left(\mu_{k} - \mu_{s}\right)}{15 \epsilon_{v}^{4}} & \text{for}\: y \leq \epsilon_{v} \\
+    \mu_{k} y & \text{for}\: y > \epsilon_v \end{cases}
 
 The following plot shows the behavior of the integrated mollifier multiplied by the smooth coefficient of friction:
 
@@ -155,11 +153,10 @@ While this approach provides a smooth transition between static and kinetic fric
 1. The product :math:`\mu(y) f_1(y)` underestimates the friction force in the static regime, which may lead to less accurate simulations of static friction.
     - This is, when :math:`\mu_k < \mu_s` and :math:`|y| \leq \epsilon_v`, :math:`\max(\mu(y) f_1(y)) < \mu_s`.
     - We could address this by scaling by :math:`\frac{\max(\mu(y) f_1(y))}{\mu_s}`, but computing the maximum is non-trivial.
-2. The combined function :math:`\mu(y) f_1(y)` is a degree 5 polynomial, which is more complex than the original mollifier :math:`f_1(y)` (degree 2). This may lead to more difficult to optimize problems. There are two options to address this:
-    a. Replacing :math:`\mu(y)` with a piecewise quadratic function would reduce the degree to 4, but it would still be more complex than the original mollifier.
-    b. Replacing :math:`\mu(y) f_1(y)` with a piecewise quadratic function that has the desired behavior. This help with 1. as well. However making it backwards compatible with the original mollifier is challenging.
+2. The combined function :math:`\mu(y) f_1(y)` is a degree 4 polynomial, which is more complex than the original mollifier :math:`f_1(y)` (degree 2). This may lead to more difficult to optimize problems.
+    - We could replace :math:`\mu(y) f_1(y)` with a piecewise quadratic function that has the desired behavior. This help with (1) as well. However making it backwards compatible with the original mollifier is challenging.
 
-If you have suggestions for improving this approach or alternative methods, please reach out on our `GitHub Discussions <https://github.com/ipc-sim/ipc-toolkit/discussions>`.
+If you have suggestions for improving this approach or alternative methods, please reach out on our `GitHub Discussions <https://github.com/ipc-sim/ipc-toolkit/discussions>`_.
 
 Future Directions
 -----------------

docs/source/_static/img/mu.png

@@ -135,15 +135,25 @@ double smooth_mu_a0(
 {
     assert(eps_a > 0);
     const double delta_mu = mu_k - mu_s;
-    if (mu_s == mu_k || y <= 0 || y >= eps_a) {
+    if (y <= 0) {
+        return 0;
+    } else if (mu_s == mu_k || y >= eps_a) {
         // If the static and kinetic friction coefficients are equal, simplify.
-        const double c = (7 / 30.) * eps_a * delta_mu;
+        const double c = (11 / 48.) * eps_a * delta_mu;
         return mu_k * tangential_adhesion_f0(y, eps_a) - c;
     } else {
         const double z = y / eps_a;
-        return y * z
-            * (z * (z * (z * (z / 3 - 1.4) + 1.5) * delta_mu - mu_s / 3)
-               + mu_s);
+        if (y < 0.5 * eps_a) {
+            return y * z
+                * (z * (z * (1 - 0.4 * z) * delta_mu - mu_s / 3.0) + mu_s);
+        } else {
+            return y * z
+                * (z
+                       * (z * (0.4 * z - 2) * delta_mu + 3 * mu_k
+                          - (10.0 / 3.0) * mu_s)
+                   - mu_k + 2 * mu_s)
+                + (3.0 / 80.0) * eps_a * delta_mu;
+        }
     }
 }
 
@@ -174,14 +184,20 @@ double smooth_mu_a2_x_minus_mu_a1_over_x3(
     const double y, const double mu_s, const double mu_k, const double eps_a)
 {
     assert(eps_a > 0);
-    if (mu_s == mu_k || y <= 0 || y >= eps_a) {
+    assert(y >= 0);
+    if (mu_s == mu_k || y >= eps_a) {
         // If the static and kinetic friction coefficients are equal, simplify.
         return mu_k * tangential_adhesion_f2_x_minus_f1_over_x3(y, eps_a);
     } else {
         const double delta_mu = mu_k - mu_s;
         const double z = 1 / eps_a;
-        return z * z
-            * (z * (z * y * (z * y * 8 - 21) + 12) * delta_mu - mu_s / y);
+        if (y < 0.5 * eps_a) {
+            return z * z * (z * (8 - 6 * z * y) * delta_mu - mu_s / y);
+        } else {
+            return z * z
+                * (z * (6 * z * y - 16) * delta_mu
+                   + (9 * mu_k - 10 * mu_s) / y);
+        }
     }
 }
 

docs/source/_static/img/mu_f1.png

@@ -15,10 +15,12 @@ double smooth_mu(
         // If the static and kinetic friction coefficients are equal, simplify.
         return mu_k;
     } else {
-        const double y_over_eps_v = std::abs(y) / eps_v;
-        return (mu_k - mu_s) * (3 - 2 * y_over_eps_v) * y_over_eps_v
-            * y_over_eps_v
-            + mu_s;
+        const double z = std::abs(y) / eps_v;
+        if (std::abs(y) < 0.5 * eps_v) {
+            return 2 * (mu_k - mu_s) * z * z + mu_s;
+        } else {
+            return -2 * (mu_k - mu_s) * (z * (z - 2) + 1) + mu_k;
+        }
     }
 }
 
@@ -30,8 +32,12 @@ double smooth_mu_derivative(
         // If the static and kinetic friction coefficients are equal, simplify.
         return 0;
     } else {
-        const double y_over_eps_v = std::abs(y) / eps_v;
-        return (mu_k - mu_s) * (6 - 6 * y_over_eps_v) * y_over_eps_v / eps_v;
+        const double z = std::abs(y) / eps_v;
+        if (std::abs(y) < 0.5 * eps_v) {
+            return 4 * (mu_k - mu_s) * z / eps_v;
+        } else {
+            return -4 * (mu_k - mu_s) * (z - 1) / eps_v;
+        }
     }
 }
 
@@ -44,11 +50,19 @@ double smooth_mu_f0(
         return mu_k * smooth_friction_f0(y, eps_v);
     } else {
         const double delta_mu = mu_k - mu_s;
-        const double c = eps_v * (17 * mu_k - 7 * mu_s) / 30;
         const double z = std::abs(y) / eps_v;
-        return y * z
-            * (z * (z * (z * (z / 3 - 1.4) + 1.5) * delta_mu - mu_s / 3) + mu_s)
-            + c;
+        if (std::abs(y) < 0.5 * eps_v) {
+            return y * z
+                * (z * (z * (1 - 0.4 * z) * delta_mu - mu_s / 3.0) + mu_s)
+                + (9.0 / 16.0) * eps_v * mu_k - (11.0 / 48.0) * eps_v * mu_s;
+        } else {
+            return y * z
+                * (z
+                       * (z * (0.4 * z - 2) * delta_mu
+                          + (3 * mu_k - (10.0 / 3.0) * mu_s))
+                   + (2 * mu_s - mu_k))
+                + 0.6 * eps_v * mu_k - (4.0 / 15.0) * eps_v * mu_s;
+        }
     }
 }
 
@@ -82,13 +96,19 @@ double smooth_mu_f2_x_minus_mu_f1_over_x3(
 {
     assert(eps_v > 0);
     if (mu_s == mu_k || std::abs(y) >= eps_v) {
-        // If the static and kinetic friction coefficients are equal, simplify.
+        // If the static and kinetic friction coefficients are equal,
+        // simplify.
         return mu_k * smooth_friction_f2_x_minus_f1_over_x3(y, eps_v);
     } else {
         const double delta_mu = mu_k - mu_s;
         const double z = 1 / eps_v;
-        return z * z
-            * (z * (z * y * (z * y * 8 - 21) + 12) * delta_mu - mu_s / y);
+        if (std::abs(y) < 0.5 * eps_v) {
+            return z * z * (z * (8 - 6 * y * z) * delta_mu - mu_s / y);
+        } else {
+            return z * z
+                * (z * (6 * y * z - 16) * delta_mu
+                   + (9 * mu_k - 10 * mu_s) / y);
+        }
     }
 }