Physics-Based Simulation

Mass Matrix and Lumping

Recall from Discretization of Weak Forms that:

$$M_{ab}={\int }_{{\Omega }^0}R(\mathbf{X},0)N_a(\mathbf{X})N_b(\mathbf{X})d\mathbf{X}.\text{\hspace{1em}(19.2.1)}$$

With the solid domain discretized into triangles $\mathcal{T}$, we have:

$$M_{ab}=\sum _{e\in \mathcal{T}}{\int }_{{\Omega }_e^0}R(\mathbf{X},0)N_a(\mathbf{X})N_b(\mathbf{X})d\mathbf{X},\text{\hspace{1em}(19.2.2)}$$

where ${\Omega }_e^0$ represents the material space of triangle $e$. Note that for linear triangle elements, since $N_i$ is nonzero only on the incident triangles of node $i$, here we only need to consider triangles with both $a$ and $b$ being their vertices.

Let us change the integration variable from $\mathbf{X}$ to $(\beta ,\gamma )$, which gives:

$$\begin{array}{rl}& {\int }_{{\Omega }_e^0}R(\mathbf{X},0)N_a(\mathbf{X})N_b(\mathbf{X})d\mathbf{X}\\ =& {\int }_0^1{\int }_0^{1-\beta }R(\beta ,\gamma ,0)N_a(\beta ,\gamma )N_b(\beta ,\gamma )|\det \left(\frac{\partial \mathbf{X}}{\partial (\beta ,\gamma )}\right)|d\gamma d\beta .\end{array}\text{\hspace{1em}(19.2.3)}$$

For simplicity, let us denote the vertices of this triangle $e$ as ${\mathbf{X}}_1$, ${\mathbf{X}}_2$, and ${\mathbf{X}}_3$, and then we have:

$$|\det \left(\frac{\partial \mathbf{X}}{\partial (\beta ,\gamma )}\right)|=|\det ([{\mathbf{X}}_2-{\mathbf{X}}_1,{\mathbf{X}}_3-{\mathbf{X}}_1])|=2A_e,$$

where $A_e$ is the area of triangle $e$. Here, $N_a$ and $N_b$ take $1-\beta -\gamma$, $\beta$, or $\gamma$ depending on the vertex indices $a$ and $b$. For example, if $a$ and $b$ correspond to the 2nd and 3rd vertices of triangle $e$, then $N_a=\beta$ and $N_b=\gamma$. Assuming uniform density, we have:

$$\begin{array}{rl}{\int }_{{\Omega }_e^0}R(\mathbf{X},0)N_a(\mathbf{X})N_b(\mathbf{X})d\mathbf{X}=& 2R{A}_e{\int }_0^1{\int }_0^{1-\beta }\beta \gamma d\gamma d\beta \\ =& 2R{A}_e{\int }_0^1{\frac{1}{2}\beta {\gamma }^2|}_{\gamma =0}^{\gamma =1-\beta }d\beta \\ =& R{A}_e{\int }_0^1\beta (1-\beta {)}^{2}d\beta \\ =& R{A}_e{\left(\frac{{\beta }^2}{2}-\frac{2{\beta }^3}{3}+\frac{{\beta }^4}{4}\right)|}_{\beta =0}^{\beta =1}=\frac{1}{12}R{A}_e.\end{array}\text{\hspace{1em}(19.2.4)}$$

With mass lumping, $M_{ab}^{\text{lump}}={\delta }_{ab}{\sum }_c{M}_{ac}$, which means:

$$M_{aa}^{\text{lump}}=\sum _{e\in \mathcal{T}}\sum _{b\in \mathcal{V}}{\int }_{{\Omega }_e^0}R(\mathbf{X},0)N_a(\mathbf{X})N_b(\mathbf{X})d\mathbf{X},\text{\hspace{1em}(19.2.5)}$$

where $\mathcal{V}$ contains all the nodes of the mesh, and all off-diagonal entries of $M^{\text{lump}}$ are $0$. Similarly, due to the locality of $N$, for each triangle element, $b$ only needs to traverse all three triangle vertices:

\begin{aligned} M_{aa}^\text{lump} & = \sum_{e \in \mathcal{T}(a)} 2 R A_e \left(\int_0^1 \int_0^{1-\beta} \beta (1- \beta - \gamma) d\gamma d\beta + \int_0^1 \int_0^{1-\beta} \beta^2 d\gamma d\beta \ &+ \int_0^1 \int_0^{1-\beta} \beta \gamma d\gamma d\beta \right) \ & = \sum_{e \in \mathcal{T}(a)} 2 R A_e \int_0^1 \beta d\gamma d\beta = \sum_{e \in \mathcal{T}(a)} 2 R A_e \int_0^1 \beta \gamma |{\gamma=0}^{\gamma =1-\beta} d\beta \ & = \sum{e \in \mathcal{T}(a)} 2 R A_e \int_0^1 \beta (1-\beta) d\beta = \sum_{e \in \mathcal{T}(a)} 2 R A_e \left. \left(\frac{\beta^2}{2} - \frac{\beta^3}{3} \right) \right|{\beta=0}^{\beta=1} \ &= \sum{e \in \mathcal{T}(a)} \frac{1}{3} R A_e, \end{aligned} \htmlId{eq:lec19:lumped_mass_matrix_sum_tri_param}{} \tag{19.2.6}

where $\mathcal{T}(a)$ denotes the set of triangles incident to node $a$. This result also explains why in Inversion-Free Elasticity when computing the mass for all the nodes, we traverse all triangles, calculate the mass of the triangle $R{A}_e$ and evenly distribute it to the three vertices. With the mass matrix computed, the momentum change and external body force terms including their energy forms are all easy to deal with.