Poisson equation with pure Neumann boundary conditions

Dear Mathematica users,

I would like to numerically solve a, as the title says, Poisson equation with pure Neumann boundary conditions

$ \nabla^2(\psi)=f$
$ \nabla(\psi)\cdot \text{n}=g$

Is it possible?

For an example I will use the demo from FEniCS project.

$ f=10\text{Exp}(-((x – 0.5)^2 + (y – 0.5)^2)/0.02)$
$ g=-\text{Sin}(5x)$

In Mathematica

f = 10*Exp[-(Power[x - 0.5, 2] + Power[y - 0.5, 2])/0.02] g = -Sin[5*x]; Needs["NDSolve`FEM`"] bmesh = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {1, 0.}, {1, 1}, {0., 1}, {0.5, 0.5}},"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,1}}]}]; mesh = ToElementMesh[bmesh, "MaxCellMeasure" -> 0.001];  m[x_, y_] =   NDSolveValue[{Laplacian[u[x, y], {x, y}] - f ==       NeumannValue[g, True] (*,DirichletCondition[u[x,y]==0, x==0.5&&y==0.5]*)}, u, {x, y} \[Element] mesh][x, y]  ddfdx[x_, y_] := Evaluate[Derivative[1, 0][m][x, y]]; ddfdy[x_, y_] := Evaluate[Derivative[0, 1][m][x, y]]; Show[ContourPlot[m[x, y], {x, y} \[Element] mesh, PlotLegends -> Automatic, Contours -> 50], VectorPlot[{ddfdx[x, y], ddfdy[x, y]}, {x, y} \[Element] mesh,    VectorColorFunction -> Hue, VectorScale -> {Small, 0.6, None}]] 

Trying to solve this in Mathematica gives a clear and understandable error

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified; the result may be off by a constant value. 

However, the result is not so clear and understandable.

enter image description here

The answer one would like to get should look something like the following image enter image description here

I tried equation elimination from this post, i.e. using

DirichletCondition[u[x, y] == 0, x == 0.5 && y == 0.5] 

to get the result. Now it looks decent if one doesn’t care about the sink which appears.. Sadly I do care as I’m interested in the gradient of $ \psi$ so such an approach leads me nowhere.

enter image description here

enter image description here

So then the question – is it possible to numerically solve Poisson equation with pure Neumann boundary conditions with Mathematica? Can anyone suggest some steps how to do this?

To add, sadly I am not a mathematician so I lack the ability to implement some routine on my own. Maybe something can be done using the weak formulation as in the example, but before trying to implement (would it actually be possible?) that I would like to learn if there is another way.