Water flux error estimator must include Neumann boundary condition
Description
The water flux error estimator computes an estimate of the conformity of the solution by integrating jumps in the solution and jumps in the gradient over grid faces.
The error estimator contains two terms:

Solution jump
\propto \gamma [[ u ]] / h_F
.At the boundary, the jump is defined by the difference between the interior solution and the Dirichlet boundary condition (zero for Neumann).

Gradient jump
\propto [[ K \nabla u ]]
.At the boundary, the jump is defined by the difference between the interior solution gradient and the Neumann boundary condition (zero for Dirichlet). This is not considered right now.
Proposal
Add the Neumann boundary condition to the "gradient jump" term in the error indicator.
How to test the implementation?
Add a unit test: For a very simple 1x1 grid, a constant initial condition and prescribed boundary conditions, compute the error indicator and verify its value.
Related issues
None