Solution (source code)

= Solution

Use the <Dirichlet Green function> with the sign convention
$$
\boxed{\Delta_{\mathbf r}G(\mathbf r,\mathbf r')=\delta(\mathbf r-\mathbf r')\text{ in }V,\qquad G(\mathbf r,\mathbf r')=0\text{ on }S.}
$$
The <normal derivative> is taken with respect to the integration variable and the outward normal of $V$. <Green's second identity> gives
$$
\int_V\left(u\Delta G-G\Delta u\right)dV
=\oint_S\left(u\partial_nG-G\partial_nu\right)dS.
$$
The <Dirac delta distribution> evaluates the first volume term at $\mathbf r'$. Substituting $\Delta u=\rho$, $u=f$ on $S$ and $G=0$ there yields
$$
u(\mathbf r')-\int_VG(\mathbf r,\mathbf r')\rho(\mathbf r)\,dV
=\oint_Sf(\mathbf r)\partial_nG(\mathbf r,\mathbf r')\,dS.
$$
Hence
$$
\boxed{u(\mathbf r')=\int_V\rho(\mathbf r)G(\mathbf r,\mathbf r')\,dV+\oint_Sf(\mathbf r)\partial_nG(\mathbf r,\mathbf r')\,dS.}
$$
One can justify applying the identity to the singular <Green function> by removing a small sphere around the source and then taking its radius to zero.