DeepAI

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

09/14/2019

### Semi-analytical calculation of the singular and hypersingular integrals for discrete Helmholtz operators in 2D BEM

Approximate solutions to elliptic partial differential equations with kn...
12/20/2019

### Fast hybrid numerical-asymptotic boundary element methods for high frequency screen and aperture problems based on least-squares collocation

We present a hybrid numerical-asymptotic (HNA) boundary element method (...
01/09/2023

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

We study in this paper new developments of the Lagrange-Galerkin method ...
01/29/2021

### The hybrid dimensional representation of permeability tensor: a reinterpretation of the discrete fracture model and its extension on nonconforming meshes

The discrete fracture model (DFM) has been widely used in the simulation...
09/17/2021

### 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...
02/17/2022

### 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...
07/22/2021

### Efficient Fast Multipole Accelerated Boundary Elements via Recursive Computation of Multipole Expansions of Integrals

In boundary element methods (BEM) in ℝ^3, matrix elements and right hand...

## 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 ,

 I=∫τ∫σ14π|x−y|φ(y)dS(y)ψ(x)dS(x),

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

 −Δu =0 inΩ, (1) u =g onΓ=∂Ω,

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,

 u(x)=∫Γu∗(x,y)∂nu(y)dS(y)−∫Γ∂n(y)u∗(x,y)g(y)dS(y),x∈Ω.

Here, is the unit normal to pointing outwards and is the fundamental solution of the Laplace operator,

 u∗(x,y)=14π|y−x|.

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

 Vt=(12I+K)gonΓ, (2)

where the layer potentials are defined by

 (Vw)(x) =∫Γu∗(x,y)w(y)dS(y), (Kw)(x) =∫Γ∂n(y)u∗(x,y)w(y)dS(y).

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.

 Γh=N⋃n=1¯τn.

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

 χτ:π→τ,χτ(x1,x2)=p1+x1(p2−p1)+x2(p3−p2),

parametrises by the reference triangle

 π={(x1,x2)∣0

We denote by

 Jτ=(p2−p1∣p3−p2)∈R3×2,gτ=√det(J⊤τJτ)

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

 S0h(π)={1},S1h(π)={1−x1,x1−x2,x2}

the set of reference functions and by

 Sph(τ)=span{φ∘χ−1τ:φ∈Sph(π)}

the local boundary element space on . We define the global space by gluing the local spaces together, i.e.

 Sph={φ:Γh→R:φ|τ∈Sph(τ) ∀τ∈Th}.

For , we moreover require that the functions are continuous.

We choose the Lagrangian basis

 φ0n(x)={1ifx∈τn,0else},φ1m(xj)={1ifj=m,0else},

where denotes the set of vertices in . Then, the ansatz

 th =N∑n=1tnφ0n∈S0h, t∈RN, gh =M∑m=1gmφ1m∈S1h, g∈RM,

for the approximate boundary data leads to the Galerkin approximation of (2): Find such that

 Vt=(12M+K)g, (3)

where the matrices and are given by

 M[n,m]=⟨φ0n,φ1m⟩,V[n,i]=⟨φ0n,Vφ0i⟩,K[n,m]=⟨φ0n,Kφ1m⟩,

with and . The brackets symbolise the usual -inner product

 ⟨u,v⟩=∫Γhu(x)v(x)dS(x).

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

 ∫Γh×Γhk(x,y)φ(y)dS(y)ψ(x)dS(x)=∑σ,τ∈Th∫σ∫τk(x,y)φ(y)dS(y)ψ(x)dS(x), (4)

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 ,

 I =∫σ∫τk(x,y)φ(y)dS(y)ψ(x)dS(x) =∫π×πgσgτk(χσ(x),χτ(y))φ(χσ(x))ψ(χτ(y))d(y)d(x),

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

 z=x−y,Z=(x+y)/2

such that

 I=∫Π∫πzq(Z+z/2,Z−z/2)dZdz

with

 Π={z=x−y:x,y∈π},πz=(π−z/2)∩(π+z/2).

The singularity of the integrand is now located at .

As shown in Figure 1, we decompose into six triangles with

 A1=(1001),A2=(0110),A3=(011−1)

and for . Thus, we obtain

 I=6∑i=1∫π∫πAizq(Z+Aiz/2,Z−Aiz/2)dZdz.

In the last step, we parametrise by using the Duffy transformation

 z(η)=(η1η1η2),Ai(η)=Aiz(η)

with Jacobian and conclude that has the representation

 I=∫(0,1)2η16∑i=1 ∫πAi(η)q(Z+Ai(η)/2,Z−Ai(η)/2)dZdη.

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

 χτ(x1,0)=χσ(x1,0),x1∈(0,1).

The kernel function is singular at this edge, i.e. at with and . The regularisation is carried out via the mappings

 I=∫(0,1)4η31η22(q(A1(η))+η35∑i=2q(Ai(η)))dη.

### 3.3 Common vertex

Let the origin in the reference domain be mapped to the common vertex, i.e.

 χτ(0,0)=χσ(0,0).

By virtue of the mappings

we obtain

 I=∫(0,1)4η31η3(q(A1(η))+q(A2(η)))dη.

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

 I=∫σ∫τ14π|y−x|dS(y)dS(x)=gτgσ4π∫π×π1|χτ(y)−χσ(x)|dydx. (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

 χσ(y)=χτ(y)=p+y1v+y2w

and the regularisation of (5) reads

 I =∫(0,1)2η16∑i=1 ∫πAi(η)q(Z+Ai(η)/2,Z−Ai(η)/2)dZdη

where the area is for . Hence, we obtain

 I=g2τ12π∫(0,1)(1|η2v+w|+1|η2w+v|+1|η2(w+v)−w|)dη. (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

 1√γ+βη2+αη22,withα>0, 4αγ−β2>0.

The anti-derivative is given by

 F(η2)=1√αln(2√α√γ+βη2+αη22+2αη2+β), (7)

see [Prudnikov1981, Section 1.2.52.8] and [Gradstein2015, Section 2.261]. Thus, the integral reduces to

 I=g2τ12π[F1(η2)+F2(η2)+F3(η2)]10,

where is with the parameters

 α1 =|v|2, β1 =2v⋅w, γ1 =|w|2, α2 =|w|2, β2 =2w⋅v, γ2 =|v|2, α3 =|w+v|2, β3 =−2(w+v)⋅w, γ3 =|w|2.

#### 4.1.2 Common edge

Let the reference mappings be given by

 χτ(y)=p+y1v+y2u,χσ(x)=p+x1v+x2w,

such that is the common edge of and starting from . Then, the integral (5) reduces to

 I =∫(0,1)4η31η22(q(A1(η))+η35∑i=2q(Ai(η)))dη (8) =gτgσ24π∫(0,1)2(1|η3(u+v)+η4w−u|+η3|η4η3(u+v)−η3u+w| +η3|η4η3(w+v)−η3w+u|+η3|η4η3u+η3(w+v)−w| +η3|η4η3(w+v)+η3u−w|)dη3dη4 =gτgσ24π5∑i=1Ii.

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 I1

Let us introduce the variables

 a=w,b=v,c=u+v.

We integrate with respect to by using (7) and obtain for the first integral

 I1=1|c|1∫0ln(|η4a+b||c|+(η4a+b)⋅c|η4a+b−c||c|+(η4a+b−c)⋅c)dη4.

 I1=1|c|ln(|a+b||c|+(a+b)⋅c|a+b−c||c|+(a+b−c)⋅c)−1|c|1∫0(h1(η4)−h0(η4))dη4, (9)

where

 h0(η4) =(η4a+b)⋅b+|η4a+b|b⋅^c|η4a+b|2+|η4a+b|(η4a+b)⋅^c, h1(η4)

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

 p=a⋅b|a|2,q2=|b|2|a|2−p2≥0

such that

 |η4a+b|=|a|√(η4+p)2+q2.

If then and the integral simplifies to

 ∫h(η4)dη4=pln(1+1/p).

Otherwise, we have and the substitution

 η4=−p+qsinh(s),dη4=qcosh(s)ds,

yields for the indefinite integral

 ∫h(η4)dη4

We use a variant of the Weierstraß substitution,

 tanh(s/2)=t,sinh(s)=2t1−t2,cosh(s)=1+t21−t2,ds=21−t2dt,

and obtain

 ∫h(η4)dη4

The integrand is now a rational function and we abbreviate it by

 11−t2β0+β1t+β2t2α0+α1t+α2t2.

We decompose it into partial fractions,

 γ11−t+γ21+t+γ3+γ4tα0+α1t+α2t2,

where

 γ1 =12β0+β1+β2α0+α1+α2, γ2 =12β0−β1+β2α0−α1+α2, γ3 =β0−(γ1+γ2)α0, γ4 =α2(γ1−γ2).

The first two terms yield

 F(t)=∫(γ11−t+γ21+t)dt=γ2ln|1+t|−γ1ln|1−t|.

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

 D=|det(a|b|^c)|2/|a|2.

For we have

 G(t)=∫γ3+γ4tα0+α1t+α2t2dt =γ42a2ln∣∣α0+α1t+α2t2∣∣ +2γ3α2−γ4α1α2√Darctan(α1+2α2t√D)

and for

 G(t)=∫γ3+γ4tα0+α1t+α2t2dt=γ4α2ln∣∣∣t+α12α2∣∣∣−2γ3α2−γ4α1α2(2α2t+α1).

We resubstitute

and set

 t0=pq+√p2+q2,t1=p+1q+√(p+1)2+q2.

Finally, we obtain

 1∫0h(η4)dη4=2q(F(t1)−F(t0)+G(t1)−G(t0)). (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 by

 H(a,b,c)=1∫0h(η4)dη4.

We 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

 a=v,b=w,c=u+v,

we have

 I2=1|c|1∫0ln(|η3a+b||c|+(η3a+b)⋅c|η3(a−c)+b||c|+(η3(a−c)+b)⋅c)dη3.

Integration by parts yields an expression almost identical to (9),

 I2=1|c|ln(|a+b||c|+(a+b)⋅c|a−c+b||c|+(a−c+b)⋅c)−1|c|1∫0(h1(η3)−h0(η3))dη3,

where and are given by

 h0(η3) =(η3a+b)⋅b+|η3a+b|b⋅^c|η3a+b|2+|η3a+b|(η3a+b)⋅^c, h1(η3) =(η3(a−c)+b)⋅b+|η3(a−c)+b|b⋅^c|η3(a−c)+b|2+|η3(a−c)+b|(η3(a−c)+b)⋅^c.

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

 χτ(y)=p+y1u1+y2u2,χσ(x)=p+x1v1+x2v2

with common vertex . Then, integration with respect to results in

 I =gτgσ12π(I1+I2).

We only consider , since is obtained by swapping and . Similar to the previous section, we introduce the variables

 a=u1+u2,b(η2)=−v1−η2v2,c=u2,

and integrate with respect to to obtain

We insert Formula (12) for the inner integral, which yields

 I1 =1|c|1∫0ln(|a+b(η2)||c|+(a+b(η2))⋅c|a−c+b(η2)||c|+(a−c+b(η2))⋅c)dη2 −1|c|1∫0(H(a−c,b(η2),c)−H(a,b(η2),c))dη2.

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

 1∫0(H(a−c,b(η2),c)−H(na,b(η2),c))dη2 ≈ n∑i=1ωi(H(a−c,b(η(i)),c)−H(a,b(η(i)),c))

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

 χτ(y)=p1+y1u1+y2u2,χσ(x)=p2+x1v1+x2v2,

and set . Analogously to the previous cases, we pull the region of integration back to by

 A:(0,1)4→π×π,A(η)=⎛⎜ ⎜ ⎜ ⎜⎝η1η1η2η3η3η4⎞⎟ ⎟ ⎟ ⎟⎠,

 I=gτgσ4π∫(0,1)4η1η3|p+η3u1+η3η4u2−η1v1−η1η2v2|dη.

Of the four iterated integrals, we compute two analytically and two by numerical quadrature, e.g.

 I≈n∑k,ℓ=1ωkωℓ∫(0,1)2η(k)η(ℓ)∣∣p+η(ℓ)u1+η4η(ℓ)u2−η(k)v1−η2η(k)v2∣∣dη2dη4,

where the two-dimensional integral is calculated analytically using (11) with

 a=η(k)v2,b=p−η(k)v1+η(ℓ)(u1+u2),c=η(ℓ)u2.

### 4.2 Double layer potential

For the double layer potential , we need to compute integrals of the form

 J=∫σ∫τ(x−y)⋅n4π|x−y|3φ(y)dS(y)dS(x)=gτgσ4π∫π×π(χσ(x)−χτ(y))⋅n|χσ(x)−χτ(y)|3φ(χτ(y))dxdy, (13)

where is the outer unit normal vector at and , i.e.

 φ(χτ(y))=a0+a1y1+a2y2

with coefficients .

#### 4.2.1 Identical elements

For identical elements , we simply have due to

 (y−x)⋅n