Defining domains with high order coordinate fields
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.
Comments (18)
-
-
reporter Yes of course.
It has to contain the full coordinate field, i.e. it could in principle replace the mesh vertices completely.
-
reporter 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.
-
reporter - edited description
-
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.
-
reporter 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...
-
Creating a coordinate field from the mesh should be easy.
-
We could generalize MeshGeometry to be a Function or a Mapping (which allows also change of coordinates for instance).
-
reporter 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.
-
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.
-
reporter - changed milestone to 1.5
Waiting for release so we can remove some deprecated code which will simplify the implementation.
-
reporter 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/
-
reporter 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?
-
reporter - removed milestone
-
reporter - removed responsible
-
reporter -
assigned issue to
- edited description
- changed milestone to 1.7
- changed version to dev
-
assigned issue to
-
reporter - changed status to resolved
-
reporter - removed milestone
Removing milestone: 1.7 (automated comment)
- Log in to comment
Cool! But shouldn't:
be:
Also does the
x
Coefficient
contain all the vertex positions or does it just contain the displacement for the higher order vertices?