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.

[ksagiyam]

@jpdean @jorgensd This is what we briefly discussed at PDESoft. Could you have a look?

[garth-wells]

It would be helpful to have a PR summary, e.g. what it sets out the achieve, what it supports, designs considered, design rationale, implementation summary, and what a MixedMesh is and how a MixedMesh is different from a Mesh.

[ksagiyam]

Updated PR summary.

Added more checks.

Fixed MixedPullback.physical_value_shape().

[ksagiyam]

Can I have another round of review on this?

[jhale]

Copied from Slack:

The implementation looks OK to me (quick look), but the tests look a bit light. In terms of the naming, I was also confused by the use of Mixed. At least when we’re talking about MixedElement, aside from it having historical meaning, my impression was that we are talking about using the Cartesian product on the elements of the list and maintaining ordering. My understanding of MixedMesh is that we are maintaining the ordering, but I didn’t get the impression it was a Cartesian product - consequently, the use of Mixed wasn’t clear to me. This is aside from the idea that Mixed Mesh has no traditional meaning in the literature - when the idea was first proposed I thought it was related to having a mesh with mixed cell type, to illustrate the confusion.

Of course these are ‘misconceptions’ on my part, but I’m just trying to get across the confusion that appeared in my own mind, so we can discuss if different naming of the same concept might be appropriate.

[jorgensd]

    cell = triangle
    elem0 = FiniteElement("Lagrange", cell, ...)
    mesh0 = Mesh(FiniteElement("Lagrange", cell, ...))
    V0 = FunctionSpace(mesh0, elem0)

    elem1 = FiniteElement("Brezzi-Douglas-Marini", cell, ...)
    mesh1 = Mesh(FiniteElement("Lagrange", cell, ...))
    V1 = FunctionSpace(mesh1, elem1)

    elem2 = FiniteElement("Discontinuous Lagrange", cell, ...)
  

This is what we use in DOLFINx at the moment.

I think the current mixed-mesh abstraction is going to be confusing when we introduce meshes with multiple cell types (ie mix of quadrilateral and triangle, or hex, yet, prism).

[ksagiyam]

My plan is to introduce MixedCell that should be shared by the MixedMesh and the MixedElement for mixed problems with multiple cell types.

We could have: single mesh + single space, single mesh + multiple spaces, and multiple meshes + multiple spaces. Current MixedFunctionSpace approach handles the last case differently inside UFL (and inside DOLFINx, I believe). Simply put, MixedMesh allows us to handle the last case just like the first two cases in UFL.

I think the problem of MixedFunctionSpace is that it behaves differently from the standard FunctionSpace, and requires special casing at high-level. For instance, we can not define Coefficient on a MixedFunctionSpace; we need to use the (special-cased) Coefficients function, which returns a list of Coefficients. We thus can not use high-level functions, such as inner and split, for coefficients on MixedFunctionSpaces. Components of an Argument should also be consistently represented by Indexed(argument, index) in UFL instead of by part, I think.

The MixedMesh approach keeps the set of basic building blocks, such as Mesh, FunctionSpace, and Coefficient, unchanged (note that MixedMesh subclasses AbstractDomain); e.g., we can do c = Coefficient(FunctionSpace(MixedMesh, MixedElement)), and directly use c in an expression. This allows us to not worry about mixedness in most parts of UFL and integrate mixed problem almost seamlessly. Indeed, in this PR, core changes only happened in GradRuleset.reference_value(), GradRuleset.reference_grad(), and pullback.apply(), which are low-level methods, and virtually no change was needed in high-level functions, and, to support mixed domain problem, Firedrake does not need to introduce new interface. The following is the basic usage in Firedrake:

mesh = ...
subm = Submesh(mesh, ...)
V0 = FunctionSpace(mesh, "CG", 1)
V1 = FunctionSpace(subm, "CG", 1)
V = V0 * V1
u = Function(V)
v = TestFunction(V)
v0, v1 = split(v)
dx1 = Measure("cell", domain=subm)
a = inner(u[0], v1) * dx1
A = assemble(a)

MixedMesh` might not be very intuitive, but I think it is the right abstraction.

[jorgensd]

For instance, we can not define Coefficient on a MixedFunctionSpace; we need to use the (special-cased) Coefficients function, which returns a list of Coefficients. We thus can not use high-level functions, such as inner and split, for coefficients on MixedFunctionSpaces. Components of an Argument should also be consistently represented by Indexed(argument, index) in UFL instead of by part, I think.

If the only thing that makes MixedMesh favorable is the possibility of creating a coefficient over all spaces, i believe we should rather make this possible with MixedFunctionSpace than adding a whole new class.

On the other hand, if everyone agrees that MixedMesh is the way to go, we should remove MixedFunctionSpace, and adapt the code in DOLFINx/FFCx to handle it. I see no benefit of having both code paths (as long as extract_blocks would work on a MixedMesh form.

and multiple meshes + multiple spaces. Current MixedFunctionSpace approach handles the last case differently inside UFL (and inside DOLFINx, I believe).

MixedFunctionSpace keeps a list of function spaces, in a similar fashion to what happens under the hood in the FunctionSpace(MixedMesh), thus I don't think the handling is too different in the two approaches.

When you say multiple meshes + multiple spaces, do you mean:

V = FunctionSpace(mesh, el0)
submesh = SubMesh(mesh, ...)
mixed_element = MixedElement([el0, el1])
Q = FunctionSpace(submesh, mixed_element)
W = MixedFunctionSpace(V, Q)

equivalent to:

    elem0 = FiniteElement("Lagrange", cell, ...)
    elem1 = FiniteElement("Brezzi-Douglas-Marini", cell, ...)
    mixed_element = MixedElement([el0, el1]) 
    elem = MixedElement([elem0, mixed_element])
    mesh0 = Mesh(FiniteElement("Lagrange", cell, ...))
    mesh1 = Mesh(FiniteElement("Lagrange", cell, ...))
    domain = MixedMesh([mesh0, mesh1])
    V = FunctionSpace(domain, elem)

or are you meaning something different?

My plan is to introduce MixedCell that should be shared by the MixedMesh and the MixedElement for mixed problems with multiple cell types.

I think we need a better name for this than "mixed", as a FunctionSpace of MixedElement on a MixedMesh with a "sub"-mesh that has a MixedCell sounds incredibly cryptic (which would be something like solve stokes equation (u, p) pair on a subdomain in your mesh, which consists of triangle and quadrilateral cells, coupled with a diffusivity coefficient living on the whole domain).

[dham]

My plan is to introduce MixedCell that should be shared by the MixedMesh and the MixedElement for mixed problems with multiple cell types.

How would this compose with a mesh which didn't have a unique cell type (because it contains triangles and quads, for example)?

[dham]

I think we should separate the name from the object. E.g. we could call it a MeshSequence or a MultiMesh.

I think the object disagreement here is caused by a long-standing difference of viewpoint about ufl.MixedFunctionSpace, which Firedrake has never used. We've always just created a FunctionSpace on a mixed element. This is partly because we distinguish between topological and geometric FunctionSpaces in order to break the inheritance cycle caused by the fact that Firedrake has first class coordinates (i.e. the mesh coordinates are a Function).

So, you can make MixedFunctionSpace a separate class, which has always been the FEniCS approach (at which point, why have a MixedElement class, you might ask) or you can use a unified FunctionSpace and push the mixed information down to the element, which has always been the Firedrake approach. In the latter case, once you have multiple domains you find yourself wanting a way of reasoning about them, which is facilitated by the class encompassing the meshes.

Basically, once you have mixed there are different levels in the dependency graph at which you can think about the objects coming together. I think this might be a case where we should just allow for both in order to let Firedrake and FEniCSx respectively just get on with it.

[dham]

My plan is to introduce MixedCell that should be shared by the MixedMesh and the MixedElement for mixed problems with multiple cell types.

How would this compose with a mesh which didn't have a unique cell type (because it contains triangles and quads, for example)?

Actually, if we follow the pattern, then the elements on a mixed topology domain are themselves a sort of mixed element. No problem.

[ksagiyam]

@jorgensd I do not think merely allowing for Coefficient(MixedFunctionSpace) will solve all issues; as long as MixedFunctionSpace does not subclass FunctionSpace and does not have ufl_domain, I think we will need to special case functions/methods in UFL for MixedFunctionSpace just as Coefficients function currently does. As I mentioned, introduction of MixedMesh allows UFL to be unaware of the mixed spaces in the most part, and it only needs to special case methods at the very bottom, i.e. GradRuleset.reference_value(), GradRuleset.reference_grad(), and pullback.apply().

I did not mean nested mixed spaces by "multiple meshes + multiple spaces". I just meant (meshA, elemA) x (meshB, elemB). Though we do not use nested mixed spaces in Firedrake, some functions such as flatten_domain_element are written with this in mind as you might have noticed. I do not think consideration of such nested mixed spaces would affect the mixed-problem design in UFL, though.

@jhale I think the use of Mixed in MixedMesh is consistent with the use of Mixed in MixedElement. If we consider assembling (i, j)-block of a matrix, we can view it as mixing mixed_mesh[i] and mixed_mesh[j] as well as mixing mixed_element[i] and mixed_element[j].

But I agree that it is confusing. As @dham suggested, MeshSequence or MultiMesh is also fine by me.

[jorgensd]

@jorgensd I do not think merely allowing for Coefficient(MixedFunctionSpace) will solve all issues; as long as MixedFunctionSpace does not subclass FunctionSpace and does not have ufl_domain, I think we will need to special case functions/methods in UFL for MixedFunctionSpace just as Coefficients function currently does. As I mentioned, introduction of MixedMesh allows UFL to be unaware of the mixed spaces in the most part, and it only needs to special case methods at the very bottom, i.e. GradRuleset.reference_value(), GradRuleset.reference_grad(), and pullback.apply().

I did not mean nested mixed spaces by "multiple meshes + multiple spaces". I just meant (meshA, elemA) x (meshB, elemB). Though we do not use nested mixed spaces in Firedrake, some functions such as flatten_domain_element are written with this in mind as you might have noticed. I do not think consideration of such nested mixed spaces would affect the mixed-problem design in UFL, though.

@jhale I think the use of Mixed in MixedMesh is consistent with the use of Mixed in MixedElement. If we consider assembling (i, j)-block of a matrix, we can view it as mixing mixed_mesh[i] and mixed_mesh[j] as well as mixing mixed_element[i] and mixed_element[j].

But I agree that it is confusing. As @dham suggested, MeshSequence or MultiMesh is also fine by me.

Let's not use MultiMesh, as that has been used in FEniCS before to mean hierarchies of meshes.

[jhale]

MeshSequence is clearer and at least to my understanding it is consistent with Python's abstract base class terminology collections.abc.Sequence.

@ksagiyam I see your analogy between MixedMesh and MixedElement at the form level, however, MixedElement was named after a widely established and accepted terminology from the finite element literature. We don't have this for MixedMesh, and so my preference is to take precise terminology from elsewhere.

[ksagiyam]

Thanks, all. Yes, let us use MeshSequence. The followings are the changes I have made according to the suggestions:

• Addressed @jorgensd 's suggestions,
• Renamed MixedMesh MeshSequence,
• Removed ufl_id from MeshSequence; I realised that my Firedrake branch broke after a recent rebase, and found that it was caused by distinguishing two MeshSequence instances made from the same set of Meshes. Thinking about this, I thought ufl_id in MeshSequence is redundant; please see updated MeshSequence.__repr__() and MeshSequence._ufl_hash_data_(). I also added a test to highlight this point.

@jhale What sort of tests do you think are missing? After this PR, I will soon make another PR that adds logics to handle restrictions (to handle cases where exterior facet integration of meshA is interior facet integration of meshB, etc). At that stage I will add one or two more examples that test facet integrations.

[ksagiyam]

@jhale @jorgensd Can we have this possibly merged?

[ksagiyam]

@jorgensd Do you have any further comments on this PR by any chance? I would like to have this PR go.

[jorgensd]

@jorgensd Do you have any further comments on this PR by any chance? I would like to have this PR go.

I’d like @garth-wells, @chrisrichardson and @michalhabera to have a look over this. If they don’t have time, I guess we should get towards merging this.

[ksagiyam]

Thanks, all.