Testing inhomogenous neumann boundary conditions, but quite general
since we can put different boundary condtions on each side.

Testing has been done for homogenous case with some fictitious source test.
See: Boundary_conditions.ipynb

two functions are written.

- ddxFaceDivBC
- faceDivBC

Need to be more tested and after that they can be merged in to main