I'm afraid that what I'm about to say is not quite you what to hear: your Neumann PDE does not have a solution for arbitrary choices of $g$.
To see why, let's compute the integral of $g$ over the surface $\Gamma$:
$$ \int_\Gamma g \ dS = \int_\Gamma \frac{\partial u}{\partial n} \ dS=\int_\Omega \nabla^2 u \ dV =\int_\Omega 0\ dV = 0.$$
Therefore, the average value of $g$ on the boundary surface $S$ must be zero. If $g$ fails to satisfy this condition, your PDE cannot be solved.
For a similar reason, your definition of the Green's function must be modified. Indeed, if we integrate $\frac{\partial G}{\partial n}$ over the surface $\Gamma$, we find that
$$ \int_\Gamma \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}dS(\vec r) = \int_\Omega \nabla_{\vec r}^2 G(\vec r,\vec r') \ dV(\vec r) = \int_V \delta(\vec r-\vec r') \ dV(\vec r) = 1.$$
So it is inconsistent to set $\frac{\partial G}{\partial n}$ equal to zero on the boundary $\Gamma$.
Fortunately, all is not lost! We can redefine the Green's function $G$ so that it satisfies
$$ \begin{cases} \nabla_{\vec r}^2 G(\vec r,\vec r') = \delta(\vec r - \vec r') &{\rm \ \ on \ \ } \Omega ,
\\
\frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}=\frac 1 A & {\rm \ \ on \ \ } \Gamma, \end{cases}
$$
where $A = \int_\Gamma dS$ is the area of the boundary surface $\Gamma$.
Now, Green's identity states that
\begin{multline}\int_\Omega \left( u(\vec r) \nabla_{\vec r}^2 G(\vec r , \vec r') - G(\vec r,\vec r') \nabla_{\vec r} u(\vec r) \right) \ dV(\vec r) \\ = \int_\Gamma \left( u(\vec r) \frac{\partial G(\vec r,\vec r')}{\partial n_{\vec r}} - G(\vec r,\vec r') \frac{\partial u(\vec r)}{\partial n_{\vec r}}\right) dS(\vec r).\end{multline}
Plugging in my new definition of $G(\vec r, \vec r')$, we see that any $u(\vec r)$ satisfying your PDE must satisfy
we can immediately see that any solution to the PDE must satisfy
$$ u(\vec r') = - \int_\Gamma G(\vec r,\vec r')g(\vec r) \ dS(\vec r) + c, \ \ \ \ \ \ (\ast )$$
where the constant $c$ is equal to $\frac 1 A \int_\Gamma u(\vec r) dS(\vec r)$, the average value of $u$ on $\Gamma$.
Having redefined the Green's function, I'll give you an explicit expression in the case where $\Omega$ is a two-dimensional circular disk of radius $1$. Here it is:
$$ G(\vec r, \vec r') = \tfrac 1 {2\pi} \ln | \vec r -\vec r' | + \tfrac 1 {2\pi} \ln |\vec r - \vec r''|,$$
where $\vec r'' = \vec r'/|\vec r'|^2$ is the image of $r'$ under an inversion about the unit circle. It is similar to the Dirichlet Green's function, expect that we have a plus sign in front of the "image" term instead of a minus sign.
I believe that if you plug this Green's function into $(\ast )$ (with the constant $c$ chosen arbitrarily), you do indeed get a solution to your PDE, provided that $g$ obeys the consistency condition $\int_\Gamma g \ dS = 0$. This is stated in these lecture notes from my university, and also in Riley, Hobson and Bence, though I would love to see a rigorous proof! I wonder if anyone knows a good reference?
Best Answer
Give two functions $v_{1}= u + \Vert f \Vert_{\infty} w$, $v_{2} = -u + \Vert f \Vert_{\infty} w$. Now, calculate $Lv_1, Lv_2$, and then use maximum principle to obtain the answer .....
Good Luck