This question has already a bounty winning answer: however, since I found it interesting, upvoted it the very first day it was posted, and found unusual difficulties in finding a nice answer, I nevertheless want to share my own researches and thoughts.
The only reference I am aware of which gives a proof of this result, is the textbook by Mikhailov [1]. Precisely, he deals with the regularity problem for the first boundary problem (Dirichlet problem) and the second boundary problem (Neumann problem) for the Poisson equation in ยง2.3, pp. 216-226. Mikhailov proves the result directly: he first solves the problem in the case of homogeneous boundary conditions, proving the following theorem:
Theorem 4 ([1], p. 217). If $f\in H^k(\Omega)$ and $\partial\Omega\in C^{k+2}$ for certain $k\ge 0$, then the generalized solutions $u(x)$ of the first and second boundary-value problems with homogeneous boundary condition for the Poisson equation belong to $H^{k+2}(\Omega)$ and satisfy (in the case of the second boundary problem it is assumed that $\int_\Omega u \mathrm{d}x=0$) the inequality
$$
\Vert u\Vert_{ H^{k+2}(\Omega)}\le C\Vert f\Vert_{ H^{k}(\Omega)}
$$
where the constant $C>0$ does not depend on f.
By using this result, he extends the regularity theorem 4 above to the case of non homogeneous boundary conditions ([2], p. 226) by reducing non homogeneous boundary conditions to homogeneous ones. He explicitly does it only for the Dirichlet problem but the same method nevertheless works for the Neumann problem: let's see this. Consider a solution $u(x)$ of the problem \eqref{eqlocalvar-Neumann} above and a function $\Phi\in H^{k+2}(\Omega)$ whose normal derivative on $\partial\Omega$ is $g\in H^{k+3/2}(\partial\Omega)$, i.e.
$$
\frac{\partial\Phi}{\partial\nu}=g\;\text{ on }\;\partial\Omega\label{1}\tag{1}
$$
Define $w=u-\Phi$: then, for every $v\in H^{1}(\Omega)$,
$$\label{GeneralizedNeumann}\tag{GN}
\begin{split}
\int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x&=\int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\
&= \int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla \Phi \cdot \nabla v \, \mathrm{d}x \\
&=\int_\Omega f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\Omega} \nabla\cdot(v\nabla\Phi) \, \mathrm{d}x + \int_{\Omega} v\Delta\Phi \, \mathrm{d}x \\
&=\int_\Omega \big[\,f + \Delta\Phi\big]v\, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x) - \int_{\partial\Omega} \nu \frac{\partial\Phi}{\partial\nu} \, \mathrm{d}\sigma(x),
\end{split}\label{2}\tag{2}
$$
and substituting \eqref{1} in \eqref{2} we get
$$
\int_{\Omega} \nabla w \cdot \nabla v \, \mathrm{d}x=\int_\Omega f_1 v \, \mathrm{d}x
$$
where $f_1=f+\Delta\Phi\in H^k(\Omega)$. This means that $w$ solves the Neumann problem \eqref{eqlocalvar-Neumann} with homogeneous boundary conditions (i.e. $g\equiv 0$), and by applying theorem 4 we have
$$
w=u-\Phi\in H^{k+2}(\Omega) \iff u=w+\Phi\in H^{k+2}(\Omega)
$$
Notes
- In Mikhailov's notation, $H^0=L^2$: also he does not use the standard notation for traces, for example expressing that $\varphi\in H^{1/2}(\partial G)$ by saying that $\varphi$ is a function in $L^2(\partial G)$ which is the trace of a function $\Phi\in H^1(G)$.
- Regarding the third boundary problem (the so called Robin problem), Mikhailov ([1], footnote *, p. 217) remarks that analysis of the regularity of solutions can be dealt as done in theorem 4 above for the solutions to the first and second boundary value problems, provided certain conditions are assumed.
[1] V. P. Mikhailov (1978), Partial differential equations, Translated from the Russian by P.C. Sinha. Revised from the 1976 Russian ed., Moscow: Mir Publishers, p. 396 MR0601389, Zbl 0388.3500.
Best Answer
As indicated in my comment, it's natural to consider $v = u$ instead, because it leads to a natural usage of the fact about $f$. Following Green's identity, we have
$$\int_{\Omega} (u \Delta u + |\nabla u|^2) \, dV = \int_{\partial \Omega} u \frac{\partial u}{\partial \vec n} d\sigma = 0.$$
Now since the first term in the integral above is nonnegative by the assumption on $f$, we can conclude that $\nabla u = 0$ almost everywhere. One can continue from here.