Physics-Based Simulation

Linear Triangle Elements

In previous discussions, we learned to calculate $\Psi$ and its derivatives with respect to $\mathbf{F}$. For simulation, however, we require $\frac{\partial \Psi }{\partial x}$ and $\frac{{\partial }^2\Psi }{\partial x^2}$. This necessitates a clear understanding of $\mathbf{F}(x)$, as it allows us to employ the chain rule to derive these derivatives with respect to $x$ effectively.

In 2D simulations, we often divide the solid domain into non-degenerate triangular elements. Assume the mapping $\mathbf{x}=\varphi (\mathbf{X})$ is linear within each triangle, thus keeping the deformation gradient $\mathbf{F}$ constant. Referencing Example 12.2.1, for a triangle defined by vertices ${\mathbf{X}}_1{\mathbf{X}}_2{\mathbf{X}}_3$, we have the equations: $${\mathbf{x}}_2-{\mathbf{x}}_1=\mathbf{F}({\mathbf{X}}_2-{\mathbf{X}}_1)\text{and}{\mathbf{x}}_3-{\mathbf{x}}_1=\mathbf{F}({\mathbf{X}}_3-{\mathbf{X}}_1),$$ where ${\mathbf{x}}_i$ denotes the world-space coordinates of the triangle vertices. This relationship leads to the expression for $\mathbf{F}$: $$\mathbf{F}=[{\mathbf{x}}_2-{\mathbf{x}}_1,{\mathbf{x}}_3-{\mathbf{x}}_1][{\mathbf{X}}_2-{\mathbf{X}}_1,{\mathbf{X}}_3-{\mathbf{X}}_1{]}^{-1}.\text{\hspace{1em}(15.1.1)}$$ Equation (15.1.1) shows that $\mathbf{F}$, derived here, maps any segment within the triangle to its world-space counterpart through linear combinations of the triangle edges ${\mathbf{X}}_2-{\mathbf{X}}_1$ and ${\mathbf{X}}_3-{\mathbf{X}}_1$. A more general and rigorous derivation of this formula will be presented in subsequent chapters.

Once $\mathbf{F}(x)$ is established, we can calculate its derivative with respect to $x$ for each triangle as follows: $$\frac{\partial [{\mathbf{F}}_{11},{\mathbf{F}}_{21},{\mathbf{F}}_{12},{\mathbf{F}}_{22}{]}^T}{\partial [{\mathbf{x}}_1^T,{\mathbf{x}}_2^T,{\mathbf{x}}_3^T{]}^T}=\left[\begin{array}{cccccc}-{\mathbf{B}}_{11}-{\mathbf{B}}_{21}& 0& {\mathbf{B}}_{11}& 0& {\mathbf{B}}_{21}& 0\\ 0& -{\mathbf{B}}_{11}-{\mathbf{B}}_{21}& 0& {\mathbf{B}}_{11}& 0& {\mathbf{B}}_{21}\\ -{\mathbf{B}}_{12}-{\mathbf{B}}_{22}& 0& {\mathbf{B}}_{12}& 0& {\mathbf{B}}_{22}& 0\\ 0& -{\mathbf{B}}_{12}-{\mathbf{B}}_{22}& 0& {\mathbf{B}}_{12}& 0& {\mathbf{B}}_{22}\end{array}\right],$$ where $\mathbf{B}=[{\mathbf{X}}_2-{\mathbf{X}}_1,{\mathbf{X}}_3-{\mathbf{X}}_1{]}^{-1}$ represents the inverse of the matrix formed by subtracting the first vertex from the second and third vertices. This matrix $\mathbf{B}$ can be precomputed at initialization along with other properties such as the volume and Lame parameters of each triangle:

Implementation 15.1.1 (Precomputation of element information, simulator.py).

# rest shape basis, volume, and lame parameters
vol = [0.0] * len(e)
IB = [np.array([[0.0, 0.0]] * 2)] * len(e)
for i in range(0, len(e)):
    TB = [x[e[i][1]] - x[e[i][0]], x[e[i][2]] - x[e[i][0]]]
    vol[i] = np.linalg.det(np.transpose(TB)) / 2
    IB[i] = np.linalg.inv(np.transpose(TB))
mu_lame = [0.5 * E / (1 + nu)] * len(e)
lam = [E * nu / ((1 + nu) * (1 - 2 * nu))] * len(e)

The Young's modulus and Poisson's ratio:

E = 1e5         # Young's modulus
nu = 0.4        # Poisson's ratio

Here, e no longer stores all edge elements as in mass-spring models but represents all triangle elements, which can be generated by modifying the meshing code as follows:

Implementation 15.1.2 (Assembling per-triangle vertex indices, square_mesh.py).

    # connect the nodes with triangle elements
    e = []
    for i in range(0, n_seg):
        for j in range(0, n_seg):
            # triangulate each cell following a symmetric pattern:
            if (i % 2)^(j % 2):
                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
                e.append([(i + 1) * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1, i * (n_seg + 1) + j + 1])
            else:
                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
                e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1, i * (n_seg + 1) + j + 1])

Triangles are arranged in a symmetric pattern and can be rendered by drawing the three edges:

Implementation 15.1.3 (Draw triangles, simulator.py).

        pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[0]]), screen_projection(x[eI[1]]))
        pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[1]]), screen_projection(x[eI[2]]))
        pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[2]]), screen_projection(x[eI[0]]))