Let $\Omega\subseteq\mathbb{R}^n$ be a bounded domain and $u\in C^0(\overline{\Omega})\cap C^2(\Omega)$ be a solution of $$\left\{\begin{matrix}-\Delta u&=&f&&\text{in }\Omega\\ u&=&0&&\text{on }\partial \Omega\end{matrix}\right.\tag{1}$$ Using Gauss's theorem, I want to conclude, that $$\int_\Omega\langle\nabla u,\nabla\varphi\rangle\;d\lambda^n=\int_\Omega f\varphi\;d\lambda^n\;\;\;\text{for all }\varphi\in C_0^1(\Omega)\;,\tag{2}$$ but I absolutely don't get it.
From $(1)$ we obtain $$-\Delta u=f\;\Rightarrow\;-\Delta u\varphi =f\varphi\;\Rightarrow\;-\int_\Omega\Delta u\varphi\;d\lambda^n=\int_\Omega f\varphi\;d\lambda^n\tag{3}$$ for all $\varphi\in C_0^1(\Omega)$. Gauss's theorem yields $$\int_\Omega\Delta u\varphi\;d\lambda^n=\int_{\partial\Omega}\langle\nabla u,\nu\rangle\;do\tag{4}\;,$$ but that doesn't seem to help at all in $(3)$ to obtain $(2)$. So, what's the trick?
Using Green's first identity what you would find is $$\int_{\Omega} \Delta u \phi = \int_{\partial \Omega} \phi \langle\nabla u, \nu \rangle - \int_{\Omega} \langle\nabla u, \nabla \phi\rangle$$ But $\phi$ is $0$ on the boundary so really you just distribute the Laplacian into two gradients, and pick up a minus sign.