Repository: ipc-sim/ipc-toolkit (language: python)

Issue:
Possibly inconsistent results for convergent formulation with respect to dhat?
**Describe the bug**
I'm using `CollisionConstraints::set_use_convergent_formulation(true)` to obtain quantities that are less tessellation-dependent and that - to my understanding - vary less with the value of `dhat`.

I am not sure if perhaps I've misunderstood the true purpose of this formulation, but I'm observing inconsistent values for different values of `dhat` in a setting in which I would expect the values to be more similar.

Consider the following example. Take a squished unit tetrahedron with vertices `(0, 0, 0)`, `(1, 0, 0)`, `(0, 1, 0)` and `(0, 0, 0.5 * dhat)`. Since the z-coordinate of the "top" vertex is proportional to `dhat`, I would expect a similar potential energy value for different values of `dhat` (noting of course that changing `dhat` for this example also changes the area of the triangles).

Here's the output I get for different values of `dhat` for this example:

```
dhat = 0.0001
  Potential: 5.26234e-09
  Scaled potential: 0.526234
dhat = 0.001
  Potential: 5.26234e-07
  Scaled potential: 0.526234
dhat = 0.01
  Potential: 5.26251e-05
  Scaled potential: 0.526251
```

Since the inputs to the barrier function are squared, I believe the "real" barrier contribution works out to `-log(d^2/dhat^2) (d^2 - dhat^2)^2`. To make this unitless, we'd have to divide by `dhat^4`, however my understanding is that the `weight` constant currently used to scale the barrier function is divided by `dhat^2` (when `dmin = 0`). The `Scaled potential` in the above output refers to dividing the potential by `dhat^2`, to make up for the missing powers.

**To Reproduce**
Here's a self-contained minimal example that I used to produce the output above:
```c++
#include <ipc/ipc.hpp>
#include <iostream>
#include <array>

int main() {
    const auto dhat_values = std::array<double, 3> { 1e-4, 1e-3, 1e-2 };

    for (const auto dhat : dhat_values) {
        Eigen::MatrixXd rest_vertices(4, 3);
        rest_vertices <<
                      0.0, 0.0, 0.0,
                1.0, 0.0, 0.0,
                0.0, 1.0, 0.0,
                0.0, 0.0, 1.0;

        Eigen::MatrixXd deformed_vertices = rest_vertices;
        deformed_vertices.row(3)(2) = 0.5 * dhat;

        Eigen::MatrixXi edges(6, 2);
        edges << 0, 1,
                0, 2,
                0, 3,
                1, 2,
                1, 3,
                2, 3;
        Eigen::MatrixXi faces(4, 3);
        faces << 0, 2, 1,
                0, 1, 3,
                0, 3, 2,
                1, 2, 3;

        ipc::CollisionMesh mesh = ipc::CollisionMesh::build_from_full_mesh(rest_vertices, edges, faces);
        ipc::CollisionConstraints constraints;
        constraints.set_use_convergent_formulation(true);
        constraints.build(mesh, deformed_vertices, dhat);

        const double potential = constraints.compute_potential(mesh, deformed_vertices, dhat);
        std::cout << "dhat = " << dhat << std::endl;
        std::cout << "  Potential: " << potential << std::endl;
        std::cout << "  Scaled potential: " << potential / (dhat * dhat) << std::endl;
    }

    return 0;
}
```

**Expected behavior**
I had expected to get the value of the "scaled potential" above when using the "convergent formulation". This might however be due to a misunderstanding on my part.

**Additional context**
This is not a big deal on my part, because I can work around it by dividing by `dhat^2`, but I thought I ought to let you know, in case this is an issue that has gone undiscovered.



Discussion / hints:
The intention of this change is to be able to express the barrier potential in the same physical units as the elasticity potential (i.e., Pa $\cdot$ m) (more on this at the bottom). With non-squared distances, we do this by writing the barrier function as

$$
b_\text{new}(d, \hat{d}) = -\kappa\hat{d}\left(\frac{d}{\hat{d}} - 1\right)^2 \log\left(\frac{d}{\hat{d}}\right).
$$

The relationship between $b_\text{old}(d, \hat{d})$ and $b_\text{new}(d, \hat{d})$ is 

$$
b_\text{new}(d, \hat{d}) = \frac{1}{\hat{d}} b_\text{old}(d, \hat{d}).
$$

In the toolkit, we only implement $b_\text{old}(d, \hat{d})$ and store this division in the `weight` for simplicity of implementation.

However, we are using distances squared so we actually want to compute

$$
b_\text{new}(d^2, \hat{d}^2) = -\kappa\hat{d}\left(\frac{d^2}{\hat{d}^2} - 1\right)^2 \log\left(\frac{d^2}{\hat{d}^2}\right).
$$

Again the relationship between $b_\text{old}(d^2, \hat{d}^2)$ and $b_\text{new}(d^2, \hat{d}^2)$ is 

$$
b_\text{new}(d, \hat{d}) = \frac{1}{\hat{d}^3} b_\text{old}(d, \hat{d}).
$$

Again, we only implement $b_\text{old}(d, \hat{d})$, so we want to store this division in the `weight` for simplicity of implementation. However, I am dividing by $\hat{d}^2$ when I should be dividing by $\hat{d}^3$ to get the correct units. I will make this change.

This will still not produce the same behavior as you expect (because we multiply the normalized barrier by $\hat{d}$), but this is the correct implementation of the algorithm we want (i.e., physically-based units for the barrier term). Alternatively, if you wanted to use the normalized barrier, you could assume the units of $\kappa$ are Pa $\cdot$ m. What I like about the intended implementation is that $\kappa$ has units of Pa and we can express $\kappa$ as a ratio of the materials Young's modulus ($0.1E$ works well in many examples). 

The intention is to treat the barriers as a thin elastic region around the mesh, and having consistent units makes it easier to pick the stiffness for this "material".

Fix patch:
diff --git a/notebooks/physical_barrier.ipynb b/notebooks/physical_barrier.ipynb
new file mode 100644
index 000000000..5e256f2b0
--- /dev/null
+++ b/notebooks/physical_barrier.ipynb
@@ -0,0 +1,238 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sympy\n",
+    "from generate_cpp_code import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dhat = sympy.Symbol('\\hat{d}')\n",
+    "d = sympy.Symbol('d')\n",
+    "dmin = sympy.Symbol('d_{\\\\text{min}}')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle - \\left(- \\hat{d} + d\\right)^{2} \\log{\\left(\\frac{d}{\\hat{d}} \\right)}$"
+      ],
+      "text/plain": [
+       "-(-\\hat{d} + d)**2*log(d/\\hat{d})"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle - \\left(-1 + \\frac{d}{\\hat{d}}\\right)^{2} \\log{\\left(\\frac{d}{\\hat{d}} \\right)}$"
+      ],
+      "text/plain": [
+       "-(-1 + d/\\hat{d})**2*log(d/\\hat{d})"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle - \\hat{d} \\left(-1 + \\frac{d}{\\hat{d}}\\right)^{2} \\log{\\left(\\frac{d}{\\hat{d}} \\right)}$"
+      ],
+      "text/plain": [
+       "-\\hat{d}*(-1 + d/\\hat{d})**2*log(d/\\hat{d})"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "def barrier(d, dhat):\n",
+    "    return -(d - dhat)**2 * sympy.log(d/dhat)\n",
+    "\n",
+    "def normalized_barrier(d, dhat):\n",
+    "    return -(d/dhat - 1)**2 * sympy.log(d/dhat)\n",
+    "\n",
+    "def physical_barrier(d, p_dhat):\n",
+    "    return dhat * normalized_barrier(d, p_dhat)\n",
+    "\n",
+    "display(barrier(d, dhat), normalized_barrier(d, dhat), physical_barrier(d, dhat))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{1}{\\hat{d}}$"
+      ],
+      "text/plain": [
+       "1/\\hat{d}"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(physical_barrier(d, dhat) / barrier(d, dhat)).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{1}{\\hat{d}^{3}}$"
+      ],
+      "text/plain": [
+       "\\hat{d}**(-3)"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(physical_barrier(d**2, dhat**2) / barrier(d**2, dhat**2)).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{1}{\\hat{d} \\left(\\hat{d} + 2 d_{\\text{min}}\\right)^{2}}$"
+      ],
+      "text/plain": [
+       "1/(\\hat{d}*(\\hat{d} + 2*d_{\\text{min}})**2)"
+      ]
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "d_param = d**2 - dmin**2\n",
+    "dhat_param = 2 * dmin * dhat + dhat**2\n",
+    "(physical_barrier(d_param, dhat_param) / barrier(d_param, dhat_param)).simplify()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "const auto t0 = d/dhat;\n",
+      "r = -std::pow(1 - t0, 2)*std::log(t0);\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(generate_code(normalized_barrier(d, dhat), \"r\").replace(\"\\hat{d}\", \"dhat\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "const auto t0 = 1.0/dhat;\n",
+      "const auto t1 = d*t0;\n",
+      "const auto t2 = 1 - t1;\n",
+      "r = t2*(2*t0*std::log(t1) - t2/d);\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(generate_code(normalized_barrier(d, dhat).diff(d), \"r\").replace(\"\\hat{d}\", \"dhat\"))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "const auto t0 = 1.0/dhat;\n",
+      "const auto t1 = d*t0;\n",
+      "const auto t2 = 1 - t1;\n",
+      "r = 4*t0*t2/d + std::pow(t2, 2)/std::pow(d, 2) - 2*std::log(t1)/std::pow(dhat, 2);\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(generate_code(normalized_barrier(d, dhat).diff(d).diff(d), \"r\").replace(\"\\hat{d}\", \"dhat\"))"
+   ]
+  },
+  {
+   "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.11.2"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/src/ipc/collisions/collision_constraints.cpp b/src/ipc/collisions/collision_constraints.cpp
index ffd961e65..92f6dbc8c 100644
--- a/src/ipc/collisions/collision_constraints.cpp
+++ b/src/ipc/collisions/collision_constraints.cpp
@@ -76,18 +76,25 @@ void CollisionConstraints::build(
 
     CollisionConstraintsBuilder::merge(storage, *this);
 
-    // This is the dhat that is used in the barrier potential (because we use
-    // squared distances).
-    const double adjusted_dhat = 2 * dmin * dhat + dhat * dhat;
-
     for (size_t ci = 0; ci < size(); ci++) {
         CollisionConstraint& constraint = (*this)[ci];
         constraint.minimum_distance = dmin;
-        if (use_convergent_formulation()) {
-            // Divide by dhat to equivalently use the "physical" barrier
-            constraint.weight /= adjusted_dhat;
+    }
+
+    if (use_convergent_formulation()) {
+        // NOTE: When using the convergent formulation we want the barrier to
+        // have units of Pa⋅m, so κ gets units of Pa and the barrier function
+        // should have units of m. See notebooks/physical_barrier.ipynb for more
+        // details.
+        const double barrier_to_physical_barrier_divisor =
+            dhat * std::pow(dhat + 2 * dmin, 2);
+
+        for (size_t ci = 0; ci < size(); ci++) {
+            CollisionConstraint& constraint = (*this)[ci];
+            constraint.weight /= barrier_to_physical_barrier_divisor;
             if (are_shape_derivatives_enabled()) {
-                constraint.weight_gradient /= adjusted_dhat;
+                constraint.weight_gradient /=
+                    barrier_to_physical_barrier_divisor;
             }
         }
     }


Test patch:
diff --git a/tests/barrier/test_barrier.cpp b/tests/barrier/test_barrier.cpp
index 98ecc01cc..a82aea7ba 100644
--- a/tests/barrier/test_barrier.cpp
+++ b/tests/barrier/test_barrier.cpp
@@ -6,11 +6,8 @@
 #include <ipc/barrier/barrier.hpp>
 
 namespace {
-double physical_barrier(const double d, const double dhat)
+double normalized_barrier(const double d, const double dhat)
 {
-    // units(d) = m and units(d̂) = m ⟹ units(b(d)) = m
-    // units(κ) = Pa ⟹ units(κ b(d)) = Pa m
-
     if (d <= 0.0) {
         return std::numeric_limits<double>::infinity();
     }
@@ -18,35 +15,40 @@ double physical_barrier(const double d, const double dhat)
         return 0;
     }
 
-    // b(d) = -d̂(d/d̂-1)²ln(d / d̂)
-    const double d_over_dhat = d / dhat;
-    const double d_over_dhat_minus_1 = d_over_dhat - 1;
-    return -dhat * d_over_dhat_minus_1 * d_over_dhat_minus_1 * log(d_over_dhat);
+    // b(d) = -(d/d̂-1)²ln(d / d̂)
+    const auto t0 = d / dhat;
+    return -std::pow(1 - t0, 2) * std::log(t0);
 }
 
-double physical_barrier_gradient(const double d, const double dhat)
+double normalized_barrier_gradient(const double d, const double dhat)
 {
     if (d <= 0.0 || d >= dhat) {
         return 0.0;
     }
-    return (d - dhat) * (-2 * d * log(d / dhat) - d + dhat) / (d * dhat);
+    const double t0 = 1.0 / dhat;
+    const double t1 = d * t0;
+    const double t2 = 1 - t1;
+    return t2 * (2 * t0 * std::log(t1) - t2 / d);
 }
 
-double physical_barrier_hessian(const double d, const double dhat)
+double normalized_barrier_hessian(const double d, const double dhat)
 {
     if (d <= 0.0 || d >= dhat) {
         return 0.0;
     }
-    const double d_minus_dhat = d - dhat;
-    return (-2 * log(d / dhat) + (d_minus_dhat * d_minus_dhat) / (d * d) - 4)
-        / dhat
-        + 4 / d;
+
+    const double t0 = 1.0 / dhat;
+    const double t1 = d * t0;
+    const double t2 = 1 - t1;
+    return 4 * t0 * t2 / d + std::pow(t2, 2) / std::pow(d, 2)
+        - 2 * std::log(t1) / std::pow(dhat, 2);
 }
 } // namespace
 
 TEST_CASE("Test barrier derivatives", "[barrier]")
 {
-    double dhat = GENERATE(range(-5, 2));
+    bool use_dist_sqr = GENERATE(false, true);
+    double dhat = GENERATE_COPY(range(use_dist_sqr ? -2 : -5, 0));
     dhat = pow(10, dhat);
 
     double d =
@@ -64,11 +66,29 @@ TEST_CASE("Test barrier derivatives", "[barrier]")
         barrier_gradient = ipc::barrier_gradient;
         barrier_hessian = ipc::barrier_hessian;
     }
+    SECTION("Normalized barrier")
+    {
+        barrier = normalized_barrier;
+        barrier_gradient = normalized_barrier_gradient;
+        barrier_hessian = normalized_barrier_hessian;
+    }
     SECTION("Barrier with physical units")
     {
-        barrier = physical_barrier;
-        barrier_gradient = physical_barrier_gradient;
-        barrier_hessian = physical_barrier_hessian;
+        barrier = [dhat](double d, double p_dhat) {
+            return dhat * normalized_barrier(d, p_dhat);
+        };
+        barrier_gradient = [dhat](double d, double p_dhat) {
+            return dhat * normalized_barrier_gradient(d, p_dhat);
+        };
+        barrier_hessian = [dhat](double d, double p_dhat) {
+            return dhat * normalized_barrier_hessian(d, p_dhat);
+        };
+    }
+
+    if (use_dist_sqr) {
+        d_vec *= d;
+        d *= d;
+        dhat *= dhat;
     }
 
     Eigen::VectorXd fgrad(1);
@@ -79,7 +99,7 @@ TEST_CASE("Test barrier derivatives", "[barrier]")
     Eigen::VectorXd grad(1);
     grad << barrier_gradient(d, dhat);
 
-    CAPTURE(dhat, d, fgrad(0), grad(0));
+    CAPTURE(dhat, d, fgrad(0), grad(0), use_dist_sqr);
     CHECK(fd::compare_gradient(fgrad, grad));
 
     // Check hessian
@@ -91,30 +111,39 @@ TEST_CASE("Test barrier derivatives", "[barrier]")
 
     grad << barrier_hessian(d, dhat);
 
-    CAPTURE(dhat, d, fgrad(0), grad(0));
+    CAPTURE(dhat, d, fgrad(0), grad(0), use_dist_sqr);
     CHECK(fd::compare_gradient(fgrad, grad));
 }
 
 TEST_CASE("Test physical barrier", "[barrier]")
 {
-    double dhat = GENERATE(range(-5, 2));
-    dhat = pow(10, dhat);
+    const bool use_dist_sqr = GENERATE(false, true);
+    const double dhat =
+        pow(10, GENERATE_COPY(range(use_dist_sqr ? -5 : -5, 0)));
 
-    double d =
-        GENERATE_COPY(take(10, random(dhat / 2, 0.9 * dhat))); // ∈ [0, d̂]
+    const double d = GENERATE_COPY(take(10, random(dhat / 2, 0.9 * dhat)));
+
+    const double p_d = (use_dist_sqr ? d : 1) * d;
+    const double p_dhat = (use_dist_sqr ? dhat : 1) * dhat;
+    const double divisor = use_dist_sqr ? std::pow(dhat, 3) : dhat;
+
+    CAPTURE(use_dist_sqr, dhat, d, divisor, p_d, p_dhat);
 
-    double b_original = ipc::barrier(d, dhat) / dhat;
-    double b_new = physical_barrier(d, dhat);
+    const double b_original = ipc::barrier(p_d, p_dhat) / divisor;
+    const double b_new = dhat * normalized_barrier(p_d, p_dhat);
 
     CHECK(b_original == Catch::Approx(b_new));
 
-    double b_original_gradient = ipc::barrier_gradient(d, dhat) / dhat;
-    double b_new_gradient = physical_barrier_gradient(d, dhat);
+    const double b_original_gradient =
+        ipc::barrier_gradient(p_d, p_dhat) / divisor;
+    const double b_new_gradient =
+        dhat * normalized_barrier_gradient(p_d, p_dhat);
 
     CHECK(b_original_gradient == Catch::Approx(b_new_gradient));
 
-    double b_original_hessian = ipc::barrier_hessian(d, dhat) / dhat;
-    double b_new_hessian = physical_barrier_hessian(d, dhat);
+    const double b_original_hessian =
+        ipc::barrier_hessian(p_d, p_dhat) / divisor;
+    const double b_new_hessian = dhat * normalized_barrier_hessian(p_d, p_dhat);
 
     CHECK(b_original_hessian == Catch::Approx(b_new_hessian));
 }
\ No newline at end of file
diff --git a/tests/friction/test_force_jacobian.cpp b/tests/friction/test_force_jacobian.cpp
index 9b8dac11d..93cfaf41a 100644
--- a/tests/friction/test_force_jacobian.cpp
+++ b/tests/friction/test_force_jacobian.cpp
@@ -258,7 +258,7 @@ TEST_CASE(
         mu = 0.5;
         dhat = 1e-3;
         // kappa = 67873353;
-        kappa = 67873353 / 10;
+        kappa = 67873353 / 10 * dhat;
         epsv_dt = 1e-4;
     }
     SECTION("square-circle-dense")