Mathematically, why does mass matrix / load vector lumping work?

Question

I know that people often replace consistent mass matrices with lumped diagonal matrices. In the past, I've also implemented a code where the load vector is assembled in a lumped fashion rather than an FEM-consistent fashion. But I've never looked into why we are allowed to do this in the first place.

What is the intuition behind lumping that allows one to apply it to mass and load vectors? What is the mathematical justification for it? In what situations is lumping not allowed / not a good approximation for mass and load vectors?

Answer 1

In the finite element method, the matrix entries and right hand side entries are defined as integrals. We can, in general, not compute these exactly and apply quadrature. But there are many quadrature formulas one could choose, and one often chooses them in a way so that (i) the error introduced by quadrature is of the same order as that due to discretization, or at least not substantially worse, and (ii) the matrix has certain properties that turn out to be convenient.

Mass lumping is an example of this working: If one chooses a particular quadrature formula (namely, the one with quadrature points located at the interpolation points of the finite element), then the resulting mass matrix happens to be diagonal. That's quite convenient for the computational implementation, and the reason why people use these quadrature formulas. It's also the reason why it "works": This particular choice of quadrature formula still has reasonably high order.

Answer 2

Diagonal matrices have obvious advantages in speeding up numerical computations, and Wolfgang Bangerth's answer is a good explanation of how to calculate a diagonal mass matrix, but it doesn't answer the OP's question "why does this work" in the sense of "why is it a good approximation to the physics you are modelling".

Conceptually, you can separate the response of an element into three parts: translational motion of a rigid body, rigid rotation about the element center of mass, and the deformation of the element.

The basic function of the element mass matrix is to represent the element KE as a quadratic form (i.e. $\frac 1 2 v^T M v$ where $v$ are the nodal velocities).

As the element size decreases, the contribution to the KE from rigid body rotation decreases faster than the contribution from translation, (for a solid element with a typical linear size of $a$, the mass is proportional to $a^3$ but the moments of inertia are proportional to $a^5$) and the contribution from the element deformation is negligible (at least for problems with small elastic strains).

Therefore, you only really need a "good" approximation to the rigid body parts of the motion, i.e. 6 DOFs, and in fact a good approximation to only the KE from rigid body translation, i.e. 3 DOFs, will converge as the element size is reduced.

The diagonal terms of the element matrix contain more than enough independent parameters to represent those 3 or 6 KE terms with sufficient accuracy. In fact for higher order elements, you can use mass diagonal mass matrices where the diagonal terms for the mid-side nodes are zero.

Note that this is a completely different situation from the element potential energy, where the contributions from rigid body translation and rotation are zero, and the only thing that matters is representing the strain energy corresponding to the element deformation. A diagonal stiffness matrix would therefore not be a feasible idea!

Answer 3

In addition to the other answers, there are scenarios in which errors in the mass matrix have no influence on the desired result.

Given a nonlinear deformation problem of the form $K(u) \space u = f(u)$ that has a unique solution $\hat{u}$, one might consider solving this as the dynamic problem $K(u) \space u + C(u) \space \dot{u} + M \space \ddot{u} = f(u)$ with mass matrix $M$ and damping matrix $C$ through some appropriate time discretization. It is easy to see that the original problem is solved when a resting state has been reached ($\dot{u} = \ddot{u} = 0$). Importantly, $M$ has no influence on this result (as long as convergence to the unique result is achieved).

The justification for mass lumping then has more to do with convergence speed rather than e.g. quadrature error¹: The ease of inverting a diagonal $M$ results in fast time step iterations with an appropriate time integration scheme (where $M^{-1}$ is the only matrix inverse).

^{1 Although reasoning about dynamic physical behavior is of course easier with a "correct" mass matrix - e.g. angular momentum may be improperly conserved by lumped mass matrices.}