differential equations - Elliptic regularity for the Neumann problem

In the case that you mentioned, we want to avoid this cut-off/difference quotients approach, since it could be hard to prove that $partial_{x_i} (xi partial_{x_j} u)$ is a valid test function.
In general, when working with regularity theory, another standard approach is to use an 'approximated problem'. However, the kind of the approximated problem, of course, depends on the PDE.
For the Neumann like problem I suggest the following approximation:

First observe that since $int_Omega f = 0$, we can easily construct a sequence ${f_n} subset C^infty(Omega)$ such that $f_n to f$ in $L^2(Omega)$ and $int_Omega f_n=0, forall n in mathbb{N}$.
Then we consider $u_n in C^infty(Omega)$ such that

$(*) -Delta u_n + frac{1}{n} u_n = f_n, mbox{ in } Omega $ and $dfrac{partial u_n}{partial nu}=0 , mbox{on }partial Omega.$

The sequence ${ u_n }$ can be obtained by the use of Theorem 2.2.2.5, p.91 and Theorem 2.5.1.1, p. 121 of Grisvard's Book p. 91.
In fact you just need to use a boostrap argument to $-Delta u_n = -frac{1}{n} u_n + f_n$.

Notice that $int_Omega u_n=nint_Omega f_n = 0$.

Now, you use $u_n$ as your test functions and obtain the following estimate:

$(**) |nabla u_n|_{L^2}^2 leq |f_n|_{L^2}^2$, $forall n in mathbb{N}$

Now you use $-Delta u_n $, as a test function to your PDE ( observe that $-Delta u_n$ is a valid test function, anyway we don't need to worry about it since the approximated equation holds everywhere).

After integrating by parts, by using $(**)$ with some standard manipulations with your boundary terms you end up with

$|D^2 u_n |_{L^2}^2 leq C(partial Omega)|f_n|_{L^2}^2$, $forall n in mathbb{N}$.
(For instance, see Grisvard's book p.132-138, in particular eq. 3.1.1.5)

The key point for the above estimation is to control the boundary elements in terms of the mean curvature of $partial Omega$.

Now, since $int_Omega u_n =0$ we conclude that
$|u_n |_{H^2}^2 leq C(partial Omega)|f_n |_{L^2}^2 $

so that

$|u_n |_{H^2}^2 leq C(partial Omega) |f|^2$, the $L^2$ norm of $f$.

In this way we obtain $uin H^2$ such that $u_n to u$ weakly in $H^2$ and strongly in $H^1$.

Observe that the latter convergence is sufficient to handle the term $dfrac{1}{n}u_n$.
Then, we can pass to the limit in equation (*) so that $u$ is a strong solution of

$-Delta u = f$ in $Omega$

$dfrac{partial u}{partial nu}=0$ on $partial Omega$

with
$|u|_{H^2} leq C(partial Omega) |f|$, the $L^2$ norm of $f$.