This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

The corresponding files can be obtained from:

• Jupyter Notebook: mass_lumping.ipynb
• Python script: mass_lumping.py

Lumping a mass matrix

Explicit dynamics simulations require the usage of lumped mass matrices i.e. diagonal mass matrices for which inversion can be done explicitly.

We show how to do this for P1 and P2 Lagrange elements. The former case can be done quite easily whereas the latter case requires a special integration scheme which is not present by default in FEniCS.

P1 elements

First method

This problem has been already tackled in the past in old forum posts, see for instance here.

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(1, 1)
V1 = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V1)
u = TrialFunction(V1)

mass_form = v*u*dx
mass_action_form = action(mass_form, Constant(1))

M_consistent = assemble(mass_form)
print("Consistent mass matrix:\n", np.array_str(M_consistent.array(), precision=3))

M_lumped = assemble(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))
print("Lumped mass matrix:\n", np.array_str(M_lumped.array(), precision=3))

Consistent mass matrix:
 [[0.083 0.042 0.042 0.   ]
 [0.042 0.167 0.083 0.042]
 [0.042 0.083 0.167 0.042]
 [0.    0.042 0.042 0.083]]
Lumped mass matrix:
 [[0.167 0.    0.    0.   ]
 [0.    0.333 0.    0.   ]
 [0.    0.    0.333 0.   ]
 [0.    0.    0.    0.167]]

For explicit dynamics simulation, the mass matrix can then be manipulated using the diagonal vector, to compute for instance \(M^{-1}w\) where \(w\) is some Function.

M_vect = assemble(mass_action_form)
w = Function(V1)
iMw = Function(V1)
iMw.vector().set_local(w.vector().get_local()/M_vect.get_local())

Second method

The above trick consists in summing all rows and affecting the corresponding mass to the corresponding diagonal degree of freedom. E.g. for element \(T\) with number of vertices \(n\), the diagonal lumped mass \(i\) is obtained as:

\[\begin{equation} \overline{M}_i = \sum_j M_{ij} = \sum_j \int_{T} N_i(x)N_j(x) \text{ dx} = \int_T N_i(x) \text{ dx} = \dfrac{|T|}{n} \end{equation}\]

for linear shape functions.

We can obtain the same result by integrating the mass form with a “vertex” scheme (i.e. quadrature points are located at the vertices of the simplex with equal weights \(\dfrac{|T|}{n}\):

\[\begin{split}\begin{align} M_{ij} &= \int_{T} N_i(x)N_j(x) \text{ dx('vertex')} \\ &= \sum_{k=1}^n \dfrac{|T|}{n} N_i(x_k)N_j(x_k) = \sum_{k=1}^n \dfrac{|T|}{n} \delta_{ik}\delta_{jk} = \dfrac{|T|}{n} \delta_{ij} \label{diagonal-mass} \tag{1} \end{align}\end{split}\]

M_lumped2 = assemble(v*u*dx(scheme="vertex", metadata={"degree":1, "representation":"quadrature"}))
print("Lumped mass matrix ('vertex' scheme):\n", np.array_str(M_lumped2.array(), precision=3))

Lumped mass matrix ('vertex' scheme):
 [[0.167 0.    0.    0.   ]
 [0.    0.333 0.    0.   ]
 [0.    0.    0.333 0.   ]
 [0.    0.    0.    0.167]]

/home/bleyerj/.local/lib/python3.6/site-packages/ffc/jitcompiler.py:234: QuadratureRepresentationDeprecationWarning:
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** representation instead.                               ***
*** ===================================================== ***
  issue_deprecation_warning()

P2 elements

The problem with the first method is that we obtain a singular matrix for P2 elements:

V2 = FunctionSpace(mesh, "CG", 2)
v = TestFunction(V2)
u = TrialFunction(V2)

mass_form = v*u*dx
mass_action_form = action(mass_form, Constant(1))

M_consistent = assemble(mass_form)
print("Consistent mass matrix:\n", np.array_str(M_consistent.array(), precision=3))

M_lumped = assemble(mass_form)
M_lumped.zero()
M_lumped.set_diagonal(assemble(mass_action_form))
print("Lumped mass matrix:\n", np.array_str(M_lumped.array(), precision=3))

Consistent mass matrix:
 [[ 0.017 -0.003 -0.003 -0.011  0.     0.     0.     0.     0.   ]
 [-0.003  0.033 -0.006  0.    -0.011  0.    -0.003 -0.011  0.   ]
 [-0.003 -0.006  0.033  0.     0.    -0.011 -0.003  0.    -0.011]
 [-0.011  0.     0.     0.178  0.044  0.044 -0.011  0.044  0.044]
 [ 0.    -0.011  0.     0.044  0.089  0.044  0.     0.     0.   ]
 [ 0.     0.    -0.011  0.044  0.044  0.089  0.     0.     0.   ]
 [ 0.    -0.003 -0.003 -0.011  0.     0.     0.017  0.     0.   ]
 [ 0.    -0.011  0.     0.044  0.     0.     0.     0.089  0.044]
 [ 0.     0.    -0.011  0.044  0.     0.     0.     0.044  0.089]]
Lumped mass matrix:
 [[0.    0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.333 0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.167 0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.167 0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.167 0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.167]]

We can however adapt the second method by considering a generalized “vertex” quadrature scheme with quadrature points located at the element degrees of freedom points i.e. at vertices and facets mid-points for P2 elements. In this case, we still have \(\eqref{diagonal-mass}\). To do this we need to implement a custom quadrature scheme called "lumped" based on the following this discussion.

This has been implemented in the file lumping_scheme.py for triangles and tetrahedra.

import lumping_scheme

M_lumped = assemble(v*u*dx(scheme="lumped", degree=2))
print("Lumped mass matrix ('lumped' scheme):\n", np.array_str(M_lumped.array(), precision=3))

Lumped mass matrix ('lumped' scheme):
 [[0.083 0.    0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.167 0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.167 0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.167 0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.083 0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.083 0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.083 0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.083 0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.083]]