Almost complete analytical integration in Galerkin BEM

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 = \int τ \int σ \frac{1}{4 π | x - y |} φ (y) d S (y) ψ (x) d S (x),

which do not vanish. In addition, the kernel function is singular at $x = y$ , so standard quadrature rules for the approximation of $I$ 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 $I$ 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 $I$ 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 $I$ 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$		$i n Ω,$		(1)
	$u$	$= g$		$o n Γ = \partial Ω,$		(1)

in a domain $Ω$ with bounded Lipschitz boundary $Γ$ . If $Ω$ is unbounded, we assume the radiation condition

The representation formula expresses the solution $u$ in terms of its boundary values only,

u (x) = \int_{Γ} u^{*} (x, y) \partial_{n} u (y) d S (y) - \int_{Γ} \partial_{n (y)} u^{*} (x, y) g (y) d S (y), x \in Ω .

Here, $n$ is the unit normal to $Γ$ pointing outwards $Ω$ and $u^{*}$ is the fundamental solution of the Laplace operator,

u^{*} (x, y) = \frac{1}{4 π | y - x |} .

The boundary value problem is hence reduced to the problem of finding the unknown Neumann trace $t = \partial_{n} u$ . To this end, we take the traces in the representation formula and insert the Dirichlet condition to obtain the boundary integral equation

V t = (\frac{1}{2} I + K) g o n Γ,

(2)

where the layer potentials are defined by

(V w) (x)

= \int_{Γ} u^{*} (x, y) w (y) d S (y),

(K w) (x)

= \int_{Γ} \partial_{n (y)} u^{*} (x, y) w (y) d S (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 $(Γ_{h}, T_{h})$ (or simply $Γ_{h}$ ) is a finite collection of non-empty and open elements $τ \subset Γ_{h}$ which satisfies:

$T_{h} = {τ_{n}}_{n = 1}^{N}$ is a triangulation of $Γ_{h}$ , i.e.

$Γ_{h} = N ⋃ n = 1 {¯ τ}_{n} .$

The intersection ${¯ τ}_{n} \cap {¯ τ}_{m}$ of two distinct elements is either empty or consists of a common vertex or edge.
Each $τ$ in $T_{h}$ is a flat triangle with vertices $p_{1}, p_{2}$ and $p_{3}$ . The reference mapping

$χ_{τ} : π \to τ, χ_{τ} (x_{1}, x_{2}) = p_{1} + x_{1} (p_{2} - p_{1}) + x_{2} (p_{3} - p_{2}),$

parametrises $τ$ by the reference triangle

$π = {(x_{1}, x_{2}) ∣ 0 < x_{1} < 1, 0 < x_{2} < x_{1}} \subset R^{2} .$

We denote by

$J_{τ} = (p_{2} - p_{1} ∣ p_{3} - p_{2}) \in R^{3 \times 2}, g_{τ} = \sqrt{det (J_{τ}^{⊤} J_{τ})}$

the Jacobian and the Gram determinant of $χ_{τ}$ respectively and assume $g_{τ} \neq 0$ .

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 $p = 0, 1$ , we denote by

S_{h}^{0} (π) = {1}, S_{h}^{1} (π) = {1 - x_{1}, x_{1} - x_{2}, x_{2}}

the set of reference functions and by

S_{h}^{p} (τ) = span {φ \circ χ_{τ}^{- 1} : φ \in S_{h}^{p} (π)}

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

S_{h}^{p} = {φ : Γ_{h} \to R : φ_{| τ} \in S_{h}^{p} (τ) \forall τ \in T_{h}} .

For $p = 1$ , we moreover require that the functions $φ$ are continuous.

We choose the Lagrangian basis

φ_{n}^{0} (x) = {\begin{matrix} 1 & i f x \in τ_{n}, 0 & e l s e \end{matrix}}, φ_{m}^{1} (x_{j}) = {\begin{matrix} 1 & i f j = m, 0 & e l s e \end{matrix}},

where ${x_{j}}_{j = 1}^{M}$ denotes the set of vertices in $Γ_{h}$ . Then, the ansatz

t_{h}

= N \sum n = 1 t_{n} φ_{n}^{0} \in S_{h}^{0},

t \in R^{N},

g_{h}

= M \sum m = 1 g_{m} φ_{m}^{1} \in S_{h}^{1},

g \in R^{M},

for the approximate boundary data leads to the Galerkin approximation of (2): Find $t \in R^{N}$ such that

V t = (\frac{1}{2} M + K) g,

(3)

where the matrices $M \in R^{N \times M}, V \in R^{N \times N}$ and $K \in R^{N \times M}$ are given by

M [n, m] = ⟨ φ_{n}^{0}, φ_{m}^{1} ⟩, V [n, i] = ⟨ φ_{n}^{0}, V φ_{i}^{0} ⟩, K [n, m] = ⟨ φ_{n}^{0}, K φ_{m}^{1} ⟩,

with $i, n = 1, \dots, N$ and $m = 1, \dots, M$ . The brackets symbolise the usual $L_{2}$ -inner product

⟨ u, v ⟩ = \int Γ_{h} u (x) v (x) d S (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 $V$ an $K$ are of the form

\int Γ_{h} \times Γ_{h} k (x, y) φ (y) d S (y) ψ (x) d S (x) = \sum σ, τ \in T_{h} \int σ \int τ k (x, y) φ (y) d S (y) ψ (x) d S (x),

(4)

where $k = u^{*}, \partial_{n (y)} u^{*}$ is the kernel function and $φ, ψ$ are trial and test functions respectively. Let $I$ be one of the summands for the non-trivial case $τ \subset supp φ$ and $σ \subset supp ψ$ . We transform back to the reference element $π$ ,

	$I$	$= \int σ \int τ k (x, y) φ (y) d S (y) ψ (x) d S (x)$
		$= \int π \times π g_{σ} g_{τ} k (χ_{σ} (x), χ_{τ} (y)) φ (χ_{σ} (x)) ψ (χ_{τ} (y)) d (y) d (x),$

and abbreviate the integrand by $q$ . Since the kernel function $k (x, y)$ is singular at $x = y$ , the integral needs to be regularised. We distinguish between four different cases: the intersection $¯ σ \cap ¯ τ$ 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

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

such that

I = \int Π \int π_{z} q (Z + z / 2, Z - z / 2) d Z d z

with

Π = {z = x - y : x, y \in π}, π_{z} = (π - z / 2) \cap (π + z / 2) .

The singularity of the integrand is now located at $z = 0$ .

Figure 1: The domain of integration $Π$ is split into six triangles $π_{i}$ .

As shown in Figure 1, we decompose $Π$ into six triangles $π_{i} = A_{i} π$ with

A_{1} = (\begin{matrix} 1 & 0 0 & 1 \end{matrix}), A_{2} = (\begin{matrix} 0 & 1 1 & 0 \end{matrix}), A_{3} = (\begin{matrix} 0 & 1 1 & - 1 \end{matrix})

and $A_{i + 3} = - A_{i}$ for $i = 1, 2, 3$ . Thus, we obtain

I = 6 \sum i = 1 \int π \int π_{A_{i} z} q (Z + A_{i} z / 2, Z - A_{i} z / 2) d Z d z .

In the last step, we parametrise $π$ by ${(0, 1)}^{2}$ using the Duffy transformation

z (η) = (\begin{matrix} η_{1} η_{1} η_{2} \end{matrix}), A_{i} (η) = A_{i} z (η)

with Jacobian $η_{1}$ and conclude that $I$ has the representation

I = \int {(0, 1)}^{2} η_{1} 6 \sum i = 1 \int π_{A_{i} (η)} q (Z + A_{i} (η) / 2, Z - A_{i} (η) / 2) d Z d η .

In Section 4, we will see that the integrand is smooth since the Jacobian $η_{1}$ cancels out the singularity of $q$ .

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

χ_{τ} (x_{1}, 0) = χ_{σ} (x_{1}, 0), x_{1} \in (0, 1) .

The kernel function is singular at this edge, i.e. at $(x, y)$ with $x_{2} = y_{2} = 0$ and $x_{1} = y_{1}$ . The regularisation is carried out via the mappings

and reads

I = \int {(0, 1)}^{4} η_{1}^{3} η_{2}^{2} (q (A_{1} (η)) + η_{3} 5 \sum i = 2 q (A_{i} (η))) 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 = \int {(0, 1)}^{4} η_{1}^{3} η_{3} (q (A_{1} (η)) + q (A_{2} (η))) 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 $V$ are of the form

I = \int σ \int τ \frac{1}{4 π | y - x |} d S (y) d S (x) = \frac{g_{τ} g_{σ}}{4 π} \int π \times π \frac{1}{| χ_{τ} (y) - χ_{σ} (x) |} d y d x .

(5)

We proceed like in Section 3 and begin with the case of identical elements.

4.1.1 Identical Elements

Let $v$ and $w$ be the edges of $τ = σ$ with starting point $p$ . Then, the triangle is parametrised by

χ_{σ} (y) = χ_{τ} (y) = p + y_{1} v + y_{2} w

and the regularisation of (5) reads

	$I$	$= \int {(0, 1)}^{2} η_{1} 6 \sum i = 1 \int π_{A_{i} (η)} q (Z + A_{i} (η) / 2, Z - A_{i} (η) / 2) d Z d η$

where the area is $∣ ∣ π_{A_{i} (η)} ∣ ∣ = {(1 - η_{1})}^{2} / 2$ for $i = 1, 2, 3$ . Hence, we obtain

I = \frac{g_{τ}^{2}}{12 π} \int (0, 1) (\frac{1}{| η_{2} v + w |} + \frac{1}{| η_{2} w + v |} + \frac{1}{| η_{2} (w + v) - w |}) d η .

(6)

Figure 2 depicts the complex continuation of the integrand for concrete values of $v$ and $w$ . 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 $f (z)$ of (6) in the complex plane.

The three terms in the integrand are of the form

\frac{1}{\sqrt{γ + β η_{2} + α η_{2}^{2}}}, w i t h α > 0, 4 α γ - β^{2} > 0.

The anti-derivative is given by

F (η_{2}) = \frac{1}{\sqrt{α}} ln (2 \sqrt{α} \sqrt{γ + β η_{2} + α η_{2}^{2}} + 2 α η_{2} + β),

(7)

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

I = \frac{g_{τ}^{2}}{12 π} {[F_{1} (η_{2}) + F_{2} (η_{2}) + F_{3} (η_{2})]}_{0}^{1},

where $F_{i}$ is $F$ with the parameters

$α_{1}$	$= {\| v \|}^{2},$	$β_{1}$	$= 2 v \cdot w,$	$γ_{1}$	$= {\| w \|}^{2},$
$α_{2}$	$= {\| w \|}^{2},$	$β_{2}$	$= 2 w \cdot v,$	$γ_{2}$	$= {\| v \|}^{2},$
$α_{3}$	$= {\| w + v \|}^{2},$	$β_{3}$	$= - 2 (w + v) \cdot w,$	$γ_{3}$	$= {\| w \|}^{2} .$

4.1.2 Common edge

Let the reference mappings be given by

χ_{τ} (y) = p + y_{1} v + y_{2} u, χ_{σ} (x) = p + x_{1} v + x_{2} w,

such that $v$ is the common edge of $σ$ and $τ$ starting from $p$ . Then, the integral (5) reduces to

$I$	$= \int {(0, 1)}^{4} η_{1}^{3} η_{2}^{2} (q (A_{1} (η)) + η_{3} 5 \sum i = 2 q (A_{i} (η))) d η$	(8)
	$= \frac{g_{τ} g_{σ}}{24 π} \int {(0, 1)}^{2} (\frac{1}{\| η_{3} (u + v) + η_{4} w - u \|} + \frac{η_{3}}{\| η_{4} η_{3} (u + v) - η_{3} u + w \|}$
	$+ \frac{η_{3}}{\| η_{4} η_{3} (w + v) - η_{3} w + u \|} + \frac{η_{3}}{\| η_{4} η_{3} u + η_{3} (w + v) - w \|}$
	$+ \frac{η_{3}}{\| η_{4} η_{3} (w + v) + η_{3} u - w \|}) d η_{3} d η_{4}$
	$= \frac{g_{τ} g_{σ}}{24 π} 5 \sum i = 1 I_{i} .$

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 ${(0, 1)}^{2}$ since the triangles do not overlap.

Figure 3: Visualisation of the integrand $f (η_{3}, η_{4})$ of (8).

First integral $I_{1}$

Let us introduce the variables

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

We integrate with respect to $η_{3}$ by using (7) and obtain for the first integral

I_{1} = \frac{1}{| c |} 1 \int 0 ln (\frac{}{| η_{4} a + b | | c | + (η_{4} a + b) \cdot c} | η_{4} a + b - c | | c | + (η_{4} a + b - c) \cdot c) d η_{4} .

Integration by parts leads to

I_{1} = \frac{1}{| c |} ln (\frac{| a + b | | c | + (a + b) \cdot c}{| a + b - c | | c | + (a + b - c) \cdot c}) - \frac{1}{| c |} 1 \int 0 (h_{1} (η_{4}) - h_{0} (η_{4})) d η_{4},

(9)

where

	$h_{0} (η_{4})$	$= \frac{(η_{4} a + b) \cdot b + \| η_{4} a + b \| b \cdot^c}{{\| η_{4} a + b \|}_{4}^{2} + \| η_{4} a + b \| (η_{4} a + b) \cdot^c},$
	$h_{1} (η_{4})$

with $^c = c / | c |$ . We note that $h_{1}$ coincides with $h_{0}$ when $b$ is replaced by $b - c$ and proceed with integrating $h = h_{0}$ . We follow the approach of [RjasanowSteinbach2007, Appendix C.2] and define

p = \frac{a \cdot b}{{| a |}^{2}}, q^{2} = \frac{{| b |}^{2}}{{| a |}^{2}} - p^{2} \geq 0

such that

| η_{4} a + b | = | a | \sqrt{{(η_{4} + p)}^{2} + q^{2}} .

If $q = 0$ then $a = p b$ and the integral simplifies to

\int h (η_{4}) d η_{4} = p ln (1 + 1 / p) .

Otherwise, we have $q > 0$ and the substitution

η_{4} = - p + q sinh (s), d η_{4} = q cosh (s) d s,

yields for the indefinite integral

	$\int h (η_{4}) d η_{4}$

We use a variant of the Weierstraß substitution,

tanh (s / 2) = t, sinh (s) = \frac{2 t}{1 - t^{2}}, cosh (s) = \frac{1 + t^{2}}{1 - t^{2}}, d s = \frac{2}{1 - t^{2}} d t,

and obtain

	$\int h (η_{4}) d η_{4}$

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

\frac{1}{1 - t^{2}} \frac{β_{0} + β_{1} t + β_{2} t^{2}}{α_{0} + α_{1} t + α_{2} t^{2}} .

We decompose it into partial fractions,

\frac{γ_{1}}{1 - t} + \frac{γ_{2}}{1 + t} + \frac{γ_{3} + γ_{4} t}{α_{0} + α_{1} t + α_{2} t^{2}},

where

	$γ_{1}$	$= \frac{1}{2} \frac{β_{0} + β_{1} + β_{2}}{α_{0} + α_{1} + α_{2}},$	$γ_{2}$	$= \frac{1}{2} \frac{β_{0} - β_{1} + β_{2}}{α_{0} - α_{1} + α_{2}},$
	$γ_{3}$	$= β_{0} - (γ_{1} + γ_{2}) α_{0},$	$γ_{4}$	$= α_{2} (γ_{1} - γ_{2}) .$

The first two terms yield

F (t) = \int (\frac{γ_{1}}{1 - t} + \frac{γ_{2}}{1 + t}) d t = γ_{2} ln | 1 + t | - γ_{1} ln | 1 - t | .

The third term depends on the discriminant $D = 4 α_{0} α_{2} - α_{1}^{2}$ of the denominator, which is non-negative due to

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

For $D > 0$ we have

	$G (t) = \int \frac{γ_{3} + γ_{4} t}{α_{0} + α_{1} t + α_{2} t^{2}} d t$	$= \frac{γ_{4}}{2 a_{2}} ln ∣ ∣ α_{0} + α_{1} t + α_{2} t^{2} ∣ ∣$
		$+ \frac{2 γ_{3} α_{2} - γ_{4} α_{1}}{α_{2} \sqrt{D}} arctan (\frac{α_{1} + 2 α_{2} t}{\sqrt{D}})$

and for $D = 0$

G (t) = \int \frac{γ_{3} + γ_{4} t}{α_{0} + α_{1} t + α_{2} t^{2}} d t = \frac{γ_{4}}{α_{2}} ln ∣ ∣ ∣ t + \frac{α_{1}}{2 α_{2}} ∣ ∣ ∣ - \frac{2 γ_{3} α_{2} - γ_{4} α_{1}}{α_{2} (2 α_{2} t + α_{1})} .

We resubstitute

and set

t_{0} = \frac{p}{q + \sqrt{p^{2} + q^{2}}}, t_{1} = \frac{p + 1}{q + \sqrt{{(p + 1)}^{2} + q^{2}}} .

Finally, we obtain

1 \int 0 h (η_{4}) d η_{4} = 2 q (F (t_{1}) - F (t_{0}) + G (t_{1}) - G (t_{0})) .

(10)

Note that the value of the integral only depends on the vectors

a, b, c

. Since it is of importance for the other cases as well, we abbreviate it by

H (a, b, c) = 1 \int 0 h (η_{4}) d η_{4} .

We conclude that $I_{1}$ can be expressed in closed form as

(11)

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

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

I_{2} = \frac{1}{| c |} 1 \int 0 ln (\frac{}{| η_{3} a + b | | c | + (η_{3} a + b) \cdot c} | η_{3} (a - c) + b | | c | + (η_{3} (a - c) + b) \cdot c) d η_{3} .

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

I_{2} = \frac{1}{| c |} ln (\frac{| a + b | | c | + (a + b) \cdot c}{| a - c + b | | c | + (a - c + b) \cdot c}) - \frac{1}{| c |} 1 \int 0 (h_{1} (η_{3}) - h_{0} (η_{3})) d η_{3},

where $h_{0}$ and $h_{1}$ are given by

	$h_{0} (η_{3})$	$= \frac{(η_{3} a + b) \cdot b + \| η_{3} a + b \| b \cdot^c}{{\| η_{3} a + b \|}_{3}^{2} + \| η_{3} a + b \| (η_{3} a + b) \cdot^c},$
	$h_{1} (η_{3})$	$= \frac{(η_{3} (a - c) + b) \cdot b + \| η_{3} (a - c) + b \| b \cdot^c}{{\| η_{3} (a - c) + b \|}_{3}^{2} + \| η_{3} (a - c) + b \| (η_{3} (a - c) + b) \cdot^c} .$

Thus, $I_{2}$ can be computed analogously to (11) by

(12)

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

$j$	$a$	$b$	$c$
$2$	$v$	$w$	$u + v$
$3$	$v$	$u$	$w + v$
$4$	$u + v + w$	$- w$	$u$
$5$	$u + v + w$	$- w$	$v + w$

Table 1: Values for

a, b, c

in (12) to compute

I_{j}

4.1.3 Common vertex

We consider the configuration

χ_{τ} (y) = p + y_{1} u_{1} + y_{2} u_{2}, χ_{σ} (x) = p + x_{1} v_{1} + x_{2} v_{2}

with common vertex $p$ . Then, integration with respect to $η_{1}$ results in

	$I$
		$= \frac{g_{τ} g_{σ}}{12 π} (I_{1} + I_{2}) .$

We only consider $I_{1}$ , since $I_{2}$ is obtained by swapping $u_{i}$ and $v_{i}$ . Similar to the previous section, we introduce the variables

a = u_{1} + u_{2}, b (η_{2}) = - v_{1} - η_{2} v_{2}, c = u_{2},

and integrate with respect to $η_{4}$ to obtain

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

	$I_{1}$	$= \frac{1}{\| c \|} 1 \int 0 ln (\frac{\| a + b (η_{2}) \| \| c \| + (a + b (η_{2})) \cdot c}{\| a - c + b (η_{2}) \| \| c \| + (a - c + b (η_{2})) \cdot c}) d η_{2}$
		$- \frac{1}{\| c \|} 1 \int 0 (H (a - c, b (η_{2}), c) - H (a, b (η_{2}), c)) d η_{2} .$

The first integral can be written in the form of $I_{1}$ from Section 4.1.2, i.e.

with $~ a = - v_{2}$ , $~ b = u_{1} + u_{2} - v_{1}$ , $~ c = u_{2}$ , 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 \int 0 (H (a - c, b (η_{2}), c) - H (n a, b (η_{2}), c)) d η_{2}$	$\approx$
		$n \sum i = 1 ω_{i} (H (a - c, b (η^{(i)}), c) - H (a, b (η^{(i)}), c))$

with weights $ω_{i} > 0$ and nodes $η^{(i)} \in [0, 1]$ .

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) = p_{1} + y_{1} u_{1} + y_{2} u_{2}, χ_{σ} (x) = p_{2} + x_{1} v_{1} + x_{2} v_{2},

and set $p = p_{1} - p_{2}$ . Analogously to the previous cases, we pull the region of integration back to ${(0, 1)}^{4}$ by

A : {(0, 1)}^{4} \to π \times π, A (η) = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} η_{1} η_{1} η_{2} η_{3} η_{3} η_{4} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠,

leading to

I = \frac{g_{τ} g_{σ}}{4 π} \int {(0, 1)}^{4} \frac{η_{1} η_{3}}{| p + η_{3} u_{1} + η_{3} η_{4} u_{2} - η_{1} v_{1} - η_{1} η_{2} v_{2} |} d η .

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

I \approx n \sum k, ℓ = 1 ω_{k} ω_{ℓ} \int {(0, 1)}^{2} \frac{η^{(k)} η^{(ℓ)}}{∣ ∣ p + η^{(ℓ)} u_{1} + η_{4} η^{(ℓ)} u_{2} - η^{(k)} v_{1} - η_{2} η^{(k)} v_{2} ∣ ∣} d η_{2} d η_{4},

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

a = η^{(k)} v_{2}, b = p - η^{(k)} v_{1} + η^{(ℓ)} (u_{1} + u_{2}), c = η^{(ℓ)} u_{2} .

4.2 Double layer potential

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

J = \int σ \int τ \frac{(x - y) \cdot n}{4 π {| x - y |}^{3}} φ (y) d S (y) d S (x) = \frac{g_{τ} g_{σ}}{4 π} \int π \times π \frac{(χ_{σ} (x) - χ_{τ} (y)) \cdot n}{{| χ_{σ} (x) - χ_{τ} (y) |}_{σ}^{3}} φ (χ_{τ} (y)) d x d y,

(13)

where $n$ is the outer unit normal vector at $τ$ and $φ \in S_{h}^{1} (τ)$ , i.e.

φ (χ_{τ} (y)) = a_{0} + a_{1} y_{1} + a_{2} y_{2}

with coefficients $a_{0}, a_{1}, a_{2} \in R$ .

4.2.1 Identical elements

For identical elements $σ = τ$ , we simply have $J = 0$ due to

(y - x) \cdot n

Almost complete analytical integration in Galerkin BEM

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

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

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

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

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

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

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

1 Introduction

2 Preliminaries

Definition 1 (Mesh).

Remark 1.

Definition 2 (Boundary element spaces).