Discontinuity of double-layer potentials I'm currently reading about solutions to boundary-value problems for Laplace's equation, and I'm a bit confused with regards to the discontinuity properties of double-layer potentials. So the text says that for a bounded domain $\Omega \subset \mathbb{R}^n$ we can look for a solution to the Dirichlet problem
\begin{align*}
\Delta u &= 0 ~~~~~x \in \Omega\\
u &=f ~~~~~x \in \partial \Omega
\end{align*}
in the form of a double-layer potential
$$u(x)=\int_{\partial \Omega} h(y) \frac{\partial \Phi}{\partial \nu_y}(x-y) \, dS_y$$
where $h \in \mathcal{C}(\partial \Omega)$.
The double-layer potential is discontinuous so that as $x$ approaches the boundary $\partial \Omega$
\begin{align*}
u_-&=\frac{h}{2} + u \\
u_+&=-\frac{h}{2} + u \\
u_- - u_+ &=h
\end{align*}
where $u_-$ is the limit as we go to the boundary from inside, and $u_+$ is the limit as we go to the boundary from the outside. 
I'm really confused. Shouldn't a harmonic function be continuous on $\overline{\Omega}\,$? If we use this method to solve electrostatics problems then the electrostatic potential will be discontinuous across the surface. However, this should not be the case.
 A: For each fixed $y\in \partial \Omega$ the function $v_{y,\nu}(x): = \frac{\partial\Phi}{\partial\nu}\Phi(x-y)$ is harmonic in $\mathbb R^n\setminus \{y\}$. This is the potential function of a unit dipole placed at $y$ and oriented along $\nu$. As such, it has a fairly complicated singularity at $y$, approaching $+\infty$, $-\infty$, and everything in between as $x\to y$. You may want to write down an explicit formula for $v_{y,\nu}(x)$ to check this.
Next, we form the integral $u(x)=\int_{\partial\Omega} h(y)\,v_y(x)\,dy$ which puts together the contributions of dipoles distributed over the surface $\partial \Omega$. The function $u$ is harmonic in $\mathbb R^n\setminus \partial\Omega$, which is an open set with at least two connected components, one of which is $\Omega$. Integration against the continuous function $h$ tones down the singularities at each $y\in\partial \Omega$, but they do not go away entirely: their residual effect is discontinuity at every $y$ such that $h(y)\ne 0$. 
There is no contradiction here, either mathematical or physical. Harmonic functions need not be continuous on the closure of their domain: this is exactly what allows us to solve boundary value problems with discontinuous boundary data. As for physical interpretation, note that we consider a surface of dipoles, not a uniformly charged surface (the latter would be described by a single layer potential).
