Physics-Based Simulation

Filter Line Search for Non-Inversion

To guarantee non-inversion just like non-interpenetration (see Filter Line Search) during the simulation, we can similarly filter the line search initial step size with a critical step size ${\alpha }^I$ that first brings the volume of any triangles to $0$. This can be obtained by solving a 1D equation per triangle: $$V({\mathbf{x}}_i+{\alpha }^I{\mathbf{p}}_i)=0,\text{\hspace{1em}(15.3.1)}$$ and taking the minimum of the solved step sizes. Here ${\mathbf{p}}_i$ is the search direction of node $i$, and in 2D, Equation (15.3.1) is equivalent to: $$\det ([{\mathbf{x}}_{21}^{\alpha },{\mathbf{x}}_{31}^{\alpha }])\equiv {\mathbf{x}}_{21,1}^{\alpha }{\mathbf{x}}_{31,2}^{\alpha }-{\mathbf{x}}_{21,2}^{\alpha }{\mathbf{x}}_{31,1}^{\alpha }=0\text{\hspace{1em}(15.3.2)}$$ with ${\mathbf{x}}_{ij}^{\alpha }={\mathbf{x}}_{ij}+{\alpha }^I{\mathbf{p}}_{ij}$ and ${\mathbf{x}}_{ij}={\mathbf{x}}_i-{\mathbf{x}}_j$, ${\mathbf{p}}_{ij}={\mathbf{p}}_i-{\mathbf{p}}_j$. Expanding Equation (15.3.2) we obtain: $$({\mathbf{x}}_{21,1}+{\alpha }^I{\mathbf{p}}_{21,1})({\mathbf{x}}_{31,2}+{\alpha }^I{\mathbf{p}}_{31,2})-({\mathbf{x}}_{21,2}+{\alpha }^I{\mathbf{p}}_{21,2})({\mathbf{x}}_{31,1}+{\alpha }^I{\mathbf{p}}_{31,1})=0,$$ which can be reorganized as a quadratic equation of ${\alpha }^I$: $$\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])({\alpha }^I{)}^2+(\det ([{\mathbf{x}}_{21},{\mathbf{p}}_{31}])+\det ([{\mathbf{p}}_{21},{\mathbf{x}}_{31}])){\alpha }^I+\det ([{\mathbf{x}}_{21},{\mathbf{x}}_{31}])=0.$$ Here, note that $\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])$ can be very tiny when the nodes do not move much or when their movement barely changes to triangle area in the current timestep, thus the equation can be degenerated into a linear one. To robustly detect this degenerate case, we cannot directly check whether $\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])$ is $0$ due to numerical errors. In fact, checking whether $\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])$ is below an epsilon is still tricky, because the scale of $\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])$ heavily depends on the scene dimension and nodal velocity during the simulation. Therefore, we use $\det ([{\mathbf{x}}_{21},{\mathbf{x}}_{31}])$ as a scaling and obtain a scaled but equivalent equation: $$\frac{\det ([{\mathbf{p}}_{21},{\mathbf{p}}_{31}])}{\det ([{\mathbf{x}}_{21},{\mathbf{x}}_{31}])}({\alpha }^I{)}^2+\frac{\det ([{\mathbf{x}}_{21},{\mathbf{p}}_{31}])+\det ([{\mathbf{p}}_{21},{\mathbf{x}}_{31}])}{\det ([{\mathbf{x}}_{21},{\mathbf{x}}_{31}])}{\alpha }^I+1=0,\text{\hspace{1em}(15.3.3)}$$ where magnitude checks can be safely performed on any coefficients with unitless thresholds.

In practice, we also need to allow some slackness so that the step size to be taken will not lead to an exactly $0$ volume. Thus, we solve ${\alpha }^I$ such that it first decreases the volume of any triangles by $90\%$, which can be realized by modifying the constant term coefficient in Equation (15.3.3) from $1$ to $0.9$:

Implementation 15.3.1 (Filter line search, NeoHookeanEnergy.py).

def init_step_size(x, e, p):
    alpha = 1
    for i in range(0, len(e)):
        x21 = x[e[i][1]] - x[e[i][0]]
        x31 = x[e[i][2]] - x[e[i][0]]
        p21 = p[e[i][1]] - p[e[i][0]]
        p31 = p[e[i][2]] - p[e[i][0]]
        detT = np.linalg.det(np.transpose([x21, x31]))
        a = np.linalg.det(np.transpose([p21, p31])) / detT
        b = (np.linalg.det(np.transpose([x21, p31])) + np.linalg.det(np.transpose([p21, x31]))) / detT
        c = 0.9  # solve for alpha that first brings the new volume to 0.1x the old volume for slackness
        critical_alpha = utils.smallest_positive_real_root_quad(a, b, c)
        if critical_alpha > 0:
            alpha = min(alpha, critical_alpha)
    return alpha

Here, if the equation does not have a positive real root, that means for this specific triangle, the step size can be taken arbitrarily large and it will not trigger inversion.

The quadratic equation can be solved as

Implementation 15.3.2 (Solve quadratic equation, utils.py).

def smallest_positive_real_root_quad(a, b, c, tol = 1e-6):
    # return negative value if no positive real root is found
    t = 0
    if abs(a) <= tol:
        if abs(b) <= tol: # f(x) = c > 0 for all x
            t = -1
        else:
            t = -c / b
    else:
        desc = b * b - 4 * a * c
        if desc > 0:
            t = (-b - math.sqrt(desc)) / (2 * a)
            if t < 0:
                t = (-b + math.sqrt(desc)) / (2 * a)
        else: # desv<0 ==> imag, f(x) > 0 for all x > 0
            t = -1
    return t

With scaled coefficients, we simply use a unitless threshold, e.g. \code{1e-6}, to check for degeneracies. If no positive real roots are found, the function simply returns $-1$.

Now as we filter the initial step size in addition to non-interpenetration:

Implementation 15.3.3 (Apply filter, time_integrator.py).

        alpha = min(BarrierEnergy.init_step_size(x, n, o, p), NeoHookeanEnergy.init_step_size(x, e, p))  # avoid interpenetration, tunneling, and inversion

and make sure all added data structures and modified functions are reflected in the time integrator, we can finally simulate the compressing square example from Moving Boundary Condition with guaranteed non-inversion (see Figure 15.3.1).