Need good reference or a proof on regularity of solution to Neumann problem Let $\Omega\subset \mathbb{R}^d $ be a bounded open subset  ($d\in \mathbb{N}$) and denote  $\partial\Omega$  its boundary which we assume to be Lipschitz.  The classical inhomogeneous Neumann problem for Laplace operator  associate to data   $f:\Omega\to\mathbb{R}$ and $g: \partial\Omega \to \mathbb{R}$ (measurable functions) consists on finding a function  $u:\Omega\to \mathbb{R}$ satisfying
\begin{equation}\label{eqlocal-Neumann}\tag{$N_1$}
-\Delta u = f  \quad\text{in}~~~ \Omega \quad\quad\text{ and } \quad\quad \frac{\partial u}{\partial \nu}=  g ~~~ \text{on}~~~ \partial \Omega.
\end{equation}

In the standard setting one usually choose $ f$ in $L^2(\Omega)$ or  in the dual space of $H^1(\Omega)$ and $g$ can be choose in the trace spaces of $H^1(\Omega)$ denote by $H^{\frac{1}{2}}(\partial\Omega)$ or in its dual $H^{-\frac{1}{2}}(\partial\Omega)$. 
Assume  $f\in L^{2}(\Omega)$ and $g \in H^{3/2}(\Omega)$ then we have the following the Green-Gauss formula
$$\label{eqgreen-Gauss}
\int_{\Omega} (-\Delta) u v \, \mathrm{d}x  =  \int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x- \int_{\partial \Omega} \gamma_{1} u \gamma_{0}v \, \mathrm{d}\sigma(x), \quad u\in H^{2}(\Omega) ~\hbox{and}~v\in H^{1}(\Omega).
$$
Henceforth, on $\partial \Omega$ we merely write $\gamma_0 v= v$ and  $ \displaystyle\gamma_1 v=\displaystyle  \frac{\partial v}{\partial \nu} $. 
Clearly from this green gauss  formula, if $u\in H^{2}(\Omega)$ and solves \eqref{eqlocal-Neumann}  then $u $ satisfies the variational problem
$$\label{eqlocalvar-Neumann}\tag{$V_1$}
   \int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x=  \int_{\partial \Omega} f v \, \mathrm{d}x + \int_{\partial \Omega}gv \, \mathrm{d}\sigma(x), \qquad \hbox{for all } ~~v\in H^{1}(\Omega).
$$
In particular if we put $v=1$ the above formulation becomes
$$\label{eqlocalcompatible-Neumann}\tag{$C_1$}
 \int_{\Omega}f\mathrm{d}x+  \int_{\partial \Omega}g\mathrm{d}\sigma(x)=0
$$
which is the is the compatibility condition. 
Viceversa we have the following global regularity result. 
Theorem
Assume $\Omega\subset \mathbb{R}^{d}$ is bounded open with $C^2$-boundary. If a  function $ u \in H^{1}(\Omega)$ is  solution to \eqref{eqlocal-Neumann} with $f\in L^{2}(\Omega)$ and $g \in H^{3/2}(\partial\Omega)$ then it belongs to $ H^{2}(\Omega).$

Question: In which book or recommendable reference can I find the proof of the above theorem`?
   I know that the of this Theorem for the corresponding Dirichlet problem has been done in the book by Brezis or by Evans. Patently, both references avoid the inhomogenous Neumann problem.

Remark
Moreover, observe that  if $u$ solves  \eqref{eqlocal-Neumann} or \eqref{eqlocalvar-Neumann}  so does $\tilde{u} = u+c$ for every $c\in\mathbb{R} $ (that is invariant under additive constant).
 A: 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.
A: A standard reference is the book of Grisvard Elliptic Problems in Nonsmooth Domains Grisvard. Chapter 2 is about regularity of elliptic equations in smooth domains. The theorem you want is Theorem 2.2.2.5. As for proofs, Grisvard mainly deals with the Dirichlet case and then in a short remark explains how to modify it for Neumann. It might not be what you want.
If you want all the proofs, take a look at Theorems 130- 136 in these lectures notes. lecture notes
