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 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:
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.
