More than one mesh tag (per integral type) for a single form

[pantolin]

Describe new/missing feature

In some problems it may be handy to define forms using different mesh tags (whose entities may intersect).
For instance, imagine that you want to build a form in which you assemble the same operator everywhere, but you have to include some extra terms for a few cells/facets (see the example below).
First possibility would be to define a single mesh tags object with two tags: for cells/facets with and without extra terms.
However, including additional terms to your form, the number of different cell/facet groups grows exponentially.

The most natural solution would be define several mesh tags and use them all in the same form (see the code below). Unfortunately, computing quantities with two (or more) MeshTags and the same IntegralType is not ideal as they cannot be combined in a single DOLFINx form.
Two (or more) forms need to be generated, which involves calling the assembler more than once and generating several (matrices, vectors) that need to be summed.

I am looking for possible alternatives:

1. Define mesh tags that contain more than one value per cell. This way a single MeshTags object would be needed and therefore 1 single form is created.
   Just looking at the code in MeshTags.h, I can’t see any particular reason why this shouldn’t work (apart from the built-in tests for checking for duplicates). I am afraid that MeshTags is used in a lot of places and extesions and breaking this “no duplicate entities” principle might blow up somewhere else.
   As mentioned in this slack discussion, another possibility would be associate every entity in the MeshTags to a std::pair of T

2. Allow more than one MeshTags to be used in the same DOLFINx form. At the end of the pipeline, the (C++) Form object simply stores (for each integral type) a container mapping tag ids to the associated entities (these ids and entities could initially come from different MeshTags objects). So again, apart from the Python code that checks that only one MeshTags object is present in the form, I don’t see any particular reason why this shouldn’t work.
   In fact, I introduced some changes in forms.py and everything worked smoothly.

Suggested user interface

from mpi4py import MPI
import ufl
from dolfinx import fem, mesh

domain = mesh.create_unit_square(MPI.COMM_WORLD, 2, 1, mesh.CellType.quadrilateral)

x = ufl.SpatialCoordinate(domain)
term_0 = 1 + x[0]
term_1 = 1 + x[1]
dx = ufl.dx(domain=domain)

# Current version
tags_0 = mesh.meshtags(domain, 2, np.array([0, 1]), np.array([0, 1]))
dx_0 = dx(subdomain_data=tags_0)

ufl_form_0 = term_0 * (dx_0(0) + dx_0(1)) + term_1 * dx_0(1)
form_0 = fem.form(ufl_form_0)

# Proposed version
tags_1 = mesh.meshtags(domain, 2, np.array([0, 1]), 1)
tags_2 = mesh.meshtags(domain, 2, np.array([1]), 2)

dx_1 = dx(subdomain_data=tags_1)
dx_2 = dx(subdomain_data=tags_2)

ufl_form_1 = term_0 * dx_1(1) + term_1 * dx_2(2)
form_1 = fem.form(ufl_form_1)  # Error