# Qforder (qfV=qfV1lump)

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?

1 Like

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:

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`