Defining domains with high order coordinate fields

[Martin Sandve Alnæs]

THIS DESCRIPTION IS DEPRECATED! But the plan remains.

Work in process, this is the syntax:

D = Domain(triangle)
assert D.coordinates() is None
P2 = VectorElement("CG", D, 2)
x = Coefficient(P2)

E = Domain(x)
assert E.coordinates() is x
V = FiniteElement("CG", E, 2)
f = Coefficient(f)

M = f*dx(E)

There will be algorithms for translating Grad(f) to (with new types) JacobianInverse(E)*LocalGrad(f).

Initially, Jacobian(E), JacobianInverse(E) and JacobianDeterminant(E) will be mapped in ffc to existing snippets J, detJ, K.

Later we can translate Jacobian(E) to LocalGrad(E.coordinates()) at an appropriate stage.

Since x is a regular Coefficient, LocalGrad(x) can easily be mapped to existing precomputed tables in FFC for evaluation of the Jacobian in quadrature points.

[Johan Hake]

Cool! But shouldn't:

P2 = FiniteElement("CG", D, 2)

be:

P2 = VectorElement("CG", D, 2)

Also does the x Coefficient contain all the vertex positions or does it just contain the displacement for the higher order vertices?

[Martin Sandve Alnæs]

Yes of course.

It has to contain the full coordinate field, i.e. it could in principle replace the mesh vertices completely.

[Martin Sandve Alnæs]

The ufl representation part is now in place. Algorithms for global to local gradients with jacobians are still left to do, plus the ffc side of things. I don't think anything has to happen on the dolfin side.

[Simone Pezzuto]

This is very interesting. I wonder how to connect DOLFIN to the generate code then. Do you have something in mind?

A possibility is to store the extra information into MeshData, with some convention on the ordering, and then use it when necessary. The functionality can be also extended to the exporters and the plot.

[Martin Sandve Alnæs]

I haven't thought too much about how to connect stuff on the C++ side of dolfin, only that the assembly will Just Work when it's done. The firedrake people have an easy time here, since they already store the coordinates of the mesh as a function.

Which operations in dolfin will need this information? I guess basically all things that use mesh coordinates...

[Anders Logg (Chalmers)]

Creating a coordinate field from the mesh should be easy.

[Simone Pezzuto]

We could generalize MeshGeometry to be a Function or a Mapping (which allows also change of coordinates for instance).

[Martin Sandve Alnæs]

There's a circular dependency issue in that a Function depends on a FunctionSpace which depends on a Mesh which depends on MeshGeometry. So making MeshGeometry into a Function would be a rather delicate operation.

[Simone Pezzuto]

What about extending the Mesh? Something like HighOrderMesh or MappedMesh or MorphedMesh which is a tuple (Mesh, Function) and it overrides all the coordinate-specific methods.

Otherwise the circular dependency can be addressed either by fixing the class design or by exploiting some metaprogramming trick (e.g. Barton-Nackman trick). But I think DOLFIN relies on dynamic polymorphism rather than static.

[Martin Sandve Alnæs]

Waiting for release so we can remove some deprecated code which will simplify the implementation.

[Martin Sandve Alnæs]

How to solve the circular dependency issue needs to be resolved for UFL and DOLFIN together before any real progress can be had on this issue. These are the related issues in dolfin:

https://bitbucket.org/fenics-project/dolfin/issue/307/ https://bitbucket.org/fenics-project/dolfin/issue/321/

[Martin Sandve Alnæs]

I think I have the solution. It is not to represent the coordinate function as an object that the mesh is created from. On the contrary, the mesh is created with an element.

# Possibly future UFL code:
cell = triangle
element = VectorElement("Lagrange", cell, 2)
mesh = Mesh(element) # Counted for unique ids just like Coefficient

# To support legacy code, whenever a cell is encountered it is interpreted as
#mesh = Mesh(VectorElement("Lagrange", cell, 1), count=-1) # Fixed count -1 like the old TestFunction design

# Like before, if you need to access x in the form, this is the notation:
x = SpatialCoordinates(mesh)

# For eventual future shape derivatives, special rules need to be implemented for the coordinate function anyway:
M = f*dx
L = derivative(M, x)

# In the form preprocessing, the SpatialCoordinates(mesh) for each mesh can still be
# handled as Coefficients as much as possible, maybe by creating Coefficient objects
# internally to let the form compiler keep it simpler.


# In the dolfin python layer, the inheritance pattern and construction can look something like this:
class dolfin.Mesh(ufl.Mesh, cpp.Mesh):
    def __init__(self, *args):
        cpp.Mesh.__init__(self, *args)
        x = self.coordinate_function() # Not actually a UFL object
        element = x.finite_element()
        ufl.Mesh.__init__(self, element, count=self.id())

Although this idea is fresh and there may be hidden issues, this design intuitively feels much better. Can anyone spot any problems?

[Martin Sandve Alnæs]

Removing milestone: 1.7 (automated comment)