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.