Physics-Based Simulation

Barrier Energy and Its Derivatives

With the point-edge distance functions implemented, we can traverse all point-edge pairs to assemble the total barrier energy and its derivatives. These will be used to solve for the search direction in the time-stepping optimization.

Since squared distances are used, here we rescale the barrier function to $$b(d^2,{\hat{d}}^2)=\{\begin{array}{ll}\frac{\kappa }{8}\hat{d}(\frac{d^2}{{\hat{d}}^2}-1)\ln \frac{d^2}{{\hat{d}}^2}& d<\hat{d}\\ 0& d\ge \hat{d},\end{array}\text{\hspace{1em}(21.3.1)}$$ so that $$\frac{{\partial }^{2}b}{\partial (\frac{d}{\hat{d}}{)}^2}({\hat{d}}^2,{\hat{d}}^2)=\kappa \hat{d}$$ still holds. Analogous to elasticity, $s=d/\hat{d}$ can be viewed as a strain measure, then the 2nd-order derivative of the energy density (per area) function $b$ w.r.t. $s$ at $s=1$ would correspond to Young's modulus times thickness $\hat{d}$, which makes $\kappa$ physically meaningful and convenient to set.

Based on Equation (20.3.1), we can derive the gradient and Hessian of the barrier potential as $$\begin{array}{rl}\nabla P_b(x)& =\frac{1}{2}\sum _{\hat{a}}{A}_{\hat{a}}\sum _{e\in \mathcal{E}-I({\mathbf{X}}_{\hat{a}})}\frac{\partial b}{\partial d}(d({\mathbf{x}}_{\hat{a}},e),{\hat{d}}^2){\nabla }_{x}d({\mathbf{x}}_{\hat{a}},e)\text{and}\\ {\nabla }^2{P}_b(x)& =\frac{1}{2}\sum _{\hat{a}}{A}_{\hat{a}}\sum _{e\in \mathcal{E}-I({\mathbf{X}}_{\hat{a}})}(\frac{{\partial }^{2}b}{\partial d^2}(d({\mathbf{x}}_{\hat{a}},e),{\hat{d}}^2){\nabla }_{x}d({\mathbf{x}}_{\hat{a}},e){\nabla }_{x}d({\mathbf{x}}_{\hat{a}},e{)}^T+\\ & \frac{\partial b}{\partial d}(d({\mathbf{x}}_{\hat{a}},e),{\hat{d}}^2){\nabla }_x^{2}d({\mathbf{x}}_{\hat{a}},e)),\end{array}$$ where $A_{\hat{a}}=\frac{\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}-1}\Vert +\Vert {\mathbf{X}}_{\hat{a}}-{\mathbf{X}}_{\hat{a}+1}\Vert }{2}$ and we omitted the superscripts and subscripts for the squared point-edge distance functions ($d({\mathbf{x}}_{\hat{a}},e)$ denotes $d_{\text{sq}}^{\text{PE}}({\mathbf{x}}_{\hat{a}},e)$ here).

The energy, gradient, and Hessian of the barrier contact potential are implemented as follows:

Implementation 21.3.1 (Barrier energy computation, BarrierEnergy.py).

    # self-contact
    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:
                    sum += 0.5 * contact_area[xI] * dhat * kappa / 8 * (s - 1) * math.log(s)

Implementation 21.3.2 (Barrier energy gradient computation, BarrierEnergy.py).

    # self-contact
    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:
                    local_grad = 0.5 * contact_area[xI] * dhat * (kappa / 8 * (math.log(s) / dhat_sqr + (s - 1) / d_sqr)) * PE.grad(x[xI], x[eI[0]], x[eI[1]])
                    g[xI] += local_grad[0:2]
                    g[eI[0]] += local_grad[2:4]
                    g[eI[1]] += local_grad[4:6]

Implementation 21.3.3 (Barrier energy Hessian computation, BarrierEnergy.py).

    # self-contact
    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:
                    d_sqr_grad = PE.grad(x[xI], x[eI[0]], x[eI[1]])
                    s = d_sqr / dhat_sqr
                    # since d_sqr is used, need to divide by 8 not 2 here for consistency to linear elasticity:
                    local_hess = 0.5 * contact_area[xI] * dhat * utils.make_PSD(kappa / (8 * d_sqr * d_sqr * dhat_sqr) * (d_sqr + dhat_sqr) * np.outer(d_sqr_grad, d_sqr_grad) \
                        + (kappa / 8 * (math.log(s) / dhat_sqr + (s - 1) / d_sqr)) * PE.hess(x[xI], x[eI[0]], x[eI[1]]))
                    index = [xI, eI[0], eI[1]]
                    for nI in range(0, 3):
                        for nJ in range(0, 3):
                            for c in range(0, 2):
                                for r in range(0, 2):
                                    IJV[0].append(index[nI] * 2 + r)
                                    IJV[1].append(index[nJ] * 2 + c)
                                    IJV[2] = np.append(IJV[2], local_hess[nI * 2 + r, nJ * 2 + c])