diff --git a/docs/source/cpp-api/adhesion.rst b/docs/source/cpp-api/adhesion.rst
index 7bfecb3c7..550fcda9a 100644
--- a/docs/source/cpp-api/adhesion.rst
+++ b/docs/source/cpp-api/adhesion.rst
@@ -27,6 +27,16 @@ Functions
- The first derivative of the tangential adhesion mollifier function divided by y.
* - :func:`tangential_adhesion_f2_x_minus_f1_over_x3`
- The second derivative of the tangential adhesion mollifier function times y minus the first derivative all divided by y cubed.
+ * - :func:`smooth_mu_a0`
+ - Compute the value of the ∫ μ(y) a₁(y) dy, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+ * - :func:`smooth_mu_a1`
+ - Compute the value of the μ(y) a₁(y), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+ * - :func:`smooth_mu_a2`
+ - Compute the value of d/dy (μ(y) a₁(y)), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+ * - :func:`smooth_mu_a1_over_x`
+ - Compute the value of the μ(y) a₁(y) / y, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+ * - :func:`smooth_mu_a2_x_minus_mu_a1_over_x3`
+ - Compute the value of the [(d/dy μ(y) a₁(y)) ⋅ y - μ(y) a₁(y)] / y³, where a₁ and a₂ are the first and second derivatives of the smooth tangential adhesion mollifier.
Normal Adhesion Potential
diff --git a/docs/source/cpp-api/friction.rst b/docs/source/cpp-api/friction.rst
index e794c8050..8dfa3d8d8 100644
--- a/docs/source/cpp-api/friction.rst
+++ b/docs/source/cpp-api/friction.rst
@@ -8,4 +8,16 @@ Smooth Mollifier
.. doxygenfunction:: smooth_friction_f1
.. doxygenfunction:: smooth_friction_f2
.. doxygenfunction:: smooth_friction_f1_over_x
-.. doxygenfunction:: smooth_friction_f2_x_minus_f1_over_x3
\ No newline at end of file
+.. doxygenfunction:: smooth_friction_f2_x_minus_f1_over_x3
+
+Smooth :math:`\mu`
+------------------
+
+.. doxygenfunction:: smooth_mu
+.. doxygenfunction:: smooth_mu_derivative
+
+.. doxygenfunction:: smooth_mu_f0
+.. doxygenfunction:: smooth_mu_f1
+.. doxygenfunction:: smooth_mu_f2
+.. doxygenfunction:: smooth_mu_f1_over_x
+.. doxygenfunction:: smooth_mu_f2_x_minus_f1_over_x3
\ No newline at end of file
diff --git a/docs/source/python-api/adhesion.rst b/docs/source/python-api/adhesion.rst
index 86392c0b8..4b090b702 100644
--- a/docs/source/python-api/adhesion.rst
+++ b/docs/source/python-api/adhesion.rst
@@ -16,4 +16,13 @@ Tangential Adhesion Potential
.. autofunction:: ipctk.tangential_adhesion_f1
.. autofunction:: ipctk.tangential_adhesion_f2
.. autofunction:: ipctk.tangential_adhesion_f1_over_x
-.. autofunction:: ipctk.tangential_adhesion_f2_x_minus_f1_over_x3
\ No newline at end of file
+.. autofunction:: ipctk.tangential_adhesion_f2_x_minus_f1_over_x3
+
+Smooth :math:`\mu`
+------------------
+
+.. autofunction:: ipctk.smooth_mu_a0
+.. autofunction:: ipctk.smooth_mu_a1
+.. autofunction:: ipctk.smooth_mu_a2
+.. autofunction:: ipctk.smooth_mu_a1_over_x
+.. autofunction:: ipctk.smooth_mu_a2_x_minus_mu_a1_over_x3
\ No newline at end of file
diff --git a/docs/source/python-api/friction.rst b/docs/source/python-api/friction.rst
index df5c1bd37..2ea1db770 100644
--- a/docs/source/python-api/friction.rst
+++ b/docs/source/python-api/friction.rst
@@ -8,4 +8,15 @@ Smooth Mollifier
.. autofunction:: ipctk.smooth_friction_f1
.. autofunction:: ipctk.smooth_friction_f2
.. autofunction:: ipctk.smooth_friction_f1_over_x
-.. autofunction:: ipctk.smooth_friction_f2_x_minus_f1_over_x3
\ No newline at end of file
+.. autofunction:: ipctk.smooth_friction_f2_x_minus_f1_over_x3
+
+Smooth :math:`\mu`
+------------------
+
+.. autofunction:: ipctk.smooth_mu
+.. autofunction:: ipctk.smooth_mu_derivative
+.. autofunction:: ipctk.smooth_mu_f0
+.. autofunction:: ipctk.smooth_mu_f1
+.. autofunction:: ipctk.smooth_mu_f2
+.. autofunction:: ipctk.smooth_mu_f1_over_x
+.. autofunction:: ipctk.smooth_mu_f2_x_minus_f1_over_x3
\ No newline at end of file
diff --git a/notebooks/separate_mu.ipynb b/notebooks/separate_mu.ipynb
new file mode 100644
index 000000000..c0ec5a66b
--- /dev/null
+++ b/notebooks/separate_mu.ipynb
@@ -0,0 +1,5447 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Smooth Static to Dynamic Coefficient of Friction"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 208,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy import *\n",
+ "import numpy as np\n",
+ "import plotly.graph_objects as go"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 209,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ " \n",
+ " \n",
+ " "
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "import plotly\n",
+ "from IPython.display import display, HTML\n",
+ "# Tomas Mazak's workaround\n",
+ "plotly.offline.init_notebook_mode()\n",
+ "display(HTML(\n",
+ " ''\n",
+ "))\n",
+ "##"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 210,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x = Symbol('x')\n",
+ "eps_v = Symbol(r'\\epsilon_v')\n",
+ "mu_s, mu_k = symbols(r'\\mu_s \\mu_k')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 211,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot(*fs, title=None, min_x=-2, max_x=2, **subs):\n",
+ " subs = {eps_v: subs.get(\n",
+ " \"eps_v\", 1e-3), mu_s: subs.get(\"mu_s\", 1), mu_k: subs.get(\"mu_k\", 0.1)}\n",
+ " xs = np.linspace(min_x*subs[eps_v], max_x*subs[eps_v], 201)\n",
+ "\n",
+ " def ys(f):\n",
+ " return np.array([f.subs({x: xi} | subs) for xi in xs], dtype=float)\n",
+ "\n",
+ " go.Figure(\n",
+ " [go.Scatter(x=xs, y=ys(f)) for f in fs],\n",
+ " layout=dict(\n",
+ " width=800, height=600, template=\"none\",\n",
+ " xaxis_title=r'$x$', title=title\n",
+ " )).show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Friction Mollifier"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 212,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle f_{0}(x) = \\begin{cases} \\left|{x}\\right| & \\text{for}\\: \\epsilon_{v} \\leq \\left|{x}\\right| \\\\\\frac{\\epsilon_{v}}{3} + \\frac{\\left|{x}\\right|^{2}}{\\epsilon_{v}} - \\frac{\\left|{x}\\right|^{3}}{3 \\epsilon_{v}^{2}} & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(f_{0}(x), Piecewise((Abs(x), \\epsilon_v <= Abs(x)), (\\epsilon_v/3 + Abs(x)**2/\\epsilon_v - Abs(x)**3/(3*\\epsilon_v**2), True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL8q499nXDhgv1gczvxqDmC/Dat4I/PIX79pHVVNEHVfv8WPMXctIV+/IQIOoUrNXr9+dOrKZ3lev9rmxvSEJV6/NlmjHqLRXb+Sy39Iv31dv+89XHLcKV2/S7A4nPnVXL+nIhXGFoJcvwOV8e8zLly/YAfOGVHaW7+8eapDboZbvxjshm2LMlu/dF5jl6jeWr/Q0D/BxYpavy1DHOviNlq/ibX4FADjWb/lJ9U+HY9Zv0KasWg6O1m/ngyOklfnWL/6fmq8dJNYv1bxRuaRP1i/smMjEK/rV78O1v85zJdXv2tI3GPpQ1e/x7q4jQbwVr8kLZW3I5xWv4CfceFASFa/3BFOC170Vb84hCo1e6BVv5T2Bl+YTFW/8GjjiLX4VL9M27+y0qRUv6lNnNzvUFS/BcB4Bg39U79iMlUwKqlTv76kMVpHVVO/GhcOhGQBU792ieqtga1Sv9L7xteeWVK/Lm6jAbwFUr+L4H8r2bFRv+dSXFX2XVG/Q8U4fxMKUb+fNxWpMLZQv/yp8dJNYlC/WBzO/GoOUL9oHVVNEHVPvyACDqFKzU6/2ubG9IQlTr+Sy39Iv31Nv0qwOJz51Uy/ApXx7zMuTL+8eapDboZLv3ReY5eo3kq/LEMc6+I2Sr/kJ9U+HY9Jv54MjpJX50i/VvFG5pE/SL8O1v85zJdHv8a6uI0G8Ea/fp9x4UBIRr84hCo1e6BFv/Bo44i1+ES/qE2c3O9QRL9gMlUwKqlDvxoXDoRkAUO/0vvG155ZQr+K4H8r2bFBv0LFOH8TCkG//Knx0k1iQL9oHVVNEHU/v9jmxvSEJT6/SLA4nPnVPL+8eapDboY7vyxDHOviNjq/nAyOklfnOL8M1v85zJc3v3yfceFASDa/8GjjiLX4NL9gMlUwKqkzv9D7xteeWTK/QMU4fxMKMb9oHVVNEHUvv0iwOJz51Sy/KEMc6+I2Kr8I1v85zJcnv/Bo44i1+CS/0PvG155ZIr9gHVVNEHUfvyBDHOviNhq/4GjjiLX4FL9gHVVNEHUPv8Bo44i1+AS/AGnjiLX49L4AAAAAAAAAAABp44i1+PQ+AGnjiLX4BD+AHVVNEHUPPwBp44i1+BQ/QEMc6+I2Gj+AHVVNEHUfP9D7xteeWSI/8GjjiLX4JD8Q1v85zJcnPzBDHOviNio/ULA4nPnVLD9wHVVNEHUvP0jFOH8TCjE/2PvG155ZMj9oMlUwKqkzP/Bo44i1+DQ/gJ9x4UBINj8Q1v85zJc3P6AMjpJX5zg/MEMc6+I2Oj/AeapDboY7P1CwOJz51Tw/4ObG9IQlPj9oHVVNEHU/P/yp8dJNYkA/RMU4fxMKQT+M4H8r2bFBP9T7xteeWUI/HBcOhGQBQz9kMlUwKqlDP6xNnNzvUEQ/9GjjiLX4RD84hCo1e6BFP4CfceFASEY/yLq4jQbwRj8Q1v85zJdHP1jxRuaRP0g/oAyOklfnSD/oJ9U+HY9JPzBDHOviNko/dF5jl6jeSj+8eapDboZLPwSV8e8zLkw/TLA4nPnVTD+Uy39Iv31NP9zmxvSEJU4/JAIOoUrNTj9sHVVNEHVPP1oczvxqDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "/Knx0k1iYD8q499nXDhgP1gczvxqDmA/Dat4I/PIXz9pHVVNEHVfP8WPMXctIV8/IQIOoUrNXj9+dOrKZ3leP9rmxvSEJV4/NlmjHqLRXT+Sy39Iv31dP+89XHLcKV0/S7A4nPnVXD+nIhXGFoJcPwOV8e8zLlw/YAfOGVHaWz+8eapDboZbPxjshm2LMls/dF5jl6jeWj/Q0D/BxYpaPy1DHOviNlo/ibX4FADjWT/lJ9U+HY9ZP0KasWg6O1k/ngyOklfnWD/6fmq8dJNYP1bxRuaRP1g/smMjEK/rVz8O1v85zJdXP2tI3GPpQ1c/x7q4jQbwVj8kLZW3I5xWP4CfceFASFY/3BFOC170VT84hCo1e6BVP5T2Bl+YTFU/8GjjiLX4VD9M27+y0qRUP6lNnNzvUFQ/BcB4Bg39Uz9iMlUwKqlTP76kMVpHVVM/GhcOhGQBUz92ieqtga1SP9L7xteeWVI/Lm6jAbwFUj+L4H8r2bFRP+dSXFX2XVE/Q8U4fxMKUT+fNxWpMLZQP/yp8dJNYlA/XhjQ2W0OUD/Q3XQdPnVPP34reT/lzU4/FurEdfMmTj/k6W8ci4BNPzb7kY/O2kw/Wu5CK+A1TD+gk5pL4pFLP1K7sEz37ko/vjWdikFNSj8y03dh46xJPwBkWC3/DUk/cLhWSrdwSD/SoIoULtVHP3TtC+iFO0c/oW7yIOGjRj+s9FUbYg5GP95PTjMre0U/hlDzxF7qRD/wxlwsH1xEP3CDosWO0EM/TFbc7M9HQz/VDyL+BMJCP1iAi1VQP0I/JngwT9S/QT+IxyhHs0NBP84+jJkPy0A/RK5yogtWQD93zOd7k8k/P/xtT5DY7j4/teFLOiscPj8+yAwy0FE9PzPCwS8MkDw/MHCa6yPXOz/OcsYdXCc7P6lqdX75gDo/XPjWxUDkOT+GvBqsdlE5P7xXcOnfyDg/nmoHNsFKOD/GlQ9KX9c3P9F5uN3+bjc/V7cxqeQRNz/17qpkVcA2P0fBU8iVejY/6M5bjOpANj90uPJomBM2P4UeSBbk8jU/t6GLTBLfNT+l4uzDZ9g1P7ehi0wS3zU/hR5IFuTyNT91uPJomBM2P+rOW4zqQDY/ScFTyJV6Nj/37qpkVcA2P1e3MankETc/0Xm43f5uNz/HlQ9KX9c3P6BqBzbBSjg/vldw6d/IOD+HvBqsdlE5P2D41sVA5Dk/rGp1fvmAOj/ScsYdXCc7PzBwmusj1zs/NsLBLwyQPD9ByAwy0FE9P7jhSzorHD4//m1PkNjuPj96zOd7k8k/P0eucqILVkA/0D6MmQ/LQD+IxyhHs0NBPyZ4ME/Uv0E/WoCLVVA/Qj/XDyL+BMJCP05W3OzPR0M/cYOixY7QQz/0xlwsH1xEP4hQ88Re6kQ/4E9OMyt7RT+s9FUbYg5GP6Nu8iDho0Y/de0L6IU7Rz/UoIoULtVHP3K4Vkq3cEg/AmRYLf8NST8403dh46xJP8I1nYpBTUo/UruwTPfuSj+gk5pL4pFLP17uQivgNUw/OPuRj87aTD/m6W8ci4BNPxjqxHXzJk4/git5P+XNTj/U3XQdPnVPP2EY0NltDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$f_0(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def f0(x):\n",
+ " \"\"\"Smooth friction mollifier\"\"\"\n",
+ " return x * x * (1 - x / (3 * eps_v)) / eps_v + eps_v / 3\n",
+ "\n",
+ "\n",
+ "sym_f0 = Piecewise(\n",
+ " (abs(x), abs(x) >= eps_v),\n",
+ " (f0(abs(x)), abs(x) < eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(\"f_{0}(x)\"), sym_f0.expand()))\n",
+ "\n",
+ "plot(sym_f0, title=r'$f_0(x)$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 213,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle f_{1}(x) = \\begin{cases} -1 & \\text{for}\\: \\epsilon_{v} < - x \\\\\\frac{2 x}{\\epsilon_{v}} + \\frac{x^{2}}{\\epsilon_{v}^{2}} & \\text{for}\\: x \\leq 0 \\\\\\frac{2 x}{\\epsilon_{v}} - \\frac{x^{2}}{\\epsilon_{v}^{2}} & \\text{for}\\: \\epsilon_{v} \\geq x \\\\1 & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(f_{1}(x), Piecewise((-1, \\epsilon_v < -x), (2*x/\\epsilon_v + x**2/\\epsilon_v**2, x <= 0), (2*x/\\epsilon_v - x**2/\\epsilon_v**2, \\epsilon_v >= x), (1, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL8q499nXDhgv1gczvxqDmC/Dat4I/PIX79pHVVNEHVfv8WPMXctIV+/IQIOoUrNXr9+dOrKZ3lev9rmxvSEJV6/NlmjHqLRXb+Sy39Iv31dv+89XHLcKV2/S7A4nPnVXL+nIhXGFoJcvwOV8e8zLly/YAfOGVHaW7+8eapDboZbvxjshm2LMlu/dF5jl6jeWr/Q0D/BxYpavy1DHOviNlq/ibX4FADjWb/lJ9U+HY9Zv0KasWg6O1m/ngyOklfnWL/6fmq8dJNYv1bxRuaRP1i/smMjEK/rV78O1v85zJdXv2tI3GPpQ1e/x7q4jQbwVr8kLZW3I5xWv4CfceFASFa/3BFOC170Vb84hCo1e6BVv5T2Bl+YTFW/8GjjiLX4VL9M27+y0qRUv6lNnNzvUFS/BcB4Bg39U79iMlUwKqlTv76kMVpHVVO/GhcOhGQBU792ieqtga1Sv9L7xteeWVK/Lm6jAbwFUr+L4H8r2bFRv+dSXFX2XVG/Q8U4fxMKUb+fNxWpMLZQv/yp8dJNYlC/WBzO/GoOUL9oHVVNEHVPvyACDqFKzU6/2ubG9IQlTr+Sy39Iv31Nv0qwOJz51Uy/ApXx7zMuTL+8eapDboZLv3ReY5eo3kq/LEMc6+I2Sr/kJ9U+HY9Jv54MjpJX50i/VvFG5pE/SL8O1v85zJdHv8a6uI0G8Ea/fp9x4UBIRr84hCo1e6BFv/Bo44i1+ES/qE2c3O9QRL9gMlUwKqlDvxoXDoRkAUO/0vvG155ZQr+K4H8r2bFBv0LFOH8TCkG//Knx0k1iQL9oHVVNEHU/v9jmxvSEJT6/SLA4nPnVPL+8eapDboY7vyxDHOviNjq/nAyOklfnOL8M1v85zJc3v3yfceFASDa/8GjjiLX4NL9gMlUwKqkzv9D7xteeWTK/QMU4fxMKMb9oHVVNEHUvv0iwOJz51Sy/KEMc6+I2Kr8I1v85zJcnv/Bo44i1+CS/0PvG155ZIr9gHVVNEHUfvyBDHOviNhq/4GjjiLX4FL9gHVVNEHUPv8Bo44i1+AS/AGnjiLX49L4AAAAAAAAAAABp44i1+PQ+AGnjiLX4BD+AHVVNEHUPPwBp44i1+BQ/QEMc6+I2Gj+AHVVNEHUfP9D7xteeWSI/8GjjiLX4JD8Q1v85zJcnPzBDHOviNio/ULA4nPnVLD9wHVVNEHUvP0jFOH8TCjE/2PvG155ZMj9oMlUwKqkzP/Bo44i1+DQ/gJ9x4UBINj8Q1v85zJc3P6AMjpJX5zg/MEMc6+I2Oj/AeapDboY7P1CwOJz51Tw/4ObG9IQlPj9oHVVNEHU/P/yp8dJNYkA/RMU4fxMKQT+M4H8r2bFBP9T7xteeWUI/HBcOhGQBQz9kMlUwKqlDP6xNnNzvUEQ/9GjjiLX4RD84hCo1e6BFP4CfceFASEY/yLq4jQbwRj8Q1v85zJdHP1jxRuaRP0g/oAyOklfnSD/oJ9U+HY9JPzBDHOviNko/dF5jl6jeSj+8eapDboZLPwSV8e8zLkw/TLA4nPnVTD+Uy39Iv31NP9zmxvSEJU4/JAIOoUrNTj9sHVVNEHVPP1oczvxqDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/AAAAAAAA8L8AAAAAAADwvwAAAAAAAPC/eJyiI7n877/ecYqO5PLvvzOAt0CC4u+/e8cpOpLL77+vR+F6FK7vv9EA3gIJiu+/5PIf0m9f77/nHafoSC7vv9mBc0aU9u6/uR6F61G47r+I9NvXgXPuv0cDeAskKO6/9UpZhjjW7b+Sy39Iv33tvx+F61G4Hu2/mnecoiO57L8Fo5I6AU3sv18HzhlR2uu/qKROQBNh67/gehSuR+HqvwiKH2PuWuq/H9JvXwfO6b8lUwWjkjrpvxoN4C2QoOi/AAAAAAAA6L/UK2UZ4ljnv5aQD3o2q+a/Ry7/If325b/rBDQRNjzlv3sUrkfheuS/+Vxtxf6y479o3nGKjuTiv8WYu5aQD+K/FYxK6gQ04b9RuB6F61Hgv/g6cM6I0t6/LHctIR/03L9GJXUCmgjbvzZFR3L5D9m/BtejcD0K17+12or9ZffUv0hQ/Bhz19K/szf4wmSq0L/4If32deDMv0a4HoXrUci/UzJVMCqpw79TIEHxY8y9v0GjAbwFErS/A9zXgXNGpL8AAAAAAAAAAAPc14FzRqQ/faMBvAUStD9yIEHxY8y9P28yVTAqqcM/YrgehetRyD8TIv32deDMP7M3+MJkqtA/SFD8GHPX0j+72or9ZffUPw3Xo3A9Ctc/PUVHcvkP2T9LJXUCmgjbPzh3LSEf9Nw/BDtwzojS3j9XuB6F61HgPxWMSuoENOE/yJi7lpAP4j9r3nGKjuTiP/xcbcX+suM/fRSuR+F65D/tBDQRNjzlP0su/yH99uU/mpAPejar5j/UK2UZ4ljnPwAAAAAAAOg/HA3gLZCg6D8oUwWjkjrpPyLSb18Hzuk/C4ofY+5a6j/kehSuR+HqP6ukTkATYes/YgfOGVHa6z8Fo5I6AU3sP5p3nKIjuew/HoXrUbge7T+Ty39Iv33tP/ZKWYY41u0/SAN4CyQo7j+K9NvXgXPuP7oehetRuO4/2YFzRpT27j/nHafoSC7vP+XyH9JvX+8/0gDeAgmK7z+uR+F6FK7vP3nHKTqSy+8/NIC3QILi7z/ecYqO5PLvP3icoiO5/O8/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$f_1(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def f1(x):\n",
+ " \"\"\"Derivative of f0\"\"\"\n",
+ " return x * (2 - x / eps_v) / eps_v\n",
+ "\n",
+ "\n",
+ "# def sign(x):\n",
+ "# return Piecewise((1, x > 0), (-1, x < 0), (0, True))\n",
+ "\n",
+ "\n",
+ "sym_f1 = Piecewise(\n",
+ " (-1, x < -eps_v),\n",
+ " (-f1(-x), x <= 0),\n",
+ " (f1(x), x <= eps_v),\n",
+ " (1, x > eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(\"f_{1}(x)\"), sym_f1.expand()))\n",
+ "\n",
+ "plot(sym_f1, title=r\"$f_1(x)$\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Smooth Coefficient of Friction\n",
+ "\n",
+ "We can use a polynomial to model a smooth transition from static to kinematic coefficient of friction."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 214,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\mu(x) = \\begin{cases} \\mu_{k} & \\text{for}\\: \\epsilon_{v} \\leq \\left|{x}\\right| \\\\\\mu_{s} + \\left(\\mu_{k} - \\mu_{s}\\right) \\left(\\frac{3 \\left|{x}\\right|^{2}}{\\epsilon_{v}^{2}} - \\frac{2 \\left|{x}\\right|^{3}}{\\epsilon_{v}^{3}}\\right) & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(\\mu(x), Piecewise((\\mu_k, \\epsilon_v <= Abs(x)), (\\mu_s + (\\mu_k - \\mu_s)*(3*Abs(x)**2/\\epsilon_v**2 - 2*Abs(x)**3/\\epsilon_v**3), True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL8q499nXDhgv1gczvxqDmC/Dat4I/PIX79pHVVNEHVfv8WPMXctIV+/IQIOoUrNXr9+dOrKZ3lev9rmxvSEJV6/NlmjHqLRXb+Sy39Iv31dv+89XHLcKV2/S7A4nPnVXL+nIhXGFoJcvwOV8e8zLly/YAfOGVHaW7+8eapDboZbvxjshm2LMlu/dF5jl6jeWr/Q0D/BxYpavy1DHOviNlq/ibX4FADjWb/lJ9U+HY9Zv0KasWg6O1m/ngyOklfnWL/6fmq8dJNYv1bxRuaRP1i/smMjEK/rV78O1v85zJdXv2tI3GPpQ1e/x7q4jQbwVr8kLZW3I5xWv4CfceFASFa/3BFOC170Vb84hCo1e6BVv5T2Bl+YTFW/8GjjiLX4VL9M27+y0qRUv6lNnNzvUFS/BcB4Bg39U79iMlUwKqlTv76kMVpHVVO/GhcOhGQBU792ieqtga1Sv9L7xteeWVK/Lm6jAbwFUr+L4H8r2bFRv+dSXFX2XVG/Q8U4fxMKUb+fNxWpMLZQv/yp8dJNYlC/WBzO/GoOUL9oHVVNEHVPvyACDqFKzU6/2ubG9IQlTr+Sy39Iv31Nv0qwOJz51Uy/ApXx7zMuTL+8eapDboZLv3ReY5eo3kq/LEMc6+I2Sr/kJ9U+HY9Jv54MjpJX50i/VvFG5pE/SL8O1v85zJdHv8a6uI0G8Ea/fp9x4UBIRr84hCo1e6BFv/Bo44i1+ES/qE2c3O9QRL9gMlUwKqlDvxoXDoRkAUO/0vvG155ZQr+K4H8r2bFBv0LFOH8TCkG//Knx0k1iQL9oHVVNEHU/v9jmxvSEJT6/SLA4nPnVPL+8eapDboY7vyxDHOviNjq/nAyOklfnOL8M1v85zJc3v3yfceFASDa/8GjjiLX4NL9gMlUwKqkzv9D7xteeWTK/QMU4fxMKMb9oHVVNEHUvv0iwOJz51Sy/KEMc6+I2Kr8I1v85zJcnv/Bo44i1+CS/0PvG155ZIr9gHVVNEHUfvyBDHOviNhq/4GjjiLX4FL9gHVVNEHUPv8Bo44i1+AS/AGnjiLX49L4AAAAAAAAAAABp44i1+PQ+AGnjiLX4BD+AHVVNEHUPPwBp44i1+BQ/QEMc6+I2Gj+AHVVNEHUfP9D7xteeWSI/8GjjiLX4JD8Q1v85zJcnPzBDHOviNio/ULA4nPnVLD9wHVVNEHUvP0jFOH8TCjE/2PvG155ZMj9oMlUwKqkzP/Bo44i1+DQ/gJ9x4UBINj8Q1v85zJc3P6AMjpJX5zg/MEMc6+I2Oj/AeapDboY7P1CwOJz51Tw/4ObG9IQlPj9oHVVNEHU/P/yp8dJNYkA/RMU4fxMKQT+M4H8r2bFBP9T7xteeWUI/HBcOhGQBQz9kMlUwKqlDP6xNnNzvUEQ/9GjjiLX4RD84hCo1e6BFP4CfceFASEY/yLq4jQbwRj8Q1v85zJdHP1jxRuaRP0g/oAyOklfnSD/oJ9U+HY9JPzBDHOviNko/dF5jl6jeSj+8eapDboZLPwSV8e8zLkw/TLA4nPnVTD+Uy39Iv31NP9zmxvSEJU4/JAIOoUrNTj9sHVVNEHVPP1oczvxqDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/YCjAZm/fuT+AEhZwKq26P+Czbigh/bs/IGmdAqrJvT8Qx7q4jQbAP5A/5fPlYME/IMw37Ajxwj/gGhzbIbTEPwDa+/lbp8Y/gLdAguLHyD+IYVSt4BLLPxiGoLSBhc0/uGnHaHgO0D/Ie8SeLGvRP0RQ/Bhz19I/UL4jdOFR1D/YnO9MDdnVPwjDFECMa9c/2AdI6vMH2T9eQj7o2azaP55JrNbTWNw/sPRGUncK3j+eGsP3WcDfPznJ6rGIvOA/m5mZmZmZ4T/8aUiBqnbiP+elUTcGU+M/3LgPivct5D9kDt1HyQblPwYSFD/G3OU/Si8PPjmv5j+z0SgTbX3nP8Zku4ysRug/DVQheUIK6T8QC7WmecfpP1P10OOcfeo/Wn7P/vYr6z+rEQvG0tHrP9Qa3gd7buw/VAWjkjoB7T+0PLQ0XIntP3ksbLwqBu4/K0Al+PB27j9R4zm2+druP3CBBMWPMe8/D4bf8v157z+yXCUOj7PvP+JwMOWN3e8/JC5bRkX37z8AAAAAAADwPyQuW0ZF9+8/4XAw5Y3d7z+xXCUOj7PvPw2G3/L9ee8/boEExY8x7z9O4zm2+druPytAJfjwdu4/eSxsvCoG7j+0PLQ0XIntP1IFo5I6Ae0/0hreB3tu7D+rEQvG0tHrP1V+z/72K+s/TPXQ45x96j8MC7WmecfpPw1UIXlCCuk/xWS7jKxG6D+v0SgTbX3nP0YvDz45r+Y/BBIUP8bc5T9gDt1HyQblP9e4D4r3LeQ/4aVRNwZT4z/8aUiBqnbiP5uZmZmZmeE/OMnqsYi84D+aGsP3WcDfP6z0RlJ3Ct4/nEms1tNY3D9UQj7o2azaP9AHSOrzB9k//MIUQIxr1z/YnO9MDdnVP0i+I3ThUdQ/RFD8GHPX0j/Ae8SeLGvRP7Rpx2h4DtA/CIagtIGFzT94YVSt4BLLP3C3QILix8g/ANr7+Vunxj/gGhzbIbTEPxjMN+wI8cI/gD/l8+VgwT8Ax7q4jQbAPyBpnQKqyb0/ALRuKCH9uz+AEhZwKq26P2AowGZv37k/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$\\mu(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def mu(x):\n",
+ " # return (mu_k - mu_s) * (0.5 * (-cos(pi * x / eps_v) + 1)) + mu_s\n",
+ " return (mu_k - mu_s) * (3 * (x / eps_v)**2 - 2 * (x / eps_v)**3) + mu_s\n",
+ "\n",
+ "\n",
+ "sym_mu = Piecewise(\n",
+ " (mu_k, abs(x) >= eps_v),\n",
+ " (mu(abs(x)), abs(x) < eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(r\"\\mu(x)\"), sym_mu))\n",
+ "\n",
+ "plot(sym_mu, title=r\"$\\mu(x)$\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Smooth Coefficient of Friction Molliifier\n",
+ "We know we want a smooth force mollifier of $\\mu(x)f_1(x)$, so we can integrate this function to get a mollifier of $\\int \\mu(x)f_1(x)dx$ that produces the desired force."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 215,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\mu(x) f_{1(x)} = \\begin{cases} - \\mu_{k} & \\text{for}\\: \\epsilon_{v} < - x \\\\\\frac{x \\left(2 + \\frac{x}{\\epsilon_{v}}\\right) \\left(\\mu_{s} + \\left(\\mu_{k} - \\mu_{s}\\right) \\left(\\frac{3 x^{2}}{\\epsilon_{v}^{2}} + \\frac{2 x^{3}}{\\epsilon_{v}^{3}}\\right)\\right)}{\\epsilon_{v}} & \\text{for}\\: x \\leq 0 \\\\\\frac{x \\left(2 - \\frac{x}{\\epsilon_{v}}\\right) \\left(\\mu_{s} + \\left(\\mu_{k} - \\mu_{s}\\right) \\left(\\frac{3 x^{2}}{\\epsilon_{v}^{2}} - \\frac{2 x^{3}}{\\epsilon_{v}^{3}}\\right)\\right)}{\\epsilon_{v}} & \\text{for}\\: \\epsilon_{v} \\geq x \\\\\\mu_{k} & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(\\mu(x) f_1(x), Piecewise((-\\mu_k, \\epsilon_v < -x), (x*(2 + x/\\epsilon_v)*(\\mu_s + (\\mu_k - \\mu_s)*(3*x**2/\\epsilon_v**2 + 2*x**3/\\epsilon_v**3))/\\epsilon_v, x <= 0), (x*(2 - x/\\epsilon_v)*(\\mu_s + (\\mu_k - \\mu_s)*(3*x**2/\\epsilon_v**2 - 2*x**3/\\epsilon_v**3))/\\epsilon_v, \\epsilon_v >= x), (\\mu_k, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL8q499nXDhgv1gczvxqDmC/Dat4I/PIX79pHVVNEHVfv8WPMXctIV+/IQIOoUrNXr9+dOrKZ3lev9rmxvSEJV6/NlmjHqLRXb+Sy39Iv31dv+89XHLcKV2/S7A4nPnVXL+nIhXGFoJcvwOV8e8zLly/YAfOGVHaW7+8eapDboZbvxjshm2LMlu/dF5jl6jeWr/Q0D/BxYpavy1DHOviNlq/ibX4FADjWb/lJ9U+HY9Zv0KasWg6O1m/ngyOklfnWL/6fmq8dJNYv1bxRuaRP1i/smMjEK/rV78O1v85zJdXv2tI3GPpQ1e/x7q4jQbwVr8kLZW3I5xWv4CfceFASFa/3BFOC170Vb84hCo1e6BVv5T2Bl+YTFW/8GjjiLX4VL9M27+y0qRUv6lNnNzvUFS/BcB4Bg39U79iMlUwKqlTv76kMVpHVVO/GhcOhGQBU792ieqtga1Sv9L7xteeWVK/Lm6jAbwFUr+L4H8r2bFRv+dSXFX2XVG/Q8U4fxMKUb+fNxWpMLZQv/yp8dJNYlC/WBzO/GoOUL9oHVVNEHVPvyACDqFKzU6/2ubG9IQlTr+Sy39Iv31Nv0qwOJz51Uy/ApXx7zMuTL+8eapDboZLv3ReY5eo3kq/LEMc6+I2Sr/kJ9U+HY9Jv54MjpJX50i/VvFG5pE/SL8O1v85zJdHv8a6uI0G8Ea/fp9x4UBIRr84hCo1e6BFv/Bo44i1+ES/qE2c3O9QRL9gMlUwKqlDvxoXDoRkAUO/0vvG155ZQr+K4H8r2bFBv0LFOH8TCkG//Knx0k1iQL9oHVVNEHU/v9jmxvSEJT6/SLA4nPnVPL+8eapDboY7vyxDHOviNjq/nAyOklfnOL8M1v85zJc3v3yfceFASDa/8GjjiLX4NL9gMlUwKqkzv9D7xteeWTK/QMU4fxMKMb9oHVVNEHUvv0iwOJz51Sy/KEMc6+I2Kr8I1v85zJcnv/Bo44i1+CS/0PvG155ZIr9gHVVNEHUfvyBDHOviNhq/4GjjiLX4FL9gHVVNEHUPv8Bo44i1+AS/AGnjiLX49L4AAAAAAAAAAABp44i1+PQ+AGnjiLX4BD+AHVVNEHUPPwBp44i1+BQ/QEMc6+I2Gj+AHVVNEHUfP9D7xteeWSI/8GjjiLX4JD8Q1v85zJcnPzBDHOviNio/ULA4nPnVLD9wHVVNEHUvP0jFOH8TCjE/2PvG155ZMj9oMlUwKqkzP/Bo44i1+DQ/gJ9x4UBINj8Q1v85zJc3P6AMjpJX5zg/MEMc6+I2Oj/AeapDboY7P1CwOJz51Tw/4ObG9IQlPj9oHVVNEHU/P/yp8dJNYkA/RMU4fxMKQT+M4H8r2bFBP9T7xteeWUI/HBcOhGQBQz9kMlUwKqlDP6xNnNzvUEQ/9GjjiLX4RD84hCo1e6BFP4CfceFASEY/yLq4jQbwRj8Q1v85zJdHP1jxRuaRP0g/oAyOklfnSD/oJ9U+HY9JPzBDHOviNko/dF5jl6jeSj+8eapDboZLPwSV8e8zLkw/TLA4nPnVTD+Uy39Iv31NP9zmxvSEJU4/JAIOoUrNTj9sHVVNEHVPP1oczvxqDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v5qZmZmZmbm/mpmZmZmZub+amZmZmZm5v6CZmZmZmbm/4DVfKcncub9qPSc2PaK6v3FCSsZV47u/jyhFH9yYvb9jlGdeDru/v9otwYjVIMG/m/rKZ/6Rwr9NdMDocizEvzYfV8V168W/1PfhICHKx78Pq6hUbMPJv9XNPrwx0su/ohXbgTTxzb/jR1c1kw3Qv+TsndFWJdG/qSzXRTU90r/tuRjCflLTv4Djo1yGYtS/EzAR+KRq1b9P+nspPGjWv7QMrh65WNe/vD1LhJc52L+qC/1rZAjZv5U4njLBwtm/Z2ZmZmZm2r/SshWtJvHav0lTIKrxYNu/9zDa5Naz27/DhKKuCOjbvztzDwnf+9u/mKgZjNrt27+79EdMp7zbvxfn2sAfZ9u/wmr4qU/s2r9RYtf2dkvav/BD66sMhNm/RbUPycGV2L+BJ7QvhIDXvzZzB4mBRNa/dnQjLCri1L+8pjgENFrTv+TAuXadrdG/RKIOk2C7z78SsjYSCdjLv/jSad0Gtce/UX0HjdRWw78WcE9XNYW9v72g5ARr/LO/h0nYletApL8AAAAAAAAAAIdJ2JXrQKQ/+KDkBGv8sz81cE9XNYW9P2x9B43UVsM/EdNp3Qa1xz8psjYSCdjLP0OiDpNgu88/5MC5dp2t0T/BpjgENFrTP3t0Iywq4tQ/OnMHiYFE1j+FJ7QvhIDXP0y1D8nBldg/9EPrqwyE2T9XYtf2dkvaP8Fq+KlP7No/GufawB9n2z+89EdMp7zbP5moGYza7ds/O3MPCd/72z/ChKKuCOjbP/Uw2uTWs9s/R1MgqvFg2z/SshWtJvHaP2hmZmZmZto/ljieMsHC2T+pC/1rZAjZP7w9S4SXOdg/tQyuHrlY1z9J+nspPGjWPxAwEfikatU/euOjXIZi1D/ruRjCflLTP6Ms10U1PdI/4eyd0VYl0T/cR1c1kw3QP5sV24E08c0/x80+vDHSyz/+qqhUbMPJP8n34SAhysc/Oh9XxXXrxT9NdMDocizEP5T6ymf+kcI/yy3BiNUgwT9ClGdeDru/P40oRR/cmL0/gUJKxlXjuz96PSc2PaK6P+A1XynJ3Lk/oJmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/mpmZmZmZuT+amZmZmZm5P5qZmZmZmbk/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$\\mu(x) f_1(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def mu_f1(x):\n",
+ " return mu(x) * f1(x)\n",
+ "\n",
+ "\n",
+ "sym_mu_f1 = Piecewise(\n",
+ " (-mu_k, x < -eps_v),\n",
+ " (-mu_f1(-x), x <= 0),\n",
+ " (mu_f1(x), x <= eps_v),\n",
+ " (mu_k, x > eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(r\"\\mu(x) f_1(x)\"), sym_mu_f1))\n",
+ "\n",
+ "plot(sym_mu_f1, title=r'$\\mu(x) f_1(x)$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The resulting force mollifier is a polynomail in $x$, so we can use sympy to get the exact integral."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 237,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{\\epsilon_{v} \\left(17 \\mu_{k} - 7 \\mu_{s}\\right)}{30}$"
+ ],
+ "text/plain": [
+ "\\epsilon_v*(17*\\mu_k - 7*\\mu_s)/30"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\int \\mu(x) f_1(x) \\mathrm{d}x = \\begin{cases} - \\epsilon_{v} \\mu_{k} + \\frac{2 \\epsilon_{v} \\mu_{s}}{3} + \\frac{\\epsilon_{v} \\left(- 7 \\mu_{k} + 7 \\mu_{s}\\right)}{5} + \\frac{\\epsilon_{v} \\left(\\mu_{k} - \\mu_{s}\\right)}{3} + \\frac{\\epsilon_{v} \\left(3 \\mu_{k} - 3 \\mu_{s}\\right)}{2} + \\mu_{k} \\left|{x}\\right| & \\text{for}\\: \\epsilon_{v} \\leq \\left|{x}\\right| \\\\\\frac{\\mu_{s} \\left|{x}\\right|^{2}}{\\epsilon_{v}} - \\frac{\\mu_{s} \\left|{x}\\right|^{3}}{3 \\epsilon_{v}^{2}} + \\frac{3 \\mu_{k} \\left|{x}\\right|^{4}}{2 \\epsilon_{v}^{3}} - \\frac{3 \\mu_{s} \\left|{x}\\right|^{4}}{2 \\epsilon_{v}^{3}} - \\frac{7 \\mu_{k} \\left|{x}\\right|^{5}}{5 \\epsilon_{v}^{4}} + \\frac{7 \\mu_{s} \\left|{x}\\right|^{5}}{5 \\epsilon_{v}^{4}} + \\frac{\\mu_{k} \\left|{x}\\right|^{6}}{3 \\epsilon_{v}^{5}} - \\frac{\\mu_{s} \\left|{x}\\right|^{6}}{3 \\epsilon_{v}^{5}} & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(\\int \\mu(x) f_1(x) \\mathrm{d}x, Piecewise((-\\epsilon_v*\\mu_k + 2*\\epsilon_v*\\mu_s/3 + \\epsilon_v*(-7*\\mu_k + 7*\\mu_s)/5 + \\epsilon_v*(\\mu_k - \\mu_s)/3 + \\epsilon_v*(3*\\mu_k - 3*\\mu_s)/2 + \\mu_k*Abs(x), \\epsilon_v <= Abs(x)), (\\mu_s*Abs(x)**2/\\epsilon_v - \\mu_s*Abs(x)**3/(3*\\epsilon_v**2) + 3*\\mu_k*Abs(x)**4/(2*\\epsilon_v**3) - 3*\\mu_s*Abs(x)**4/(2*\\epsilon_v**3) - 7*\\mu_k*Abs(x)**5/(5*\\epsilon_v**4) + 7*\\mu_s*Abs(x)**5/(5*\\epsilon_v**4) + \\mu_k*Abs(x)**6/(3*\\epsilon_v**5) - \\mu_s*Abs(x)**6/(3*\\epsilon_v**5), True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL8q499nXDhgv1gczvxqDmC/Dat4I/PIX79pHVVNEHVfv8WPMXctIV+/IQIOoUrNXr9+dOrKZ3lev9rmxvSEJV6/NlmjHqLRXb+Sy39Iv31dv+89XHLcKV2/S7A4nPnVXL+nIhXGFoJcvwOV8e8zLly/YAfOGVHaW7+8eapDboZbvxjshm2LMlu/dF5jl6jeWr/Q0D/BxYpavy1DHOviNlq/ibX4FADjWb/lJ9U+HY9Zv0KasWg6O1m/ngyOklfnWL/6fmq8dJNYv1bxRuaRP1i/smMjEK/rV78O1v85zJdXv2tI3GPpQ1e/x7q4jQbwVr8kLZW3I5xWv4CfceFASFa/3BFOC170Vb84hCo1e6BVv5T2Bl+YTFW/8GjjiLX4VL9M27+y0qRUv6lNnNzvUFS/BcB4Bg39U79iMlUwKqlTv76kMVpHVVO/GhcOhGQBU792ieqtga1Sv9L7xteeWVK/Lm6jAbwFUr+L4H8r2bFRv+dSXFX2XVG/Q8U4fxMKUb+fNxWpMLZQv/yp8dJNYlC/WBzO/GoOUL9oHVVNEHVPvyACDqFKzU6/2ubG9IQlTr+Sy39Iv31Nv0qwOJz51Uy/ApXx7zMuTL+8eapDboZLv3ReY5eo3kq/LEMc6+I2Sr/kJ9U+HY9Jv54MjpJX50i/VvFG5pE/SL8O1v85zJdHv8a6uI0G8Ea/fp9x4UBIRr84hCo1e6BFv/Bo44i1+ES/qE2c3O9QRL9gMlUwKqlDvxoXDoRkAUO/0vvG155ZQr+K4H8r2bFBv0LFOH8TCkG//Knx0k1iQL9oHVVNEHU/v9jmxvSEJT6/SLA4nPnVPL+8eapDboY7vyxDHOviNjq/nAyOklfnOL8M1v85zJc3v3yfceFASDa/8GjjiLX4NL9gMlUwKqkzv9D7xteeWTK/QMU4fxMKMb9oHVVNEHUvv0iwOJz51Sy/KEMc6+I2Kr8I1v85zJcnv/Bo44i1+CS/0PvG155ZIr9gHVVNEHUfvyBDHOviNhq/4GjjiLX4FL9gHVVNEHUPv8Bo44i1+AS/AGnjiLX49L4AAAAAAAAAAABp44i1+PQ+AGnjiLX4BD+AHVVNEHUPPwBp44i1+BQ/QEMc6+I2Gj+AHVVNEHUfP9D7xteeWSI/8GjjiLX4JD8Q1v85zJcnPzBDHOviNio/ULA4nPnVLD9wHVVNEHUvP0jFOH8TCjE/2PvG155ZMj9oMlUwKqkzP/Bo44i1+DQ/gJ9x4UBINj8Q1v85zJc3P6AMjpJX5zg/MEMc6+I2Oj/AeapDboY7P1CwOJz51Tw/4ObG9IQlPj9oHVVNEHU/P/yp8dJNYkA/RMU4fxMKQT+M4H8r2bFBP9T7xteeWUI/HBcOhGQBQz9kMlUwKqlDP6xNnNzvUEQ/9GjjiLX4RD84hCo1e6BFP4CfceFASEY/yLq4jQbwRj8Q1v85zJdHP1jxRuaRP0g/oAyOklfnSD/oJ9U+HY9JPzBDHOviNko/dF5jl6jeSj+8eapDboZLPwSV8e8zLkw/TLA4nPnVTD+Uy39Iv31NP9zmxvSEJU4/JAIOoUrNTj9sHVVNEHVPP1oczvxqDlA//Knx0k1iUD+gNxWpMLZQP0TFOH8TClE/6FJcVfZdUT+M4H8r2bFRPzBuowG8BVI/1PvG155ZUj94ieqtga1SPxoXDoRkAVM/vqQxWkdVUz9iMlUwKqlTPwbAeAYN/VM/qk2c3O9QVD9O27+y0qRUP/Jo44i1+FQ/lvYGX5hMVT86hCo1e6BVP9wRTgte9FU/gJ9x4UBIVj8kLZW3I5xWP8i6uI0G8FY/bEjcY+lDVz8Q1v85zJdXP7RjIxCv61c/WPFG5pE/WD/6fmq8dJNYP54MjpJX51g/QpqxaDo7WT/mJ9U+HY9ZP4q1+BQA41k/LkMc6+I2Wj/S0D/BxYpaP3ZeY5eo3lo/GuyGbYsyWz+8eapDboZbP2AHzhlR2ls/BJXx7zMuXD+oIhXGFoJcP0ywOJz51Vw/8D1cctwpXT+Uy39Iv31dPzhZox6i0V0/3ObG9IQlXj9+dOrKZ3lePyICDqFKzV4/xo8xdy0hXz9sHVVNEHVfPwyreCPzyF8/WBzO/GoOYD8q499nXDhgP/yp8dJNYmA/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "Lq7LA2uvOD9UQooU3Y04P3jWSCVPbDg/nmoHNsFKOD/C/sVGMyk4P+eShFelBzg/DCdDaBfmNz8xuwF5icQ3P1ZPwIn7ojc/euN+mm2BNz+gdz2r3183P8QL/LtRPjc/6p+6zMMcNz8ONHndNfs2PzPIN+6n2TY/WFz2/hm4Nj998LQPjJY2P6KEcyD+dDY/xhgyMXBTNj/srPBB4jE2PxBBr1JUEDY/NtVtY8buNT9aaSx0OM01P4D96oSqqzU/pJGplRyKNT/JJWimjmg1P+65JrcARzU/Ek7lx3IlNT844qPY5AM1P1x2YulW4jQ/ggoh+sjAND+mnt8KO580P8wynhutfTQ/8MZcLB9cND8VWxs9kTo0Pzrv2U0DGTQ/XoOYXnX3Mz+EF1dv59UzP6irFYBZtDM/zT/UkMuSMz/y05KhPXEzPxdoUbKvTzM/PPwPwyEuMz9hkM7TkwwzP4YkjeQF6zI/qrhL9XfJMj/QTAoG6qcyP/TgyBZchjI/GXWHJ85kMj8+CUY4QEMyP2OdBEmyITI/qYmG4wYAMj8TzVherd0xP5t+U0H/uTE/IAhSJ1+UMT+MBjJ8OWwxPwQOnysFQTE/plGrQEQSMT9ZLzV2hN8wP1KfGbhfqDA/aoczlXxsMD8+8ieijiswP5RU/pmtyi8/K3EVQ0czLz/QjjD9opAuP4bR4neV4i0/KIFQgg8pLT+2JWIaHmQsPwZrkV3qkys/ZMxPW7m4Kj+TBwfJ69IpP45Xs5f94ig/1XYXa4XpJz+jaYryM+cmP2kPXyPT3CU/OHzlVEXLJD+VGQY+hLMjPyiPdtSfliI/8HKIDb11IT8XwZGAFFIgP9Q23tXhWR4/KqE/K1sPHD+9Vfc4accZPxMNBS76hBc/9F+LRRJLFT+bEtaexxwTP5zvk9g9/RA/cmOI3kLfDT9V+q69Re4JPyvQCA/jLQY/6CSfq2ykAj9cNB23OLD+Pi5Ay2UQnvg+tADU6Swe8z5zIRb3w3fsPsdFZ8rgA+Q+8D7Zv6vp2T7VBSnt0nTNPi+F3fQGbro+qeKdJn6mmj4AAAAAAAAAAKninSZ+ppo+zoXd9AZuuj4PBint0nTNPj0/2b+r6dk++EVnyuAD5D6rIRb3w3fsPrQA1OksHvM+LkDLZRCe+D5wNB23OLD+PvIkn6tspAI/N9AID+MtBj9g+q69Re4JP4ljiN5C3w0/qu+T2D39ED+pEtaexxwTP/Rfi0USSxU/Gg0FLvqEFz/EVfc4accZPzKhPytbDxw/2jbe1eFZHj8YwZGAFFIgP/dyiA29dSE/L4921J+WIj+VGQY+hLMjPzh85VRFyyQ/bA9fI9PcJT+paYryM+cmP9l2F2uF6Sc/jVezl/3iKD+ZBwfJ69IpP2rMT1u5uCo/B2uRXeqTKz+2JWIaHmQsPyqBUIIPKS0/jNHid5XiLT/YjjD9opAuPypxFUNHMy8/mFT+ma3KLz9A8ieijiswP2uHM5V8bDA/Up8ZuF+oMD9ZLzV2hN8wP6NRq0BEEjE/BA6fKwVBMT+KBjJ8OWwxPyEIUidflDE/mH5TQf+5MT8TzVherd0xP62JhuMGADI/Y50ESbIhMj8+CUY4QEMyPxp1hyfOZDI/9eDIFlyGMj/QTAoG6qcyP6u4S/V3yTI/hiSN5AXrMj9ikM7TkwwzPzz8D8MhLjM/F2hRsq9PMz/y05KhPXEzP84/1JDLkjM/qasVgFm0Mz+EF1dv59UzP2CDmF519zM/Ou/ZTQMZND8WWxs9kTo0P/DGXCwfXDQ/zDKeG619ND+mnt8KO580P4IKIfrIwDQ/XXZi6VbiND844qPY5AM1PxRO5cdyJTU/7rkmtwBHNT/JJWimjmg1P6SRqZUcijU/gP3qhKqrNT9aaSx0OM01PzbVbWPG7jU/EUGvUlQQNj/srPBB4jE2P8gYMjFwUzY/ooRzIP50Nj998LQPjJY2P1hc9v4ZuDY/NMg37qfZNj8ONHndNfs2P+qfuszDHDc/xQv8u1E+Nz+gdz2r3183P3zjfpptgTc/Vk/AifuiNz8xuwF5icQ3PwwnQ2gX5jc/6JKEV6UHOD/E/sVGMyk4P51qBzbBSjg/eNZIJU9sOD9UQooU3Y04Py6uywNrrzg/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$\\int \\mu(x) f_1(x) \\mathrm{d}x$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "mu_f0 = integrate(mu_f1(x), x)\n",
+ "\n",
+ "# Constant of integration s.t. mu_f0(eps_v) = mu_k * eps_v\n",
+ "c = mu_k * eps_v - mu_f0.subs(x, eps_v)\n",
+ "display(c.simplify())\n",
+ "mu_f0 = mu_f0\n",
+ "\n",
+ "mu_f0.simplify()\n",
+ "\n",
+ "sym_mu_f0 = Piecewise(\n",
+ " (mu_k * abs(x) - c, abs(x) >= eps_v),\n",
+ " (mu_f0.expand().subs(x, abs(x)), abs(x) < eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(r\"\\int \\mu(x) f_1(x) \\mathrm{d}x\"), sym_mu_f0))\n",
+ "\n",
+ "plot(sym_mu_f0, title=r'$\\int \\mu(x) f_1(x) \\mathrm{d}x$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 217,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\operatorname{Poly}{\\left( \\frac{\\mu_{k} - \\mu_{s}}{3 \\epsilon_{v}^{5}} x^{6} + \\frac{- 7 \\mu_{k} + 7 \\mu_{s}}{5 \\epsilon_{v}^{4}} x^{5} + \\frac{3 \\mu_{k} - 3 \\mu_{s}}{2 \\epsilon_{v}^{3}} x^{4} - \\frac{\\mu_{s}}{3 \\epsilon_{v}^{2}} x^{3} + \\frac{\\mu_{s}}{\\epsilon_{v}} x^{2}, x, domain=\\mathbb{Z}\\left(\\mu_{k}, \\mu_{s}, \\epsilon_{v}\\right) \\right)}$"
+ ],
+ "text/plain": [
+ "Poly((\\mu_k - \\mu_s)/(3*\\epsilon_v**5)*x**6 + (-7*\\mu_k + 7*\\mu_s)/(5*\\epsilon_v**4)*x**5 + (3*\\mu_k - 3*\\mu_s)/(2*\\epsilon_v**3)*x**4 - \\mu_s/(3*\\epsilon_v**2)*x**3 + \\mu_s/\\epsilon_v*x**2, x, domain='ZZ(\\mu_k,\\mu_s,\\epsilon_v)')"
+ ]
+ },
+ "execution_count": 217,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "poly(mu_f0, x)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 218,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\operatorname{Poly}{\\left( \\frac{\\mu_{k} - \\mu_{s}}{3 \\epsilon_{v}^{5}} x^{6} + \\frac{- 7 \\mu_{k} + 7 \\mu_{s}}{5 \\epsilon_{v}^{4}} x^{5} + \\frac{3 \\mu_{k} - 3 \\mu_{s}}{2 \\epsilon_{v}^{3}} x^{4} - \\frac{\\mu_{s}}{3 \\epsilon_{v}^{2}} x^{3} + \\frac{\\mu_{s}}{\\epsilon_{v}} x^{2} + \\frac{17 \\epsilon_{v} \\mu_{k}}{30} - \\frac{7 \\epsilon_{v} \\mu_{s}}{30}, x, domain=\\mathbb{Z}\\left(\\mu_{k}, \\mu_{s}, \\epsilon_{v}\\right) \\right)}$"
+ ],
+ "text/plain": [
+ "Poly((\\mu_k - \\mu_s)/(3*\\epsilon_v**5)*x**6 + (-7*\\mu_k + 7*\\mu_s)/(5*\\epsilon_v**4)*x**5 + (3*\\mu_k - 3*\\mu_s)/(2*\\epsilon_v**3)*x**4 - \\mu_s/(3*\\epsilon_v**2)*x**3 + \\mu_s/\\epsilon_v*x**2 + 17*\\epsilon_v*\\mu_k/30 - 7*\\epsilon_v*\\mu_s/30, x, domain='ZZ(\\mu_k,\\mu_s,\\epsilon_v)')"
+ ]
+ },
+ "execution_count": 218,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "def mu_f0(x, eps_v, mu_s, mu_k):\n",
+ " \"\"\"Integral of mu(x) f1(x).\"\"\"\n",
+ " d_mu = mu_k - mu_s\n",
+ " x_over_eps_v = x / eps_v\n",
+ " return x * x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * d_mu / 3 - 1.4 * d_mu\n",
+ " ) + 1.5 * d_mu\n",
+ " ) - mu_s/3\n",
+ " ) + mu_s) + eps_v * (17 * mu_k - 7 * mu_s) / 30\n",
+ "\n",
+ "\n",
+ "poly(nsimplify(mu_f0(x, eps_v, mu_s, mu_k)), x)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Checking my work"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 219,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{\\epsilon_{v}^{6} \\left(17 \\mu_{k} - 7 \\mu_{s}\\right) + 10 x^{2} \\left(3 \\epsilon_{v}^{4} \\mu_{s} - x \\left(\\epsilon_{v}^{3} \\mu_{s} - x \\left(4.5 \\epsilon_{v}^{2} \\left(\\mu_{k} - \\mu_{s}\\right) + x \\left(4.2 \\epsilon_{v} \\left(- \\mu_{k} + \\mu_{s}\\right) + x \\left(\\mu_{k} - \\mu_{s}\\right)\\right)\\right)\\right)\\right)}{30 \\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "(\\epsilon_v**6*(17*\\mu_k - 7*\\mu_s) + 10*x**2*(3*\\epsilon_v**4*\\mu_s - x*(\\epsilon_v**3*\\mu_s - x*(4.5*\\epsilon_v**2*(\\mu_k - \\mu_s) + x*(4.2*\\epsilon_v*(-\\mu_k + \\mu_s) + x*(\\mu_k - \\mu_s))))))/(30*\\epsilon_v**5)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{\\epsilon_{v} \\mu_{k}}{3} + \\frac{\\mu_{k} x^{2}}{\\epsilon_{v}} - \\frac{\\mu_{k} x^{3}}{3 \\epsilon_{v}^{2}}$"
+ ],
+ "text/plain": [
+ "\\epsilon_v*\\mu_k/3 + \\mu_k*x**2/\\epsilon_v - \\mu_k*x**3/(3*\\epsilon_v**2)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "f = mu_f0(x, eps_v, mu_s, mu_k).simplify()\n",
+ "display(f)\n",
+ "display(f.subs(mu_s, mu_k).simplify().expand())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 220,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{6 \\epsilon_{v}^{4} \\mu_{s} - 3 \\epsilon_{v}^{3} \\mu_{s} x + 18 \\epsilon_{v}^{2} \\mu_{k} x^{2} - 18 \\epsilon_{v}^{2} \\mu_{s} x^{2} - 21 \\epsilon_{v} \\mu_{k} x^{3} + 21 \\epsilon_{v} \\mu_{s} x^{3} + 6 \\mu_{k} x^{4} - 6 \\mu_{s} x^{4}}{3 \\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "(6*\\epsilon_v**4*\\mu_s - 3*\\epsilon_v**3*\\mu_s*x + 18*\\epsilon_v**2*\\mu_k*x**2 - 18*\\epsilon_v**2*\\mu_s*x**2 - 21*\\epsilon_v*\\mu_k*x**3 + 21*\\epsilon_v*\\mu_s*x**3 + 6*\\mu_k*x**4 - 6*\\mu_s*x**4)/(3*\\epsilon_v**5)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{2 \\mu_{k}}{\\epsilon_{v}} - \\frac{\\mu_{k} x}{\\epsilon_{v}^{2}}$"
+ ],
+ "text/plain": [
+ "2*\\mu_k/\\epsilon_v - \\mu_k*x/\\epsilon_v**2"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "f1_over_x = nsimplify((f.diff(x) / x).simplify())\n",
+ "display(f1_over_x)\n",
+ "display(f1_over_x.subs(mu_s, mu_k).simplify().expand())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 221,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{- 30 \\epsilon_{v}^{3} \\mu_{s} + 360 \\epsilon_{v}^{2} \\mu_{k} x - 360 \\epsilon_{v}^{2} \\mu_{s} x - 630 \\epsilon_{v} \\mu_{k} x^{2} + 630 \\epsilon_{v} \\mu_{s} x^{2} + 240 \\mu_{k} x^{3} - 240 \\mu_{s} x^{3}}{30 \\epsilon_{v}^{5} x}$"
+ ],
+ "text/plain": [
+ "(-30*\\epsilon_v**3*\\mu_s + 360*\\epsilon_v**2*\\mu_k*x - 360*\\epsilon_v**2*\\mu_s*x - 630*\\epsilon_v*\\mu_k*x**2 + 630*\\epsilon_v*\\mu_s*x**2 + 240*\\mu_k*x**3 - 240*\\mu_s*x**3)/(30*\\epsilon_v**5*x)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle - \\frac{\\mu_{k}}{\\epsilon_{v}^{2} x}$"
+ ],
+ "text/plain": [
+ "-\\mu_k/(\\epsilon_v**2*x)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "smooth_friction_f2_x_minus_f1_over_x3 = nsimplify(\n",
+ " ((f.diff(x).diff(x) * x - f.diff(x)) / x**3).simplify())\n",
+ "display(smooth_friction_f2_x_minus_f1_over_x3)\n",
+ "display(smooth_friction_f2_x_minus_f1_over_x3.subs(\n",
+ " mu_s, mu_k).simplify().expand())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Safe Division"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 222,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{- \\epsilon_{v}^{3} \\mu_{s} - x^{2} \\left(3 \\epsilon_{v} - 2 x\\right) \\left(\\mu_{k} - \\mu_{s}\\right) + 6 x \\left(\\epsilon_{v} - x\\right) \\left(2 \\epsilon_{v} - x\\right) \\left(\\mu_{k} - \\mu_{s}\\right)}{\\epsilon_{v}^{5} x}$"
+ ],
+ "text/plain": [
+ "(-\\epsilon_v**3*\\mu_s - x**2*(3*\\epsilon_v - 2*x)*(\\mu_k - \\mu_s) + 6*x*(\\epsilon_v - x)*(2*\\epsilon_v - x)*(\\mu_k - \\mu_s))/(\\epsilon_v**5*x)"
+ ]
+ },
+ "execution_count": 222,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "smooth_mu_f2_x_minus_f1_over_x3 = \\\n",
+ " (((mu(x) * f1(x)).diff(x) * x - mu(x) * f1(x)) / x**3).simplify()\n",
+ "smooth_mu_f2_x_minus_f1_over_x3"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 223,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle - \\frac{\\mu_{s}}{\\epsilon_{v}^{2} x} + \\frac{12 \\mu_{k}}{\\epsilon_{v}^{3}} - \\frac{12 \\mu_{s}}{\\epsilon_{v}^{3}} - \\frac{21 \\mu_{k} x}{\\epsilon_{v}^{4}} + \\frac{21 \\mu_{s} x}{\\epsilon_{v}^{4}} + \\frac{8 \\mu_{k} x^{2}}{\\epsilon_{v}^{5}} - \\frac{8 \\mu_{s} x^{2}}{\\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "-\\mu_s/(\\epsilon_v**2*x) + 12*\\mu_k/\\epsilon_v**3 - 12*\\mu_s/\\epsilon_v**3 - 21*\\mu_k*x/\\epsilon_v**4 + 21*\\mu_s*x/\\epsilon_v**4 + 8*\\mu_k*x**2/\\epsilon_v**5 - 8*\\mu_s*x**2/\\epsilon_v**5"
+ ]
+ },
+ "execution_count": 223,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "smooth_mu_f2_x_minus_f1_over_x3.expand()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 224,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle 8 \\mu_{k} x^{2} z^{5} - 21 \\mu_{k} x z^{4} + 12 \\mu_{k} z^{3} - 8 \\mu_{s} x^{2} z^{5} + 21 \\mu_{s} x z^{4} - 12 \\mu_{s} z^{3} - \\frac{\\mu_{s} z^{2}}{x}$"
+ ],
+ "text/plain": [
+ "8*\\mu_k*x**2*z**5 - 21*\\mu_k*x*z**4 + 12*\\mu_k*z**3 - 8*\\mu_s*x**2*z**5 + 21*\\mu_s*x*z**4 - 12*\\mu_s*z**3 - \\mu_s*z**2/x"
+ ]
+ },
+ "execution_count": 224,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "smooth_mu_f2_x_minus_f1_over_x3.expand().subs({1/eps_v: Symbol('z')})"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Tangential Adhesion"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 225,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle a_{0}(x) = \\begin{cases} 0 & \\text{for}\\: x \\leq 0 \\\\\\frac{x^{2}}{\\epsilon_{v}} - \\frac{x^{3}}{3 \\epsilon_{v}^{2}} & \\text{for}\\: \\left|{x}\\right| < 2 \\epsilon_{v} \\\\\\frac{4 \\epsilon_{v}}{3} & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(a_{0}(x), Piecewise((0, x <= 0), (x**2/\\epsilon_v - x**3/(3*\\epsilon_v**2), Abs(x) < 2*\\epsilon_v), (4*\\epsilon_v/3, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iYL92cRsN4C1gv99xio7k8l+/0gDeAgmKX7/FjzF3LSFfv7gehetRuF6/rK3YX3ZPXr+fPCzUmuZdv5LLf0i/fV2/hlrTvOMUXb956SYxCKxcv2x4eqUsQ1y/YAfOGVHaW79TliGOdXFbv0YldQKaCFu/OrTIdr6fWr8tQxzr4jZavyDSb18Hzlm/E2HD0ytlWb8H8BZIUPxYv/p+arx0k1i/7Q2+MJkqWL/gnBGlvcFXv9QrZRniWFe/x7q4jQbwVr+6SQwCK4dWv67YX3ZPHla/oWez6nO1Vb+U9gZfmExVv4iFWtO841S/exSuR+F6VL9uowG8BRJUv2IyVTAqqVO/VcGopE5AU79IUPwYc9dSvzzfT42XblK/Lm6jAbwFUr8i/fZ14JxRvxaMSuoENFG/CBueXinLUL/8qfHSTWJQv95xio7k8k+/xI8xdy0hT7+srdhfdk9Ov5LLf0i/fU2/eOkmMQisTL9gB84ZUdpLv0YldQKaCEu/LEMc6+I2Sr8UYcPTK2VJv/p+arx0k0i/4JwRpb3BR7/GuriNBvBGv67YX3ZPHka/lPYGX5hMRb96FK5H4XpEv2IyVTAqqUO/SFD8GHPXQr8ubqMBvAVCvxaMSuoENEG//Knx0k1iQL/EjzF3LSE/v5DLf0i/fT2/YAfOGVHaO78sQxzr4jY6v/h+arx0kzi/yLq4jQbwNr+U9gZfmEw1v2AyVTAqqTO/LG6jAbwFMr/8qfHSTWIwv5DLf0i/fS2/KEMc6+I2Kr/IuriNBvAmv2AyVTAqqSO/+Knx0k1iIL8wQxzr4jYav2AyVTAqqRO/IEMc6+I2Cr8AQxzr4jb6vgAAAAAAAAAAAEMc6+I2+j5AQxzr4jYKP2AyVTAqqRM/QEMc6+I2Gj8AqvHSTWIgP2AyVTAqqSM/0Lq4jQbwJj8wQxzr4jYqP5DLf0i/fS0/AKrx0k1iMD8wbqMBvAUyP2AyVTAqqTM/mPYGX5hMNT/IuriNBvA2P/h+arx0kzg/MEMc6+I2Oj9gB84ZUdo7P5DLf0i/fT0/yI8xdy0hPz/8qfHSTWJAPxiMSuoENEE/MG6jAbwFQj9IUPwYc9dCP2QyVTAqqUM/fBSuR+F6RD+U9gZfmExFP7DYX3ZPHkY/yLq4jQbwRj/gnBGlvcFHP/x+arx0k0g/FGHD0ytlST8sQxzr4jZKP0gldQKaCEs/YAfOGVHaSz946SYxCKxMP5TLf0i/fU0/rK3YX3ZPTj/EjzF3LSFPP+Bxio7k8k8//Knx0k1iUD8KG55eKctQPxaMSuoENFE/Iv32deCcUT8wbqMBvAVSPzzfT42XblI/SFD8GHPXUj9WwaikTkBTP2IyVTAqqVM/bqMBvAUSVD98FK5H4XpUP4iFWtO841Q/lPYGX5hMVT+iZ7Pqc7VVP67YX3ZPHlY/ukkMAiuHVj/IuriNBvBWP9QrZRniWFc/4pwRpb3BVz/uDb4wmSpYP/p+arx0k1g/CPAWSFD8WD8UYcPTK2VZPyDSb18Hzlk/LkMc6+I2Wj86tMh2vp9aP0YldQKaCFs/VJYhjnVxWz9gB84ZUdpbP2x4eqUsQ1w/eukmMQisXD+GWtO84xRdP5LLf0i/fV0/oDws1JrmXT+srdhfdk9eP7gehetRuF4/xo8xdy0hXz/UAN4CCYpfP+Bxio7k8l8/dnEbDeAtYD/8qfHSTWJgP4Lix5i7lmA/CBueXinLYD+QU3Qkl/9gPxaMSuoENGE/nMQgsHJoYT8i/fZ14JxhP6g1zTtO0WE/MG6jAbwFYj+2pnnHKTpiPzzfT42XbmI/whcmUwWjYj9IUPwYc9diP86I0t7gC2M/VsGopE5AYz/c+X5qvHRjP2IyVTAqqWM/6Gor9pfdYz9uowG8BRJkP/Tb14FzRmQ/fBSuR+F6ZD8CTYQNT69kP4iFWtO842Q/Dr4wmSoYZT+U9gZfmExlPxwv3SQGgWU/omez6nO1ZT8ooImw4ellP67YX3ZPHmY/NBE2PL1SZj+6SQwCK4dmP0KC4seYu2Y/yLq4jQbwZj9O845TdCRnP9QrZRniWGc/WmQ730+NZz/gnBGlvcFnP2jV52or9mc/7g2+MJkqaD90RpT2Bl9oP/p+arx0k2g/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAhoPhSfjLpD7Hnt8KO5/EPrFwWYXNANc+ftTbjMBF5D7chaQpFWfvPlgLs9DOafY+rOg/fkk9/j7VP9SQy5IDP706jh8qjQg/srdlrY4JDj+BRDm+4wESP6RAZmfROxU/l7lFc3ewGD99mOMCvV0cP1LjpZvEICA/sxbFmOEsIj/uWlUJKVIkP66k3H2OjyY/n+jghgXkKD9PG+i0gU4rP3MxeJj2zS0/y48L4aswMD8xbSVhzIMxP0OrzJTW3zI/R0REREREND+RMs83j7A1P3pwsDcxJDc/RfgqDKSeOD9IxIF9YR86P93O91PjpTs/RRLQV6MxPT/YiE1RG8I+P3aWWYRiK0A/Y/whI433QD9fcyFpSsVBP5R4eTpXlEI/JIlLe3BkQz88IrkPUzVEPwbB49u7BkU/puLsw2fYRT9IBParE6pGPw6jIHh8e0c/JjyODF9MSD+9TGBNeBxJP+1RuB6F60k/6ci3ZEK5Sj/ZLoADbYVLP98AM9/BT0w/Krzx2/0XTT/h3d3d3d1NPybjGMkeoU4/KEnEgX1hTz+KxgB2Ww9QPwIW+fVDbFA/lNHbslbHUD/Wt7keciBRP1qHo6t0d1E/tf6pyzzMUT953N3wqB5SPzzfT42XblI/lMUQE+e7Uj8QTjH0dQZTP0g3wqIiTlM/0T/UkMuSUz88JngwT9RTPx6pvvOLElQ/Doe4TGBNVD+bfnatqoRUP15OCYhJuFQ/6bSBThvoVD/QcPBy/hNVP6dAZmfRO1U/A+PznXJfVT94FqqIwH5VP5yZmZmZmVU//SrTQtyvVT81iWf2ZsFVP9VyZyYYzlU/dqbjRM7VVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/peLsw2fYVT+l4uzDZ9hVP6Xi7MNn2FU/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$a_0(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def a0_part1(x):\n",
+ " \"\"\"Smooth adhesion mollifier\"\"\"\n",
+ " return x * x / eps_v * (1 - x / (3 * eps_v))\n",
+ "\n",
+ "\n",
+ "def a0_part2(x):\n",
+ " \"\"\"Smooth adhesion mollifier\"\"\"\n",
+ " return 4 * eps_v / 3\n",
+ "\n",
+ "\n",
+ "sym_a0 = Piecewise(\n",
+ " (0, x <= 0),\n",
+ " (a0_part1(x), abs(x) < 2 * eps_v),\n",
+ " (a0_part2(x), abs(x) >= 2 * eps_v),\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(\"a_{0}(x)\"), sym_a0.expand()))\n",
+ "\n",
+ "plot(sym_a0, title=r'$a_0(x)$', x_offset=-1, max_x=3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 226,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle a_{1}(x) = \\begin{cases} 0 & \\text{for}\\: x \\leq 0 \\\\\\frac{2 x}{\\epsilon_{v}} - \\frac{x^{2}}{\\epsilon_{v}^{2}} & \\text{for}\\: x < 2 \\epsilon_{v} \\\\0 & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(a_{1}(x), Piecewise((0, x <= 0), (2*x/\\epsilon_v - x**2/\\epsilon_v**2, x < 2*\\epsilon_v), (0, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iUL9YHM78ag5Qv2kdVU0QdU+/IQIOoUrNTr/a5sb0hCVOv5LLf0i/fU2/S7A4nPnVTL8DlfHvMy5Mv7x5qkNuhku/dF5jl6jeSr8tQxzr4jZKv+Un1T4dj0m/ngyOklfnSL9W8UbmkT9Ivw7W/znMl0e/x7q4jQbwRr+An3HhQEhGvziEKjV7oEW/8GjjiLX4RL+pTZzc71BEv2IyVTAqqUO/GhcOhGQBQ7/S+8bXnllCv4vgfyvZsUG/Q8U4fxMKQb/8qfHSTWJAv2gdVU0QdT+/2ubG9IQlPr9KsDic+dU8v7x5qkNuhju/LEMc6+I2Or+eDI6SV+c4vw7W/znMlze/fp9x4UBINr/waOOItfg0v2AyVTAqqTO/0vvG155ZMr9CxTh/Ewoxv2gdVU0QdS+/SLA4nPnVLL8sQxzr4jYqvwzW/znMlye/8GjjiLX4JL/Q+8bXnlkiv2gdVU0QdR+/KEMc6+I2Gr/waOOItfgUv2AdVU0QdQ+/4GjjiLX4BL/AaOOItfj0vgAAAAAAAAAAAGnjiLX49D4AaeOItfgEP4AdVU0QdQ8/8GjjiLX4FD8wQxzr4jYaP3AdVU0QdR8/2PvG155ZIj/waOOItfgkPxDW/znMlyc/MEMc6+I2Kj9QsDic+dUsP2gdVU0QdS8/RMU4fxMKMT/U+8bXnlkyP2QyVTAqqTM/9GjjiLX4ND+An3HhQEg2PxDW/znMlzc/oAyOklfnOD8wQxzr4jY6P7x5qkNuhjs/TLA4nPnVPD/c5sb0hCU+P2wdVU0QdT8//Knx0k1iQD9ExTh/EwpBP4zgfyvZsUE/1PvG155ZQj8aFw6EZAFDP2IyVTAqqUM/qk2c3O9QRD/yaOOItfhEPzqEKjV7oEU/gJ9x4UBIRj/IuriNBvBGPxDW/znMl0c/WPFG5pE/SD+eDI6SV+dIP+Yn1T4dj0k/LkMc6+I2Sj92XmOXqN5KP7x5qkNuhks/BJXx7zMuTD9MsDic+dVMP5TLf0i/fU0/3ObG9IQlTj8iAg6hSs1OP2wdVU0QdU8/WBzO/GoOUD/8qfHSTWJQP6A3FakwtlA/RMU4fxMKUT/oUlxV9l1RP4zgfyvZsVE/MG6jAbwFUj/U+8bXnllSP3aJ6q2BrVI/GhcOhGQBUz++pDFaR1VTP2IyVTAqqVM/BsB4Bg39Uz+qTZzc71BUP07bv7LSpFQ/8mjjiLX4VD+W9gZfmExVPziEKjV7oFU/3BFOC170VT+An3HhQEhWPyQtlbcjnFY/yLq4jQbwVj9sSNxj6UNXPxDW/znMl1c/tGMjEK/rVz9W8UbmkT9YP/p+arx0k1g/ngyOklfnWD9CmrFoOjtZP+Yn1T4dj1k/irX4FADjWT8uQxzr4jZaP9LQP8HFilo/dl5jl6jeWj8Y7IZtizJbP7x5qkNuhls/YAfOGVHaWz8ElfHvMy5cP6giFcYWglw/TLA4nPnVXD/wPVxy3CldP5TLf0i/fV0/NlmjHqLRXT/a5sb0hCVeP3506spneV4/IgIOoUrNXj/GjzF3LSFfP2odVU0QdV8/Dqt4I/PIXz9ZHM78ag5gPyvj32dcOGA//Knx0k1iYD/OcAM+P4xgP6A3FakwtmA/cv4mFCLgYD9ExTh/EwphPxaMSuoENGE/6FJcVfZdYT+6GW7A54dhP4vgfyvZsWE/XaeRlsrbYT8vbqMBvAViPwE1tWytL2I/0/vG155ZYj+lwthCkINiP3eJ6q2BrWI/SVD8GHPXYj8bFw6EZAFjP+zdH+9VK2M/vqQxWkdVYz+Qa0PFOH9jP2IyVTAqqWM/NPlmmxvTYz8GwHgGDf1jP9iGinH+JmQ/qk2c3O9QZD97FK5H4XpkP03bv7LSpGQ/H6LRHcTOZD/xaOOItfhkP8Mv9fOmImU/lfYGX5hMZT9nvRjKiXZlPzmEKjV7oGU/C0s8oGzKZT/cEU4LXvRlP67YX3ZPHmY/gJ9x4UBIZj9SZoNMMnJmPyQtlbcjnGY/9vOmIhXGZj/IuriNBvBmP5qByvj3GWc/bEjcY+lDZz89D+7O2m1nPw/W/znMl2c/4ZwRpb3BZz+0YyMQr+tnP4QqNXugFWg/VvFG5pE/aD8ouFhRg2loP/p+arx0k2g/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA9zXgXNGpD99owG8BRK0P3IgQfFjzL0/YDJVMCqpwz9VuB6F61HIPwUi/fZ14Mw/uTf4wmSq0D9IUPwYc9fSP7vaiv1l99Q/DdejcD0K1z89RUdy+Q/ZP0YldQKaCNs/MnctIR/03D/+OnDOiNLeP1S4HoXrUeA/F4xK6gQ04T/ImLuWkA/iP2vecYqO5OI//Fxtxf6y4z99FK5H4XrkP+sENBE2POU/Si7/If325T+YkA96NqvmP9UrZRniWOc/AAAAAAAA6D8cDeAtkKDoPyhTBaOSOuk/ItJvXwfO6T8Iih9j7lrqP+J6FK5H4eo/qaROQBNh6z9gB84ZUdrrPwijkjoBTew/mnecoiO57D8ehetRuB7tP5PLf0i/fe0/9kpZhjjW7T9HA3gLJCjuP4j029eBc+4/uB6F61G47j/YgXNGlPbuP+cdp+hILu8/5fIf0m9f7z/SAN4CCYrvP65H4XoUru8/eccpOpLL7z81gLdAguLvP95xio7k8u8/eJyiI7n87z8AAAAAAADwP3icoiO5/O8/3nGKjuTy7z80gLdAguLvP3nHKTqSy+8/rkfhehSu7z/RAN4CCYrvP+XyH9JvX+8/5x2n6Egu7z/XgXNGlPbuP7gehetRuO4/h/Tb14Fz7j9GA3gLJCjuP/RKWYY41u0/kct/SL997T8dhetRuB7tP5p3nKIjuew/BaOSOgFN7D9fB84ZUdrrP6ikTkATYes/4HoUrkfh6j8Iih9j7lrqPx3Sb18Hzuk/I1MFo5I66T8cDeAtkKDoPwAAAAAAAOg/1CtlGeJY5z+WkA96NqvmP0Yu/yH99uU/5wQ0ETY85T94FK5H4XrkP/dcbcX+suM/Zt5xio7k4j/JmLuWkA/iPxOMSuoENOE/ULgehetR4D/2OnDOiNLePyx3LSEf9Nw/QCV1ApoI2z8yRUdy+Q/ZP/zWo3A9Ctc/uNqK/WX31D9GUPwYc9fSP7I3+MJkqtA/+CH99nXgzD9JuB6F61HIP0kyVTAqqcM/JyBB8WPMvT82owG8BRK0P3zb14FzRqQ/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$a_1(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def a1(x):\n",
+ " \"\"\"Derivative of a0\"\"\"\n",
+ " return (x / eps_v) * (2 - x / eps_v)\n",
+ "\n",
+ "\n",
+ "# def sign(x):\n",
+ "# return Piecewise((1, x > 0), (-1, x < 0), (0, True))\n",
+ "\n",
+ "\n",
+ "sym_a1 = Piecewise(\n",
+ " (0, x <= 0),\n",
+ " (a1(x), x < 2 * eps_v),\n",
+ " (0, x >= 2 * eps_v),\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(\"a_{1}(x)\"), sym_a1.expand()))\n",
+ "\n",
+ "plot(sym_a1, title=r\"$a_1(x)$\", min_x=-1, max_x=3)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Smooth Coefficient of Adhesion Mollifier"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 227,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\mu(x) a_{1(x)} = \\begin{cases} 0 & \\text{for}\\: x \\leq 0 \\\\\\frac{x \\left(2 - \\frac{x}{\\epsilon_{v}}\\right) \\left(\\mu_{s} + \\left(\\mu_{k} - \\mu_{s}\\right) \\left(\\frac{3 x^{2}}{\\epsilon_{v}^{2}} - \\frac{2 x^{3}}{\\epsilon_{v}^{3}}\\right)\\right)}{\\epsilon_{v}} & \\text{for}\\: \\epsilon_{v} \\geq x \\\\\\frac{\\mu_{k} x \\left(2 - \\frac{x}{\\epsilon_{v}}\\right)}{\\epsilon_{v}} & \\text{for}\\: x \\leq 2 \\epsilon_{v} \\\\0 & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(\\mu(x) a_1(x), Piecewise((0, x <= 0), (x*(2 - x/\\epsilon_v)*(\\mu_s + (\\mu_k - \\mu_s)*(3*x**2/\\epsilon_v**2 - 2*x**3/\\epsilon_v**3))/\\epsilon_v, \\epsilon_v >= x), (\\mu_k*x*(2 - x/\\epsilon_v)/\\epsilon_v, x <= 2*\\epsilon_v), (0, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iUL9YHM78ag5Qv2kdVU0QdU+/IQIOoUrNTr/a5sb0hCVOv5LLf0i/fU2/S7A4nPnVTL8DlfHvMy5Mv7x5qkNuhku/dF5jl6jeSr8tQxzr4jZKv+Un1T4dj0m/ngyOklfnSL9W8UbmkT9Ivw7W/znMl0e/x7q4jQbwRr+An3HhQEhGvziEKjV7oEW/8GjjiLX4RL+pTZzc71BEv2IyVTAqqUO/GhcOhGQBQ7/S+8bXnllCv4vgfyvZsUG/Q8U4fxMKQb/8qfHSTWJAv2gdVU0QdT+/2ubG9IQlPr9KsDic+dU8v7x5qkNuhju/LEMc6+I2Or+eDI6SV+c4vw7W/znMlze/fp9x4UBINr/waOOItfg0v2AyVTAqqTO/0vvG155ZMr9CxTh/Ewoxv2gdVU0QdS+/SLA4nPnVLL8sQxzr4jYqvwzW/znMlye/8GjjiLX4JL/Q+8bXnlkiv2gdVU0QdR+/KEMc6+I2Gr/waOOItfgUv2AdVU0QdQ+/4GjjiLX4BL/AaOOItfj0vgAAAAAAAAAAAGnjiLX49D4AaeOItfgEP4AdVU0QdQ8/8GjjiLX4FD8wQxzr4jYaP3AdVU0QdR8/2PvG155ZIj/waOOItfgkPxDW/znMlyc/MEMc6+I2Kj9QsDic+dUsP2gdVU0QdS8/RMU4fxMKMT/U+8bXnlkyP2QyVTAqqTM/9GjjiLX4ND+An3HhQEg2PxDW/znMlzc/oAyOklfnOD8wQxzr4jY6P7x5qkNuhjs/TLA4nPnVPD/c5sb0hCU+P2wdVU0QdT8//Knx0k1iQD9ExTh/EwpBP4zgfyvZsUE/1PvG155ZQj8aFw6EZAFDP2IyVTAqqUM/qk2c3O9QRD/yaOOItfhEPzqEKjV7oEU/gJ9x4UBIRj/IuriNBvBGPxDW/znMl0c/WPFG5pE/SD+eDI6SV+dIP+Yn1T4dj0k/LkMc6+I2Sj92XmOXqN5KP7x5qkNuhks/BJXx7zMuTD9MsDic+dVMP5TLf0i/fU0/3ObG9IQlTj8iAg6hSs1OP2wdVU0QdU8/WBzO/GoOUD/8qfHSTWJQP6A3FakwtlA/RMU4fxMKUT/oUlxV9l1RP4zgfyvZsVE/MG6jAbwFUj/U+8bXnllSP3aJ6q2BrVI/GhcOhGQBUz++pDFaR1VTP2IyVTAqqVM/BsB4Bg39Uz+qTZzc71BUP07bv7LSpFQ/8mjjiLX4VD+W9gZfmExVPziEKjV7oFU/3BFOC170VT+An3HhQEhWPyQtlbcjnFY/yLq4jQbwVj9sSNxj6UNXPxDW/znMl1c/tGMjEK/rVz9W8UbmkT9YP/p+arx0k1g/ngyOklfnWD9CmrFoOjtZP+Yn1T4dj1k/irX4FADjWT8uQxzr4jZaP9LQP8HFilo/dl5jl6jeWj8Y7IZtizJbP7x5qkNuhls/YAfOGVHaWz8ElfHvMy5cP6giFcYWglw/TLA4nPnVXD/wPVxy3CldP5TLf0i/fV0/NlmjHqLRXT/a5sb0hCVeP3506spneV4/IgIOoUrNXj/GjzF3LSFfP2odVU0QdV8/Dqt4I/PIXz9ZHM78ag5gPyvj32dcOGA//Knx0k1iYD/OcAM+P4xgP6A3FakwtmA/cv4mFCLgYD9ExTh/EwphPxaMSuoENGE/6FJcVfZdYT+6GW7A54dhP4vgfyvZsWE/XaeRlsrbYT8vbqMBvAViPwE1tWytL2I/0/vG155ZYj+lwthCkINiP3eJ6q2BrWI/SVD8GHPXYj8bFw6EZAFjP+zdH+9VK2M/vqQxWkdVYz+Qa0PFOH9jP2IyVTAqqWM/NPlmmxvTYz8GwHgGDf1jP9iGinH+JmQ/qk2c3O9QZD97FK5H4XpkP03bv7LSpGQ/H6LRHcTOZD/xaOOItfhkP8Mv9fOmImU/lfYGX5hMZT9nvRjKiXZlPzmEKjV7oGU/C0s8oGzKZT/cEU4LXvRlP67YX3ZPHmY/gJ9x4UBIZj9SZoNMMnJmPyQtlbcjnGY/9vOmIhXGZj/IuriNBvBmP5qByvj3GWc/bEjcY+lDZz89D+7O2m1nPw/W/znMl2c/4ZwRpb3BZz+0YyMQr+tnP4QqNXugFWg/VvFG5pE/aD8ouFhRg2loP/p+arx0k2g/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAh0nYletApD/4oOQEa/yzPzVwT1c1hb0/Xn0HjdRWwz8F02ndBrXHPx2yNhIJ2Ms/TaIOk2C7zz/kwLl2na3RP8GmOAQ0WtM/e3QjLCri1D86cweJgUTWP4AntC+EgNc/SbUPycGV2D/zQ+urDITZP1Ni1/Z2S9o/w2r4qU/s2j8a59rAH2fbP7z0R0ynvNs/magZjNrt2z87cw8J3/vbP8SEoq4I6Ns/9zDa5Naz2z9IUyCq8WDbP9CyFa0m8do/aGZmZmZm2j+WOJ4ywcLZP6kL/WtkCNk/vD1LhJc52D+0DK4euVjXP0z6eyk8aNY/ETAR+KRq1T9846NchmLUP+q5GMJ+UtM/oyzXRTU90j/h7J3RViXRP9xHVzWTDdA/mxXbgTTxzT/VzT68MdLLP/yqqFRsw8k/z/fhICHKxz8yH1fFdevFP010wOhyLMQ/lPrKZ/6Rwj/LLcGI1SDBP0KUZ14Ou78/jShFH9yYvT+CQkrGVeO7P3o9JzY9oro/4DVfKcncuT+gmZmZmZm5Py196IL6lrk/5SfVPh2PuT/CmV/NAYK5P8fShy6ob7k/8dJNYhBYuT9BmrFoOju5P7cos0EmGbk/U35S7dPxuD8Sm49rQ8W4P/l+arx0k7g/BSrj32dcuD84nPnVHCC4P5DVrZ6T3rc/Dtb/OcyXtz+xne+nxku3P3ssfeiC+rY/aoKo+wCktj9/n3HhQEi2P7qD2JlC57U/Gi/dJAaBtT+goX+CixW1P0nbv7LSpLQ/HNydtdsutD8WpBmLprOzPzMzMzMzM7M/donqrYGtsj/fpj/7kSKyP2uLMhtkkrE/HzfDDfj8sD/5qfHSTWKwP/LHe9XKhK8/PMpPqn06rj/cWl8ktOWsP7h5qkNuhqs/5iYxCKwcqj9fYvNxbaioPyMs8YCyKac/M4QqNXugpT+Pap+OxwykPzDfT42XbqI/LeI7MevFoD/V5sb0hCWeP+kljdE6qpo/lIHK+PcZlz/V+X5qvHSTP0EdVU0QdY8/H4CaWrbWhz8rHM78ag6AP8ri32dcOHA/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$\\mu(x) a_1(x)$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "def mu_a1(x):\n",
+ " return mu(x) * a1(x)\n",
+ "\n",
+ "\n",
+ "sym_mu_a1 = Piecewise(\n",
+ " (0, x <= 0),\n",
+ " (mu_a1(x), x <= eps_v),\n",
+ " (mu_k * a1(x), x <= 2 * eps_v),\n",
+ " (0, x > 2 * eps_v)\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(r\"\\mu(x) a_1(x)\"), sym_mu_a1))\n",
+ "\n",
+ "plot(sym_mu_a1, title=r'$\\mu(x) a_1(x)$', min_x=-1, max_x=3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 228,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\int \\mu(x) a_1(x) \\mathrm{d}x = \\begin{cases} 0 & \\text{for}\\: x \\leq 0 \\\\\\frac{\\mu_{s} x^{2}}{\\epsilon_{v}} - \\frac{\\mu_{s} x^{3}}{3 \\epsilon_{v}^{2}} + \\frac{x^{4} \\left(3 \\mu_{k} - 3 \\mu_{s}\\right)}{2 \\epsilon_{v}^{3}} + \\frac{x^{5} \\left(- 7 \\mu_{k} + 7 \\mu_{s}\\right)}{5 \\epsilon_{v}^{4}} + \\frac{x^{6} \\left(\\mu_{k} - \\mu_{s}\\right)}{3 \\epsilon_{v}^{5}} & \\text{for}\\: \\epsilon_{v} \\geq x \\\\\\frac{7 \\epsilon_{v}^{3} \\left(- \\mu_{k} + \\mu_{s}\\right) + 10 \\mu_{k} x^{2} \\left(3 \\epsilon_{v} - x\\right)}{30 \\epsilon_{v}^{2}} & \\text{for}\\: x \\leq 2 \\epsilon_{v} \\\\\\frac{\\epsilon_{v} \\left(33 \\mu_{k} + 7 \\mu_{s}\\right)}{30} & \\text{otherwise} \\end{cases}$"
+ ],
+ "text/plain": [
+ "Eq(\\int \\mu(x) a_1(x) \\mathrm{d}x, Piecewise((0, x <= 0), (\\mu_s*x**2/\\epsilon_v - \\mu_s*x**3/(3*\\epsilon_v**2) + x**4*(3*\\mu_k - 3*\\mu_s)/(2*\\epsilon_v**3) + x**5*(-7*\\mu_k + 7*\\mu_s)/(5*\\epsilon_v**4) + x**6*(\\mu_k - \\mu_s)/(3*\\epsilon_v**5), \\epsilon_v >= x), ((7*\\epsilon_v**3*(-\\mu_k + \\mu_s) + 10*\\mu_k*x**2*(3*\\epsilon_v - x))/(30*\\epsilon_v**2), x <= 2*\\epsilon_v), (\\epsilon_v*(33*\\mu_k + 7*\\mu_s)/30, True)))"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "application/vnd.plotly.v1+json": {
+ "config": {
+ "plotlyServerURL": "https://plot.ly"
+ },
+ "data": [
+ {
+ "type": "scatter",
+ "x": {
+ "bdata": "/Knx0k1iUL9YHM78ag5Qv2kdVU0QdU+/IQIOoUrNTr/a5sb0hCVOv5LLf0i/fU2/S7A4nPnVTL8DlfHvMy5Mv7x5qkNuhku/dF5jl6jeSr8tQxzr4jZKv+Un1T4dj0m/ngyOklfnSL9W8UbmkT9Ivw7W/znMl0e/x7q4jQbwRr+An3HhQEhGvziEKjV7oEW/8GjjiLX4RL+pTZzc71BEv2IyVTAqqUO/GhcOhGQBQ7/S+8bXnllCv4vgfyvZsUG/Q8U4fxMKQb/8qfHSTWJAv2gdVU0QdT+/2ubG9IQlPr9KsDic+dU8v7x5qkNuhju/LEMc6+I2Or+eDI6SV+c4vw7W/znMlze/fp9x4UBINr/waOOItfg0v2AyVTAqqTO/0vvG155ZMr9CxTh/Ewoxv2gdVU0QdS+/SLA4nPnVLL8sQxzr4jYqvwzW/znMlye/8GjjiLX4JL/Q+8bXnlkiv2gdVU0QdR+/KEMc6+I2Gr/waOOItfgUv2AdVU0QdQ+/4GjjiLX4BL/AaOOItfj0vgAAAAAAAAAAAGnjiLX49D4AaeOItfgEP4AdVU0QdQ8/8GjjiLX4FD8wQxzr4jYaP3AdVU0QdR8/2PvG155ZIj/waOOItfgkPxDW/znMlyc/MEMc6+I2Kj9QsDic+dUsP2gdVU0QdS8/RMU4fxMKMT/U+8bXnlkyP2QyVTAqqTM/9GjjiLX4ND+An3HhQEg2PxDW/znMlzc/oAyOklfnOD8wQxzr4jY6P7x5qkNuhjs/TLA4nPnVPD/c5sb0hCU+P2wdVU0QdT8//Knx0k1iQD9ExTh/EwpBP4zgfyvZsUE/1PvG155ZQj8aFw6EZAFDP2IyVTAqqUM/qk2c3O9QRD/yaOOItfhEPzqEKjV7oEU/gJ9x4UBIRj/IuriNBvBGPxDW/znMl0c/WPFG5pE/SD+eDI6SV+dIP+Yn1T4dj0k/LkMc6+I2Sj92XmOXqN5KP7x5qkNuhks/BJXx7zMuTD9MsDic+dVMP5TLf0i/fU0/3ObG9IQlTj8iAg6hSs1OP2wdVU0QdU8/WBzO/GoOUD/8qfHSTWJQP6A3FakwtlA/RMU4fxMKUT/oUlxV9l1RP4zgfyvZsVE/MG6jAbwFUj/U+8bXnllSP3aJ6q2BrVI/GhcOhGQBUz++pDFaR1VTP2IyVTAqqVM/BsB4Bg39Uz+qTZzc71BUP07bv7LSpFQ/8mjjiLX4VD+W9gZfmExVPziEKjV7oFU/3BFOC170VT+An3HhQEhWPyQtlbcjnFY/yLq4jQbwVj9sSNxj6UNXPxDW/znMl1c/tGMjEK/rVz9W8UbmkT9YP/p+arx0k1g/ngyOklfnWD9CmrFoOjtZP+Yn1T4dj1k/irX4FADjWT8uQxzr4jZaP9LQP8HFilo/dl5jl6jeWj8Y7IZtizJbP7x5qkNuhls/YAfOGVHaWz8ElfHvMy5cP6giFcYWglw/TLA4nPnVXD/wPVxy3CldP5TLf0i/fV0/NlmjHqLRXT/a5sb0hCVeP3506spneV4/IgIOoUrNXj/GjzF3LSFfP2odVU0QdV8/Dqt4I/PIXz9ZHM78ag5gPyvj32dcOGA//Knx0k1iYD/OcAM+P4xgP6A3FakwtmA/cv4mFCLgYD9ExTh/EwphPxaMSuoENGE/6FJcVfZdYT+6GW7A54dhP4vgfyvZsWE/XaeRlsrbYT8vbqMBvAViPwE1tWytL2I/0/vG155ZYj+lwthCkINiP3eJ6q2BrWI/SVD8GHPXYj8bFw6EZAFjP+zdH+9VK2M/vqQxWkdVYz+Qa0PFOH9jP2IyVTAqqWM/NPlmmxvTYz8GwHgGDf1jP9iGinH+JmQ/qk2c3O9QZD97FK5H4XpkP03bv7LSpGQ/H6LRHcTOZD/xaOOItfhkP8Mv9fOmImU/lfYGX5hMZT9nvRjKiXZlPzmEKjV7oGU/C0s8oGzKZT/cEU4LXvRlP67YX3ZPHmY/gJ9x4UBIZj9SZoNMMnJmPyQtlbcjnGY/9vOmIhXGZj/IuriNBvBmP5qByvj3GWc/bEjcY+lDZz89D+7O2m1nPw/W/znMl2c/4ZwRpb3BZz+0YyMQr+tnP4QqNXugFWg/VvFG5pE/aD8ouFhRg2loP/p+arx0k2g/",
+ "dtype": "f8"
+ },
+ "y": {
+ "bdata": "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAqeKdJn6mmj7Ohd30Bm66Pg8GKe3SdM0+Fz/Zv6vp2T7gRWfK4APkPo4hFvfDd+w+xADU6Swe8z4uQMtlEJ74PnA0Hbc4sP4+8iSfq2ykAj830AgP4y0GP1X6rr1F7gk/fGOI3kLfDT+i75PYPf0QP6ES1p7HHBM/+1+LRRJLFT8aDQUu+oQXP8RV9zhpxxk/NKE/K1sPHD/YNt7V4VkePxfBkYAUUiA/9HKIDb11IT8sj3bUn5YiP5YZBj6EsyM/OHzlVEXLJD9tD18j09wlP6dpivIz5yY/2HYXa4XpJz+OV7OX/eIoP5QHB8nr0ik/Z8xPW7m4Kj8Na5Fd6pMrP7QlYhoeZCw/LoFQgg8pLT+M0eJ3leItP9SOMP2ikC4/KHEVQ0czLz+UVP6ZrcovPz7yJ6KOKzA/aoczlXxsMD9Qnxm4X6gwP1cvNXaE3zA/pFGrQEQSMT8CDp8rBUExP44GMnw5bDE/HwhSJ1+UMT+YflNB/7kxPw/NWF6t3TE/pYmG4wYAMj9fnQRJsiEyP9gKEhM/QzI/1IHn/cRkMj/lC00qPYYyP5SyCrmgpzI/an/oyujIMj/0e66ADuoyP7ixJPsKCzM/RSoTW9crMz8h70HBbEwzP9gJeU7EbDM/9IOAI9eMMz/+ZiBhnqwzP4K8ICgTzDM/CI5JmS7rMz8b5WLV6Qk0P0TLNP09KDQ/D0qHMSRGND8FayKTlWM0P7A3zkKLgDQ/m7lSYf6cND9O+ncP6Lg0P1YDBm5B1DQ/Ot7EnQPvND+GlHy/Jwk1P8Qv9fOmIjU/fbn2W3o7NT88O0kYm1M1P4q+tEkCazU/8kwBEamBNT//7/aOiJc1PzmxXeSZrDU/LJr9MdbANT9gtJ6YNtQ1P2AJCTm05jU/tqIENEj4NT/uiVmq6wg2P47Iz7yXGDY/JGgvjEUnNj84ckA57jQ2P1XwyuSKQTY/BOyWrxRNNj/Qbmy6hFc2P0KCEybUYDY/5i9UE/xoNj9Egfai9W82P+Z/wvW5dTY/WTWALEJ6Nj8kq/dnh302P9Lq8MiCfzY/7v0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/7f0zcC2ANj/t/TNwLYA2P+39M3AtgDY/",
+ "dtype": "f8"
+ }
+ }
+ ],
+ "layout": {
+ "height": 600,
+ "template": {
+ "data": {
+ "scatter": [
+ {
+ "type": "scatter"
+ }
+ ]
+ }
+ },
+ "title": {
+ "text": "$\\int \\mu(x) a_1(x) \\mathrm{d}x$"
+ },
+ "width": 800,
+ "xaxis": {
+ "title": {
+ "text": "$x$"
+ }
+ }
+ }
+ },
+ "text/html": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "mu_a0 = integrate(mu_a1(x), x)\n",
+ "\n",
+ "# Constant of integration s.t. mu_a0(eps_v) = mu_k * a0_part1(eps_v)\n",
+ "c = (mu_k * a0_part1(eps_v) - mu_a0.subs(x, eps_v)).simplify()\n",
+ "mu_a0 = mu_a0\n",
+ "\n",
+ "mu_a0.simplify()\n",
+ "\n",
+ "sym_mu_a0 = Piecewise(\n",
+ " (0, x <= 0),\n",
+ " (mu_a0, x <= eps_v),\n",
+ " ((mu_k * a0_part1(x) - c).simplify(), x <= 2 * eps_v),\n",
+ " ((mu_k * a0_part2(x) - c).simplify(), x > 2 * eps_v),\n",
+ ")\n",
+ "\n",
+ "display(Eq(Symbol(r\"\\int \\mu(x) a_1(x) \\mathrm{d}x\"), sym_mu_a0))\n",
+ "\n",
+ "plot(sym_mu_a0, title=r'$\\int \\mu(x) a_1(x) \\mathrm{d}x$', min_x=-1, max_x=3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 229,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{\\mu_{s} x^{2}}{\\epsilon_{v}} - \\frac{\\mu_{s} x^{3}}{3 \\epsilon_{v}^{2}} + \\frac{x^{4} \\left(3 \\mu_{k} - 3 \\mu_{s}\\right)}{2 \\epsilon_{v}^{3}} + \\frac{x^{5} \\left(- 7 \\mu_{k} + 7 \\mu_{s}\\right)}{5 \\epsilon_{v}^{4}} + \\frac{x^{6} \\left(\\mu_{k} - \\mu_{s}\\right)}{3 \\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "\\mu_s*x**2/\\epsilon_v - \\mu_s*x**3/(3*\\epsilon_v**2) + x**4*(3*\\mu_k - 3*\\mu_s)/(2*\\epsilon_v**3) + x**5*(-7*\\mu_k + 7*\\mu_s)/(5*\\epsilon_v**4) + x**6*(\\mu_k - \\mu_s)/(3*\\epsilon_v**5)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# print(pycode(mu_a0.subs({eps_v: \"eps_v\", mu_s: \"mu_s\", mu_k: \"mu_k\"})))\n",
+ "display(mu_a0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 230,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def mu_a0(x, eps_v, mu_s, mu_k):\n",
+ " \"\"\"Integral of mu(x) a1(x).\"\"\"\n",
+ " d_mu = mu_k - mu_s\n",
+ " x_over_eps_v = x / eps_v\n",
+ " c = 7 * eps_v * d_mu / 30\n",
+ " return x * x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * (\n",
+ " x_over_eps_v * d_mu / 3 - 7 * d_mu / 5\n",
+ " ) + 3 * d_mu / 2\n",
+ " ) - mu_s/3\n",
+ " ) + mu_s) + c"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 231,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{7 \\epsilon_{v} \\mu_{k}}{30} - \\frac{7 \\epsilon_{v} \\mu_{s}}{30} + \\frac{\\mu_{s} x^{2}}{\\epsilon_{v}} - \\frac{\\mu_{s} x^{3}}{3 \\epsilon_{v}^{2}} + \\frac{3 \\mu_{k} x^{4}}{2 \\epsilon_{v}^{3}} - \\frac{3 \\mu_{s} x^{4}}{2 \\epsilon_{v}^{3}} - \\frac{7 \\mu_{k} x^{5}}{5 \\epsilon_{v}^{4}} + \\frac{7 \\mu_{s} x^{5}}{5 \\epsilon_{v}^{4}} + \\frac{\\mu_{k} x^{6}}{3 \\epsilon_{v}^{5}} - \\frac{\\mu_{s} x^{6}}{3 \\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "7*\\epsilon_v*\\mu_k/30 - 7*\\epsilon_v*\\mu_s/30 + \\mu_s*x**2/\\epsilon_v - \\mu_s*x**3/(3*\\epsilon_v**2) + 3*\\mu_k*x**4/(2*\\epsilon_v**3) - 3*\\mu_s*x**4/(2*\\epsilon_v**3) - 7*\\mu_k*x**5/(5*\\epsilon_v**4) + 7*\\mu_s*x**5/(5*\\epsilon_v**4) + \\mu_k*x**6/(3*\\epsilon_v**5) - \\mu_s*x**6/(3*\\epsilon_v**5)"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "display(a := mu_a0(x, eps_v, mu_s, mu_k).expand())"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 232,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{2 \\mu_{s} x}{\\epsilon_{v}} - \\frac{\\mu_{s} x^{2}}{\\epsilon_{v}^{2}} + \\frac{6 \\mu_{k} x^{3}}{\\epsilon_{v}^{3}} - \\frac{6 \\mu_{s} x^{3}}{\\epsilon_{v}^{3}} - \\frac{7 \\mu_{k} x^{4}}{\\epsilon_{v}^{4}} + \\frac{7 \\mu_{s} x^{4}}{\\epsilon_{v}^{4}} + \\frac{2 \\mu_{k} x^{5}}{\\epsilon_{v}^{5}} - \\frac{2 \\mu_{s} x^{5}}{\\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "2*\\mu_s*x/\\epsilon_v - \\mu_s*x**2/\\epsilon_v**2 + 6*\\mu_k*x**3/\\epsilon_v**3 - 6*\\mu_s*x**3/\\epsilon_v**3 - 7*\\mu_k*x**4/\\epsilon_v**4 + 7*\\mu_s*x**4/\\epsilon_v**4 + 2*\\mu_k*x**5/\\epsilon_v**5 - 2*\\mu_s*x**5/\\epsilon_v**5"
+ ]
+ },
+ "execution_count": 232,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "a.diff(x)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 233,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\frac{2 \\mu_{s} x}{\\epsilon_{v}} - \\frac{\\mu_{s} x^{2}}{\\epsilon_{v}^{2}} + \\frac{6 \\mu_{k} x^{3}}{\\epsilon_{v}^{3}} - \\frac{6 \\mu_{s} x^{3}}{\\epsilon_{v}^{3}} - \\frac{7 \\mu_{k} x^{4}}{\\epsilon_{v}^{4}} + \\frac{7 \\mu_{s} x^{4}}{\\epsilon_{v}^{4}} + \\frac{2 \\mu_{k} x^{5}}{\\epsilon_{v}^{5}} - \\frac{2 \\mu_{s} x^{5}}{\\epsilon_{v}^{5}}$"
+ ],
+ "text/plain": [
+ "2*\\mu_s*x/\\epsilon_v - \\mu_s*x**2/\\epsilon_v**2 + 6*\\mu_k*x**3/\\epsilon_v**3 - 6*\\mu_s*x**3/\\epsilon_v**3 - 7*\\mu_k*x**4/\\epsilon_v**4 + 7*\\mu_s*x**4/\\epsilon_v**4 + 2*\\mu_k*x**5/\\epsilon_v**5 - 2*\\mu_s*x**5/\\epsilon_v**5"
+ ]
+ },
+ "execution_count": 233,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "(mu(x) * a1(x)).expand()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 234,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle 8 \\mu_{k} x^{2} z^{5} - 21 \\mu_{k} x z^{4} + 12 \\mu_{k} z^{3} - 8 \\mu_{s} x^{2} z^{5} + 21 \\mu_{s} x z^{4} - 12 \\mu_{s} z^{3} - \\frac{\\mu_{s} z^{2}}{x}$"
+ ],
+ "text/plain": [
+ "8*\\mu_k*x**2*z**5 - 21*\\mu_k*x*z**4 + 12*\\mu_k*z**3 - 8*\\mu_s*x**2*z**5 + 21*\\mu_s*x*z**4 - 12*\\mu_s*z**3 - \\mu_s*z**2/x"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/latex": [
+ "$\\displaystyle \\operatorname{Poly}{\\left( \\left(8 \\mu_{k} x^{2} - 8 \\mu_{s} x^{2}\\right) z^{5} + \\left(- 21 \\mu_{k} x + 21 \\mu_{s} x\\right) z^{4} + \\left(12 \\mu_{k} - 12 \\mu_{s}\\right) z^{3} - \\frac{\\mu_{s}}{x} z^{2}, z, domain=\\mathbb{Z}\\left(x, \\mu_{k}, \\mu_{s}\\right) \\right)}$"
+ ],
+ "text/plain": [
+ "Poly((8*\\mu_k*x**2 - 8*\\mu_s*x**2)*z**5 + (-21*\\mu_k*x + 21*\\mu_s*x)*z**4 + (12*\\mu_k - 12*\\mu_s)*z**3 - \\mu_s/x*z**2, z, domain='ZZ(x,\\mu_k,\\mu_s)')"
+ ]
+ },
+ "execution_count": 234,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "a2 = (((mu(x) * a1(x)).diff(x) * x - mu(x) * a1(x)) / x**3).simplify()\n",
+ "display(a2.expand().subs({1/eps_v: Symbol('z')}))\n",
+ "Poly(a2.expand().subs({1/eps_v: Symbol('z')}), Symbol('z'))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.13.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/python/src/adhesion/adhesion.cpp b/python/src/adhesion/adhesion.cpp
index cf6152f19..aeb319269 100644
--- a/python/src/adhesion/adhesion.cpp
+++ b/python/src/adhesion/adhesion.cpp
@@ -142,4 +142,102 @@ void define_adhesion(py::module_& m)
The second derivative of the tangential adhesion mollifier function times y minus the first derivative all divided by y cubed.
)ipc_Qu8mg5v7",
"y"_a, "eps_a"_a);
+
+ m.def(
+ "smooth_mu_a0", &smooth_mu_a0,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the ∫ μ(y) a₁(y) dy, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+
+ Note:
+ The `a0`/`a1` are unrelated to the `a0`/`a1` in the normal adhesion.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static adhesion.
+ mu_k: Coefficient of kinetic adhesion.
+ eps_a: Velocity threshold below which static adhesion force is applied.
+
+ Returns:
+ The value of the integral at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_a"_a);
+
+ m.def(
+ "smooth_mu_a1", &smooth_mu_a1,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the μ(y) a₁(y), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+
+ Note:
+ The `a1` is unrelated to the `a1` in the normal adhesion.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static adhesion.
+ mu_k: Coefficient of kinetic adhesion.
+ eps_a: Velocity threshold below which static adhesion force is applied.
+
+ Returns:
+ The value of the product at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_a"_a);
+
+ m.def(
+ "smooth_mu_a2", &smooth_mu_a2,
+ R"ipc_Qu8mg5v7(
+ Compute the value of d/dy (μ(y) a₁(y)), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+
+ Note:
+ The `a1`/`a2` are unrelated to the `a1`/`a2` in the normal adhesion.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static adhesion.
+ mu_k: Coefficient of kinetic adhesion.
+ eps_a: Velocity threshold below which static adhesion force is applied.
+
+ Returns:
+ The value of the derivative at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_a"_a);
+
+ m.def(
+ "smooth_mu_a1_over_x", &smooth_mu_a1_over_x,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the μ(y) a₁(y) / y, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+
+ Notes:
+ The `x` in the function name refers to the parameter `y`.
+ The `a1` is unrelated to the `a1` in the normal adhesion.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static adhesion.
+ mu_k: Coefficient of kinetic adhesion.
+ eps_a: Velocity threshold below which static adhesion force is applied.
+
+ Returns:
+ The value of the product at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_a"_a);
+
+ m.def(
+ "smooth_mu_a2_x_minus_mu_a1_over_x3",
+ &smooth_mu_a2_x_minus_mu_a1_over_x3,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the [(d/dy μ(y) a₁(y)) ⋅ y - μ(y) a₁(y)] / y³, where a₁ and a₂ are the first and second derivatives of the smooth tangential adhesion mollifier.
+
+ Notes:
+ The `x` in the function name refers to the parameter `y`.
+ The `a1`/`a2` are unrelated to the `a1`/`a2` in the normal adhesion.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static adhesion.
+ mu_k: Coefficient of kinetic adhesion.
+ eps_a: Velocity threshold below which static adhesion force is applied.
+
+ Returns:
+ The value of the expression at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_a"_a);
}
diff --git a/python/src/bindings.cpp b/python/src/bindings.cpp
index 90413d2a3..3e5422dec 100644
--- a/python/src/bindings.cpp
+++ b/python/src/bindings.cpp
@@ -86,6 +86,7 @@ PYBIND11_MODULE(ipctk, m)
// friction
define_smooth_friction_mollifier(m);
+ define_smooth_mu(m);
// implicits
define_plane_implicit(m);
diff --git a/python/src/collisions/tangential/tangential_collision.cpp b/python/src/collisions/tangential/tangential_collision.cpp
index 9e1a24076..0292088ab 100644
--- a/python/src/collisions/tangential/tangential_collision.cpp
+++ b/python/src/collisions/tangential/tangential_collision.cpp
@@ -122,8 +122,11 @@ void define_tangential_collision(py::module_& m)
&TangentialCollision::normal_force_magnitude,
"Normal force magnitude")
.def_readwrite(
- "mu", &TangentialCollision::mu,
- "Ratio between normal and tangential forces (e.g., friction coefficient)")
+ "mu_s", &TangentialCollision::mu_s,
+ "Ratio between normal and static tangential forces (e.g., friction coefficient)")
+ .def_readwrite(
+ "mu_k", &TangentialCollision::mu_k,
+ "Ratio between normal and kinetic tangential forces (e.g., friction coefficient)")
.def_readwrite("weight", &TangentialCollision::weight, "Weight")
.def_property(
"weight_gradient",
diff --git a/python/src/collisions/tangential/tangential_collisions.cpp b/python/src/collisions/tangential/tangential_collisions.cpp
index e30c33f1f..15aa015b9 100644
--- a/python/src/collisions/tangential/tangential_collisions.cpp
+++ b/python/src/collisions/tangential/tangential_collisions.cpp
@@ -16,6 +16,14 @@ void define_tangential_collisions(py::module_& m)
double>(&TangentialCollisions::build),
"mesh"_a, "vertices"_a, "collisions"_a, "normal_potential"_a,
"normal_stiffness"_a, "mu"_a)
+ .def(
+ "build",
+ py::overload_cast<
+ const CollisionMesh&, Eigen::ConstRef,
+ const NormalCollisions&, const NormalPotential&, double, double,
+ double>(&TangentialCollisions::build),
+ "mesh"_a, "vertices"_a, "collisions"_a, "normal_potential"_a,
+ "normal_stiffness"_a, "mu_s"_a, "mu_k"_a)
.def(
"build",
[](TangentialCollisions& self, const CollisionMesh& mesh,
@@ -23,23 +31,25 @@ void define_tangential_collisions(py::module_& m)
const NormalCollisions& collisions,
const NormalPotential& normal_potential,
const double normal_stiffness,
- Eigen::ConstRef mus) {
+ Eigen::ConstRef mu_s,
+ Eigen::ConstRef mu_k) {
self.build(
mesh, vertices, collisions, normal_potential,
- normal_stiffness, mus);
+ normal_stiffness, mu_s, mu_k);
},
"mesh"_a, "vertices"_a, "collisions"_a, "normal_potential"_a,
- "normal_stiffness"_a, "mus"_a)
+ "normal_stiffness"_a, "mu_s"_a, "mu_k"_a)
.def(
"build",
py::overload_cast<
const CollisionMesh&, Eigen::ConstRef,
const NormalCollisions&, const NormalPotential&, const double,
Eigen::ConstRef,
+ Eigen::ConstRef,
const std::function&>(
&TangentialCollisions::build),
"mesh"_a, "vertices"_a, "collisions"_a, "normal_potential"_a,
- "normal_stiffness"_a, "mus"_a, "blend_mu"_a)
+ "normal_stiffness"_a, "mu_s"_a, "mu_k"_a, "blend_mu"_a)
.def(
"__len__", &TangentialCollisions::size,
"Get the number of friction collisions.")
diff --git a/python/src/friction/CMakeLists.txt b/python/src/friction/CMakeLists.txt
index c81476d0e..c3d8d6d7b 100644
--- a/python/src/friction/CMakeLists.txt
+++ b/python/src/friction/CMakeLists.txt
@@ -1,5 +1,6 @@
set(SOURCES
smooth_friction_mollifier.cpp
+ smooth_mu.cpp
)
target_sources(ipctk PRIVATE ${SOURCES})
\ No newline at end of file
diff --git a/python/src/friction/bindings.hpp b/python/src/friction/bindings.hpp
index 82ee13fa8..13a318c26 100644
--- a/python/src/friction/bindings.hpp
+++ b/python/src/friction/bindings.hpp
@@ -2,4 +2,5 @@
#include
-void define_smooth_friction_mollifier(py::module_& m);
\ No newline at end of file
+void define_smooth_friction_mollifier(py::module_& m);
+void define_smooth_mu(py::module_& m);
\ No newline at end of file
diff --git a/python/src/friction/smooth_mu.cpp b/python/src/friction/smooth_mu.cpp
new file mode 100644
index 000000000..ada7654e9
--- /dev/null
+++ b/python/src/friction/smooth_mu.cpp
@@ -0,0 +1,127 @@
+#include
+
+#include
+
+using namespace ipc;
+
+void define_smooth_mu(py::module_& m)
+{
+ m.def(
+ "smooth_mu", &smooth_mu,
+ R"ipc_Qu8mg5v7(
+ Smooth coefficient from static to kinetic friction.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the μ at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_derivative", &smooth_mu_derivative,
+ R"ipc_Qu8mg5v7(
+ Compute the derivative of the smooth coefficient from static to kinetic friction.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the derivative at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_f0", &smooth_mu_f0,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth friction mollifier.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the integral at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_f1", &smooth_mu_f1,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the μ(y) f₁(y), where f₁ is the first derivative of the smooth friction mollifier.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the product at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_f2", &smooth_mu_f2,
+ R"ipc_Qu8mg5v7(
+ Compute the value of d/dy (μ(y) f₁(y)), where f₁ is the first derivative of the smooth friction mollifier.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the derivative at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_f1_over_x", &smooth_mu_f1_over_x,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the μ(y) f₁(y) / y, where f₁ is the first derivative of the smooth friction mollifier.
+
+ Note:
+ The `x` in the function name refers to the parameter `y`.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the product at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+
+ m.def(
+ "smooth_mu_f2_x_minus_mu_f1_over_x3",
+ &smooth_mu_f2_x_minus_mu_f1_over_x3,
+ R"ipc_Qu8mg5v7(
+ Compute the value of the [(d/dy μ(y) f₁(y)) ⋅ y - μ(y) f₁(y)] / y³, where f₁ and f₂ are the first and second derivatives of the smooth friction mollifier.
+
+ Note:
+ The `x` in the function name refers to the parameter `y`.
+
+ Parameters:
+ y: The tangential relative speed.
+ mu_s: Coefficient of static friction.
+ mu_k: Coefficient of kinetic friction.
+ eps_v: Velocity threshold below which static friction force is applied.
+
+ Returns:
+ The value of the expression at y.
+ )ipc_Qu8mg5v7",
+ "y"_a, "mu_s"_a, "mu_k"_a, "eps_v"_a);
+}
diff --git a/src/ipc/adhesion/adhesion.cpp b/src/ipc/adhesion/adhesion.cpp
index 2a6d39f12..336d7236a 100644
--- a/src/ipc/adhesion/adhesion.cpp
+++ b/src/ipc/adhesion/adhesion.cpp
@@ -2,6 +2,8 @@
#include "adhesion.hpp"
+#include
+
#include
namespace ipc {
@@ -101,7 +103,7 @@ double tangential_adhesion_f2(const double y, const double eps_a)
return 0;
}
- return (-y / eps_a + 1) * 2 / eps_a; // -2y/ϵ² + 2/ϵ
+ return (2 - 2 * y / eps_a) / eps_a; // -2y/ϵ² + 2/ϵ
}
double tangential_adhesion_f1_over_x(const double y, const double eps_a)
@@ -125,4 +127,62 @@ tangential_adhesion_f2_x_minus_f1_over_x3(const double y, const double eps_a)
return -1 / (y * eps_a * eps_a);
}
+// ~~ Smooth μ variants ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// Here a0, a1, and a2 refer to the mollifier functions above.
+
+double smooth_mu_a0(
+ const double y, const double mu_s, const double mu_k, const double eps_a)
+{
+ assert(eps_a > 0);
+ const double delta_mu = mu_k - mu_s;
+ if (mu_s == mu_k || y <= 0 || y >= eps_a) {
+ // If the static and kinetic friction coefficients are equal, simplify.
+ const double c = (7 / 30.) * 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);
+ }
+}
+
+double smooth_mu_a1(
+ const double y, const double mu_s, const double mu_k, const double eps_a)
+{
+ return smooth_mu(y, mu_s, mu_k, eps_a) * tangential_adhesion_f1(y, eps_a);
+}
+
+double smooth_mu_a2(
+ const double y, const double mu_s, const double mu_k, const double eps_a)
+{
+ return smooth_mu_derivative(y, mu_s, mu_k, eps_a)
+ * tangential_adhesion_f1(y, eps_a)
+ + smooth_mu(y, mu_s, mu_k, eps_a) * tangential_adhesion_f2(y, eps_a);
+}
+
+double smooth_mu_a1_over_x(
+ const double y, const double mu_s, const double mu_k, const double eps_a)
+{
+ // This is a known formulation: μ(y) f₁(y) / y
+ // where we use the robust division by y to avoid division by zero.
+ return smooth_mu(y, mu_s, mu_k, eps_a)
+ * tangential_adhesion_f1_over_x(y, eps_a);
+}
+
+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) {
+ // 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);
+ }
+}
+
} // namespace ipc
diff --git a/src/ipc/adhesion/adhesion.hpp b/src/ipc/adhesion/adhesion.hpp
index 1d1d75164..8f0b8269b 100644
--- a/src/ipc/adhesion/adhesion.hpp
+++ b/src/ipc/adhesion/adhesion.hpp
@@ -80,6 +80,61 @@ double tangential_adhesion_f1_over_x(const double y, const double eps_a);
double
tangential_adhesion_f2_x_minus_f1_over_x3(const double y, const double eps_a);
+// ~~ Smooth μ variants ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+// NOTE: Here a0, a1, and a2 refer to the mollifier functions above.
+
+/// @brief Compute the value of the ∫ μ(y) a₁(y) dy, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+/// @note The `a0`/`a1` are unrelated to the `a0`/`a1` in the normal adhesion.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static adhesion.
+/// @param mu_k Coefficient of kinetic adhesion.
+/// @param eps_a Velocity threshold below which static adhesion force is applied.
+/// @return The value of the integral at y.
+double smooth_mu_a0(
+ const double y, const double mu_s, const double mu_k, const double eps_a);
+
+/// @brief Compute the value of the μ(y) a₁(y), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+/// @note The `a1` is unrelated to the `a1` in the normal adhesion.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static adhesion.
+/// @param mu_k Coefficient of kinetic adhesion.
+/// @param eps_a Velocity threshold below which static adhesion force is applied.
+/// @return The value of the product at y.
+double smooth_mu_a1(
+ const double y, const double mu_s, const double mu_k, const double eps_a);
+
+/// @brief Compute the value of d/dy (μ(y) a₁(y)), where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+/// @note The `a1`/`a2` are unrelated to the `a1`/`a2` in the normal adhesion.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static adhesion.
+/// @param mu_k Coefficient of kinetic adhesion.
+/// @param eps_a Velocity threshold below which static adhesion force is applied.
+/// @return The value of the derivative at y.
+double smooth_mu_a2(
+ const double y, const double mu_s, const double mu_k, const double eps_a);
+
+/// @brief Compute the value of the μ(y) a₁(y) / y, where a₁ is the first derivative of the smooth tangential adhesion mollifier.
+/// @note The `x` in the function name refers to the parameter `y`.
+/// @note The `a1` is unrelated to the `a1` in the normal adhesion.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static adhesion.
+/// @param mu_k Coefficient of kinetic adhesion.
+/// @param eps_a Velocity threshold below which static adhesion force is applied.
+/// @return The value of the product at y.
+double smooth_mu_a1_over_x(
+ const double y, const double mu_s, const double mu_k, const double eps_a);
+
+/// @brief Compute the value of the [(d/dy μ(y) a₁(y)) ⋅ y - μ(y) a₁(y)] / y³, where a₁ and a₂ are the first and second derivatives of the smooth tangential adhesion mollifier.
+/// @note The `x` in the function name refers to the parameter `y`.
+/// @note The `a1`/`a2` are unrelated to the `a1`/`a2` in the normal adhesion.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static adhesion.
+/// @param mu_k Coefficient of kinetic adhesion.
+/// @param eps_a Velocity threshold below which static adhesion force is applied.
+/// @return The value of the expression at y.
+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);
+
/// @}
} // namespace ipc
diff --git a/src/ipc/collisions/tangential/tangential_collision.hpp b/src/ipc/collisions/tangential/tangential_collision.hpp
index 4a0d40577..793eea1c9 100644
--- a/src/ipc/collisions/tangential/tangential_collision.hpp
+++ b/src/ipc/collisions/tangential/tangential_collision.hpp
@@ -85,8 +85,11 @@ class TangentialCollision : virtual public CollisionStencil {
/// @brief Normal force magnitude
double normal_force_magnitude;
- /// @brief Ratio between normal and tangential forces (e.g., friction coefficient)
- double mu;
+ /// @brief Ratio between normal and static tangential forces (e.g., friction coefficient)
+ double mu_s;
+
+ /// @brief Ratio between normal and kinetic tangential forces (e.g., friction coefficient)
+ double mu_k;
/// @brief Weight
double weight = 1;
diff --git a/src/ipc/collisions/tangential/tangential_collisions.cpp b/src/ipc/collisions/tangential/tangential_collisions.cpp
index 8b1b0a839..466b32e88 100644
--- a/src/ipc/collisions/tangential/tangential_collisions.cpp
+++ b/src/ipc/collisions/tangential/tangential_collisions.cpp
@@ -17,10 +17,12 @@ void TangentialCollisions::build(
const NormalCollisions& collisions,
const NormalPotential& normal_potential,
const double normal_stiffness,
- Eigen::ConstRef mus,
+ Eigen::ConstRef mu_s,
+ Eigen::ConstRef mu_k,
const std::function& blend_mu)
{
- assert(mus.size() == vertices.rows());
+ assert(mu_s.size() == vertices.rows());
+ assert(mu_k.size() == vertices.rows());
const Eigen::MatrixXi& edges = mesh.edges();
const Eigen::MatrixXi& faces = mesh.faces();
@@ -40,7 +42,8 @@ void TangentialCollisions::build(
normal_stiffness);
const auto& [v0i, v1i, _, __] = FC_vv.back().vertex_ids(edges, faces);
- FC_vv.back().mu = blend_mu(mus(v0i), mus(v1i));
+ FC_vv.back().mu_s = blend_mu(mu_s(v0i), mu_s(v1i));
+ FC_vv.back().mu_k = blend_mu(mu_k(v0i), mu_k(v1i));
}
FC_ev.reserve(C_ev.size());
@@ -51,8 +54,11 @@ void TangentialCollisions::build(
const auto& [vi, e0i, e1i, _] = FC_ev.back().vertex_ids(edges, faces);
const double edge_mu =
- (mus(e1i) - mus(e0i)) * FC_ev.back().closest_point[0] + mus(e0i);
- FC_ev.back().mu = blend_mu(edge_mu, mus(vi));
+ (mu_s(e1i) - mu_s(e0i)) * FC_ev.back().closest_point[0] + mu_s(e0i);
+ FC_ev.back().mu_s = blend_mu(edge_mu, mu_s(vi));
+ const double edge_mu_k =
+ (mu_k(e1i) - mu_k(e0i)) * FC_ev.back().closest_point[0] + mu_k(e0i);
+ FC_ev.back().mu_k = blend_mu(edge_mu_k, mu_k(vi));
}
FC_ee.reserve(C_ee.size());
@@ -72,11 +78,21 @@ void TangentialCollisions::build(
c_ee, c_ee.dof(vertices, edges, faces), normal_potential,
normal_stiffness);
- double ea_mu =
- (mus(ea1i) - mus(ea0i)) * FC_ee.back().closest_point[0] + mus(ea0i);
- double eb_mu =
- (mus(eb1i) - mus(eb0i)) * FC_ee.back().closest_point[1] + mus(eb0i);
- FC_ee.back().mu = blend_mu(ea_mu, eb_mu);
+ double ea_mu_s =
+ (mu_s(ea1i) - mu_s(ea0i)) * FC_ee.back().closest_point[0]
+ + mu_s(ea0i);
+ double eb_mu_s =
+ (mu_s(eb1i) - mu_s(eb0i)) * FC_ee.back().closest_point[1]
+ + mu_s(eb0i);
+ FC_ee.back().mu_s = blend_mu(ea_mu_s, eb_mu_s);
+
+ double ea_mu_k =
+ (mu_k(ea1i) - mu_k(ea0i)) * FC_ee.back().closest_point[0]
+ + mu_k(ea0i);
+ double eb_mu_k =
+ (mu_k(eb1i) - mu_k(eb0i)) * FC_ee.back().closest_point[1]
+ + mu_k(eb0i);
+ FC_ee.back().mu_k = blend_mu(ea_mu_k, eb_mu_k);
}
FC_fv.reserve(C_fv.size());
@@ -86,10 +102,15 @@ void TangentialCollisions::build(
normal_stiffness);
const auto& [vi, f0i, f1i, f2i] = FC_fv.back().vertex_ids(edges, faces);
- double face_mu = mus(f0i)
- + FC_fv.back().closest_point[0] * (mus(f1i) - mus(f0i))
- + FC_fv.back().closest_point[1] * (mus(f2i) - mus(f0i));
- FC_fv.back().mu = blend_mu(face_mu, mus(vi));
+ double face_mu_s = mu_s(f0i)
+ + FC_fv.back().closest_point[0] * (mu_s(f1i) - mu_s(f0i))
+ + FC_fv.back().closest_point[1] * (mu_s(f2i) - mu_s(f0i));
+ FC_fv.back().mu_s = blend_mu(face_mu_s, mu_s(vi));
+
+ double face_mu_k = mu_k(f0i)
+ + FC_fv.back().closest_point[0] * (mu_k(f1i) - mu_k(f0i))
+ + FC_fv.back().closest_point[1] * (mu_k(f2i) - mu_k(f0i));
+ FC_fv.back().mu_k = blend_mu(face_mu_k, mu_k(vi));
}
}
diff --git a/src/ipc/collisions/tangential/tangential_collisions.hpp b/src/ipc/collisions/tangential/tangential_collisions.hpp
index ffb2b6be7..12de1a669 100644
--- a/src/ipc/collisions/tangential/tangential_collisions.hpp
+++ b/src/ipc/collisions/tangential/tangential_collisions.hpp
@@ -29,10 +29,25 @@ class TangentialCollisions {
const NormalPotential& normal_potential,
double normal_stiffness,
double mu)
+ {
+ this->build(
+ mesh, vertices, collisions, normal_potential, normal_stiffness, mu,
+ mu);
+ }
+
+ void build(
+ const CollisionMesh& mesh,
+ Eigen::ConstRef vertices,
+ const NormalCollisions& collisions,
+ const NormalPotential& normal_potential,
+ double normal_stiffness,
+ double mu_s,
+ double mu_k)
{
this->build(
mesh, vertices, collisions, normal_potential, normal_stiffness,
- Eigen::VectorXd::Constant(vertices.rows(), mu));
+ Eigen::VectorXd::Constant(vertices.rows(), mu_s),
+ Eigen::VectorXd::Constant(vertices.rows(), mu_k));
}
void build(
@@ -41,7 +56,8 @@ class TangentialCollisions {
const NormalCollisions& collisions,
const NormalPotential& normal_potential,
const double normal_stiffness,
- Eigen::ConstRef mus,
+ Eigen::ConstRef mu_k,
+ Eigen::ConstRef mu_s,
const std::function& blend_mu =
default_blend_mu);
diff --git a/src/ipc/friction/CMakeLists.txt b/src/ipc/friction/CMakeLists.txt
index c48f84f73..07cfd2205 100644
--- a/src/ipc/friction/CMakeLists.txt
+++ b/src/ipc/friction/CMakeLists.txt
@@ -1,6 +1,8 @@
set(SOURCES
smooth_friction_mollifier.cpp
smooth_friction_mollifier.hpp
+ smooth_mu.cpp
+ smooth_mu.hpp
)
target_sources(ipc_toolkit PRIVATE ${SOURCES})
\ No newline at end of file
diff --git a/src/ipc/friction/smooth_friction_mollifier.cpp b/src/ipc/friction/smooth_friction_mollifier.cpp
index e1a5da578..cc534d63b 100644
--- a/src/ipc/friction/smooth_friction_mollifier.cpp
+++ b/src/ipc/friction/smooth_friction_mollifier.cpp
@@ -20,7 +20,9 @@ double smooth_friction_f1(const double y, const double eps_v)
if (std::abs(y) >= eps_v) {
return 1;
}
- return y * (2 - y / eps_v) / eps_v;
+
+ const double y_over_eps_v = y / eps_v;
+ return y_over_eps_v * (2 - y_over_eps_v);
}
double smooth_friction_f2(const double y, const double eps_v)
diff --git a/src/ipc/friction/smooth_friction_mollifier.hpp b/src/ipc/friction/smooth_friction_mollifier.hpp
index e245432f8..a37104dae 100644
--- a/src/ipc/friction/smooth_friction_mollifier.hpp
+++ b/src/ipc/friction/smooth_friction_mollifier.hpp
@@ -55,6 +55,7 @@ double smooth_friction_f2(const double y, const double eps_v);
/// \end{cases}
/// \f\]
///
+/// @note The `x` in the function name refers to the parameter `y`.
/// @param y The tangential relative speed.
/// @param eps_v Velocity threshold below which static friction force is applied.
/// @return The value of the derivative of smooth_friction_f0 divided by y.
@@ -69,6 +70,7 @@ double smooth_friction_f1_over_x(const double y, const double eps_v);
/// \end{cases}
/// \f\]
///
+/// @note The `x` in the function name refers to the parameter `y`.
/// @param y The tangential relative speed.
/// @param eps_v Velocity threshold below which static friction force is applied.
/// @return The derivative of f1 times y minus f1 all divided by y cubed.
diff --git a/src/ipc/friction/smooth_mu.cpp b/src/ipc/friction/smooth_mu.cpp
new file mode 100644
index 000000000..8416aa9c2
--- /dev/null
+++ b/src/ipc/friction/smooth_mu.cpp
@@ -0,0 +1,95 @@
+#include "smooth_mu.hpp"
+
+#include
+
+#include
+#include
+
+namespace ipc {
+
+double smooth_mu(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ assert(eps_v > 0);
+ if (mu_s == mu_k || std::abs(y) >= eps_v) {
+ // 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;
+ }
+}
+
+double smooth_mu_derivative(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ assert(eps_v > 0);
+ if (mu_s == mu_k || std::abs(y) >= eps_v) {
+ // 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;
+ }
+}
+
+double smooth_mu_f0(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ assert(eps_v > 0);
+ if (mu_s == mu_k || std::abs(y) >= eps_v) {
+ // If the static and kinetic friction coefficients are equal, simplify.
+ 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;
+ }
+}
+
+double smooth_mu_f1(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ // This is a known formulation: μ(y) f₁(y)
+ return smooth_mu(y, mu_s, mu_k, eps_v) * smooth_friction_f1(y, eps_v);
+}
+
+double smooth_mu_f2(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ // Apply the chain rule:
+ return smooth_mu_derivative(y, mu_s, mu_k, eps_v)
+ * smooth_friction_f1(y, eps_v)
+ + smooth_mu(y, mu_s, mu_k, eps_v) * smooth_friction_f2(y, eps_v);
+}
+
+double smooth_mu_f1_over_x(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ // This is a known formulation: μ(y) f₁(y) / y
+ // where we use the robust division by y to avoid division by zero.
+ return smooth_mu(y, mu_s, mu_k, eps_v)
+ * smooth_friction_f1_over_x(y, eps_v);
+}
+
+double smooth_mu_f2_x_minus_mu_f1_over_x3(
+ const double y, const double mu_s, const double mu_k, const double eps_v)
+{
+ assert(eps_v > 0);
+ if (mu_s == mu_k || std::abs(y) >= eps_v) {
+ // 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);
+ }
+}
+
+} // namespace ipc
\ No newline at end of file
diff --git a/src/ipc/friction/smooth_mu.hpp b/src/ipc/friction/smooth_mu.hpp
new file mode 100644
index 000000000..1c7f7fd89
--- /dev/null
+++ b/src/ipc/friction/smooth_mu.hpp
@@ -0,0 +1,70 @@
+#pragma once
+
+namespace ipc {
+
+/// @brief Smooth coefficient from static to kinetic friction.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the μ at y.
+double smooth_mu(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the derivative of the smooth coefficient from static to kinetic friction.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the derivative at y.
+double smooth_mu_derivative(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth friction mollifier.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the integral at y.
+double smooth_mu_f0(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the value of the μ(y) f₁(y), where f₁ is the first derivative of the smooth friction mollifier.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the product at y.
+double smooth_mu_f1(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the value of d/dy (μ(y) f₁(y)), where f₁ is the first derivative of the smooth friction mollifier.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the derivative at y.
+double smooth_mu_f2(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the value of the μ(y) f₁(y) / y, where f₁ is the first derivative of the smooth friction mollifier.
+/// @note The `x` in the function name refers to the parameter `y`.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the product at y.
+double smooth_mu_f1_over_x(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+/// @brief Compute the value of the [(d/dy μ(y) f₁(y)) ⋅ y - μ(y) f₁(y)] / y³, where f₁ and f₂ are the first and second derivatives of the smooth friction mollifier.
+/// @note The `x` in the function name refers to the parameter `y`.
+/// @param y The tangential relative speed.
+/// @param mu_s Coefficient of static friction.
+/// @param mu_k Coefficient of kinetic friction.
+/// @param eps_v Velocity threshold below which static friction force is applied.
+/// @return The value of the expression at y.
+double smooth_mu_f2_x_minus_mu_f1_over_x3(
+ const double y, const double mu_s, const double mu_k, const double eps_v);
+
+} // namespace ipc
\ No newline at end of file
diff --git a/src/ipc/potentials/friction_potential.cpp b/src/ipc/potentials/friction_potential.cpp
index 747eddb4b..28a4a3f39 100644
--- a/src/ipc/potentials/friction_potential.cpp
+++ b/src/ipc/potentials/friction_potential.cpp
@@ -1,5 +1,7 @@
#include "friction_potential.hpp"
+#include
+
namespace ipc {
FrictionPotential::FrictionPotential(const double eps_v) : Super()
@@ -7,19 +9,22 @@ FrictionPotential::FrictionPotential(const double eps_v) : Super()
set_eps_v(eps_v);
}
-double FrictionPotential::f0(const double x) const
+double FrictionPotential::mu_f0(
+ const double x, const double mu_s, const double mu_k) const
{
- return smooth_friction_f0(x, eps_v());
+ return smooth_mu_f0(x, mu_s, mu_k, eps_v());
}
-double FrictionPotential::f1_over_x(const double x) const
+double FrictionPotential::mu_f1_over_x(
+ const double x, const double mu_s, const double mu_k) const
{
- return smooth_friction_f1_over_x(x, eps_v());
+ return smooth_mu_f1_over_x(x, mu_s, mu_k, eps_v());
}
-double FrictionPotential::f2_x_minus_f1_over_x3(const double x) const
+double FrictionPotential::mu_f2_x_minus_mu_f1_over_x3(
+ const double x, const double mu_s, const double mu_k) const
{
- return smooth_friction_f2_x_minus_f1_over_x3(x, eps_v());
+ return smooth_mu_f2_x_minus_mu_f1_over_x3(x, mu_s, mu_k, eps_v());
}
} // namespace ipc
\ No newline at end of file
diff --git a/src/ipc/potentials/friction_potential.hpp b/src/ipc/potentials/friction_potential.hpp
index c9aabdd44..cf0f1afea 100644
--- a/src/ipc/potentials/friction_potential.hpp
+++ b/src/ipc/potentials/friction_potential.hpp
@@ -25,10 +25,33 @@ class FrictionPotential : public TangentialPotential {
}
protected:
- double f0(const double x) const override;
- double f1_over_x(const double x) const override;
- double f2_x_minus_f1_over_x3(const double x) const override;
+ /// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the integral at x.
+ double
+ mu_f0(const double x, const double mu_s, const double mu_k) const override;
+ /// @brief Compute the value of the [μ(y) f₁(y)] / x, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the product at x.
+ double mu_f1_over_x(
+ const double x, const double mu_s, const double mu_k) const override;
+
+ /// @brief Compute the value of [(d/dx (μ(y) f₁(y))) x - μ(y) f₁(y)] / x³, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the derivative at x.
+ double mu_f2_x_minus_mu_f1_over_x3(
+ const double x, const double mu_s, const double mu_k) const override;
+
+ /// @brief Check if the potential is dynamic.
+ /// @param speed The tangential relative speed.
+ /// @return True if the potential is dynamic, false otherwise.
bool is_dynamic(const double speed) const override
{
return speed > eps_v();
diff --git a/src/ipc/potentials/tangential_adhesion_potential.cpp b/src/ipc/potentials/tangential_adhesion_potential.cpp
index 852689b82..8f70f0536 100644
--- a/src/ipc/potentials/tangential_adhesion_potential.cpp
+++ b/src/ipc/potentials/tangential_adhesion_potential.cpp
@@ -10,19 +10,22 @@ TangentialAdhesionPotential::TangentialAdhesionPotential(const double eps_a)
set_eps_a(eps_a);
}
-double TangentialAdhesionPotential::f0(const double x) const
+double TangentialAdhesionPotential::mu_f0(
+ const double x, const double mu_s, const double mu_k) const
{
- return tangential_adhesion_f0(x, eps_a());
+ return smooth_mu_a0(x, mu_s, mu_k, eps_a());
}
-double TangentialAdhesionPotential::f1_over_x(const double x) const
+double TangentialAdhesionPotential::mu_f1_over_x(
+ const double x, const double mu_s, const double mu_k) const
{
- return tangential_adhesion_f1_over_x(x, eps_a());
+ return smooth_mu_a1_over_x(x, mu_s, mu_k, eps_a());
}
-double TangentialAdhesionPotential::f2_x_minus_f1_over_x3(const double x) const
+double TangentialAdhesionPotential::mu_f2_x_minus_mu_f1_over_x3(
+ const double x, const double mu_s, const double mu_k) const
{
- return tangential_adhesion_f2_x_minus_f1_over_x3(x, eps_a());
+ return smooth_mu_a2_x_minus_mu_a1_over_x3(x, mu_s, mu_k, eps_a());
}
} // namespace ipc
\ No newline at end of file
diff --git a/src/ipc/potentials/tangential_adhesion_potential.hpp b/src/ipc/potentials/tangential_adhesion_potential.hpp
index e21f2c336..c3290e4d4 100644
--- a/src/ipc/potentials/tangential_adhesion_potential.hpp
+++ b/src/ipc/potentials/tangential_adhesion_potential.hpp
@@ -25,10 +25,33 @@ class TangentialAdhesionPotential : public TangentialPotential {
}
protected:
- double f0(const double x) const override;
- double f1_over_x(const double x) const override;
- double f2_x_minus_f1_over_x3(const double x) const override;
+ /// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static adhesion.
+ /// @param mu_k Coefficient of kinetic adhesion.
+ /// @return The value of the integral at x.
+ double
+ mu_f0(const double x, const double mu_s, const double mu_k) const override;
+ /// @brief Compute the value of the [μ(y) f₁(y)] / x, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static adhesion.
+ /// @param mu_k Coefficient of kinetic adhesion.
+ /// @return The value of the product at x.
+ double mu_f1_over_x(
+ const double x, const double mu_s, const double mu_k) const override;
+
+ /// @brief Compute the value of [(d/dx (μ(y) f₁(y))) x - μ(y) f₁(y)] / x³, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static adhesion.
+ /// @param mu_k Coefficient of kinetic adhesion.
+ /// @return The value of the derivative at x.
+ double mu_f2_x_minus_mu_f1_over_x3(
+ const double x, const double mu_s, const double mu_k) const override;
+
+ /// @brief Check if the potential is dynamic.
+ /// @param speed The tangential relative speed.
+ /// @return True if the potential is dynamic, false otherwise.
bool is_dynamic(const double speed) const override
{
return speed > eps_a();
diff --git a/src/ipc/potentials/tangential_potential.cpp b/src/ipc/potentials/tangential_potential.cpp
index 18ab0c03d..be34c9e84 100644
--- a/src/ipc/potentials/tangential_potential.cpp
+++ b/src/ipc/potentials/tangential_potential.cpp
@@ -140,22 +140,33 @@ double TangentialPotential::operator()(
const TangentialCollision& collision,
Eigen::ConstRef velocities) const
{
- // μ N(xᵗ) f₀(‖u‖) (where u = T(xᵗ)ᵀv)
+ //
+ // μₛ = μₖ:
+ // μ N(xᵗ) f₀(‖u‖) (where u = T(xᵗ)ᵀv)
+ //
+ // μₛ != μₖ:
+ // N(xᵗ) g₀(‖u‖) (where g₀(x) = ∫ μ(x) f₀(x) dx)
+ //
// Compute u = PᵀΓv
const VectorMax2d u = collision.tangent_basis.transpose()
* collision.relative_velocity(velocities);
- return collision.weight * collision.mu * collision.normal_force_magnitude
- * f0(u.norm());
+ return collision.weight * collision.normal_force_magnitude
+ * mu_f0(u.norm(), collision.mu_s, collision.mu_k);
}
VectorMax12d TangentialPotential::gradient(
const TangentialCollision& collision,
Eigen::ConstRef velocities) const
{
- // ∇ₓ μ N(xᵗ) f₀(‖u‖) (where u = T(xᵗ)ᵀv)
- // = μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u
+ //
+ // μₛ = μₖ:
+ // ∇ᵥ μ N(xᵗ) f₀(‖u‖) = μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u
+ //
+ // μₛ != μₖ:
+ // ∇ᵥ N(xᵗ) g₀(‖u‖)/‖u‖ T(xᵗ) u = N(xᵗ) g₁(‖u‖)/‖u‖ T(xᵗ) u
+ //
// Compute u = PᵀΓv
const VectorMax2d u = collision.tangent_basis.transpose()
@@ -166,14 +177,14 @@ VectorMax12d TangentialPotential::gradient(
collision.relative_velocity_matrix().transpose()
* collision.tangent_basis;
- // Compute f₁(‖u‖)/‖u‖
- const double f1_over_norm_u = f1_over_x(u.norm());
+ // Compute μ(‖u‖) f₁(‖u‖)/‖u‖
+ const double mu_f1_over_norm_u =
+ mu_f1_over_x(u.norm(), collision.mu_s, collision.mu_k);
- // μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u ∈ ℝⁿ
- // (n×2)(2×1) = (n×1)
+ // μ(‖u‖) N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u ∈ (n×2)(2×1) = (n×1)
return T
- * ((collision.weight * collision.mu * collision.normal_force_magnitude
- * f1_over_norm_u)
+ * ((collision.weight * collision.normal_force_magnitude
+ * mu_f1_over_norm_u)
* u);
}
@@ -182,9 +193,15 @@ MatrixMax12d TangentialPotential::hessian(
Eigen::ConstRef velocities,
const PSDProjectionMethod project_hessian_to_psd) const
{
- // ∇ₓ μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u (where u = T(xᵗ)ᵀ v)
- // = μ N T [(f₂(‖u‖)‖u‖ − f₁(‖u‖))/‖u‖³ uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
- // = μ N T [f₂(‖u‖) uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
+ //
+ // μₛ = μₖ:
+ // ∇ᵥ μ N(xᵗ) f₁(‖u‖)/‖u‖ T(xᵗ) u
+ // = μ N T [(f₁'(‖u‖)‖u‖ − f₁(‖u‖))/‖u‖³ uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
+ // = μ N T [f₂(‖u‖) uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
+ //
+ // μₛ != μₖ:
+ // ∇ᵥ N(xᵗ) g₁(‖u‖)/‖u‖ T(xᵗ) u = N T [g₂(‖u‖) uuᵀ + g₁(‖u‖)/‖u‖ I] Tᵀ
+ //
// Compute u = PᵀΓv
const VectorMax2d u = collision.tangent_basis.transpose()
@@ -198,12 +215,12 @@ MatrixMax12d TangentialPotential::hessian(
// Compute ‖u‖
const double norm_u = u.norm();
- // Compute f₁(‖u‖)/‖u‖
- const double f1_over_norm_u = f1_over_x(norm_u);
+ // Compute μ(‖u‖) f₁(‖u‖)/‖u‖
+ const double mu_f1_over_norm_u =
+ mu_f1_over_x(norm_u, collision.mu_s, collision.mu_k);
- // Compute μ N(xᵗ)
- const double scale =
- collision.weight * collision.mu * collision.normal_force_magnitude;
+ // Compute N(xᵗ)
+ const double scale = collision.weight * collision.normal_force_magnitude;
MatrixMax12d hess;
if (is_dynamic(norm_u)) {
@@ -222,7 +239,7 @@ MatrixMax12d TangentialPotential::hessian(
// I - uuᵀ/‖u‖² = ūūᵀ / ‖u‖² (where ū⋅u = 0)
const Eigen::Vector2d u_perp(-u[1], u[0]);
hess = // grouped to reduce number of operations
- (T * ((scale * f1_over_norm_u / (norm_u * norm_u)) * u_perp))
+ (T * ((scale * mu_f1_over_norm_u / (norm_u * norm_u)) * u_perp))
* (u_perp.transpose() * T.transpose());
}
} else if (norm_u == 0) {
@@ -232,15 +249,16 @@ MatrixMax12d TangentialPotential::hessian(
if (project_hessian_to_psd != PSDProjectionMethod::NONE && scale <= 0) {
hess.setZero(collision.ndof(), collision.ndof()); // -PSD = NSD ⟹ 0
} else {
- hess = scale * f1_over_norm_u * T * T.transpose();
+ hess = scale * mu_f1_over_norm_u * T * T.transpose();
}
} else {
// ∇²D(v) = μ N T [f₂(‖u‖) uuᵀ + f₁(‖u‖)/‖u‖ I] Tᵀ
// ⟹ only need to project the inner 2x2 matrix to PSD
- const double f2 = f2_x_minus_f1_over_x3(norm_u);
+ const double f2 =
+ mu_f2_x_minus_mu_f1_over_x3(norm_u, collision.mu_s, collision.mu_k);
MatrixMax2d inner_hess = f2 * u * u.transpose();
- inner_hess.diagonal().array() += f1_over_norm_u;
+ inner_hess.diagonal().array() += mu_f1_over_norm_u;
inner_hess *= scale; // NOTE: negative scaling will be projected out
inner_hess = project_to_psd(inner_hess, project_hessian_to_psd);
@@ -297,13 +315,14 @@ VectorMax12d TangentialPotential::force(
// Compute τ = PᵀΓv
const VectorMax2d tau = T.transpose() * velocities;
- // Compute f₁(‖τ‖)/‖τ‖
- const double f1_over_norm_tau = f1_over_x(tau.norm());
+ // Compute μ(‖τ‖) f₁(‖τ‖)/‖τ‖
+ const double mu_s = no_mu ? 1.0 : collision.mu_s;
+ const double mu_k = no_mu ? 1.0 : collision.mu_k;
+ const double mu_f1_over_norm_tau = mu_f1_over_x(tau.norm(), mu_s, mu_k);
// F = -μ N f₁(‖τ‖)/‖τ‖ T τ
// NOTE: no_mu -> leave mu out of this function (i.e., assuming mu = 1)
- return -collision.weight * (no_mu ? 1.0 : collision.mu) * N
- * f1_over_norm_tau * T * tau;
+ return -collision.weight * N * mu_f1_over_norm_tau * T * tau;
}
MatrixMax12d TangentialPotential::force_jacobian(
@@ -414,20 +433,22 @@ MatrixMax12d TangentialPotential::force_jacobian(
jac_tau = T.transpose(); // Tᵀ ∇ᵥv = Tᵀ
}
- // Compute f₁(‖τ‖)/‖τ‖
+ // Compute μ f₁(‖τ‖)/‖τ‖
const double tau_norm = tau.norm();
- const double f1_over_norm_tau = f1_over_x(tau_norm);
+ const double mu_f1_over_norm_tau =
+ mu_f1_over_x(tau_norm, collision.mu_s, collision.mu_k);
- // Compute ∇(f₁(‖τ‖)/‖τ‖)
- VectorMax12d grad_f1_over_norm_tau;
+ // Compute ∇(μ f₁(‖τ‖)/‖τ‖)
+ VectorMax12d grad_mu_f1_over_norm_tau;
if (tau_norm == 0) {
// lim_{x→0} f₂(x)x² = 0
- grad_f1_over_norm_tau.setZero(n);
+ grad_mu_f1_over_norm_tau.setZero(n);
} else {
// ∇ (f₁(‖τ‖)/‖τ‖) = (f₂(‖τ‖)‖τ‖ - f₁(‖τ‖)) / ‖τ‖³ τᵀ ∇τ
- double f2 = f2_x_minus_f1_over_x3(tau_norm);
+ double f2 = mu_f2_x_minus_mu_f1_over_x3(
+ tau_norm, collision.mu_s, collision.mu_k);
assert(std::isfinite(f2));
- grad_f1_over_norm_tau = f2 * tau.transpose() * jac_tau;
+ grad_mu_f1_over_norm_tau = f2 * tau.transpose() * jac_tau;
}
// Premultiplied values
@@ -439,15 +460,15 @@ MatrixMax12d TangentialPotential::force_jacobian(
// = -μ f₁(‖τ‖)/‖τ‖ (T τ) [∇N]ᵀ
if (need_jac_N_or_T) {
- J = f1_over_norm_tau * T_times_tau * grad_N.transpose();
+ J = mu_f1_over_norm_tau * T_times_tau * grad_N.transpose();
}
- // + -μ N T τ [∇(f₁(‖τ‖)/‖τ‖)]
- J += N * T_times_tau * grad_f1_over_norm_tau.transpose();
+ // + -N T τ [∇(μ f₁(‖τ‖)/‖τ‖)]
+ J += N * T_times_tau * grad_mu_f1_over_norm_tau.transpose();
// + -μ N f₁(‖τ‖)/‖τ‖ [∇T] τ
if (need_jac_N_or_T) {
- const VectorMax2d scaled_tau = N * f1_over_norm_tau * tau;
+ const VectorMax2d scaled_tau = N * mu_f1_over_norm_tau * tau;
for (int i = 0; i < n; i++) {
// ∂J/∂xᵢ = ∂T/∂xᵢ * τ
J.col(i) += jac_T.middleRows(i * n, n) * scaled_tau;
@@ -455,14 +476,14 @@ MatrixMax12d TangentialPotential::force_jacobian(
}
// + -μ N f₁(‖τ‖)/‖τ‖ T [∇τ]
- J += N * f1_over_norm_tau * T * jac_tau;
+ J += N * mu_f1_over_norm_tau * T * jac_tau;
// NOTE: ∇ₓw(x) is not local to the collision pair (i.e., it involves more
- // than the 4 collisioning vertices), so we do not have enough information
+ // than the 4 colliding vertices), so we do not have enough information
// here to compute the gradient. Instead this should be handled outside of
- // the function. For a simple multiplicitive model (∑ᵢ wᵢ Fᵢ) this can be
+ // the function. For a simple multiplicative model (∑ᵢ wᵢ Fᵢ) this can be
// done easily.
- J *= -collision.weight * collision.mu;
+ J *= -collision.weight;
return J;
}
diff --git a/src/ipc/potentials/tangential_potential.hpp b/src/ipc/potentials/tangential_potential.hpp
index 379590a9c..b09f7cccb 100644
--- a/src/ipc/potentials/tangential_potential.hpp
+++ b/src/ipc/potentials/tangential_potential.hpp
@@ -141,9 +141,33 @@ class TangentialPotential : public Potential {
const double dmin = 0) const;
protected:
- virtual double f0(const double x) const = 0;
- virtual double f1_over_x(const double x) const = 0;
- virtual double f2_x_minus_f1_over_x3(const double x) const = 0;
+ /// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the integral at x.
+ virtual double
+ mu_f0(const double x, const double mu_s, const double mu_k) const = 0;
+
+ /// @brief Compute the value of the [μ(y) f₁(y)] / x, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the product at x.
+ virtual double mu_f1_over_x(
+ const double x, const double mu_s, const double mu_k) const = 0;
+
+ /// @brief Compute the value of [(d/dx (μ(y) f₁(y))) x - μ(y) f₁(y)] / x³, where f₁ is the first derivative of the smooth mollifier.
+ /// @param x The tangential relative speed.
+ /// @param mu_s Coefficient of static friction.
+ /// @param mu_k Coefficient of kinetic friction.
+ /// @return The value of the derivative at x.
+ virtual double mu_f2_x_minus_mu_f1_over_x3(
+ const double x, const double mu_s, const double mu_k) const = 0;
+
+ /// @brief Check if the speed is dynamic.
+ /// @param speed The tangential relative speed.
+ /// @return True if the speed is dynamic, false otherwise.
virtual bool is_dynamic(const double speed) const = 0;
};
diff --git a/tests/src/tests/adhesion/test_adhesion.cpp b/tests/src/tests/adhesion/test_adhesion.cpp
index 39d5a2667..064e396b3 100644
--- a/tests/src/tests/adhesion/test_adhesion.cpp
+++ b/tests/src/tests/adhesion/test_adhesion.cpp
@@ -139,4 +139,67 @@ TEST_CASE("Tangential adhesion mollifier", "[adhesion][tangential]")
f2 /= x;
CHECK(f2 == Catch::Approx(fd_f2[0]).margin(MARGIN).epsilon(EPSILON));
+}
+
+TEST_CASE("Tangential smooth mu adhesion mollifier", "[adhesion][tangential]")
+{
+ static constexpr double EPSILON = 1e-4;
+ static constexpr double MARGIN = 1e-6;
+ // NOTE: Use h=1e-10 for finite difference because min eps_a=1e-8
+ static constexpr double H = 1e-10;
+ const double mu_s = GENERATE(range(0.0, 1.0, 0.1));
+ const double mu_k = GENERATE(range(0.0, 1.0, 0.1));
+ const double eps_a = std::pow(10, GENERATE(range(-8, 0, 1)));
+ const double x = std::pow(10, GENERATE(range(-8, 0, 1)));
+
+ if (x == 1e-8 && eps_a == 1e-8) {
+ return;
+ }
+
+ CAPTURE(x, eps_a, x / eps_a, mu_s, mu_k);
+
+ Eigen::Matrix X;
+ X << x;
+
+ // Check gradient
+
+ Eigen::VectorXd fd_a1(1);
+ fd::finite_gradient(
+ X,
+ [&](const Eigen::VectorXd& _X) {
+ return smooth_mu_a0(_X[0], mu_s, mu_k, eps_a);
+ },
+ fd_a1, fd::AccuracyOrder::SECOND, H);
+
+ CHECK(
+ smooth_mu_a1(x, mu_s, mu_k, eps_a)
+ == Catch::Approx(fd_a1[0]).margin(MARGIN).epsilon(EPSILON));
+
+ CHECK(
+ smooth_mu_a1_over_x(x, mu_s, mu_k, eps_a) * x
+ == Catch::Approx(fd_a1[0]).margin(MARGIN).epsilon(EPSILON));
+
+ // Check hessian
+ if (x == eps_a) {
+ return;
+ }
+
+ Eigen::VectorXd fd_a2(1);
+ fd::finite_gradient(
+ X,
+ [&](const Eigen::VectorXd& _X) {
+ return smooth_mu_a1(_X[0], mu_s, mu_k, eps_a);
+ },
+ fd_a2, fd::AccuracyOrder::SECOND, H);
+
+ CHECK(
+ smooth_mu_a2(x, mu_s, mu_k, eps_a)
+ == Catch::Approx(fd_a2[0]).margin(MARGIN).epsilon(EPSILON));
+
+ double a2 = smooth_mu_a2_x_minus_mu_a1_over_x3(x, mu_s, mu_k, eps_a);
+ a2 *= x * x * x;
+ a2 += smooth_mu_a1(x, mu_s, mu_k, eps_a);
+ a2 /= x;
+
+ CHECK(a2 == Catch::Approx(fd_a2[0]).margin(MARGIN).epsilon(EPSILON));
}
\ No newline at end of file
diff --git a/tests/src/tests/friction/CMakeLists.txt b/tests/src/tests/friction/CMakeLists.txt
index 96bb01b64..50245f75f 100644
--- a/tests/src/tests/friction/CMakeLists.txt
+++ b/tests/src/tests/friction/CMakeLists.txt
@@ -3,6 +3,7 @@ set(SOURCES
test_force_jacobian.cpp
test_friction.cpp
test_smooth_friction_mollifier.cpp
+ test_smooth_mu.cpp
# Benchmarks
diff --git a/tests/src/tests/friction/test_friction.cpp b/tests/src/tests/friction/test_friction.cpp
index 41979f00b..8a82f57dd 100644
--- a/tests/src/tests/friction/test_friction.cpp
+++ b/tests/src/tests/friction/test_friction.cpp
@@ -147,7 +147,8 @@ bool read_ipc_friction_data(
E, F, mmcvids, lambda, coords, bases, tangential_collisions);
tests::mmcvids_to_collisions(E, F, mmcvids, collisions);
for (int i = 0; i < tangential_collisions.size(); i++) {
- tangential_collisions[i].mu = mu;
+ tangential_collisions[i].mu_s = mu;
+ tangential_collisions[i].mu_k = mu;
}
return true;
@@ -226,7 +227,8 @@ TEST_CASE(
CHECK(
collision.normal_force_magnitude
== Catch::Approx(expected_collision.normal_force_magnitude));
- CHECK(collision.mu == Catch::Approx(expected_collision.mu));
+ CHECK(collision.mu_s == Catch::Approx(expected_collision.mu_s));
+ CHECK(collision.mu_k == Catch::Approx(expected_collision.mu_k));
CHECK(
collision.vertex_ids(E, F) == expected_collision.vertex_ids(E, F));
diff --git a/tests/src/tests/friction/test_smooth_mu.cpp b/tests/src/tests/friction/test_smooth_mu.cpp
new file mode 100644
index 000000000..ad8321e2e
--- /dev/null
+++ b/tests/src/tests/friction/test_smooth_mu.cpp
@@ -0,0 +1,85 @@
+#include
+#include
+#include
+
+#include
+
+#include
+
+using namespace ipc;
+
+TEST_CASE("Smooth mu", "[friction][mollifier][mu]")
+{
+
+ static constexpr double EPSILON = 1e-4;
+ static constexpr double MARGIN = 1e-6;
+ // NOTE: Use h=1e-10 for finite difference because min eps_v=1e-8
+ static constexpr double H = 1e-10;
+ const double mu_s = GENERATE(range(0.0, 1.0, 0.1));
+ const double mu_k = GENERATE(range(0.0, 1.0, 0.1));
+ const double eps_v = std::pow(10, GENERATE(range(-8, 0, 1)));
+ const double x = std::pow(10, GENERATE(range(-8, 0, 1)));
+
+ if (x == 1e-8 && eps_v == 1e-8) {
+ return;
+ }
+
+ CAPTURE(x, eps_v, x / eps_v, mu_s, mu_k);
+
+ Eigen::Matrix X;
+ X << x;
+
+ // Check gradient
+
+ Eigen::VectorXd fd_f1(1);
+ fd::finite_gradient(
+ X,
+ [&](const Eigen::VectorXd& _X) {
+ return smooth_mu_f0(_X[0], mu_s, mu_k, eps_v);
+ },
+ fd_f1, fd::AccuracyOrder::SECOND, H);
+
+ CHECK(
+ smooth_mu_f1(x, mu_s, mu_k, eps_v)
+ == Catch::Approx(fd_f1[0]).margin(MARGIN).epsilon(EPSILON));
+
+ CHECK(
+ smooth_mu_f1_over_x(x, mu_s, mu_k, eps_v) * x
+ == Catch::Approx(fd_f1[0]).margin(MARGIN).epsilon(EPSILON));
+
+ // Check hessian
+ if (x == eps_v) {
+ return;
+ }
+
+ Eigen::VectorXd fd_mu1(1);
+ fd::finite_gradient(
+ X,
+ [&](const Eigen::VectorXd& _X) {
+ return smooth_mu(_X[0], mu_s, mu_k, eps_v);
+ },
+ fd_mu1, fd::AccuracyOrder::SECOND, H);
+
+ CHECK(
+ smooth_mu_derivative(x, mu_s, mu_k, eps_v)
+ == Catch::Approx(fd_mu1[0]).margin(MARGIN).epsilon(EPSILON));
+
+ Eigen::VectorXd fd_f2(1);
+ fd::finite_gradient(
+ X,
+ [&](const Eigen::VectorXd& _X) {
+ return smooth_mu_f1(_X[0], mu_s, mu_k, eps_v);
+ },
+ fd_f2, fd::AccuracyOrder::SECOND, H);
+
+ CHECK(
+ smooth_mu_f2(x, mu_s, mu_k, eps_v)
+ == Catch::Approx(fd_f2[0]).margin(MARGIN).epsilon(EPSILON));
+
+ double f2 = smooth_mu_f2_x_minus_mu_f1_over_x3(x, mu_s, mu_k, eps_v);
+ f2 *= x * x * x;
+ f2 += smooth_mu_f1(x, mu_s, mu_k, eps_v);
+ f2 /= x;
+
+ CHECK(f2 == Catch::Approx(fd_f2[0]).margin(MARGIN).epsilon(EPSILON));
+}
\ No newline at end of file