Matrix Free Solver for DTR with Neumann BCs

[Luca Greco]

Hi! I am trying to develop a generic matrix free solver for DTR problems, so far everything seems to be working fine except for dishomogeneus Neumann BCs. Attached is an example equation I am trying to solve. 

The issue I am facing is related to the solution. In order to assess to correctness of the implementation I am confronting the solution of the same problem with a matrix based solver, inspecting both of the solutions in paraview highlights an error in the upper bound of the solution. Moreover, is seems that the solution is qualitatively correct, two screenshots are provided.

The code used I also attached, it is just an adaptation of step-37 with some ideas from step-59 for calculation of boundary integrals with FEFaceEvaluation.

Thanks in advance!