PDE Laplace equation. Integral representation form and Green function Let $\Omega$ be a domain in $\mathbb{R}^{d}$ and assume that for any $y \in \Omega$ there is a function $h_{y} \in C^{2}(\overline{\Omega})$ such that
 \begin{equation}
 \label{eq8.1}
  \begin{cases}
   \Delta h_{y}(x) = 0 \text{ in } \Omega \\
   h_{y}(x) = E(x,y) \text{ on } \partial \Omega
  \end{cases}
 \end{equation}
Where $E$ denotes the fundamental solution to $\Delta$ in dimension $d$. Under those assumptions we define the Green function $G(x,y) = E(x,y) - h_{y}(x)$
Let $\Omega$ be a bounded domain such that one can define the Green function $G$. Then for any $u \in C^{2}(\overline{\Omega})$ and any $y \in \Omega$, we have
$$ u(y) = - \int_{\Omega} G \Delta u dx - \int_{\partial \Omega} \partial_v G(x,y) u(x) d\sigma(x) $$
Where $v$ is the outer normal of $\partial \Omega$. I don't see how the above integral representation is supposed to directly follow from the definition of the Green function and the fact that if $E$ is a fundamental solution of $\Delta$ we have:
$$
u(y) = - \int_{\Omega} E \Delta u dx + \int_{\partial \Omega} E(x,y) \partial_vu(x) - \partial_vE(x,y) u(x) d\sigma(x)
$$
 A: It follows from the divergence theorem applied to the vector field $h_y\nabla u - \nabla h_y u$. Then, the we get the terms term $\partial_v u h_y = \partial_vE$ and $\partial_v h_y u$ which cancels the term from $\partial_v G u$, so we recover the original integral representation for $u$. 
A: First of all, remember what $E$ is. In $3$ or more dimensions, it is
$$E(x,y) = \frac{1}{|x-y|^{n-2}}$$
If you take the second derivative and sum you obtain $$\Delta E(x,y)  = \sum_{i=1}^{n}\frac{\partial^2 E(x,y)}{\partial x_i^2} = 0 \mbox{ for $x\neq y$}$$
this is the same as saying $\Delta E(x,y) = \delta_x(y)$, where $\delta_x$ function is the function such that $$\delta_x = \delta(y-x)$$  Where $\delta$ is the Dirac Delta Function:
$$\delta(x) = \begin{cases} 
      +\infty & x= 0 \\
      0 & x\neq 0 \\
   \end{cases}$$
This isn't a rigorous function definition, because $\infty$ in just one point is meaningless. However it can be defined rigorously. 
Well, the delta funciton has the nice property that 
$$\int_{-\infty}^{\infty} f(x)\delta(x-a) \ dx = f(a)$$
The intuition is that if you have such function $E$ with the property that its laplacian is a delta function, then 
$$\Delta_x E(x,y) = \delta_x(y-x) \implies \int_{-\infty}^{\infty}f(x)\Delta E(x,y) = \int_{\infty}^{\infty} f(x)\delta_x(y-x) = f(x)$$
See that you have $f(x)$ on one side and our magical function $E$ on the other. We're gonna do that for more dimensions, but rigorously.
Consider the second Green identity:
$$\int_{\Omega} u\Delta w - w\Delta u \ dy = \int_{\partial \Omega}\left(u\frac{\partial w}{\partial n}-w\frac{\partial u}{\partial n}\right)$$
So by doing $u=$ and $w=E$ we have:
$$\int_{\Omega} u\Delta E - E\Delta u \ dy = \int_{\partial \Omega}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right)$$
we cannot integrate $E$ on the entire $\Omega$ because it is not defined when $x=y$, so we need to separate this integral in the region $\Omega_\epsilon$, which is a region formed by $\Omega$ minus a ball around the point of indefinition like this:
$$\Omega_\epsilon = \Omega-B_x(\epsilon)$$
So if we separate $\Omega = \Omega_\epsilon \cup  B_x(\epsilon)$, we know that the integral on $\Omega_\epsilon$ is $0$ because there's no indefinition there, and we hav to deal with the integral on $\Omega_\epsilon$. Note that $\partial\Omega_\epsilon = \partial \Omega \cup \partial B_x(\epsilon)$, so we end up with:
$$\int_{\Omega_\epsilon} u\Delta E - E\Delta u \ dy =\int_{\Omega_\epsilon}-E\Delta u  \ dy \implies \\ \int_{\partial \Omega_\epsilon}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right) = \int_{\partial \Omega}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right)+  \int_{\partial B_x(\epsilon)}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right)$$
First note that 
$$\lim_{\epsilon\to 0}\int_{\Omega_\epsilon}-E\Delta u  \ dy = \int_{\Omega}-E\Delta u  \ dy$$
since it does not depend on $\epsilon$ and by an additional unexplained argument here, $E$ is integrable on $x=y$, so we can write that integral because it is finite.
Now by another argument I'll not make here, 
$$\lim_{\epsilon\to 0} \left(\int_{\partial \Omega}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right)+  \int_{\partial B_x(\epsilon)}\left(u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n}\right)\right) = \\ \int_{\partial \Omega} u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n} dy + u(x)$$
so we end up with
$$\int_{\partial \Omega} u\frac{\partial E}{\partial n}-E\frac{\partial u}{\partial n} dy + u(x) = \int_{\Omega}-E\Delta u  \ dy \tag{1}$$
Note that it only solves half of the problem. We don't know for example how does $\frac{\partial u}{\partial n}$ looks like on $\partial \Omega$. This is mitigated by introducing the corrector function $h_x(y)$ that you talked about. 
The unexplained parts and the corrector function details can be found here
