MeshView aware assembler

I'm dumping some notes here on what assembler needs to handle for future meshview implementation, to get the discussion started again.

Some details here are likely to be missing or maybe wrong, focus is on the high level concepts and feasibility. All code is (obviously) pseudocode.

# (1) Loop over groups of relevant cells, facets or
# other entities needed for one tabulate_tensor
# (example: entities = (cell,) for cell_integral)
# (example: entities = (facet, cell) for exterior_facet_integral)
# (example: entities = (facet, cell0, cell1) for interior_facet_integral)
# where facet = (tdim-1, facet_index)
# where cell = (tdim, cell_index)
for entities in integration_context:

    # (2) Get mesh coordinate dofs for all mesh cells involved.
    # (Meshviews do not define their own coordinates)
    for (d, cell_index) in entities if d == tdim:
            xdofs += mesh.coordinate_dofs(cell_index)

    # (3) Loop over all entities and all unique dofmaps on each entity:
    for (d, e) in entitites:
        # (3a) dofmaps_on_entity_dim is built before the loop over entities
        # and holds basic dofmaps defined only on one mesh or meshview
        for dofmap in dofmaps_on_entity_dim[d]:
             # (3a) Examples:
             # for a dofmap on a mesh or cell-meshview, d == tdim and e is a cell number
             # for a dofmap on a facet-meshview, d == tdim-1 and e is a facet number
             dofs += dofmap.cell_dofs(e)

    # (4) Restrict all functions to relevant dofs
    for each function f:
        # (4a) If a function can live in a mixed function space
        # with components living on different meshviews,
        # this iterates over the basic dofmaps for separate components
        for dofmap in function space of f:
            # (4b) Need to get e.g. function on both cell0 and cell1 for dS
            for (d, e) in entities if d == dofmap.tdim():
                # (4c) Reusing dofmap.cell_dofs(e) from (3) is probably more efficient
                w += f.restrict(dofmap.cell_dofs(e))

    # (5) Compute element tensor, the main difference here
    # is that we must agree with generated code on the
    # ordering of xdofs and w in the new situations.
    integral.tabulate_tensor(A, xdofs, w, [facet, ...])

    # (6) Insert A using dofmaps from (3)

Here I'm making a few assumptions:

A MeshView has a parent Mesh, a dim <= mesh.tdim, and a mapping of entity indices to/from the Mesh
A DofMap only lives on a single Mesh or MeshView
Functions that have one component on the boundary and one on the cells will need one DofMap for each
A DofMap on facets of a Mesh will take facet indices to its dm.cell_dofs(index)
A form still only involves one mesh, even if multiple MeshViews are involved (i.e. not covering multimesh situation here)

There are a lot of other design and implementation questions hidden in the details, some of them are:

Should the DofMap class cover both the Mesh and MeshView case?

Should there be a MixedFunctionSpace class with multiple DofMaps?

Should there be a MixedFunction class with one MixedFunctionSpace or multiple FunctionSpace?

etc.

Another central choice is that in (6) we need to know which dofmaps to use. The new thing is that test and trial functions may have components not active everywhere, and this varies per integral. Do we let each integral object communicate which dofmaps (blocks) to use for this A, or do we compile forms separately for each block (split forms on test/trial function components based on meshview id).

Comments (2)