Not yet implemented:

- Vectorization of quadrature loop
- Anisotropic amount of quadrature points per direction

CI (no -ffast-math) produces a slightly relaxed error.

The check fails in recent CI, but I think it should be correct
anyway. Also, I dont think that error bounds below single precision
accuracy are good for numdiffs.

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.

I also made the order of the terms more consistent and it is now much
easier to compare the skeletons and the boundaries.