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 parentMesh
, a dim <= mesh.tdim, and a mapping of entity indices to/from theMesh
-
A
DofMap
only lives on a singleMesh
orMeshView
-
Functions
that have one component on the boundary and one on the cells will need oneDofMap
for each -
A
DofMap
on facets of aMesh
will take facet indices to itsdm.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)
-
reporter -
- removed milestone
- Log in to comment
Lets discuss this more after the release, just dumping my thoughts here so we can get started. Some things need support in ufl/ffc, their design will depend a bit on other design choices.