# Inconsistent behaviour of vertex quadrature scheme / mass-lumping for variational form with more than 2 functions?

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

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

(f,g)_h = \int_\Omega \mathcal{I}_h ( f(x) \cdot g(x) ) dx = \sum_{z\in\mathcal{N}_h} f(z) \cdot g(z) \int_\Omega \phi_z dx.

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

(\phi_{z_i}, \phi_{z_j})_h = \sum_{z\in\mathcal{N}_h} \phi_{z_i}(z) \cdot \phi_{z_j}(z) \int_\Omega \phi_z dx = \delta_{z_i, z_j} \int_\Omega \phi_{z_i} dx .

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

(f,A v)_h = \int_\Omega \mathcal{I}_h ( f(x) \cdot A(x) v(x) ) dx = \sum_{z\in\mathcal{N}_h} f(z) \cdot A(z) v(z) \int_\Omega \phi_z dx.

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

(\phi_{z_i},a \phi_{z_j})_h = \sum_{z\in\mathcal{N}_h} \phi_{z_i} (z) a(z) \phi_{z_j}(z) \int_\Omega \phi_z dx = \delta_{z_i, z_j} a(z) \int_\Omega \phi_{z_i} dx.

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 a simple 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:
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

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:

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

You are misunderstanding

as the matrix A couples test functions from different dofs (i.e. those in the vector space at the same block.

You can see this band with a slightly modified version of your code:

from dolfinx.fem import functionspace
import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx.mesh import create_unit_square, CellType
from dolfinx.fem import assemble_matrix
from dolfinx.fem import Function, form

from ufl import TrialFunction, TestFunction, Measure, dot, as_tensor, inner

"""
Sanity test for mass-lumping
"""
version = dolfinx.__version__
print("Using dolfinx version ", version)

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

n = 1
domain = create_unit_square(MPI.COMM_WORLD, n, n, CellType.quadrilateral)

dim = 2

V = functionspace(domain, ("Lagrange", 1, (dim,)))
W = functionspace(domain, ("Lagrange", 1, (dim, dim)))

u = TrialFunction(V)
u0 = Function(V)
v = TestFunction(V)

f_I = as_tensor(((1, 0), (0, 1)))
f_M = as_tensor(((1, 1), (1, 1)))

eval_list = [
("dxL", inner(u, dot(f_I, v))*dxL),
("dxLf", inner(u, dot(f_M, v))*dxL)

]
for description, eq in eval_list:
X = assemble_matrix(form(eq))
X.scatter_reverse()
C = X.to_dense()
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)

Yielding

sing dolfinx version  0.8.0.0
---
dxL
assembled matrix with shape (8, 8)
[[0.25 0.   0.   0.   0.   0.   0.   0.  ]
[0.   0.25 0.   0.   0.   0.   0.   0.  ]
[0.   0.   0.25 0.   0.   0.   0.   0.  ]
[0.   0.   0.   0.25 0.   0.   0.   0.  ]
[0.   0.   0.   0.   0.25 0.   0.   0.  ]
[0.   0.   0.   0.   0.   0.25 0.   0.  ]
[0.   0.   0.   0.   0.   0.   0.25 0.  ]
[0.   0.   0.   0.   0.   0.   0.   0.25]]
diagonals only True
---
dxLf
assembled matrix with shape (8, 8)
[[0.25 0.25 0.   0.   0.   0.   0.   0.  ]
[0.25 0.25 0.   0.   0.   0.   0.   0.  ]
[0.   0.   0.25 0.25 0.   0.   0.   0.  ]
[0.   0.   0.25 0.25 0.   0.   0.   0.  ]
[0.   0.   0.   0.   0.25 0.25 0.   0.  ]
[0.   0.   0.   0.   0.25 0.25 0.   0.  ]
[0.   0.   0.   0.   0.   0.   0.25 0.25]
[0.   0.   0.   0.   0.   0.   0.25 0.25]]
diagonals only False
the following indices are affected
[[0 1 2 3 4 5 6 7]
[1 0 3 2 5 4 7 6]]

There are other ways of lumping the matrix (I don’t have a nice example at hand, but it could be done combining DOLFINx kernels and numba).

dP was used in Legacy FEniCS, and has not been implemented in DOLFINx at the time of writing.

Thank you so much for the detailed reply!

That makes a lot of sense and is reassuring since my misunderstanding only lies in my interpretation of the discrete inner product. Also, there is fortunately no difference between the mathematical definition of the discrete inner product and its implementation in FEniCSx.