Thoughts After First Implementation

Merging the hypercube refinement and mesh extrusion generation in dune-codegen using the common interface worked reasonably well. The current blockstructured tests for the Poisson problem as well as the 1D linear transport problem for the mesh extrusion could be reproduced. I didn't test the Stokes problem, since this requires additional thought for the datalayout. Also, the vectorization of the blockstructured code does not work, since the common version uses indirection for the vector accessing. Besides these issues, the places where I had to write specialization for each case are minimal. Most of those are within the geometry implementation, which was expected. Addionally, the computation of the dof indices had to be specialized in small parts. More on these divergences is written below.

Divergences in the Implementation

Obviously the geometry implementation of the hypercube refinement and the mesh extrusion differ from each other. In addition to the geometry_dofs, the computations of the reference normals and the facet jacobian (transformation from dim - 1 reference element into real reference element) are separated. Maybe this could be unified by implementing those functions as part of a reference element class for each different type of elements, s.t. only the geometry_dofs need to be implemented separatly.
For the macro element interior facet integrals the index of the neighbouring cell is explicitly assumed to be below wrt to the last index (extrusion index). This obviously won't work in the hypercube refinement case. To fix this the decision, which kernel to generate, needs to be reworked.
There are still additional boundary predicates needed for the boundary integrals in the hypercube case. This is mostly due to the incorrect boundary sets.
Since there is no interface for the datalayout, I explicitly assume the datalayout in each case.
The order of indices gets reversed to compute the linear indices in the hypercube case, see below for the reasoning.
Integrals for the top and bottom boundary part are added in the extrusion case.
The generation of macro element inner face integrals is only done for the extrusion case, so there might be some issues in the hypercube case.

Additions to the Interface

The boundary and interior functions were implemented. These implentations are strictly confined to the isl set and therefore not that usefull, the below.
The simplex entity set is generalized s.t. it can start from a value other than 0. This helped for the definition of the interior set.

Changes / Improvements on the Interface

The access to some sort of 'boundary' or 'interior' entity set is quite helpfull for the code generation. But these boundary or interior sets are not the boundary or interior of an ISL set, since these contain entities that are not needed for the code generation or exclude entities that are. For example consider a 5x5 hypercube refinement with the following entity set for vertical faces
```
{[i, j]: 0<=i<6 and 0<=j<5}
```
The interior of the set in the strict sense would be
```
{[i, j] : 1<=i<5 and 1<=j<4}
```
This excludes the vertical faces i,j with 1<=i<5 and j=0,4, which are needed to generate the loop over all interior faces. Thus the documentation of these functions needs to be changed and the naming might also be adjusted. Additionally, for me it is currently not clear how to compute the correct boundary and interior sets.
The order of indices with respect to the geometric direction is not clear. During the computation of geometry_dofs in the hypercube refinement case the index of a cell entity is translated into a coordinate in the macro element. This needs to relate each index to a direction, which I assumed to be just the order of the indices, i.e. first index is x1, second index is x2 and so on. But the linearization for tensor product sets uses another order of the indices, in that case the fastest changing index is the last index, i.e. the reversed order from before. This results in diverged geometric coordinates and linearized cell indices. Maybe this could be fixed through a datalayout description.
This is more of an implementation note, but the geometry_dofs implementation for a structured topology needs access to geometry_dofs of base topology, if more complicated constructs like the extrusion of a blockstructured grid should be considered. In my implementation of this I've added a base_geometry member to the MeshGeometry class and used a free standing function to get the geometry implementation for the base topology. Maybe this could be done in a smarter way.
The hypercube refinement doesn't use the base topology at all, this means that a blockstructured refinement of an extruded grid is not possible. Although, to be honest, I don't know how that should look like.

Changes / Improvements to the dune-codegen Implementation

The entity set should be determined outside of the visitor, most likely within kernel_generation.py. This would allow to use the entity sets to set the local id of the entities for the next kernel, whereas currently this is set explicitly. This could also help with identifing neighbour cells for macro element interior facets.
The refinement should be passed as part of the MeshGeometry and not by some metadata. I assume that using the ufl mesh in the generation would entail additional changes to the generation pipeline.

wence- / mesh-structure