[BUG]: Problem when assembling different integrals for different tags

[pantolin]

Summarize the issue

I encountered this problem that I suspect it is a bug when calling the different kernels for different cell tags.

By setting the same quadrature degree for both measures, or setting a and b to the same value, everything works as expected.

How to reproduce the bug

Just run the MWE below with dolfinx v0.9.0

Minimal Example (Python)

from mpi4py import MPI

import numpy as np
import ufl
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import assemble_scalar, form
from dolfinx.mesh import create_unit_square, meshtags

msh = create_unit_square(MPI.COMM_WORLD, 1, 1, CellType.triangle)


# Setting different quadrature degrees and values for a and b, makes it fail
cells = np.array([0, 1])
tags = np.array([0, 1])
cell_tags = meshtags(msh, msh.topology.dim, cells, tags)
dx = ufl.dx(domain=msh, subdomain_data=cell_tags)
dx0 = dx(subdomain_id=(0, 1), degree=2)
dx1 = dx(subdomain_id=(0), degree=3)

a, b = 1.0, 2.0
ufl_form_0 = a * dx0
ufl_form_1 = b * dx1
ufl_form = ufl_form_0 + ufl_form_1

expected = 1.0 * a + 0.5 * b

computed_0 = assemble_scalar(form(ufl_form_0)) + assemble_scalar(form(ufl_form_1))
print(f"Two assemblies -- Expected: {expected}. Computed: {computed_0}")

computed_1 = assemble_scalar(form(ufl_form))
print(f"One assembly   -- Expected: {expected}. Computed: {computed_1}")

Output (Python)

Two assemblies -- Expected: 2.0. Computed: 2.0
One assembly   -- Expected: 2.0. Computed: 2.5

Version

0.9.0

DOLFINx git commit

No response

Installation

conda on osx arm64

Additional information

No response

[jorgensd]

I think the issue is this around these lines ( https://github.com/FEniCS/dolfinx/blame/f1c5f307c70f414872a090e349f1d88211f54b58/python/dolfinx/fem/forms.py#L369), as it seems like the order of the integrals here are different that the one in the generated code.

FFCx orders the b integral prior to the two a integrals, while the referenced code goes through it as
a (id 0), b (id 0), a (id 1)

[jorgensd]

So I think we can simplify the Python code to:

        # Extract subdomain ids from ufcx_form
        integral_offsets = [ufcx_form.form_integral_offsets[i] for i in range(4)]
        for i in range(3):
            integral_type = IntegralType(i)
            for j in range(integral_offsets[i], integral_offsets[i + 1]):
                subdomain_ids[integral_type.name].append(
                    ufcx_form.form_integral_ids[j]
                )

However, we would still have to rewrite some C++ logic, as it doesn't use the integral index (

dolfinx/cpp/dolfinx/fem/utils.h

Line 531 in f1c5f30

integrals.insert({{IntegralType::cell, id, form_idx},

).