# How to formulate lumped mass matrix in FEM

When solving time dependent PDE's using the finite element method, for example say the heat equation, if we use explicit time stepping then we have to solve a linear system because of the mass matrix. For example if we stick with the heat equation example,

$\frac{\partial{u}}{\partial{t}} = c\nabla{}^{2}u$

then using forward Euler we get

$M(\frac{u^{n+1}-u^{n}}{dt}) = -cKu^{n}$

and thus even though we are using an explicit time stepping scheme we still have to solve a linear system. This is obviously a major problem since the primary advantage of using explicit schemes is to NOT have to solve a linear system. I have read that a common way to get around this problem is to instead use a "lumped" mass matrix which transforms the regular (consistent?) mass matrix into a diagonal matrix and thus makes inversion trivial. Upon doing a google search however I am still not entirely sure how this lumped mass matrix is created. For example looking at he paper NUMERICAL EXPERIMENTS ON MASS LUMPING FOR THE ADVECTION-DIFFUSION EQUATION by Edson Wendland Harry and Edmar Schulz they create their lumped mass matrix by simply summing all coefficients onto the diagonal. So for example if our original consistent mass matrix was:

\begin{pmatrix} 4 & 2 & 1 & 2 \\ 2 & 4 & 2 & 1 \\ 1 & 2 & 4 & 2 \\ 2 & 1 & 2 & 4\end{pmatrix}

then the lumped mass matrix would be:

\begin{pmatrix} 9 & 0 & 0 & 0 \\ 0 & 9 & 0 & 0 \\ 0 & 0 & 9 & 0 \\ 0 & 0 & 0 & 9\end{pmatrix}

My question then is: Is this the correct way to form the lumped mass matrix? What disadvantages exist when using the lumped mass matrix instead of the full consistent mass matrix in terms of accuracy? The authors of the paper I mentioned actually suggested not using the lumped mass matrix, although it seemed they were only using an implicit time stepping scheme which I thought was odd given that the primary reason to use such matrices is for explicit methods.

Note: I would never use forward Euler to solve the heat equation, that was just an example. Also if it matters my problem is solving the Navier Stokes equations where the nonlinear term is treated explicitly and the diffusion term is treated implicitly.

Thanks

• As far as I understand it, the mass matrix is constant with respect to time, so you only need to perform one LU factorization, and then the remaining work would be $O(n^{2})$ operations on different right-hand sides, which doesn't seem onerous, because you'd be doing something similar with matrix-vector multiplies anyway. Commented May 21, 2015 at 19:18
• Yes I could do that if I was using a direct solver, but if I am using PCG or some other iterative solver I don't think that would help Commented May 22, 2015 at 16:42
• I personally don't trust mass lumping mathematically. Computationally, it doesn't give you any advantage unless you aim for explicit time stepping, in which case the diagonal mass matrix makes much easier to solve. If you're using an implicit time stepping method, you don't gain any sparsity in the matrix. I think you only gain error at that point by not using a consistent matrix.
– Paul
Commented Apr 7, 2017 at 21:35
• I’m surprised no one has mentioned the method of Fried and Markus (1975) for quadrilaterals, which uses nodes at Lobatto points to avoid a loss of truncation error. Not an issue until you get to cubics, but precludes serendipity elements. The idea has been extended to triangles, but requires a special basis and quadrature. Commented Jan 29, 2020 at 15:10

I do not think that there is a definite answer to this, because it might change from one topic to other (and also depends on the type of elements you are using). There are some recent papers talking about that, as well [2]. So, it is not a closed discussion. Furthermore, you can have different inertial components (at least in mechanics), when you have elements with kinematic constraints as beams or shells.

Zienkiewicz (See [1], section 16.2.4) discuss three methods for lumping the mass matrix

1. The row sum method $$M^{(\text{lumped})}_{ii} = \sum_j M_{ij}$$

2. Diagonal scaling $$M^{(\text{lumped})}_{ii} = c M_{ii}$$ with $c$ adjusted to satisfy $\sum_j M^{(\text{lumped})}_{jj} = \int_\Omega \rho d\Omega$.

3. Evaluation of $M$ using a quadrature involving only the nodal points and thus automatically yielding a diagonal matrix for standard element shape function in wich $N_i=0$ for $x=x_j$, $i\neq j$.

Not all the methods work in all the cases, for example, the row sum method does not work for 8-node serendipity elements since it would lead to negative masses.

I have used method 2 with the factor being the total mass of the element ($M_{\text{tot}}$) divided by the trace of the matrix ($\mathrm{Tr}(M)$), i.e.

$$M^{(\text{lumped})}_{ii} = \frac{M_{\text{tot}}}{\mathrm{Tr}(M)} M_{ii} \quad \text{(no summation on i)} \enspace .$$

I have also used method 3 with the so-called Spectral Element Methods with Lobatto nodes (using these locations as nodes and integration points), that automatically lead to diagonal matrices.

From [1], you can see this figure describing some of the methods for some element types

## References

[1] Zhu, J., Z. R. L. Taylor, and O. C. Zienkiewicz. "The finite element method: its basis and fundamentals." (2005): 54-102.

[2] Felippa, Carlos A., Qiong Guo, and K. C. Park. "Mass matrix templates: General description and 1d examples." Archives of Computational Methods in Engineering 22.1 (2015): 1-65.