Mass matrix lumping

Question

Dear all,

I have been looking in the documentation but I haven't seen it mentioned a single time so I was wondering: is it possible to ask FEniCS to lump the mass matrix?

Thanks in advance!
Vincent

Answer 1

This is not possible in the form language. But you can assemble the action of a matrix and get a vector. Then you use this to set the diagonal of the matrix.

mass_form = v*u*dx()
mass_action_form = action(mass_form, Constant(1))
M_lumped = assemble(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))

Comments

Thank you for the answer!

I am a little confused though: you are saying that I need to assemble the action of the matrix, that is to use a matrix-free method, right? Then I do not really understand what you mean by setting the diagonal of the matrix (since we need not assemble the matrix for a matrix-free method...)

---

I updated the answer so it now present the mass_action_form. This you can do what ever you want with. You can also exchange the Constant(1) with any scalar form such as a Function, which then when assembled will yield your matrix free evaluation.

However if you want the lumped matrix, which you actually asked for you need to assemble a vector and set the diagonal.

---

Thanks! It helps a lot! :)

---

Be advised, however, that with parabolic and/or higher order continuous Lagrange interpolation lumping by using action (i.e. summing the elements of the consistent mass matrix rows) can lead lead to null or negative diagonal terms. For example, with parabolic shape functions on a triangle you get a singular mass matrix.
A possible alternative is to get the diagonal of the consistent mass matrix, and scale it up to get the correct mass. See e.g. "The finite element method" by Zienkiewicz & Taylor", vol 2, chap. 16.