introduce MixedMesh class

[ksagiyam]

First half of https://github.com/FEniCS/ufl/pull/264

Introduce MixedMesh class for multi-mesh problems.

Note that extract_domains() now uses more robust sort_domains() inside, so it might behave slightly differently.

Edit 12-09-2024:

MixedMesh class represents a collection of meshes (e.g., submeshes) that, along with a MixedElement, can represent a mixed function space defined across multiple domains. The motivation is to allow for treating Arguments and Coefficients on a mixed function space defined across multiple domains just like those on a mixed function space defined on a single mesh.

Specifically, the following becomes possible (see tests for more):

    cell = triangle
    mesh0 = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
    mesh1 = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
    domain = MixedMesh([mesh0, mesh1])
    elem0 = FiniteElement("Lagrange", cell, 1, (), identity_pullback, H1)
    elem1 = FiniteElement("Lagrange", cell, 2, (), identity_pullback, H1)
    elem = MixedElement([elem0, elem1])
    V = FunctionSpace(domain, elem)
    v = TestFunction(V)
    v0, v1 = split(v)

For now, one can only perform cell integrations when MixedMeshes are used. This is because, e.g., an interior facet integration on a submesh may either be interior or exterior facet integration on the parent mesh, and we need a logic to apply default restrictions on coefficients defined on the participating meshes. This is the second half of https://github.com/FEniCS/ufl/pull/264.

Also, currently, all component meshes must have the same cell type (and thus the same topological dimension) -- we are to remove this limitation in the future.

Core changes:

• GradRuleSet.{reference_value, reference_grad} work component-wise (component of the mixed space) if the FunctionSpace is defined on a MixedMesh, so that each component is associated with a component of the MixedMesh, say domain_i (JacobianInverse(domain_i) is then well defined).
• extract_arguments_and_coefficients is now called extract_terminals_with_domain, and it now also collects GeometricQuantys, so that we can correctly handle, say SpatialCoordinate(mesh1) and SpatialCoordinate(mesh2), in problem solving environments.