Qforder (qfV=qfV1lump)

Mik · January 13, 2021, 10:27am

Dear All,

I have a question regarding the qforder, more specifically on qfV=qfV1lump.

When I implement qfV=qfV1lump instead of default value, my program (P1 space) runs almost 3 to 5 times faster, with improved accuracy.

If I am not wrong, is it doing lumping of the mass matrix?
I wonder if the lumping can be as simple as adding qfV=qfV1lump to my varf?

Thanks in advance for your answer.

ivan · November 7, 2021, 9:19am

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