Support dolfinx::fem::Expression with co-dim 0 submeshes

[jorgensd]

Describe new/missing feature

With the co-dim 0 support in assembly (and interpolation post #3177 ), it would be natural to extend dolfinx::fem::Expression to handle test/trial-functions and coefficients from different sub-meshes, linked by entity_maps.

There are two different ways this could be supported, either by passing in entity_maps at Expression construction
or through eval (which would be closer to what is done in #3111 and with #2730).
For either approach, we need to designate a "base"-mesh that we map all entities from.

Suggested user interface

mesh = ...
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
u = dolfinx.fem.Function(V)
mesh0, map0, _, _ = dolfinx.mesh.create_submesh(mesh, ..., i)
mesh1, map1, _, _ = dolfinx.mesh.create_submesh(mesh, ..., j)
u0 = dolfinx.fem.functionspace(mesh, ("DG", 2, (gdim, )))
u1 = dolfinx.fem.functionspace(mesh, ("N1curl", 1))
# Interface 1, would work nicely with interpolation
expr = dolfinx.fem.Expression(ufl.inner(u0, u1)*u, V.element.interpolation_points(),
                              domain=mesh, entity_maps={mesh0:map0_inverse, mesh1:map1_inverse})
# Interface 2, would require additional information in interpolation wrapper regarding entity maps.
expr_2 = dolfinx.fem.Expression(ufl.inner(u0, u1)*u, V.element.interpolation_points()
expr_2.eval(mesh, np.arange(num_cells_in_mesh), entity_maps={mesh0:map0_inverse, mesh1:map1_inverse})

u_new = dolfinx.fem.Function(V)
# Interpolation with interface 1
u_new.interpolate(expr)

# Interpolation with interface 2
u_new.interpolate(expr_2, domain=mesh, entity_maps={mesh0:map0_inverse, mesh1:map1_inverse})