Physics-Based Simulation

Mass-Spring Potential Energy

In this case study, we focus exclusively on incorporating the mass-spring elasticity potential into our system. The concept of mass-spring elasticity is akin to treating each edge of the mesh as if it were a spring. This approach is inspired by Hooke's Law, allowing us to formulate the potential energy on edge $e$ as follows:

$$P_e(x)=l^2\frac{1}{2}k{\left(\frac{\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2}{l^2}-1\right)}^2,\text{\hspace{1em}(4.3.1)}$$

Here, ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ represent the current positions of the two endpoints of the edge. The variable $l$ denotes the original length of the edge, and $k$ is a parameter controlling the spring's stiffness. Notably, when the distance between the two endpoints $\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert$ equals the original length $l$, the potential energy $P_e(x)$ attains its global minimum value of $0$, indicating no force is exerted.

An important aspect of this formulation is the inclusion of $l^2$ at the beginning. This is analogous to integrating the spring energy across the solid and choosing edges as quadrature points. This integration helps maintain a consistent relationship between the stiffness behavior and the parameter $k$, regardless of mesh resolution variations.

Another deviation from standard spring energy formulations is our avoidance of the square root operation. We directly use $\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2$, making our model polynomial in nature. This simplification yields more streamlined expressions for the gradient and Hessian:

$$\frac{\partial P_e}{\partial {\mathbf{x}}_1}(x)=-\frac{\partial P_e}{\partial {\mathbf{x}}_2}(x)=2k\left(\frac{\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2}{l^2}-1\right)({\mathbf{x}}_1-{\mathbf{x}}_2),$$

$$\begin{array}{rl}\frac{{\partial }^2{P}_e}{\partial {\mathbf{x}}_1^2}(x)& =\frac{{\partial }^2{P}_e}{\partial {\mathbf{x}}_2^2}(x)=-\frac{{\partial }^2{P}_e}{\partial {\mathbf{x}}_1{\mathbf{x}}_2}(x)=-\frac{{\partial }^2{P}_e}{\partial {\mathbf{x}}_2{\mathbf{x}}_1}(x)\\ & =\frac{4k}{l^2}({\mathbf{x}}_1-{\mathbf{x}}_2)({\mathbf{x}}_1-{\mathbf{x}}_2{)}^T+2k\left(\frac{\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2}{l^2}-1\right)\mathbf{I}\\ & =\frac{2k}{l^2}\left(2({\mathbf{x}}_1-{\mathbf{x}}_2)({\mathbf{x}}_1-{\mathbf{x}}_2{)}^T+(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }^2-l^2)\mathbf{I}\right).\end{array}$$

The total potential energy of the system, denoted as $P(x)$, can be derived by summing the potential energy across all edges. This is calculated using Equation (4.3.1). Thus, the total potential energy is expressed as: $$P(x)=\sum _e{P}_e(x)$$ where the summation is taken over all edges in the mesh.

Implementation 4.3.1 (MassSpringEnergy.py).

import numpy as np
import utils

def val(x, e, l2, k):
    sum = 0.0
    for i in range(0, len(e)):
        diff = x[e[i][0]] - x[e[i][1]]
        sum += l2[i] * 0.5 * k[i] * (diff.dot(diff) / l2[i] - 1) ** 2
    return sum

def grad(x, e, l2, k):
    g = np.array([[0.0, 0.0]] * len(x))
    for i in range(0, len(e)):
        diff = x[e[i][0]] - x[e[i][1]]
        g_diff = 2 * k[i] * (diff.dot(diff) / l2[i] - 1) * diff
        g[e[i][0]] += g_diff
        g[e[i][1]] -= g_diff
    return g

def hess(x, e, l2, k):
    IJV = [[0] * (len(e) * 16), [0] * (len(e) * 16), np.array([0.0] * (len(e) * 16))]
    for i in range(0, len(e)):
        diff = x[e[i][0]] - x[e[i][1]]
        H_diff = 2 * k[i] / l2[i] * (2 * np.outer(diff, diff) + (diff.dot(diff) - l2[i]) * np.identity(2))
        H_local = utils.make_PSD(np.block([[H_diff, -H_diff], [-H_diff, H_diff]]))
        # add to global matrix
        for nI in range(0, 2):
            for nJ in range(0, 2):
                indStart = i * 16 + (nI * 2 + nJ) * 4
                for r in range(0, 2):
                    for c in range(0, 2):
                        IJV[0][indStart + r * 2 + c] = e[i][nI] * 2 + r
                        IJV[1][indStart + r * 2 + c] = e[i][nJ] * 2 + c
                        IJV[2][indStart + r * 2 + c] = H_local[nI * 2 + r, nJ * 2 + c]
    return IJV

In dealing with the Hessian matrix of the mass-spring energy, a key consideration is its non-symmetric positive definite (SPD) nature. To address this, a specific modification is employed: we neutralize the negative eigenvalues of the local Hessian corresponding to each edge. This is done prior to incorporating these local Hessians into the global matrix. The process involves setting negative eigenvalues to zero, thus ensuring that the resulting global Hessian matrix adheres more closely to the desired SPD properties. This modification is crucial for Newton's method.

Implementation 4.3.2 (Positive Semi-Definite Projection).

import numpy as np
import numpy.linalg as LA

def make_PSD(hess):
    [lam, V] = LA.eigh(hess)    # Eigen decomposition on symmetric matrix
    # set all negative Eigenvalues to 0
    for i in range(0, len(lam)):
        lam[i] = max(0, lam[i])
    return np.matmul(np.matmul(V, np.diag(lam)), np.transpose(V))