Case Study: Square Drop* - Physics-Based Simulation

Case Study: Square Drop

To conclude, let's consider a case study where we simulate a square dropped onto a fixed planar ground. Building on our previous mass-spring model for an elastic square, we augment a barrier potential into its Incremental Potential and apply the filter line search scheme to manage the contact between the square's degrees of freedom (DOFs) and the ground.

The excutable Python project for this section can be found at https://github.com/phys-sim-book/solid-sim-tutorial under the 3_contact folder. MUDA GPU implementations can be found at https://github.com/phys-sim-book/solid-sim-tutorial-gpu under the simulators/3_contact folder.

If we further limit the planar ground to be horizontal, e.g. at $y=y_0$, its signed distance function can be made even simpler than Equation (7.1.1): $d (x) = x_{y} - y_{0}, \nabla d (x) = [01], \nabla^{2} d (x) = 0 . (8.3.1)$ Combining it with Equation (7.2.4) and Equation (7.2.5), we can conveniently implement the gradient and Hessian computation for the barrier potential of this horizontal ground:

Implementation 8.3.1 (Barrier energy value, gradient, and Hessian, BarrierEnergy.py).

import math
import numpy as np

dhat = 0.01
kappa = 1e5

def val(x, y_ground, contact_area):
    sum = 0.0
    for i in range(0, len(x)):
        d = x[i][1] - y_ground
        if d < dhat:
            s = d / dhat
            sum += contact_area[i] * dhat * kappa / 2 * (s - 1) * math.log(s)
    return sum

def grad(x, y_ground, contact_area):
    g = np.array([[0.0, 0.0]] * len(x))
    for i in range(0, len(x)):
        d = x[i][1] - y_ground
        if d < dhat:
            s = d / dhat
            g[i][1] = contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d))
    return g

def hess(x, y_ground, contact_area):
    IJV = [[0] * len(x), [0] * len(x), np.array([0.0] * len(x))]
    for i in range(0, len(x)):
        IJV[0][i] = i * 2 + 1
        IJV[1][i] = i * 2 + 1
        d = x[i][1] - y_ground
        if d < dhat:
            IJV[2][i] = contact_area[i] * dhat * kappa / (2 * d * d * dhat) * (d + dhat)
        else:
            IJV[2][i] = 0.0
    return IJV

For the filter line search, with the position in the last iteration $\mathbf{x}$ and a search direction $\mathbf{p}$ of a specific node, the signed distance function is simply \[ d(\mathbf{x} + \alpha \mathbf{p}) = \mathbf{x}_y + \alpha \mathbf{p}_y - y_0, \] where $\alpha$ is the step size, and there is only one positive real root $\alpha = (y_0 - \mathbf{x}_y) / \mathbf{p}_y$ when $\mathbf{p}_y < 0$ since $\mathbf{x}_y > y_0$ (no interpenetration up to current iteration). Taking the minimum of the positive real root per node then gives us the step size upper bound $\alpha_C$ defined in Equation (8.2.1):

Implementation 8.3.2 (Ground CCD, BarrierEnergy.py).

def init_step_size(x, y_ground, p):
    alpha = 1
    for i in range(0, len(x)):
        if p[i][1] < 0:
            alpha = min(alpha, 0.9 * (y_ground - x[i][1]) / p[i][1])
    return alpha

Here we scale the upper bound by $0.9\times$ so that exact touching configurations with $d=0$ and $b = \infty$ (floating-point number overflow) can be avoided.

Then once we make sure the step size upper bound is used to initialize the line search

Implementation 8.3.3 (Filter line search, time_integrator.py).

        # filter line search
        alpha = BarrierEnergy.init_step_size(x, y_ground, p)  # avoid interpenetration and tunneling
        while IP_val(x + alpha * p, e, x_tilde, m, l2, k, y_ground, contact_area, h) > E_last:
            alpha /= 2

and that the contact area weights for all nodes are calculated

Implementation 8.3.4 (Contact area, simulator.py).

contact_area = [side_len / n_seg] * len(x)     # perimeter split to each node

and passed to our simulator, we can simulate the square drop with mass-spring stiffness k=2e4 and time step size h=0.01 as shown in Figure 8.3.1.

**Figure 8.3.1.** A mass-spring elastic square is dropped onto the ground with $0$ initial velocity under gravity. Here we show the frames when the square is: just dropped, first touching the ground, compressed to the maximum in this simulation, and becoming static.

Remark 8.3.1 (Contact Layer Integration). Since in practice, contact forces are only exerted on the boundary of the solids, the barrier potential should be integrated only on the boundary as well. This also explains why in our case study the contact area weight per node is simply calculated as the diameter of the square evenly distributed onto each boundary node. However, as mass-spring elasticity cannot guarantee that all interior nodes will stay inside the boundary of the solid, we simply apply the barrier potential to all nodal DOFs of the square.