Physics-Based Simulation

Precomputing Normal and Tangent Information

To make the temporally discretized friction force integrable, we must explicitly discretize certain normal and tangent information. This information only needs to be calculated once at the beginning of each time step, as it will remain constant during each optimization.

First, we need to calculate ${\lambda }^n$ for each point-edge pair using ${\mathbf{x}}^n$. Recall that we used squared distances as input for the barrier functions, so ${\lambda }^n$ should be calculated using the chain rule as follows:

$$\begin{array}{rl}{\lambda }_{\hat{a},e}^n& =\frac{1}{2}{A}_{\hat{a}}\left(-\frac{\partial b(d_{\text{sq}}^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e),{\hat{d}}^2)}{\partial d^{\text{PE}}}\right)=\frac{1}{2}{A}_{\hat{a}}\left(-\frac{\partial b(d_{\text{sq}}^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e),{\hat{d}}^2)}{\partial d_{\text{sq}}^{\text{PE}}}\frac{\partial d_{\text{sq}}^{\text{PE}}}{\partial d^{\text{PE}}}\right)\\ & =\frac{1}{2}{A}_{\hat{a}}\left(-\frac{\partial b(d_{\text{sq}}^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e),{\hat{d}}^2)}{\partial d_{\text{sq}}^{\text{PE}}}\right)2d^{\text{PE}}.\end{array}$$

According to the scaled barrier function taking squared distance as input (Equation (21.3.1)), we can derive

$$\frac{\partial b(d_{\text{sq}},{\hat{d}}^2)}{\partial d_{\text{sq}}}=\{\begin{array}{ll}\frac{\kappa }{8}\hat{d}\left(\frac{1}{{\hat{d}}^2}\ln \frac{d_{\text{sq}}}{{\hat{d}}^2}+\frac{1}{d_{\text{sq}}}\left(\frac{d_{\text{sq}}}{{\hat{d}}^2}-1\right)\right)& \text{if }d<\hat{d};\\ 0& \text{if }d\ge \hat{d}.\end{array}$$

Remark 22.2.1. The set of point-edge pairs for friction in our semi-implicit friction setting is fixed in each time step and is different from the set of normal contact pairs. The set for friction only contains those pairs with $d^{\text{PE}}({\mathbf{x}}_{\hat{a}}^n,e)<\hat{d}$, and this does not change with the optimization variable ${\mathbf{x}}^{n+1}$ in the current time step.

Now for the tangent information, the key is to keep the normal and the barycentric coordinate of the closest point on the edge constant. For the $k$-th point-edge pair, if we denote the node indices for the point and edge as $p$, $e_0$, and $e_1$, then we can write the relative sliding velocity as

$${\mathbf{v}}_k=(\mathbf{I}-\mathbf{n}{\mathbf{n}}^T)({\mathbf{v}}_p-((1-r){\mathbf{v}}_{e_0}+r{\mathbf{v}}_{e_1})),$$

where $r={\mathrm{argmin}}_c\Vert {\mathbf{x}}_p-((1-c){\mathbf{x}}_{e_0}+c{\mathbf{x}}_{e_1})\Vert$ is the barycentric coordinate and $\mathbf{n}=({\mathbf{x}}_p-((1-r){\mathbf{x}}_{e_0}+r{\mathbf{x}}_{e_1}))/\Vert {\mathbf{x}}_p-((1-r){\mathbf{x}}_{e_0}+r{\mathbf{x}}_{e_1})\Vert$ is the normal of the edge. Here we see that $r$ and $\mathbf{n}$ are both dependent on $\mathbf{x}$, so directly integrating ${\mathbf{v}}_k$ is nontrivial. By calculating $\mathbf{n}$ and $r$ using ${\mathbf{x}}^n$, we obtain the semi-implicit relative sliding velocity

$${\bar{\mathbf{v}}}_k=(\mathbf{I}-{\mathbf{n}}^n({\mathbf{n}}^n{)}^T)({\mathbf{v}}_p-((1-r^n){\mathbf{v}}_{e_0}+r^n{\mathbf{v}}_{e_1})),$$

and now only the velocities are dependent on ${\mathbf{x}}^{n+1}$, which makes both integration and differentiation straightforward. If we denote ${\hat{\mathbf{v}}}_k={\mathbf{v}}_p-((1-r^n){\mathbf{v}}_{e_0}+r^n{\mathbf{v}}_{e_1})$, we obtain

$$\frac{\partial {\bar{\mathbf{v}}}_k}{\partial {\hat{\mathbf{v}}}_k}=(\mathbf{I}-{\mathbf{n}}^n({\mathbf{n}}^n{)}^T)\text{and}\frac{\partial {\hat{\mathbf{v}}}_k}{\partial [{\mathbf{x}}_p^T,{\mathbf{x}}_{e_0}^T,{\mathbf{x}}_{e_1}^T{]}^T}=\frac{1}{\hat{h}}\left[\begin{array}{ccc}\mathbf{I}& (r^n-1)\mathbf{I}& -r^n\mathbf{I}\end{array}\right].$$

Code

Next, let's look at the code. Implementation 22.2.1 calculates the barycentric coordinate of the closest point and the normal given point-edge nodal positions. The idea is to orthogonally project ${\mathbf{x}}_p$ onto the edge.

Implementation 22.2.1 (Calculating contact point and normal, PointEdgeDistance.py).

# compute normal and the parameterization of the closest point on the edge
def tangent(p, e0, e1):
    e = e1 - e0
    ratio = np.dot(e, p - e0) / np.dot(e, e)
    if ratio < 0:    # point(p)-point(e0) expression
        n = p - e0
    elif ratio > 1:  # point(p)-point(e1) expression
        n = p - e1
    else:            # point(p)-line(e0e1) expression
        n = p - ((1 - ratio) * e0 + ratio * e1)
    return [n / np.linalg.norm(n), ratio]

Then, Implementation 22.2.2 traverses all non-incident point-edge pairs with a distance smaller than $\hat{d}$, calculates $\lambda$, and calls the above function to calculate $\mathbf{n}$ and $r$.

As in Frictional Contact, these lines of code are executed at the beginning of each time step in time_integrator.py, and the information for each friction pair is stored and passed to the energy, gradient, and Hessian computation functions that we will discuss next.

Implementation 22.2.2 (Semi-implicit friction precomputation, BarrierEnergy.py).

    # self-contact
    mu_lambda_self = []
    dhat_sqr = dhat * dhat
    for xI in bp:
        for eI in be:
            if xI != eI[0] and xI != eI[1]: # do not consider a point and its incident edge
                d_sqr = PE.val(x[xI], x[eI[0]], x[eI[1]])
                if d_sqr < dhat_sqr:
                    s = d_sqr / dhat_sqr
                    # since d_sqr is used, need to divide by 8 not 2 here for consistency to linear elasticity
                    # also, lambda = -\partial b / \partial d = -(\partial b / \partial d^2) * (\partial d^2 / \partial d)
                    mu_lam = mu * -0.5 * contact_area[xI] * dhat * (kappa / 8 * (math.log(s) / dhat_sqr + (s - 1) / d_sqr)) * 2 * math.sqrt(d_sqr)
                    [n, r] = PE.tangent(x[xI], x[eI[0]], x[eI[1]]) # normal and closest point parameterization on the edge
                    mu_lambda_self.append([xI, eI[0], eI[1], mu_lam, n, r])