Mass Lumping in dolfinx

[maxre]

While porting my code from dolfin to dolfinx, it looks like the way to implement mass-lumping does not work anymore.

Back in dolfin one would have written:

dxL = dx(metadata={'quadrature_degree': 1, 'quadrature_scheme': 'vertex'})

In my MWE however the printed mass matrix is not diagonal as it should be, but does not even take the provided metadata into account…
MWE:

 import numpy as np
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import VectorFunctionSpace, Function
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.triangle)
V = VectorFunctionSpace(domain, ("Lagrange", 1), dim =2)
from dolfinx import fem
import ufl
u = ufl.TrialFunction(V)
u0 = Function(V)
v = ufl.TestFunction(V)

# Mass Lumping
dxL = ufl.dx(metadata={'quadrature_degree': 1, 'quadrature_scheme': 'vertex'}) 
Eqa = ufl.inner(u, v)*ufl.dx
Eqb = ufl.inner(u,v)*dxL
Xa= assemble_matrix(fem.form(Eqa))
Xb= assemble_matrix(fem.form(Eqb))
Xa.assemble()
Xb.assemble()
Xa.convert("dense")
Ca = Xa.getDenseArray()
Xb.convert("dense")
Cb = Xa.getDenseArray()
print("w/o mass lumping \n",Ca)
print("w/ mass lumping \n",Cb)

Output:

w/o mass lumping 
 [[0.00390625 0.         0.00065104 ... 0.         0.         0.        ]
 [0.         0.00390625 0.         ... 0.         0.         0.        ]
 [0.00065104 0.         0.00130208 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.00390625 0.         0.00065104]
 [0.         0.         0.         ... 0.         0.00130208 0.        ]
 [0.         0.         0.         ... 0.00065104 0.         0.00130208]]
w/ mass lumping 
 [[0.00390625 0.         0.00065104 ... 0.         0.         0.        ]
 [0.         0.00390625 0.         ... 0.         0.         0.        ]
 [0.00065104 0.         0.00130208 ... 0.         0.         0.        ]
 ...
 [0.         0.         0.         ... 0.00390625 0.         0.00065104]
 [0.         0.         0.         ... 0.         0.00130208 0.        ]
 [0.         0.         0.         ... 0.00065104 0.         0.00130208]]

What am I missing?

[dokken]

@maxre

Xb.convert("dense")
Cb = Xa.getDenseArray()

You have a typo here.
This should be

Xb.convert("dense")
Cb = Xb.getDenseArray()

[maxre]

Thanks for noticing the typo! What a rookie mistake on my side…
Also, I was able to solve the problem myself in the meantime based on the following topic:

Set integration rule to closed/vertex

I want to set the integration rule to closed for all terms of my problem and not only the mass matrix. I have a MWE on the unit interval with P1 elements and a term f = x from dolfinx import mesh, fem import ufl from mpi4py import MPI nx = 2 domain = mesh.create_unit_interval(MPI.COMM_WORLD, nx) V = fem.FunctionSpace(domain, ("CG", 1)) x = ufl.SpatialCoordinate( domain ) function = x[0] v = ufl.TestFunction(V) b = function*v*ufl.dx linear_form = fem.form(b) b = fem.petsc.create_vector(lin…

The correct way would be:
dxL = ufl.Measure("dx", metadata = {"quadrature_rule": "vertex"})
I am not fully sure about further meta-settings. But for now it works in my application.

The whole MWE that compares the methods.

import numpy as np
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.fem import VectorFunctionSpace, Function 

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.triangle)
V = VectorFunctionSpace(domain, ("Lagrange", 1), dim =2)
from dolfinx import fem
import numpy
import ufl
u = ufl.TrialFunction(V)
u0 = Function(V)
v = ufl.TestFunction(V)
# Mass Lumping
dxL = ufl.Measure("dx", metadata = {"quadrature_rule": "vertex"})
Eqa = ufl.inner(u, v)*ufl.dx
Eqb = ufl.inner(u,v)*dxL

Xa= assemble_matrix(fem.form(Eqa))
Xb= assemble_matrix(fem.form(Eqb))
Xa.assemble()
Xb.assemble()
Xa.convert("dense")
Ca = Xa.getDenseArray()
Xb.convert("dense")
Cb = Xb.getDenseArray()
print("w/o mass lumping \n",Ca)
print("w/ mass lumping \n",Cb)
print("difference: ", np.max(np.abs(Ca-Cb)))
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

print("only values on the diagonal: ",only_diagonal)

[Holzi]

Hey there,
I also wanted to implement mass lumping and found this discussion. While this works nicely in 2D, it seems to make issues in 1D. In this MWE

mesh_1d = df.mesh.create_unit_interval(MPI.COMM_WORLD,10)
V_1d = df.fem.FunctionSpace(mesh_1d,("Lagrange",1))

u = ufl.TrialFunction(V_1d)
v = ufl.TestFunction(V_1d)

lumped_integration_dx = ufl.Measure("dx", metadata={"quadrature_rule": "vertex"})
mass_form = df.fem.form(ufl.inner(u,v)* lumped_integration_dx)
MassMatrix = df.fem.petsc.assemble_matrix(mass_form, bcs=[])
MassMatrix.assemble()

I get the following traceback (running dolfinx 0.5.2 on Ubuntu22)

Traceback (most recent call last):
  File "/home/user/Workspace/domainsplitting/src/utilities.py", line 166, in <module>
    MWE_mass_lumping()
  File "/home/user/Workspace/domainsplitting/src/utilities.py", line 118, in MWE_mass_lumping
    mass_form = df.fem.form(ufl.inner(u,v)* lumped_integration_dx)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 166, in form
    return _create_form(form)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 161, in _create_form
    return _form(form)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 135, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 168, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/jit.py", line 232, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, parameters=parameters)
  File "/usr/lib/python3/dist-packages/ffcx/compiler.py", line 107, in compile_ufl_objects
    code = generate_code(ir, parameters)
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/codegeneration.py", line 50, in generate_code
    code_integrals = [integral_generator(integral_ir, parameters) for integral_ir in ir.integrals]
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/codegeneration.py", line 50, in <listcomp>
    code_integrals = [integral_generator(integral_ir, parameters) for integral_ir in ir.integrals]
  File "/usr/lib/python3/dist-packages/ffcx/codegeneration/integrals.py", line 27, in generator
    logger.info(f"--- type: {ir.integral_type}")
AttributeError: 'numpy.ndarray' object has no attribute 'integral_type'

Are you aware of any work arounds for that instead of manually summing up the rows (this should be equivalent in 1D)?

[dokken]

I’ve fixed this in: Fix interval integrals with vertex rules by jorgensd · Pull Request #562 · FEniCS/ffcx · GitHub