# Mass Lumping in dolfinx

**URL:** <https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916>\
**Category:** dolfinx\
**Created:** [April 12, 2023, 4:14pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916 "2023-04-12T16:14:59Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [April 12, 2023, 4:14pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916/1 "2023-04-12T16:14:59Z")

</div>

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:

```auto
 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:

```auto
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?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 12, 2023, 4:58pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916/2 "2023-04-12T16:58:13Z")

</div>

> [@maxre](#):
>
> ```auto
> Xb.convert("dense")
> Cb = Xa.getDenseArray()
> 
> ```

You have a typo here.  
This should be

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

```

---

<div class="post-metadata">

**Author:** ![maxre](https://avatars.discourse-cdn.com/v4/letter/m/f08c70/32.png) [@maxre](https://fenicsproject.discourse.group/u/maxre)\
**Post date:** [April 12, 2023, 6:09pm UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916/3 "2023-04-12T18:09:19Z")

</div>

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](https://fenicsproject.discourse.group/t/set-integration-rule-to-closed-vertex/10695):
>
> 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.

```auto
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)

```

---

<div class="post-metadata">

**Author:** ![Holzi](https://avatars.discourse-cdn.com/v4/letter/h/5f9b8f/32.png) [@Holzi](https://fenicsproject.discourse.group/u/Holzi)\
**Post date:** [June 12, 2023, 7:43am UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916/4 "2023-06-12T07:43:04Z")

</div>

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

```python
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)

```auto
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)?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 12, 2023, 8:34am UTC](https://fenicsproject.discourse.group/t/mass-lumping-in-dolfinx/10916/5 "2023-06-12T08:34:42Z")

</div>

I’ve fixed this in: [Fix interval integrals with vertex rules by jorgensd · Pull Request #562 · FEniCS/ffcx · GitHub](https://github.com/FEniCS/ffcx/pull/562)