I'm trying to combine two ways of looking at the Laplacian $\Delta$ on $\mathbb R^n$ (and on other domains).

Firstly, it is well known that this operator is essentially self-adjoint on $C_c^\infty(\mathbb R^n)$. Secondly, I know that for $f,g \in C^\infty_c(\mathbb R^n)$ it holds that $\langle u, \Delta u \rangle = \langle \nabla u, \nabla v \rangle $, where the inner-product is the $L^2$ one.

Now the domain $\mathcal D (\Delta) \subset L^2$ on which the Laplacian is self-adjoint must contain $C_c^\infty (\mathbb R^n)$, but I can't seem to find much more information on this. By the previous identity I feel like $\Delta u \in L^2$ should imply that $u \in H^1(\mathbb R^n)$, but I can't seem to find a way to prove this. I know it is possible to prove that $\mathcal D(\Delta) \subset H^2(\mathbb R^n$) via elliptic regularity etc, but it seems to me like there "should" be an elementary argument for showing $\mathcal D(\Delta) \subset H^1$ (preferably via the non-Fourier-transform characterisation of weak derivatives) but I can't seem to find it anywhere.

Edit: To clarify my question: I am looking for a simple proof of $\Delta u = f \in L^2$ implies $u \in H^1$ where $\Delta u$ is defined as the self-adjoint extension of the Laplacian $\Delta$ on $C^\infty_c$, preferrably one which also works on domains other than $\mathbb R^n$.