How to formulate lumped mass matrix in FEM

Question

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

Answer 1

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.