Passing from Classical Formulation to Weak Forulation on PDEs. Integration by parts in $n$ dimensions? I am reading about the Variational Method for solving PDEs and ODEs. The process to pass from a classical solution to a weak solution is pretty clear when we are dealing with ODEs, by using Integration by parts. For example, for the problem: 
$$ \left\{\begin{matrix}
-u(t)''+u(t)=f(t) & t \in [a,b] \\ 
 u(a)=u(b)=0 \ \  
\end{matrix}\right.
\quad \quad \quad \quad \boldsymbol{(^{*})} 
$$
We start by picking a classical solution $u$ and multiply on both sides of the equation by a smooth compactly supported function $\varphi$ we obtain:
$$ -u'(b)\varphi(b)+u'(a)\varphi(a) + \int_a^{b}u'\varphi' + \int_a^{b}u\varphi = \int_a^{b}f\varphi$$
i.e.:
$$\int_a^{b}u'\varphi' + \int_a^{b}u\varphi = \int_a^{b}f\varphi$$
As both $\varphi(b)=\varphi(a)=0$. This is the weak (variational) formulation, obtained from the classical formulation.

Now I want to do the same for a PDE

Let's take for example the Laplacian: $$\left\{\begin{matrix}
- \Delta u + u=f & in \ \Omega  \quad \quad \quad \quad \quad \quad \boldsymbol{(^{*})}\\
\ \ \ \quad  \quad \  u=0 & \ in \ \partial \Omega  \quad  \quad  \quad  \quad  \quad   \quad  \quad
\end{matrix}\right.$$
If we repeat the process by using $u \in C^2(\overline{\Omega})$, what my book says is that we conclude: $$\int_{\Omega} \sum_{i=1}^n u_{x_i}v_{x_i} + \int_{\Omega}uv= \int_{\Omega}fv \quad \forall v \in C^1_c(\Omega) \quad  \boldsymbol{(^{**})}$$
But it does not justify why. Can someone clarify me how to do this? I understand, due to the fact that for the partial derivatives $(uv)_{x_i}=u_{x_i}v+v_{x_i}u$, we can also do some kind of integration by parts like this (Let's suppose $n=2$ and we take the first term): 
$$\int_{\Omega}  u_{x_1x_1}v+ u_{x_2x_2}v \ dx = \int_{\alpha_1}\int_{\alpha_2}u_{x_1x_1}v+ u_{x_2x_2}v \ dx_1 dx_2=$$ $$ \int_{\alpha_1}( u_{x_1}v_{x_1}  ]_{\alpha2} - \int_{\alpha_2} u_{x_1}v_{x_1}\ dx_1)  \  dx_2  + \int_{\alpha_2} u_{x_2x_2}v\ dx_1 dx_2= ...$$
Then, if I were able to (I mean rigorously, I know it must be due to the compact support) prove the term $ u_{x_1}v_{x_1}  ]_{\alpha2} =0$, then I would get the desired equation (Doing the same thing for the term $\int_{\alpha_2} u_{x_2x_2}v\ dx_1 dx_2$ ).
So my doubts are: 


*

*$\alpha_1, \alpha_2$ are put there as placeholders, I don't really know how to rigorously afirm we integrate along the curve that bounds the domain $\Omega$... 

*Connected with 1, I don't really know how to confirm $ u_{x_1}v_{x_1}  ]_{\alpha2} =0$. I suppose rigorously afirming 1  and using the fact that $v$ is compactly supported would help.

*Is the general process correct? I haven't been able to find any information on integrating by parts in multiple integrals.
Aditional Information: Functional Analysis, Sobolev Spaces and Partial Differential Equations, by Haim Brezis affirms the integration by parts formula holds for weak derivatives (One dimension) and uses it to solve ODEs. But it does not provide the same formula (The one I try to apply here) for weak partial derivatives, although it does prove the formula to partially derivate the product of two funtions. 
The process being carried to pass from a classical solution to a weak solution is also different. In ODEs, it is carried out by multiplying by $\varphi \in H^1_0(I)$ and using the formula of Integration by Parts for weak derivatives. In PDEs, it is carried out by multiplying by $v \in C^1_c(\Omega)$ , then (Using the process described here?) and then argumenting by density to conclude the result for $H^1_0(\Omega)$.
 A: D1X, this fact has been said here What is precisely the definition of Elliptic Partial Differential Equation? 
We remember that if $\Omega \subset \mathbb{R}^n$ is a limited open with boundary $\partial \Omega$ of class $C^1$, the following are true Green's identity (with $u,v \in C^2(\Omega)\cap C^1(\overline{\Omega})$)
\begin{align*}
\displaystyle \int_\Omega v(x) \Delta u(x) dx + \int_{\Omega} \nabla u(x) \cdot \nabla v(x) dx &= \int_{\partial \Omega} v \partial_\nu u d \sigma 
\end{align*}
for simplicity consider the following problem
\begin{align*}
\displaystyle \begin{cases} -\Delta u =f & \Omega \\ u=0 & \partial \Omega  \end{cases}
\end{align*}
where in general $f \in L^2(\Omega)$, buf $f \in C(\overline{\Omega})$, and $u \in C^2(\Omega) \cap C(\overline{\Omega})$ is a classical solution, the problem above is true pointwise. Define
\begin{align*}
\displaystyle (\nabla u , \nabla v)_{L^2}:= \sum_{j=1}^{n} (\partial_j u , \partial_j v)_{L^2}
\end{align*}
Assuming $u \in C^2(\Omega) \cap C(\overline{\Omega})$ such that $u=0$ on $\partial \Omega$, multiplied by the $L^2$-scalar product Poisson equation by a test function $\varphi \in \mathcal{D}(\Omega)=C_{c}^\infty(\Omega)$ , we get
\begin{align*}
\displaystyle (\varphi, -\Delta u)_{L^2} = (\varphi , f)_{L^2} 
\end{align*}
I apply the Green's identity with $v=\varphi$, i.e. since $\mathrm{supp}(\varphi) \subset \Omega$ we have $\varphi(x)=0$ on $\partial \Omega$ and the integral term on the boundary vanishes, then
\begin{align*}
\displaystyle (\nabla \varphi, \nabla u)_{L^2}=(\varphi , f)_{L^2} , \forall \varphi \in \mathcal{D}(\Omega)
\end{align*}
Since $H_{0}^1(\Omega)$ is the closure of $\mathcal{D}(\Omega)$ in the $H^1$-norm we have that the above is true $\forall v \in H_{0}^1(\Omega)$, since therei is $\lbrace \varphi_k \rbrace \subset \mathcal{D}(\Omega)$ such that $\varphi_k \rightarrow v$ in $H^1(\Omega)$ i.e. $\varphi_k \rightarrow v$ and $\partial_j \varphi_k \rightarrow \partial_j v$ in $L^2(\Omega)$, and here we go from classic formulation to weak:
If $f \in L^2(\Omega)$ is fixed, a function $u \in H_{0}^1(\Omega)$ is called a weak solution for the Dirichlet problem if 
\begin{align*}
\displaystyle (\nabla u, \nabla v)_{L^2}=(v , f)_{L^2} , \forall v \in H_{0}^1(\Omega)
\end{align*}
