# Double layer potential operator at a surface, for coplanar source and observation points

### Background

I'm implementing a Laplace equation solver with the boundary element method for an internal Neumann problem (i.e. bounded domain, with Neumann boundary conditions). The solution is only unique up to a constant, so I'm using the deflation method [e.g. this and this], to get a solution.

I am using a triangular mesh for the object's surface, using piecewise constant expansion and testing functions to obtain the matrix equation. The object can be assumed to be a cuboid (i.e. corners, flat surfaces).

I have a conceptual question on numerically computing the double-layer potential operator, which involves an inner integral over the source triangle (for which I use an analytical technique), and an outer integral over the testing triangle (for which I use numerical quadrature). When the source and test triangles overlap, I am using the Cauchy principal value concept to extract out the residue term which accounts for the strong singularity. My question is about the other, off-diagonal matrix entries, when source and observation points are on different triangles.

### Question

I am observing serious numerical issues, where I can get close to the correct solution, but if I change the order of numerical integration for the outer integral (i.e. number of quadrature points), the solution degrades drastically, even if I increase the order. I am wondering if the reason is related to the following:

The double-layer operator involves the normal derivative of the Green's function, $$\partial G/\partial n$$, where $$G = 1/r$$, and $$r=|\vec{r} - \vec{r}\,'|$$ is the distance between the source $$\vec{r}\,'$$ and observer $$\vec{r}$$. I am confused about whether I can use \begin{align} \frac{\partial G}{\partial n} = \hat{n}\cdot\nabla G\tag{1} \end{align} in the computation. When the source and observer triangles are not the same, $$r > 0$$, so the function G should be differentiable, therefore I should be able to use (1) in my computation. But \begin{align} \nabla G = -(\vec{r} - \vec{r}\,') \frac{1}{r^3} \tag{2}. \end{align} which means that whenever $$\vec{r}\,'$$ and $$\vec{r}$$ are coplanar, $$\hat{n}\cdot\nabla G = 0$$ according to (2). But when the source and observer are on the same triangle, that term is suddenly non-zero. Is this correct?

The above doesn't seem right to me physically. The double-layer potential should represent a layer of tiny dipoles along the surface of the object, and I don't understand how those dipoles can produce an extremely strong field locally, and an exactly zero field everywhere else on that plane. I'm imagining a tiny dipole, which should certainly produce a field along its bisecting plane, as would all the other tiny dipoles on that same plane. Moreover, it's easy to see that those fields due to the tiny dipoles should add up on that flat surface. Either this physical interpretation of the double-layer potential is wrong, or (1) must not be valid right at the surface - which is it?