Hey Mik,
What type of problem are you solving? I’m aware of some papers that show the benefits of mass-lumping can be problem dependent:
- Hinton, E., Rock, T., & Zienkiewicz, O. C. (1976). A note on mass lumping and related processes in the finite element method. Earthquake Engineering & Structural Dynamics, 4(3), 245–249. https://doi.org/10.1002/eqe.4290040305
- Wendland, E., & Schulz, H. E. (2005). Numerical experiments on mass lumping for the advection-diffusion equation. Revista Minerva, 2(2), 227-233. http://www.fipai.org.br/Minerva%2002(02)%2012.pdf
I’ve made a small example to test how to confirm you get a lumped mass matrix:
mesh Th = square(2, 2);
plot(Th,wait=true);
fespace Vh(Th, P1);
Vh u, v;
cout << "ndof = " << Vh.ndof << endl;
// P1 mass matrix
varf vp1 (u, v) = int2d(Th,qft=qf1pTlump) (u*v);
matrix A = vp1(Vh, Vh);
cout << A << endl;
// Extract COO matrix
int nnz = A.nnz;
int[int] I(nnz), J(nnz);
real[int] Avals(nnz);
[I,J,Avals] = A;
// Copy the values into a dense array for nicer output
int n = Vh.ndof;
real[int,int] B(n,n);
B = 0;
for [i, Ai : Avals] {
B(I(i),J(i)) = Ai;
}
cout << B << endl;
The output I get for B
is:
9 9
0.08333333333 0 0 0 0 0 0 0 0
0 0.125 0 0 0 0 0 0 0
0 0 0.04166666667 0 0 0 0 0 0
0 0 0 0.125 0 0 0 0 0
0 0 0 0 0.25 0 0 0 0
0 0 0 0 0 0.125 0 0 0
0 0 0 0 0 0 0.04166666667 0 0
0 0 0 0 0 0 0 0.125 0
0 0 0 0 0 0 0 0 0.08333333333
which is a diagonal matrix, like you would expect from lumping. If you look at the sum of the diagonal values: 1/4 + 4 × 1/8 + 2 × 1/12 + 2 × 1/24 = 1/4 + 1/2 + 1/4 = 1.
I suspect that with P2 elements you can pull off mass lumping by using the following variational form:
// P2 mass matrix
varf vp2 (u, v) = int2d(Th,qft=qf1pTlump)(0.5*u*v)
+ int2d(Th,qft=qf2pT)(0.5*u*v);
These results seem to match the following resource:
Lumping a mass matrix — Numerical tours of continuum mechanics using FEniCS master documentation
The reason I’ve added the 0.5
coefficients in the vp2
form is since, both quadrature rules qf1pTlump
and qf2pT
use |T_k|/3
as the weight values of the quadrature. Since P2 elements use 6 nodal values (3 corners + 3 edges) the weight should be |T_k|/6