Feature/improve box facetcoupling

This "tpfa-like" flux approximation on interior Neumann boundaries significantly reduces the condition number of the linear system and removes the problem of not being able to reduce the residual up to machine precision. It now runs stable and converges as expected.

This is wip because it has to wait for !1431 (merged). The diff will be much better to review after that is merged.