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=n\int_\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 $u\in 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$.
Best Answer
Dorian, aren't you messing up things a little? Surely you can expand any $L^2$ function in a series of eigenfunctions for the elliptic operator, but please notice that this simple fact already requires quite a detailed theory of elliptic operators on bounded domains. In order to prove the existence of eigenfunctions you must be able to solve the equation $Lu=\lambda u$, and if you want to use your expansion for regularity results, you need to study the properties of the eigenfunctions, their growth etc.
But actually you are right, in a sense. It is indeed possible to prove existence and regularity of solution at the interior of the domain by using essentially the same methods as on the whole space. However, if you want to control the properties of the solution at the boundary, then this requires new tools. As a minimum, in the lucky situation of a smooth boundary, you can reduce to the case of a half space, but no less than that. If you are not convinced, think of the fact that some results cease to be true if you drop the assumption that the boundary is Lipschitz or satisfies some suitable cone condition.
As you suspect, it is also possible to do frequency space analysis much in the same way as on the whole space, but I would not call this easy. There is a beautiful set of notes "Lectures on semiclassical analysis" available on the web, by Evans and Zworsky, see Theorem 3.17 there (they prove interior Schauder estimates using Paley-Littlewood, apparently following a suggestion of H.Smith). I repeat: this is interior regularity, the behaviour at the boundary is substantially more difficult.