Draft: Resolve "flux reconstruction and transport issues in simplex grid"

What does this MR do?

This MR aims to solve #196 by refactoring the Flux reconstruction machinery.

• One of the main failures of the current machinery is that the Pk elements used as RT-skeleton do not match the value in the normal direction of the complete RT element. This is problematic for cubes.
• Another problem is that because it needs the RT-skeleton and -volume elements, the availability of flux reconstruction depends on whether these elements are provided. With the refactor, we are only limited by the complete RT elements implemented in dune-localfunctions.
• The current approach is very sensible to break if there are internal changes in PDELab because of the direct usage of the engines and assemblers.

The idea here is to use standard PDELab utilities to build the flux reconstruction by prescription. That is, on the one hand, accumulate \int_T \sigma_h^\star \cdot r for the internal degrees of freedom of the RT element and \int_F(\sigma_h^\star \cdot n_F) (q\cdot n_F) for the degrees of freedom in the facets where \sigma_h^\star is an RT element and $r/q$ corresponds to the RT element corresponding to its interior/exterior degrees of freedom. On the other hand, compute the residual of the already computed local operator but using RT elements as trial functions. To do so, the value of the RT element corresponds \sigma_h^\star to the gradient of the original trial function \nabla_h w_h, and, on the facets, values of the original trial function w_h correspond to the out-flux of the RT element, \sigma_h^\star \cdot n. Once mass matrix and residual vector are assembled, is just a matter of solving the system Ax=b to get the values of the degrees of freedom corresponding to the reconstructed flux for RT-elements.

• To make these calculations, it is imperative to identify which degree of freedom in the RT element corresponds to which sub entity.
• The mass matrix will have blocks (per entity) that are linear independent to each other. If we can additionally identify these blocks, calculating direct inverses is pretty quick.

Is there something that needs to be double checked?

Fill this in

Can this MR be accepted?

• Implemented ...
  • Identify DOFs (degrees of freedom) with local keys for each entity
  • Create a local operator that makes l2 norm on interior DOFs and l2 norm on the normal direction for DOFs on the intersection
  • Update FV Richards operator to apply flux reconstruction using only RT elements (without RT-skeleton and -volume)
  • Update DG Richards operator to apply flux reconstruction using only RT elements (without RT-skeleton and -volume)
  • Update FV Transport operator to apply flux reconstruction using only RT elements (without RT-skeleton and -volume)
  • Update DG Transport operator to apply flux reconstruction using only RT elements (without RT-skeleton and -volume)
  • (required) Cache local keys in local operators. This is very expensive!
  • (optional) cache evaluation of the RT local basis
  • Use blocked mass matrix to apply direct solvers on each element/facet
  • (optional) Cache mass matrix when the grid view has no changed between updates
• Added/Updated tests:
  • Check that ensures that RT local keys finder make sense
  • ...
• Added/Updated documentation
• Pipeline passing
• Squash option set
• Delete branch option set
• Added entry to CHANGELOG.md

Assignee: If the Squash option is set, check/update the commit message right before merging!

Related issues

Closes #196