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

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?

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.

• Awesome answer, as always. I would be also very interested in your opinion on the second part of the question, when lumping is not allowed/bad approximation, if anything comes to mind. – Anton Menshov Jun 5 at 1:37
• @AntonMenshov: It would seem that it would be difficult (maybe impossible?) to get a good approximation via lumping for higher order elements, since (e.g. diagonal) lumping in that case would to be equivalent to a lower order quadrature applied to higher order polynomials. – Paul Jun 5 at 3:16
• @WolfgangBangerth: I think I understand now. So, it’s like using newton-cotes rules for integration instead of gaussian quadrature. Since each lagrange interpolation functions have unit values at one specific node, migrating the quadrature points to the nodes results in only the diagonal terms becoming nonzero (at least, for linear elements). – Paul Jun 5 at 4:05
• @Paul, I have used mass lumping for higher order elements (up to 14) and it works fine. It is customary for the Spectral Element Method, where you use the nodes as quadrature points. In that case, the order of approximation of $2n - 3$. – nicoguaro Jun 5 at 11:44
• The important point is that for higher order elements, one needs to define the "lumped" mass matrix via particular quadrature formulas. The original form -- from which the term "lumped" originates -- added the off-diagonal entries to the diagonal, but that only works if they are all positive. If you apply Gaussian quadrature, this is true for lowest-order elements, but not for higher-order elements. – Wolfgang Bangerth Jun 6 at 13:18

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 error1: 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.

• The idea is correct in principle, but using a second-order in time formulation is a poor choice because it leads to an oscillating behavior of the system that is only damped by the $C\dot u$ term. A better choice is to only use $C\dot u$ and forego the $M\ddot u$ term altogether. But the same argument applies for the $C\dot u$ term, and in fact in many cases $C$ can also lead to a mass matrix. – Wolfgang Bangerth Jun 6 at 13:20
• I should say "mass matrix" in quotes here, since the term "mass matrix" of course originally came from the second order term where $\ddot u$ is generally multiplied by the density which, when integrated over an element, yields a "mass". – Wolfgang Bangerth Jun 6 at 13:21

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!