I stumbled upon some behavior of the vertex quadrature which is unintuitive to me…

I use the vertex quadrature scheme to apply **mass-lumping** by calling

```
Measure("dx", metadata = {"quadrature_rule": "vertex"})
```

Mathematically that means, that I want to replace the standard L^2 product by a discrete inner product, usually defined for two (or more) continuous functions f,g such that

where \phi_z are the Lagrangian shape functions and \mathcal{N}_h the nodes of the mesh. This leads to a diagonal mass matrix since

It further also works when **evaluating more than two functions** by choosing e.g. g= A v where A is a matrix function:

Also in this case the resulting matrix should be diagonal since for a scalar function a:

This seems to be also consistent with the description in the ffcx package. There it says:

The vertex scheme, i.e.,

averaging the function value in the

> vertices and multiplying with the simplex volume, is only of

order 1 and inferior to other generic schemes in terms of

error reduction. Equation systems generated with the vertex

scheme have some properties that other schemes lack, e.g., the

mass matrix is asimple diagonal matrix. This may be

prescribed in certain cases.

However in practice when using more than two functions, the matrix is not diagonal any more. Consider the following MWE:

```
import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx.mesh import create_unit_square, CellType
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import Function, form, VectorFunctionSpace
from ufl import TrialFunction, TestFunction, Measure, inner, dx, dP,dot
# dolfinx 0.7.3 or 0.6.0
"""
Sanity test for mass-lumping
"""
version = dolfinx.__version__
print("Using dolfinx version ", version)
n = 2
domain = create_unit_square(MPI.COMM_WORLD, n, n, CellType.triangle)
dim =2
if float(version[0:3]) >= 0.7:
from dolfinx.fem import functionspace, ElementMetaData
V = functionspace(domain, ElementMetaData("Lagrange", 1, shape =(dim,)))
W = functionspace(domain, ElementMetaData("Lagrange", 1, shape =(dim,dim)))
else:
from dolfinx.fem import VectorFunctionSpace, TensorFunctionSpace
V = VectorFunctionSpace(domain, ("Lagrange", 1), dim =dim)
W = TensorFunctionSpace(domain, ("Lagrange", 1), shape =(dim,dim))
u = TrialFunction(V)
u0 = Function(V)
v = TestFunction(V)
f = Function(W)
f.x.array[:] = np.random.rand(len(f.x.array[:]))
def only_diagonal(Cb):
"""
checks if a matrix has only diagonal entries
"""
only_diagonal=True
for i in range(len(Cb)):
for j in range(len(Cb[0])):
if i!=j:
if Cb[i][j] != 0.0:
only_diagonal=False
return only_diagonal
def non_conforming_indices(C):
"""
return list of indices of a matrix C which are non-diagonal and non-zero
"""
# indexes where C is non zero
non_zero_indexes = np.transpose((np.absolute(C)>0.0).nonzero())
# bool array to sort out diagonals
non_diag_subindex = non_zero_indexes.T[0] != non_zero_indexes.T[1]
# resulting indices
res = (non_zero_indexes)[non_diag_subindex]
return res
#NOTE - For vertex quadrature rule order will always be reduced to 1
dxL = Measure("dx", metadata = {"quadrature_rule": "vertex"})
eval_list = [ ("L2", inner(u, v)*dx),
("dxL", inner(u,v)*dxL),
("L2f", inner(u, dot(f,v))*dx),
("dxLf", inner(u,dot(f,v))*dxL) ] #,
# ("dPf", inner(u,dot(f,v))*dP),
# ("dP", inner(u,v)*dP)
# ]
for description,eq in eval_list:
X= assemble_matrix(form(eq))
X.assemble()
X.convert("dense")
C = X.getDenseArray()
np.savetxt("matrices/matrix-"+description+".out", np.around(C, decimals=2, out=None), delimiter=',', fmt='%f')
print("---")
print(description)
print("assembled matrix with shape", np.shape(C), "\n", C)
diag_bool = only_diagonal(C)
print("diagonals only", diag_bool)
if diag_bool == False:
print("the following indices are affected \n", non_conforming_indices(C).T)
```

This leads me to the questions:

- ) Is this a misunderstanding on my side or a mistake?
- ) Is there any way to enforce a diagonal matrix?
- ) What is the point of the Measure “dP”? That one only throws an error.