Use new flux BC to create a nonlinear pressure field.
Allow for higher order RT and BDM elements.

Working:
2D: RT0 x P0, RT1 x P1, BDM1 x P0, BDM2 x P1
3D: RT0 x P0

Not working:
2D: RT2 x P2 - solver does not converge.
3D: RT1 x P1 - fluxes in z direction are not continuous. Rest looks good.
3D: BDM1 x P0 - no interpolation provided by dune-localfunctions.