I have seen an example written in Fenics where it was possible to get the lumped mass matrix. I have tried to use it in Fenicsx but I have obtained an error when I use constant and as_backend_type

u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(u_, v_) * ufl.dx
aV = ufl.action(a, ufl.Constant(1))
dV = assemble(mass_action_form)
as_backend_type(dV).vec().reciprocal()

Please provide a minimal working example (i.e. a code that can reproduce the error, that people can copy-paste.
The code above is not DOLFINx code, so I cannot tell you where you have gone wrong.

This would give you the PETSc.Vector in legacy dolfin

I would use a dolfinx.fem.Function, as done in the previous example, as the âConstantâ cannot be used within action. You could of course do something along the lines of

from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.fem import FunctionSpace
import ufl
from ufl import dx, action
from petsc4py import PETSc
domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 2, mesh.CellType.triangle)
V = FunctionSpace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx
option = 2
if option == 1:
mass_action_form = ufl.replace(mass_form, {u:fem.Constant(domain, PETSc.ScalarType(1))})
elif option == 2:
Q = FunctionSpace(domain, ("DG", 0))
one = fem.Function(Q)
one.x.array[:] = 1
mass_action_form = action(mass_form, one)
else:
raise ValueError("Invalid option")
M_lumped = assemble_matrix(fem.form(mass_form))
M_lumped.zeroEntries()
M_lumped.setDiagonal(assemble_vector(fem.form(mass_action_form)))
print(M_lumped[:, :])

Please note that I had to make several changes to your code, as most of it didnât really make sense.

Thank you very much for providing the method, it did indeed work. I knew very little about ufl.replace before, so I checked its instructions and now have a better understanding of it.

Thank you for your prompt response, but it may be a bit difficult for me to understand this.

I find myself very confused about the relationship between different variables and data structures used in FEniCSx, particularly how they interact and are managed in the backend.

Could you please provide some insights or resources that could help clarifying?

Okay, I will take the time to read and understand your article and tutorial.

Also, why is DG-0 used here? I guess DG-0 may indicate that the entire domain is a constant of 1, while CG-1 only has a value of 1 at each node. However, it seems that the final lumped mass matrix obtained is the same.

It just has fewer degrees of freedom, and you are guaranteed that the basis value is 1 over the whole domain. It should also be the case for P1, but not necessarily for higher order polynomials.

I have tried it. CG or DG elements of any order can obtain the same lumped mass matrix. This means that the number of columns in the matrix of mass_form is not equal to the number of rows in the vector of one, while the sum of elements in each row of mass_form can still be correctly solved and assembled to the diagonal after action. I am confused about this.

I am still thinking about whether it would be better to directly operate on the matrix assembled by petsc to solve the lumped mass matrix instead of using some functions in UFL. For example:

...
Assuming M_consistent is already PETSc.Mat
...
M_lumped = M_consistent.copy()
M_lumped.setOption(PETSc.Mat.Option.NEW_NONZERO_ALLOCATION_ERR, False)
M_lumped.zeroEntries()
rstart, rend = M_consistent.getOwnershipRange()
for i in range(rstart, rend):
row_values = M_consistent.getRow(i)[1]
row_sum = np.sum(row_values)
M_lumped.setValue(i, i, row_sum, addv=False)
M_lumped.assemble()
print(M_lumped[:, :])

Will this approach be less efficient than using UFL in a parallel environment? Can you provide some suggestionsďź

This bit can be improved by using dolfinx.fem.petsc.create_matrix instead of using assemble. One could even create a custom sparisty pattern here with diagonal only.

Consider the following as an illustration (using DOLFINx v0.8.0):

from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.fem import functionspace
from dolfinx.cpp.la import SparsityPattern, petsc
import ufl
from ufl import dx, action
from petsc4py import PETSc
import numpy as np
from time import perf_counter
N = 200
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N, mesh.CellType.triangle)
V = functionspace(domain, ("Lagrange", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
mass_form = v*u*dx
option = 2
if option == 1:
mass_action_form = ufl.replace(mass_form, {u:fem.Constant(domain, PETSc.ScalarType(1))})
elif option == 2:
Q = functionspace(domain, ("DG", 0))
one = fem.Function(Q)
one.x.array[:] = 1
mass_action_form = action(mass_form, one)
else:
raise ValueError("Invalid option")
def create_lumped_matrix():
pattern = SparsityPattern(
domain.comm,
[V.dofmap.index_map, V.dofmap.index_map],
[V.dofmap.index_map_bs, V.dofmap.index_map_bs],
)
all_blocks = V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts
pattern.insert_diagonal(np.arange(all_blocks, dtype=np.int32))
pattern.finalize()
M_lumped = petsc.create_matrix(domain.comm, pattern)
M_lumped.zeroEntries()
diagonal = assemble_vector(fem.form(mass_action_form))
M_lumped.setDiagonal(diagonal)
M_lumped.assemble()
return M_lumped
start = perf_counter()
M_L = create_lumped_matrix()
end = perf_counter()
start_m = perf_counter()
M = assemble_matrix(fem.form(mass_form))
M.assemble()
end_m = perf_counter()
print(f"{end-start:.3e} s for lumped mass matrix")
print(f"{end_m-start_m:.3e} s for full mass matrix")

yielding

5.941e-03 s for lumped mass matrix
1.871e-02 s for full mass matrix

for N=200

and

3.031e-02 s for lumped mass matrix
1.264e-01 s for full mass matrix

for N=500.
Note that this is just comparing assembling the original matrix and the fully lumped version. We havenât even started the expensive operations you sketched above for making the lumped system: