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 RTskeleton 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 RTskeleton 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
dunelocalfunctions
.  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 outflux 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 RTelements.
 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 RTskeleton and volume) 
Update DG Richards operator to apply flux reconstruction using only RT elements (without RTskeleton and volume) 
Update FV Transport operator to apply flux reconstruction using only RT elements (without RTskeleton and volume) 
Update DG Transport operator to apply flux reconstruction using only RT elements (without RTskeleton 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