Bootstrapping

[miklos1]

A mesh has coordinates which is a function. A function has a function space. A function space has a mesh. A mesh has coordinates which is a function. A function has a function space. A function space has a mesh. A mesh has coordinates which is a function. A function has a function space. A function space has a mesh. A mesh has coordinates which is a function. ...

While this circular dependency is very nice for certain reasons (coordinates are no special), our implementation in Firedrake barely works, and it actually does not work in some cases (periodic meshes with fancy function spaces, higher-order meshes). This patch aims to change that by using robust design that works.

Design principles

The idea is to break the cyclic dependency under the hood, by making a distinction between topological (without coordinates) and geometric (with coordinates) versions of mesh, function space, and function.

A mesh topology is self-contained. One can define coordinateless function space on a mesh topology, and then a coordinateless function on a coordinateless function space. A coordinateless function is a perfect candidate for describing the coordinate field; in fact, a coordinateless function has all the data to define a mesh geometry. Then a function space is a coordinateless function space with a mesh geometry, and a function is a coordinateless function with a mesh geometry. (A function can also be interpreted as a pair of coordinateless functions, one for the coordinates and one for the actual values.)

- But then the coordinate field of a mesh is special again?!

It is not that special (it can still be of any function space), and it is all under the hood. Using the "official" API, a mesh still returns a regular function as its coordinates. By creating a function which references the the actual coordinate field (which is coordinateless function) both as values and coordinates, we are even memory efficient.

In fact, mesh.coordinates.function_space().mesh() is mesh still holds.

- Martin may not like it.

This is about refactoring Firedrake, and Martin does not care about that. In fact, these changes make it easier to cooperate with current and future UFL, since it is mostly parallel to UFL's design. ('current' and 'future' in our reference frame.)

API changes

I argue for the removal of the "changing mesh coordinates by assignment" approach, in favour of a referentially transparent solution. While the new infrastructure makes it easier to change mesh coordinates in proper way, it would require more work to avoid the invalidation of existing function spaces and functions on that mesh. (Note: this is not about changing the coordinate values, but about changing the coordinate function, generally with a function which has a different function space.)

As far as I know, there are currently two uses of changing coordinates:

1. Creating a mesh from an arbitrary coordinate field

There is a new API to achieve this. If f is a Function, then Mesh(f) returns a mesh geometry whose coordinate field is the values of f. No data is duplicated.

2. Changing the coordinates of an existing function

Currently, there is no public API to do this. However, thanks to the flexibility of the design, if f is a Function and mesh is a mesh geometry, then

Function(functionspace.WithGeometry(f.function_space(), mesh), val=f.topological)

will return a new Function object which refers to f's data values and has mesh as its geometry. Again, no data duplication.

I can add a public API if we agree what that should be.

Required changes outside Firedrake

Spurious UFL check in MixedElement

If a ufl.MixedElement is reconstructed over another domain, it checks that the value shape of the subelements remains the same. This does not hold, if one of the subelements is vector-valued (e.g. RT, BDM, etc.) and the geometric dimension changes.

UFL cell reconstruction does not work for OuterProductCells

Who would have thought that a cell name and geometric dimension is not enough to create a cell? (Or what is a ufl.Cell('OuterProductCell', 3)?)

ufl.OuterProductCell equality is too strict

Currently we use the following ufl.OuterProductCells:

OuterProductCell(Cell("interval"), Cell("interval"))
OuterProductCell(Cell("interval", 2), Cell("interval"))
OuterProductCell(Cell("interval", 2), Cell("interval"), gdim=3)
OuterProductCell(Cell("interval", 3), Cell("interval"))
OuterProductCell(Cell("triangle"), Cell("interval"))
OuterProductCell(Cell("triangle", 3), Cell("interval"))
OuterProductCell(Cell("quadrilateral"), Cell("interval"))
OuterProductCell(Cell("quadrilateral", 3), Cell("interval"))

After this patch, we use the following ufl.OuterProductCells:

OuterProductCell(Cell("interval"), Cell("interval"))
OuterProductCell(Cell("interval"), Cell("interval"), gdim=3)
OuterProductCell(Cell("triangle"), Cell("interval"))
OuterProductCell(Cell("quadrilateral"), Cell("interval"))

We should either loosen equality in UFL, or simplify ffc/codesnippets.py, firedrake/io.py, and who knows what else...

[dham]

I don't understand the choice of use cases. Isn't the primary use case that a mesh exists and we want to replace its coordinates (and have this apply automagically to all of the FunctionSpaces and therefore Functions which live on this mesh?)

[miklos1]

@dham: We might do that, but mutation is evil, and therefore to be avoided unless necessary. Considering the implementation of that, it would not be a great fuss for FunctionSpace, but it would be for Function, primarily because a Function is a ufl.Coefficient, which has to be initialised with a UFL element. The re-initialisation of the ufl.Coefficient would require inappropriate intimacy from Function towards ufl.Coefficient, and even more so to update automagically.

[dham]

OK. I can buy that. However this implies that someone needs to add some spectacularly clear documentation to the manual to explain how this works. (E.g. That changing the value of the coordinates is allowed but that changing the function space creates a new mesh).

[wence-]

I believe this breaks functionality that some people in Leeds are using. Why do all function spaces need to be defined on the same mesh?

[miklos1]

Currently every function space class is topological, and WithGeometry adds a mesh geometry. Therefore MixedFunctionSpace is topological, so it has topological function spaces as subspaces. A WithGeometry can surround a MixedFunctionSpace to add a mesh geometry.

I didn't know Firedrake can solve problems on different meshes. Can we add a test that exercises the feature that the Leeds people use?

[wence-]

Wait, what? A mesh doesn't have a cell?

[miklos1]

It does, but ufl_cell() now triggers init().

[wence-]

Superseded by #636, but I think this is now fundamentally sound.