Matrix inversion at code generation time does only work to
a very limited extent (up to n=4). We can instead assemble
the tensor in C++ and invert it there (e.g. using Dune::FieldMatrix)

Calculate the evaluation and jacobians of geometry mappings of cells
using sum factorization. By avoiding all geometry mapping of intersections
we can avoid the face-twist problem and just make the grid consistens (not
yet done).

TODO:
- Make grid consistent
- Vectorization
- 3D

This will be needed to replace nonlinearities in Richards code
with a cubic spline evaluation without bloating dune-perftool with
the implementation details of such procedure.

Also provides a test, that defines a custom square function.

In ULF the jump is defined as jump(v) = v('+') - v('-'). Accordingly
the outer unit normal has positve restriction and ('+') refers to self
und ('-') to neighbor. Unfortunatly you can also find '+' and '-'
swapped in the literature.

This reverts commit 5081753a.

Needs monkey patching of the preprocessing to not
transfer to reference element...

It makes the code generation workflow more complicated without
contributing anything substantial

Instationary stuff not fixed so far.

The visitor state regarding indices needs to be correctly reset
for child iteration in these cases.