Log In Sign Up

Almost complete analytical integration in Galerkin BEM

In this work, semi-analytical formulae for the numerical evaluation of surface integrals occurring in Galerkin boundary element methods (BEM) in 3D are derived. The integrals appear as the entries of BEM matrices and are formed over pairs of surface triangles. Since the integrands become singular if the triangles have non-empty intersection, the transformation presented by Sauter and Schwab is used to remove the singularities. It is shown that the resulting integrals admit analytical formulae if the triangles are identical or share a common edge. Moreover, the four-dimensional integrals are reduced to one- or two-dimensional integrals for triangle pairs with common vertices or disjoint triangles respectively. The efficiency and accuracy of the formulae is demonstrated in numerical experiments.


New error estimates of Lagrange-Galerkin methods for the advection equation

We study in this paper new developments of the Lagrange-Galerkin method ...

On Overcoming the Transverse Boundary Error of the SU/PG Scheme for Moving Conductor Problems

Conductor moving in magnetic field is quite common in electrical equipme...

Two continuous (4, 5) pairs of explicit 9-stage Runge-Kutta methods

An 11-dimensional family of embedded (4, 5) pairs of explicit 9-stage Ru...

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 low-rank 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 four-dimensional, 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 finite-part 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 one-dimensional integral for triangles with a common vertex.

2 Preliminaries

We consider the numerical solution of the Laplace problem


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


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 non-empty and open elements which satisfies:

  1. is a triangulation of , i.e.

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

  2. 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 non-standard 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 piece-wise constant and piece-wise 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


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


where is the kernel function and are trial and test functions respectively. Let be one of the summands for the non-trivial 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

  1. the whole element,

  2. exactly one edge,

  3. exactly one point,

  4. 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


The singularity of the integrand is now located at .

Figure 1: The domain of integration is split into six triangles .

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


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


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.

Figure 2: Visualisation of the integrand of (6) in the complex plane.

The three terms in the integrand are of the form

The anti-derivative is given by


see [Prudnikov1981, Section] 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


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.

Figure 3: Visualisation of the integrand of (8).

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



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,


The first two terms yield

The third term depends on the discriminant of the denominator, which is non-negative due to

For we have

and for

We resubstitute

and set

Finally, we obtain


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 by

We conclude that can be expressed in closed form as


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


Because this applies to the other integrals as well, we only list the parameters , and in Table 1.

Table 1: Values for in (12) to compute .

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 Far-field

Although the far-field 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 two-dimensional 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


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