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.

So far, it had to restore indexing information on the visitor
for the second branch. This can be avoided by index pushdown,
just as we do with the sum of two tensors.

It spares us some hassle with function manglers.

Make anisotropic quadrature order possible for non DG examples.