Status and design ideas for parameterized geometries

Next steps towards parameterized geometries

At least for the foreseeable future we'll stick to coordinate elements where the dofs are point values. In particular we'll use only continuous Lagrange, so the MeshGeometry only needs a degree() property added.

I suggest then a model of MeshGeometry that is a small generalization of the 'array of vertex coordinates' model it is today: just make it an array of point coordinates, where the vertices are a subset of the points.

Here's an attempt at making a task list for the next refactoring steps towards the endgoal:

Task A:

Distinguish between vertices and other points in MeshGeometry interface. (Vertices is a subset of points, and for affine mesh they're the same.)

1) Add functions and use these instead of x() and size() where it's clear that the vertices are what's intended:

#!c++
      num_vertices()
      vertex_coordinates(vertex_index)

2) Add functions and use these where access to x() and size() should include all points:

#!c++
      num_points()
      point_coordinates(std::size_t point_index)

3) Add access to points which are defined on entities

#!c++
      std::size_t index = mesh.geometry().get_entity_index(dim, order, entity_index);
      // index is the same as entity_index for vertices (dim = 0)
      Point p = mesh.geometry().point(index);

4) Add function to Cell

#!c++
      void Cell::get_coordinate_dofs(vector<double>& coordinate_dofs) const;

4) Deprecate x(), size(), etc. with message telling about the new functions. These old functions can map to vertex_coordinates etc.

5) Mesh has some interface for accessing coordinates without going through geometry(), similar changes are needed there.

This can be done without breaking anything and will help give us a more seamless transition.

Task B:

Make it possible to initialize mesh with degree > 1.

1) Extend MeshEditor::open to take a degree argument, and pass it on to make MeshGeometry::degree() return this instead of 1:

#!c++
  MeshEditor editor;
  editor.open(mesh, tdim, gdim, degree);

  // All cells must be initialised before entities (below)
  editor.init_cells(n_cells);
  editor.add_cells(i, i0, i1, i2, i3);

  // Vertices can be initialised before or after entities (below)
  editor.init_vertices(n_vertices);
  editor.add_vertex(i, point);

  // For degree > 1, initialise Edges (and possibly Faces too) so that points can be added to them
  editor.init_entities();
  ...
  // Add a point on an entity. For entities with multiple points, "order" is the indexing on that entity 
  editor.add_entity_point(entity_dimension, order, entity_index, point); 
  ...
  editor.close();
  assert(mesh.geometry().degree() == degree);

Numbering for P1,P2 is trivial, for P3 is fairly easy (two points on edge is easy to order), but for P4 or higher in 3D it might be best to use e.g. ufc::dofmap::tabulate_dof_coordinates. Is this a likely use scenario?

Task C:

Uflacs updates.

1) If necessary, fix computation of x and J which may be slightly broken for degree > 1 because there are no working tests atm.

Task D:

UFC updates.

1) Renaming vertex_coordinates to coordinate_dofs.

2) Add ufc::domain class (see issue).

Task E etc:

1) I/O - some work already done for XDMF/HDF5 output - needs cleaning up

2) Read from file (which formats should be supported?)

3) plotting - low priority

4) interpolation?

5) collision detection / point-in-cell with curved boundaries - could be tricky

Concepts

Coordinate element The spatial coordinates of the geometric domain are represented by a mapping from reference coordinates on each cell, represented in an arbitrary finite element basis. This leads to a new property of the Mesh: the coordinate element. The value size of the vector valued coordinate element is the geometric dimension. A minimal initial development can be adding a degree to dolfin Mesh and only support Lagrange.

Relation to isoparametric elements The "iso" in "isoparametric elements" refers to that the same element is used for representation of geometry and function. This feature is more generic, as the representation of geometry is parameterized by any finite element independently of the functions in a form. What's a good name for this? "Parametric geometry", "finite element parameterized geometry", "parameterized geometry"?

Relation to NURBS Generalizing this feature to NURBS etc later involves a number of other changes to mesh, dofmaps, functions, etc. The introduction of "coordinate dofs" as a generalization of "vertex coordinates" is one step that will probably ease the implementation of NURBS geometries later. Apart from that lets keep NURBS out of the current work.

Coordinate dofs The more general concept "coordinate dofs" will replace "vertex coordinates" in UFC and related code in FFC/DOLFIN. For Lagrange degree 1 coordinate dofs and vertex coordinates coincide, possibly apart from conventions on interleaving or blockwise component dofs.

Numbering of coordinate dofs Numbering coordinate dofs the same way as other finite element dofs for each scalar component is crucial, as any other choice will lead to indirection and/or unnecessary duplication of precomputed finite element tables in the generated kernels. Organizing dofs as one block for each component (dofs for x ..., dofs for y ..., dofs for z ...) seems most convenient in the generated kernels for two reasons: it makes the coordinate function and fem functions work the same, and it will allow a dense inner product

x[i] = sum_j xdofs[i*fe_dim + j]*FE[q][j]

instead of

x[i] = sum_j xdofs[i + j*gdim]*FE[q][j]

which has better data locality.

Code generation

Uflacs already supports generation of tabulate_tensor for geometries described by arbitrary vector valued finite elements. This only needs to be adapted to iterative changes in UFL and UFC.

UFL specification

Currently a ufl.Domain can have a ufl.Coefficient to describe the coordinates. This will change so a ufl.Domain only provides the FiniteElement the coordinate field is parameterized by over a reference cell.

Tying this in with dolfin requires adding degree or element to dolfin::Mesh or dolfin::MeshGeometry. How the coordinate dofs are stored is orthogonal to all concerns in UFL/FFC/UFLACS.

DOLFIN/UFL integration in Python

It is crucial that the C++ dolfin::Mesh contains all information about the coordinate element such that the UFL FiniteElement can be constructed for a dolfin.Mesh that was initially created on the C++ side. The dolfin.Mesh currently 'has a' ufl.Domain, but will later 'be a' ufl.Mesh, mirroring how a dolfin.Function 'is a' ufl.Coefficient.

Assembly, DOLFIN/UFC integration in C++

Assembly has been possible with uflacs generated code for a year through a workaround where the mesh vertex coordinates are ignored by the generated kernel and a Function specifies the coordinate dofs instead. The plan is to change this behaviour and let DOLFIN pass coordinate dofs in place of vertex coordinates as a generalization of todays implementation.

The key point in the Assembler is to get the coordinate dofs from the Mesh and pass them on to ufc instead of the vertex coordinates. Looking at the code in Assembler.cpp, the vertex_coordinates are fetched via this function:

#!C++
class Cell {
    // FIXME: This function is part of a UFC transition
    /// Get cell vertex coordinates
    void get_vertex_coordinates(std::vector<double>& coordinates) const
    {
      const std::size_t gdim = _mesh->geometry().dim();
      const std::size_t num_vertices = this->num_vertices();
      const unsigned int* vertices = this->entities(0);
      coordinates.resize(num_vertices*gdim);
      for (std::size_t i = 0; i < num_vertices; i++)
        for (std::size_t j = 0; j < gdim; j++)
          coordinates[i*gdim + j] = _mesh->geometry().x(vertices[i])[j];
    }
}

which instead needs to get the coordinate dofs for a cell. See note on numbering of coordinate dofs above.

DOLFIN Mesh class

The main difficulty on the DOLFIN side is that any place that uses the MeshGeometry coordinates currently assumes an affine mesh. Key question: should the .x() methods of MeshGeometry still refer to vertex coordinates, with a new function .coordinate_dofs() added? Or should .x() be changed to return coordinate dofs? Adding a new function MeshGeometry::coordinate_dofs() seems a good choice as an initial step, allowing fixing users of .x() independently. Global numbering also enters the picture here.

It seems there is agreement that returning a pointer to the vertex storage is a bad idea and should be removed. An alternative interface to editing coordinates must be provided in place. For parameterized geometries, we want some way to copy the coordinates to/from (and increment by) a Function in a compatible FunctionSpace. We probably also need an interface to set a single vertex coordinate. These editing interfaces all need to invalidate the BBTree that the Mesh contains, which would solve some other issues there.

DOLFIN Mesh file support

Simone has python prototype code for reading quadratic (and cubic?) gmsh geometries into dolfin Functions.

Chris already has prototyped storing of quadratic meshes to XDMF files (see XDMFFile::write_quadratic)

Issues of consistent global and local numbering.

Less essential parts (but still important)

Plot support Mesh and Function. Can be done by using VTK quadratic elements and/or refining the mesh and interpolating before plotting.

Refinement Topology is the same, but must use the geometric parameterization to determine geometric coordinates and not the coordinate dofs directly. Update to CellEntity::midpoint()

Other? List places where .x() is used here? MeshSmoothing, MeshTransformation, most CellTypes, Cell, Vertex, SubDomain, LocalMeshData, ALE, GenericBoundingBoxTree, PointSource, Expression

Most of these can easily be changed to accommodate changes to underlying data.

Dependencies

See also FunctionSpace_in_UFL.