1 Introduction
Whenever unbounded domains appear in the modelling of physical problems, boundary element methods (BEM) present a particularly effective tool for the numerical simulation. Instead of discretising the underlying boundary value problem directly, BEM operate on the corresponding integral equations posed on the boundary. Hence, shape functions are defined on the surface and not in the volume, which results in less degrees of freedoms overall. BEM are therefore applied in numerous fields of science, ranging from computational elasticity to electromagnetic and acoustic scattering.
However, the price one pays is the occurrence of dense matrices that are expensive to calculate. A Galerkin approximation of the Laplace equation with trial and test functions and requires the calculation of integrals over surface triangles and ,
which do not vanish. In addition, the kernel function is singular at , so standard quadrature rules for the approximation of may perform poorly. Whereas the problem of fully populated matrices can be solved by hierarchical lowrank approximation [Bebendorf2008], several algorithms for the efficient calculation of are available in the literature. The general strategy is to use coordinate transformations that render the integral suitable for numerical integration. Approaches based on polar or Duffy coordinates [Duffy1982, MousaviSukumar2010] yield regular integrals for a wide selection of kernel functions [SchwabWendland1992], which can be approximated by product quadrature rules [Sauter1998]. However, since the integrals are essentially fourdimensional, numerical quadrature is expensive. In order to reduce the computational effort, analytical integration can be carried out for specific kernels to obtain lower dimensional integrals [RjasanowSteinbach2007, Taylor2003]. If the integral is only defined in weak sense as a finitepart or Hadamard integral, then integration by parts is often the appropriate solution [DautrayLions1990].
The main contribution of this article is the derivation of analytical formulae for the complete integration of for singular cases. Our approach is based on the regularisation method by S. Erichsen and S. A. Sauter [Sauter1998] which removes the singularities by applying a variant of the Duffy transformation. We show that the resulting representation admits closed formulae of for identical triangles as well as triangles with a common edge and reduce it to a onedimensional integral for triangles with a common vertex.
2 Preliminaries
We consider the numerical solution of the Laplace problem
(1)  
in a domain with bounded Lipschitz boundary . If is unbounded, we assume the radiation condition
The representation formula expresses the solution in terms of its boundary values only,
Here, is the unit normal to pointing outwards and is the fundamental solution of the Laplace operator,
The boundary value problem is hence reduced to the problem of finding the unknown Neumann trace . To this end, we take the traces in the representation formula and insert the Dirichlet condition to obtain the boundary integral equation
(2) 
where the layer potentials are defined by
Neumann or mixed boundary conditions can be treated similarly. We refer to [McLean2000] for more details.
For the numerical solution of (2) with BEM, we discretise the boundary with finite elements.
Definition 1 (Mesh).
A mesh (or simply ) is a finite collection of nonempty and open elements which satisfies:

is a triangulation of , i.e.
The intersection of two distinct elements is either empty or consists of a common vertex or edge.

Each in is a flat triangle with vertices and . The reference mapping
parametrises by the reference triangle
We denote by
the Jacobian and the Gram determinant of respectively and assume .
Remark 1.
Certainly, the particular choice of the reference element is not important for our approach. The reason why we use a nonstandard nonetheless lies in the fact that the presentation in Section 3 becomes simpler and it is in accordance with the literature referenced there.
On the triangular mesh, we define piecewise constant and piecewise linear ansatz functions.
Definition 2 (Boundary element spaces).
For , we denote by
the set of reference functions and by
the local boundary element space on . We define the global space by gluing the local spaces together, i.e.
For , we moreover require that the functions are continuous.
We choose the Lagrangian basis
where denotes the set of vertices in . Then, the ansatz
for the approximate boundary data leads to the Galerkin approximation of (2): Find such that
(3) 
where the matrices and are given by
with and . The brackets symbolise the usual inner product
The boundary integral equation is now reduced to a system of linear equations, which can be solved efficiently with direct or iterative methods.
3 Integral Regularisation
The entries of an are of the form
(4) 
where is the kernel function and are trial and test functions respectively. Let be one of the summands for the nontrivial case and . We transform back to the reference element ,
and abbreviate the integrand by . Since the kernel function is singular at , the integral needs to be regularised. We distinguish between four different cases: the intersection may consist either of

the whole element,

exactly one edge,

exactly one point,

be empty.
In the following, we summarise the regularisation introduced in [Sauter1998]. For the most part, we adhere to the version of [SauterSchwab2011, Chapter 5].
3.1 Identical elements
For identical elements , we substitute
such that
with
The singularity of the integrand is now located at .
As shown in Figure 1, we decompose into six triangles with
and for . Thus, we obtain
In the last step, we parametrise by using the Duffy transformation
with Jacobian and conclude that has the representation
In Section 4, we will see that the integrand is smooth since the Jacobian cancels out the singularity of .
3.2 Common edge
If the two triangles intersect at exactly one edge, we can proceed similarly to the first case. Let and be chosen in such a way that the common edge is parametrised by
The kernel function is singular at this edge, i.e. at with and . The regularisation is carried out via the mappings
and reads
3.3 Common vertex
Let the origin in the reference domain be mapped to the common vertex, i.e.
By virtue of the mappings
we obtain
In summary, the regularisation yields integral representations with smooth integrands on the unit cube. In this form, the integral can be approximated efficiently by quadrature rules and the quadrature error decays exponentially with the quadrature order.
4 Calculation of integrals
Instead of applying quadrature rules directly, we calculate parts of the regularised integrals analytically.
4.1 Single layer potential
With the discretisation provided in Section 2, the entries of the single layer potential are of the form
(5) 
We proceed like in Section 3 and begin with the case of identical elements.
4.1.1 Identical Elements
Let and be the edges of with starting point . Then, the triangle is parametrised by
and the regularisation of (5) reads
where the area is for . Hence, we obtain
(6) 
Figure 2 depicts the complex continuation of the integrand for concrete values of and . It is smooth on the real axis, since the edges are linearly independent, but has poles and branch cuts in the complex domain.
The three terms in the integrand are of the form
The antiderivative is given by
(7) 
see [Prudnikov1981, Section 1.2.52.8] and [Gradstein2015, Section 2.261]. Thus, the integral reduces to
where is with the parameters
4.1.2 Common edge
Let the reference mappings be given by
such that is the common edge of and starting from . Then, the integral (5) reduces to
(8)  
In comparison to the previous case, the integrand is not necessarily smooth in the real domain. When the two triangles lie in the same plane, it has poles as seen in Figure 3. However, they only occur outside of since the triangles do not overlap.
First integral
Let us introduce the variables
We integrate with respect to by using (7) and obtain for the first integral
Integration by parts leads to
(9) 
where
with . We note that coincides with when is replaced by and proceed with integrating . We follow the approach of [RjasanowSteinbach2007, Appendix C.2] and define
such that
If then and the integral simplifies to
Otherwise, we have and the substitution
yields for the indefinite integral
We use a variant of the Weierstraß substitution,
and obtain
The integrand is now a rational function and we abbreviate it by
We decompose it into partial fractions,
where
The first two terms yield
The third term depends on the discriminant of the denominator, which is nonnegative due to
For we have
and for
We resubstitute
and set
Finally, we obtain
(10) 
Note that the value of the integral only depends on the vectors
. Since it is of importance for the other cases as well, we abbreviate it byWe conclude that can be expressed in closed form as
(11) 
with , , .
Remaining integrals
For the remaining integrals, we integrate with respect to the fourth variable firstly. With
we have
Integration by parts yields an expression almost identical to (9),
where and are given by
Thus, can be computed analogously to (11) by
(12) 
Because this applies to the other integrals as well, we only list the parameters , and in Table 1.
4.1.3 Common vertex
We consider the configuration
with common vertex . Then, integration with respect to results in
We only consider , since is obtained by swapping and . Similar to the previous section, we introduce the variables
and integrate with respect to to obtain
We insert Formula (12) for the inner integral, which yields
The first integral can be written in the form of from Section 4.1.2, i.e.
with , , , and its value is hence given by (11). Because it is not possible to integrate the remaining integral analytically, we approximate it numerically with a quadrature rule
with weights and nodes .
4.1.4 Farfield
Although the farfield does not constitute a singular case, the analytical formulae are still applicable. Let the elements be given by
and set . Analogously to the previous cases, we pull the region of integration back to by
leading to
Of the four iterated integrals, we compute two analytically and two by numerical quadrature, e.g.
where the twodimensional integral is calculated analytically using (11) with
4.2 Double layer potential
For the double layer potential , we need to compute integrals of the form
(13) 
where is the outer unit normal vector at and , i.e.
with coefficients .
4.2.1 Identical elements
For identical elements , we simply have due to