Physics-Based Simulation

Point-Edge Distance

Next, we calculate the point-edge distance and its derivatives. These will be used to solve for the contact forces. For a node $\mathbf{p}$ and an edge ${\mathbf{e}}_0{\mathbf{e}}_1$, their squared distance is defined as

$$d_{\text{sq}}^{\text{PE}}(\mathbf{p},{\mathbf{e}}_0,{\mathbf{e}}_1)=\underset{\lambda }{\min }\Vert \mathbf{p}-((1-\lambda ){\mathbf{e}}_0+\lambda {\mathbf{e}}_1){\Vert }^2\text{s.t.}\lambda \in [0,1],$$

which is the closest squared distance between $\mathbf{p}$ and any point on ${\mathbf{e}}_0{\mathbf{e}}_1$.

Remark 21.2.1 (Distance Calculation Optimization). Here, we use the squared unsigned distances for evaluating the contact energies. This approach avoids taking square roots, which can complicate the expression of the derivatives and increase numerical rounding errors during computation. Additionally, unsigned distances can be directly extended for codimensional pairs, such as point-point pairs, which are useful when simulating particle contacts in 2D. They also do not suffer from locking issues, as signed distances do, when there are large displacements.

Fortunately, $d_{\text{sq}}^{\text{PE}}(\mathbf{p},{\mathbf{e}}_0,{\mathbf{e}}_1)$ is a piece-wise smooth function w.r.t. the DOFs: $$d_{\text{sq}}^{\text{PE}}(\mathbf{p},{\mathbf{e}}_0,{\mathbf{e}}_1)=\{\begin{array}{l}\Vert \mathbf{p}-{\mathbf{e}}_0{\Vert }^2\text{if}({\mathbf{e}}_1-{\mathbf{e}}_0)\cdot (\mathbf{p}-{\mathbf{e}}_0)<0,\\ \Vert \mathbf{p}-{\mathbf{e}}_1{\Vert }^2\text{if}({\mathbf{e}}_1-{\mathbf{e}}_0)\cdot (\mathbf{p}-{\mathbf{e}}_0)>\Vert {\mathbf{e}}_1-{\mathbf{e}}_0{\Vert }^2,\\ \frac{1}{\Vert {\mathbf{e}}_1-{\mathbf{e}}_0{\Vert }^2}(\det ([\mathbf{p}-{\mathbf{e}}_0,{\mathbf{e}}_1-{\mathbf{e}}_0]){)}^2\text{otherwise,}\end{array}\text{\hspace{1em}(21.2.1)}$$ where the smooth expression can be determined by checking whether the node is inside the orthogonal span of the edge. Given these smooth expressions, we can differentiate each of them and obtain the derivatives of the point-edge distance function. The implementations are as follows:

Implementation 21.2.1 (Point-Edge distance calculation (Hessian omitted), PointEdgeDistance.py).

import numpy as np

import distance.PointPointDistance as PP
import distance.PointLineDistance as PL

def val(p, e0, e1):
    e = e1 - e0
    ratio = np.dot(e, p - e0) / np.dot(e, e)
    if ratio < 0:    # point(p)-point(e0) expression
        return PP.val(p, e0)
    elif ratio > 1:  # point(p)-point(e1) expression
        return PP.val(p, e1)
    else:            # point(p)-line(e0e1) expression
        return PL.val(p, e0, e1)

def grad(p, e0, e1):
    e = e1 - e0
    ratio = np.dot(e, p - e0) / np.dot(e, e)
    if ratio < 0:    # point(p)-point(e0) expression
        g_PP = PP.grad(p, e0)
        return np.reshape([g_PP[0:2], g_PP[2:4], np.array([0.0, 0.0])], (1, 6))[0]
    elif ratio > 1:  # point(p)-point(e1) expression
        g_PP = PP.grad(p, e1)
        return np.reshape([g_PP[0:2], np.array([0.0, 0.0]), g_PP[2:4]], (1, 6))[0]
    else:            # point(p)-line(e0e1) expression
        return PL.grad(p, e0, e1)

It can be verified that the point-edge distance function is $C^1$-continuous everywhere, including at the interfaces between different segments. For the point-point case, we have:

Implementation 21.2.2 (Point-Point distance calculation, PointPointDistance.py).

import numpy as np

def val(p0, p1):
    e = p0 - p1
    return np.dot(e, e)

def grad(p0, p1):
    e = p0 - p1
    return np.reshape([2 * e, -2 * e], (1, 4))[0]

def hess(p0, p1):
    H = np.array([[0.0] * 4] * 4)
    H[0, 0] = H[1, 1] = H[2, 2] = H[3, 3] = 2
    H[0, 2] = H[1, 3] = H[2, 0] = H[3, 1] = -2
    return H

For the point-line case, the distance evaluations can be implemented as follows, and the derivatives can be obtained using symbolic differentiation tools.

Implementation 21.2.3 (Point-Line distance calculation (Hessian omitted), PointLineDistance.py).

import numpy as np

def val(p, e0, e1):
    e = e1 - e0
    numerator = e[1] * p[0] - e[0] * p[1] + e1[0] * e0[1] - e1[1] * e0[0]
    return numerator * numerator / np.dot(e, e)

def grad(p, e0, e1):
    g = np.array([0.0] * 6)
    t13 = -e1[0] + e0[0]
    t14 = -e1[1] + e0[1]
    t23 = 1.0 / (t13 * t13 + t14 * t14)
    t25 = ((e0[0] * e1[1] + -(e0[1] * e1[0])) + t14 * p[0]) + -(t13 * p[1])
    t24 = t23 * t23
    t26 = t25 * t25
    t27 = (e0[0] * 2.0 + -(e1[0] * 2.0)) * t24 * t26
    t26 *= (e0[1] * 2.0 + -(e1[1] * 2.0)) * t24
    g[0] = t14 * t23 * t25 * 2.0
    g[1] = t13 * t23 * t25 * -2.0
    t24 = t23 * t25
    g[2] = -t27 - t24 * (-e1[1] + p[1]) * 2.0
    g[3] = -t26 + t24 * (-e1[0] + p[0]) * 2.0
    g[4] = t27 + t24 * (p[1] - e0[1]) * 2.0
    g[5] = t26 - t24 * (p[0] - e0[0]) * 2.0
    return g