Exterior Mixed Boundary Value Problem

−1

2I+K_κ



γ^0,extu, s



∂Ω

=⟨V_κg_N, s⟩_∂Ω for all s∈H^−1/2(∂Ω), and

⟨D_κγ^0,extu, t⟩_∂Ω =



−1

2I −K_κ^∗

 gN, t



∂Ω

for allt∈H^1/2(∂Ω), (3.52) respectively.

The arguments for the unique solvability of the boundary integral equation (3.51) and the corresponding variational problem (3.52) for κ² not coinciding with an eigenvalue of (3.36) are the same as in Section 3.5.2.

Theorem 3.30. If u∈H_loc¹ (Ω^ext, ∆+κ²) is a solution to the exterior Neumann problem (3.48), then the Dirichlet trace γ^0,extu satisfies the boundary integral equation (3.51) and u has the representation (3.49).

Conversely, if γ^0,extu satisfies the boundary integral equation (3.51), then the represen-tation formula (3.49) defines a solution u ∈ H_loc¹ (Ω^ext, ∆+κ²) to the interior Neumann problem (3.48).

Note that Theorems 3.25 and 3.30 ensure that for all κ∈R+ it holds



−1

2I−K_κ^∗



g_N∈ImDκ

and thus one of the methods mentioned at the end of the previous section can be used to compute a solution to (3.51) withκ² coinciding with an eigenvalue of (3.36).

3.5.6 Exterior Mixed Boundary Value Problem

Finally, let us consider the exterior mixed boundary value problem



















∆u+κ²u= 0 inΩ^ext, γ^0,extu=g_D onΓ_D, γ^1,extu=g_N onΓ_N,











∇u(x), x

∥x∥



−iκu(x)









=O

 1

∥x∥²



for ∥x∥ → ∞

(3.53)

with an unbounded domain Ω^ext := R³ \Ω and boundary conditions gD ∈ H^1/2(ΓD), g_N∈H^−1/2(Γ_N). In the smooth case the solution to (3.53) is given by

u(x) =−



ΓD

γ^1,extu(y)v_κ(x,y) ds_y−



ΓN

g_N(y)v_κ(x,y) ds_y +



ΓD

gD(y) ∂vκ

∂n_y(x,y) dsy+



ΓN

γ^0,extu(y) ∂vκ

∂n_y(x,y) dsy for x∈Ω^ext with unknown Cauchy data γ^0,extu|_ΓN and γ^1,extu|_ΓD. Similarly as for the interior mixed problem we have the boundary integral equations

(V_κγ^1,extu)(x) =−1

2γ^0,extu(x) + (K_κγ^0,extu)(x) for x∈Γ_D⊂∂Ω, (D_κγ^0,extu)(x) =−1

2γ^1,extu(x)−(K_κ^∗γ^1,extu)(x) for x∈Γ_N⊂∂Ω.

(3.54)

Inserting the functions t, sdefined by (3.39) into (3.54) yields (V_κs)(x)−(K_κt)(x) =−1

2g˜_D(x) + (K_κg˜_D)(x)−(V_κg˜_N)(x) for x∈Γ_D, (K_κ^∗s)(x) + (D_κt)(x) =−1

2g˜_N(x)−(K_κ^∗g˜_N)(x)−(D_κg˜_D)(x) for x∈Γ_N

(3.55)

and the equivalent variational formulation

a(s, t, q, r) =F(q, r) for all q∈H^−1/2(ΓD), r ∈H^1/2(ΓN) (3.56) with unknown functions s∈H^−1/2(Γ_D),t∈H^1/2(Γ_N), the bilinear form

a(s, t, q, r) :=⟨V_κs, q⟩_ΓD− ⟨K_κt, q⟩_ΓD+⟨K_κ^∗s, r⟩_ΓN+⟨D_κt, r⟩_ΓN and the right-hand side

F(q, r) :=



−1

2I+Kκ



˜ gD, q



ΓD

− ⟨V_κg˜N, q⟩_ΓD+



−1

2I−K_κ^∗



˜ gN, r



ΓN

− ⟨D_κg˜D, r⟩_ΓN. Theorem 3.31. If u∈H_loc¹ (Ω^ext, ∆+κ²) is a solution to the mixed problem (3.53), then the functions s, t satisfy the system of boundary integral equations (3.55) and u has the representation

u(x) =−(V_κ(s+ ˜g_N))(x) + (W_κ(t+ ˜g_D))(x) for x∈Ω^ext. (3.57) Conversely, if s, t satisfy the system of boundary integral equations (3.55), then the representation formula (3.57) defines a solution u ∈ H_loc¹ (Ω^ext, ∆+κ²) to the interior mixed problem (3.53).

4 Discretization and Numerical Realization

In this section we will describe the discretization of the boundary integral equations covered in the previous part. We will discuss the collocation method as well as the more complicated Galerkin scheme. The solution will be sought in approximating function spaces introduced below.

We restrict ourselves to triangular meshes approximating the corresponding boundary

∂Ω, i.e.,

∂Ω ≈

E



k=1

τ_k

and assume that neighbouring triangular elements either share a whole edge or a sin-gle vertex. Each element τ_k ⊂ R³ with nodes x^k1,x^k2,x^k3 can be described via the parametrization

R_k(ξ) :=x^k1+R_kξ for ξ∈ˆτ , (4.1) whereR_k denotes the matrix

R_k :=

x^k2−x^k1 x^k3 −x^k1

=





x^k₁² −x^k₁¹ x^k₁³ −x^k₁¹ x^k₂² −x^k₂¹ x^k₂³ −x^k₂¹ x^k₃² −x^k₃¹ x^k₃³ −x^k₃¹



∈R^3×2 (4.2) and τˆ⊂R² represents the reference triangle (see Figure 4.1a)

ˆ

τ :={ξ∈R²: 0< ξ1 <1,0< ξ2<1−ξ1}.

Remark 4.1. Note that a function f defined on τ_k can be identified with a function fˆ defined on τˆand vice versa, i.e.,

f(x) =f(R_k(ξ)) =: ˆf(ξ) for x∈τ_k,ξ ∈τ .ˆ

In the following text we use the symbol ∂Ω to denote the discretized boundary.

4.1 Piecewise Constant Basis Functions

For every elementτ_k we define the functionψ_k (see Figure 4.2a) as ψk(x) :=

1 for x∈τ_k, 0 otherwise.

Following Remark 4.1, the functionψ_k can alternatively be defined via the reference func-tion

ψ(ξ) :=ˆ

1 for ξ∈τ ,ˆ 0 otherwise.

50 4 Discretization and Numerical Realization

ξ1

ξ₂

ˆ τ

1 1

(a) Reference triangle.

τ_k x^k1

x^k3

x^k2 x1

x₃

x₂

(b) General triangle.

Figure 4.1: Triangular mesh elements.

τ_k ψk(x)

(a) Piecewise constant basis function.

x^k ϕk(x)

(b) Piecewise affine basis function.

Figure 4.2: Basis functions.

Furthermore, we define the linear spaceT_ψ(∂Ω) as Tψ(∂Ω) := span{ψ_k}^E_k=1,

whereE denotes the number of elements. Obviously,dimT_ψ(∂Ω) =Eand every complex-valued function g_ψ ∈T_ψ(∂Ω)can be represented as

g_ψ =

E



k=1

g_kψ_k

and thus it can be identified with a vectorg:= [g₁, . . . , g_E]^T∈C^E.

A smooth enough function gcan be approximated by a linear combination of piecewise constant functions in several ways. The simplest approach is to set

g_k=g(x^k∗)

withx^k∗ denoting the midpoint ofτk. Alternatively, we can set

However, since we deal with functions in Lebesgue and Sobolev spaces, where we do not distinguish between functions differing on a measure zero set, the above mentioned approxi-mations may prove impossible. Hence, we introduce the projectionP_ψ:L²(∂Ω)→T_ψ(∂Ω) minimizing the error

∥g−P_ψg∥_L2(∂Ω), i.e.,

P_ψg:= arg min

gψ∈T_ψ(∂Ω)

∥g−g_ψ∥_L2(∂Ω).

To find the coefficients gk for a real-valued function g we introduce the function f:R^E →Ras

For the partial derivatives off we get

∂f

Moreover, for the Hessian of f we have H_f[j, i] = ∂²f

∂g_j∂g_i = 2⟨ψ_i, ψ_j⟩,

which is the double Gram matrix corresponding to the basis ofT_ψ(∂Ω). Hence, the Hessian is positive definite and to find the minimum of f we set _∂g^∂f

j = 0 to obtain the system of

52 4 Discretization and Numerical Realization

where∆k denotes the surface ofτk, we get for the piecewise constant approximation g^k= 1

∆k



τ_k

g(x) ds_x.

Let us now consider the case when g= Reg+ i Img is complex-valued. We thus search for coefficientsg_k= Reg_k+ i Img_k. For theL² error norm we get

and therefore, the problem is reduced to searching for appropriate coefficients of the two real-valued functions Reg andImg as was described in the previous paragraph.

The spaceTψ(∂Ω)will be used for the approximation of the Neumann data, i.e., of the normal derivatives on ∂Ω.

4.2 Continuous Piecewise Affine Basis Functions

LetN denote the number of all nodes of a given discretization. Then we define the family of functions {ϕ_k}^N_k=1 continuous over the whole discretized boundary (see Figure 4.2b) as

ϕ_k(x) := piecewise affine otherwise.

Furthermore, we define the linear spaceT_ϕ(∂Ω)as T_ϕ(∂Ω) := span{ϕ_k}^N_k=1.

Obviously, dimT_ϕ(∂Ω) = N. A restriction of a function ϕ_k to an element τ_ℓ ⊂ suppϕ_k can be identified with one of the functions ϕˆ1,ϕˆ2,ϕˆ3 defined on the reference triangle as

ˆ

ϕ1(ξ) := 1−ξ1−ξ2, ϕˆ2(ξ) :=ξ1, ϕˆ3(ξ) :=ξ2 for ξ∈τ .ˆ

k\τ_k x^k1 x^k2 x^k3 local indices

1 1 2 3

global indices

2 2 4 3

3 4 5 3

... ... ... ...

τ₁ τ3

τ2

x³

x⁴ x²

x¹ x⁵

Figure 4.3: Connectivity table.

We use the connectivity table (see Figure 4.3) to determine the mapping between global and local indices of nodes building the triangular mesh. For example, the point x⁴ is the second vertex of the triangleτ2, i.e., x⁴=x²². Furthermore, forϕ4|_τ2 we have

ϕ4|_τ2(x) =ϕ4|_τ2(R₄(ξ)) = ˆϕ2(ξ) =ξ1 for x∈τ2,ξ∈ˆτ .

As in the case of piecewise constant basis functions we can identify a function T_ϕ(∂Ω)∋g_ϕ=

N



k=1

g_kϕ_k

with a vector [g₁, . . . , g_N]^T ∈ C^N. For the approximation of a smooth enough functiong we can either use an interpolation, i.e.,

g_k=g(x^k).

For a function g ∈L²(∂Ω) it is more appropriate to define the projection P_ϕ: L²(∂Ω)→ Tϕ(∂Ω) as

Pϕg:= arg min

gϕ∈Tϕ(∂Ω)

∥g−gϕ∥_L2(∂Ω).

The coefficients for real-valued functions can be found in the same way as for the piecewise constant functions as the solution to the system of linear equations

N



k=1

g_k⟨ϕ_k, ϕ_j⟩_L2(∂Ω)=⟨g, ϕ_j⟩_L2(∂Ω) for j∈ {1, . . . , N}.

In the case of complex-valued functions we can separately search for the real parts and the imaginary parts of the coefficients.

The space T_ϕ(∂Ω) will be used for the approximation of the Dirichlet data, i.e., the values of the solution on∂Ω.

4.3 Discretized Boundary Integral Equations

In the following sections we describe the discretization of the boundary integral equations derived in Section 3.5. Although the collocation method will be mentioned, we will prefer

54 4 Discretization and Numerical Realization

the well studied Galerkin scheme based on the corresponding variational formulations. To approximate the Dirichlet and Neumann boundary data we use continuous piecewise affine functions and piecewise constant functions, respectively.

4.3.1 Interior Dirichlet Boundary Value Problem

The solution to the interior Dirichlet boundary value problem (3.24) is given by the rep-resentation formula (3.25) with unknown Neumann data g_N := γ^1,intu. To compute the missing values we use the Fredholm integral equation of the first kind

(V_κg_N)(x) = 1

2g_D(x) + (K_κg_D)(x) forx∈∂Ω or equivalently



∂Ω

g_N(y)v_κ(x,y) ds_y= 1

2g_D(x) +



∂Ω

g_D(y)∂v_κ

∂ny

(x,y) ds_y forx∈∂Ω (4.3) with the fundamental solution vκ and the normal derivative

∂v_κ

∂n_y(x,y) = 1 4π

e^{iκ∥x−y∥}

∥x−y∥³(1−iκ∥x−y∥)⟨x−y,n(y)⟩.

The corresponding variational formulation reads

⟨V_κg_N, s⟩_∂Ω =

1

2I+K_κ

 g_D, s



∂Ω

for all s∈H^−1/2(∂Ω) or



∂Ω

s(x)



∂Ω

gN(y)vκ(x,y) dsydsx

= 1 2



∂Ω

s(x)g_D(x) ds_x+



∂Ω

s(x)



∂Ω

g_D(y)∂v_κ

∂n_y(x,y) ds_yds_x. We seek the solution in the approximate form

g_N≈g_N,h :=

E



ℓ=1

g_ℓ^Nψ_ℓ∈T_ψ(∂Ω) (4.4)

which can be identified with a vectorg^N∈C^E. Furthermore, we use theL² projection to approximate the given Dirichlet data

gD≈gD,h:=

N



j=1

g_j^Dϕj ∈Tϕ(∂Ω), (4.5)

which corresponds to a vector g^D ∈ C^N. The solution gN,h is given by the Galerkin For the left-hand side of (4.6) we obtain

E The right-hand side of (4.6) yields

N

56 4 Discretization and Numerical Realization

The Galerkin equations (4.6) can thus be rewritten into the system of linear equations Vκ,hg^N=

1

2Mh+Kκ,h

 g^D.

The approximate solution to (3.24) can be obtained using the discretized representation formula

u(x)≈u_h(x) :=

E



ℓ=1

g_ℓ^N



τℓ

v_κ(x,y) ds_y−

N



j=1

g^D_j



∂Ω

ϕ_j(y)∂v_κ

∂n_y(x,y) ds_y. (4.11) Another approach to the discretization of the boundary integral equations is the col-location method. Because the relation (4.3) is valid for all x ∈ ∂Ω, we can insert the midpointsx^k∗ into (4.3) to obtain the system of linear equations

V_κ,hg^N=

1

2M_h+K_κ,h

 g^D

with

C^E×E ∋V_κ,h[k, ℓ] := 1 4π



τℓ

e^{iκ∥xk∗−y∥}

∥x^k∗−y∥ds_y, R^E×N ∋Mh[k, j] :=ϕj(x^k∗),

C^E×N ∋K_κ,h[k, j] :=



∂Ω

ϕ_j(y) 1 4π

e^{iκ∥xk∗−y∥}

∥x^k∗−y∥³(1−iκ∥x^k∗−y∥)⟨x^k∗−y,n(y)⟩ds_y. Contrary to the Galerkin scheme, the system matrix V_κ,h arising from the collocation method is in general non-symmetric and the stability of the method is still an open problem.

Although the collocation scheme can be derived for other boundary integral equations in the same way, in the following sections we prefer the Galerkin discretization.

4.3.2 Interior Neumann Boundary Value Problem

The solution to the interior Neumann boundary value problem (3.32) is given by the representation formula (3.33) or by its discretized version (4.11). To compute the missing Dirichlet datag_D:=γ^0,intu we use the hypersingular equation

(DκgD)(x) = 1

2gN(x)−(K_κ^∗gN)(x) forx∈∂Ω or

−γ^1,int



∂Ω

g_D(y)∂v_κ

∂ny

(x,y) ds_y = 1

2g_N(x)−



∂Ω

g_N(y)∂v_κ

∂nx

(x,y) ds_y for x∈∂Ω

with

The corresponding variational problem reads

⟨D_κgD, t⟩_∂Ω = Inserting the approximations of the Cauchy data (4.4), (4.5) into (4.12) we obtain the Galerkin equations Using Theorem 3.14, the left-hand side of (4.13) yields

N

For the right-hand side of (4.13) we obtain

E

58 4 Discretization and Numerical Realization

Γ_D

Γ_N

E_D (Dirichlet elements) E_N (Neumann elements) ND (Dirichlet nodes) N_N(Neumann nodes)

Figure 4.4: Triangulation of∂Ω for the mixed boundary value problem.

with matricesM_h andK_κ,hdefined by (4.8) and (4.9), respectively. The Galerkin equations (4.13) can now be represented as

D_κ,hg^D=

1

2M^T_h −K^T_κ,h

 g^N.

Note that in the system of linear equations above the transpositions are not conjugate transpositions and only involve reordering of the original matrices.

4.3.3 Interior Mixed Boundary Value Problem

For the solution to the interior mixed boundary value problem (3.37) we consider the symmetric approach given by the variational formulation

a(s, t, q, r) =F(q, r) for all q∈H^−1/2(ΓD), r ∈H^1/2(ΓN) (4.15) with unknown functions s∈H^−1/2(Γ_D),t∈H^1/2(Γ_N), the bilinear form

a(s, t, q, r) :=⟨V_κs, q⟩_ΓD − ⟨K_κt, q⟩_ΓD+⟨K_κ^∗s, r⟩_ΓN+⟨D_κt, r⟩_ΓN, the right-hand side

F(q, r) :=

1

2I+K_κ



˜ g_D, q



ΓD

− ⟨V_κ˜g_N, q⟩_ΓD+

1

2I−K_κ^∗



˜ g_N, r



ΓN

− ⟨D_κ˜g_D, r⟩_ΓN and some fixed prolongations of the Cauchy datag˜D,g˜N.

Before we proceed, we divide the boundary elements into two groups, E_D denoting the set of all elements with the Dirichlet boundary condition and EN denoting the set of Neumann elements. Similarly, we divide the set of boundary nodes into sets ND and NN

(see Figure 4.4). Note that in the following text we also use the symbols E_D,E_N,N_D and NN to denote the cardinality of the corresponding sets.

Note that for smooth enough functions t, swe have t|_ΓD = 0 and s|_ΓN = 0. Thus, we seek the unknown functions in the approximate forms

t≈t_h := 

i∈NN

t_iϕ_i ∈T_ϕ(Γ_N), s≈s_h := 

j∈ED

s_jψ_j ∈T_ψ(Γ_D) (4.16)

represented by vectors t ∈ C^NN and s ∈ C^ED, respectively. Moreover, we define the

N = 0 with Γ_N⁻ denoting the Neumann part of the boundary without elements adjacent to Γ_D and g˜_N,h|_ΓD = 0. Inserting the approximate representations (4.16), (4.17) into the variational formulation (4.15) we obtain the system of Galerkin equations

a(s_h, t_h, ψ_k, ϕ_ℓ) =F(ψ_k, ϕ_ℓ) for allk∈E_D, ℓ∈N_N. (4.18) For the left-hand side of (4.18) we obtain

⟨V_κs, ψ_k⟩_ΓD = 

Note that the matrices from the previous formulae are submatrices of (4.7), (4.9) and (4.14) with dimensions

V_κ,h ∈C^ED×ED, K_κ,h ∈C^ED×NN, D_κ,h∈C^NN×NN. For the right-hand side of (4.18) we have

⟨˜gD, ψk⟩_ΓD = 

60 4 Discretization and Numerical Realization

with Γ_D⁺ denoting the Dirichlet part of the boundary together with Neumann elements adjacent to Γ_D. Again, the matrices from the above relations are submatrices of (4.8), (4.9), (4.7) and (4.14) with dimensions

M¯h ∈R^ED×ND, K¯κ,h∈C^ED×ND, V¯κ,h∈C^ED×EN,

¯¯

Mh ∈R^EN×NN, K¯¯κ,h∈C^EN×NN, D¯κ,h ∈C^NN×ND. The system of Galerkin equations can now be rewritten as



Note that similarly as for the interior Neumann problem the transpositions do not involve the complex conjugation.

4.3.4 Exterior Dirichlet Boundary Value Problem

In the following three sections we describe the solution to exterior boundary value problems.

Since the discretization techniques are identical to those mentioned above, the discussion will be more succinct.

The approximate solution to the exterior Dirichlet boundary value problem (3.43) is given by the discretized representation formula (3.44), i.e.,

u(x)≈u_h(x) :=− with unknown Neumann data g_N,h. To obtain the missing data we use the variational formulation corresponding to the Fredholm integral equation (3.45), i.e.,

⟨V_κg_N, s⟩_∂Ω = Inserting the approximations of the Cauchy data (4.4), (4.5) into (4.20) we obtain the system of Galerkin equations

E

and the related system of linear equations V_κ,hg^N=

with the matrices given by (4.7), (4.8) and (4.9).

4.3.5 Exterior Neumann Boundary Value Problem

The solution to the exterior Neumann boundary value problem (3.48) is given by the representation formula (3.49) and its discretized variant (4.19). To compute the missing Dirichlet data we use the hypersingular boundary integral equation (3.51) and the related variational problem Substituting the approximate boundary data g_D,h,g_N,hinto (4.21) we obtain the Galerkin equations

62 4 Discretization and Numerical Realization

and the system of linear equations Dκ,hg^D=



−1

2M^T_h −K^T_κ,h

 g^N,

given by the matrices (4.14), (4.8) and (4.9). Again, the transpositions do not involve the complex conjugation.

4.3.6 Exterior Mixed Boundary Value Problem

Lastly, let us consider the exterior mixed boundary value problem (3.53). Similarly as in the case of the interior mixed problem we use the symmetric approach given by the variational formulation

a(s, t, q, r) =F(q, r) for all q∈H^−1/2(Γ_D), r ∈H^1/2(Γ_N) (4.22) with unknown functions s∈H^−1/2(Γ_D),t∈H^1/2(Γ_N), the bilinear form

a(s, t, q, r) :=⟨V_κs, q⟩_ΓD− ⟨K_κt, q⟩_ΓD+⟨K_κ^∗s, r⟩_ΓN+⟨D_κt, r⟩_ΓN the right-hand side

F(q, r) :=



−1

2I+K_κ



˜ g_D, q



Γ_D

− ⟨V_κ˜g_N, q⟩_ΓD+



−1

2I−K_κ^∗



˜ g_N, r



Γ_N

− ⟨D_κ˜g_D, r⟩_ΓN and some fixed prolongations of the Cauchy data˜gD,g˜N. Inserting the discretized functions s_h, t_hfrom (4.16) and the prolongationsg˜_D,g˜_Nfrom (4.17) into (4.22) we obtain the system of Galerkin equations

a(s_h, t_h, ψ_k, ϕ_ℓ) =F(ψ_k, ϕ_ℓ) for all k∈E_D, ℓ∈N_N and the related system of linear equations



V_κ,h −K_κ,h K^T_κ,h D_κ,h

 s t



=

 −V¯_κ,h −¹₂M¯_h+ ¯K_κ,h

−¹₂M¯¯^T_h −K¯¯^T_κ,h −D¯_κ,h

 g^N g^D



with the matrices defined in Section 4.3.3. The transpositions in the above given system do not involve the complex conjugation of the elements.

4.4 Integration over Boundary Elements

There are several approaches to the computation of the double integrals arising from the Galerkin discretization of the boundary integral equations. One way is to use a numerical quadrature for both integrals. It should be, however, noted, that numerical integration introduces an additional error and one has to be cautious when dealing with singularities in the integrands. As a workaround to this problem the Duffy transformation proposed in

[8] (see also [16]) can be used. However, we will not describe this technique and will prefer analytic integration when possible.

Due to the complexity of the integrands it seems impossible to compute both integrals analytically. However, it is possible to integrate the inner integral analytically and for the remaining one use a suitable quadrature scheme. In this section we describe a combination of analytic and numerical integration as proposed in [17].

4.4.1 Numerical Integration

In the following text we assume the discretization of∂Ω given by

∂Ω≈

E



k=1

τ_k,

whereτ_kare triangles inR³. Recall from (4.1) and (4.2), that for every elementτ we have the parametrisation

R(ξ) :=x¹+Rξ =x¹+

x²−x¹ ξ₁+

x³−x¹

ξ₂ forξ ∈τ ,ˆ whereτˆdenotes the reference triangle (see Figure 4.1a)

ˆ

τ :={ξ∈R²: 0< ξ₁ <1,0< ξ₂<1−ξ₁}.

For an arbitrary function f defined on τ we thus have



τ

f(x) ds_x =



ˆ τ

f(R(ξ))



x²−x¹

×

x³−x¹

 ds_ξ= 2∆_τ



ˆ τ

f(ξ) dsˆ _ξ withfˆ(ξ) :=f(R(ξ))and ∆τ denoting the surface of the triangleτ, i.e.,

∆_τ =



τ

1 ds_x.

Hence, it is only sufficient to develop a quadrature scheme corresponding to the refer-ence triangle τˆ. We approximate the integration as



ˆ τ

fˆ(ξ) ds_ξ≈ 1 2

M



k=1

ω_kfˆ(ξ^k) with unknown parameters ω_k and ξ^k∈τˆ for k∈ {1, . . . , M}.

In numerical examples presented in the last section we use a 7-point scheme, which is exact for polynomials up to order 5, i.e., for functions in the linear space

P5(ˆτ) := span

 ξ^a₁ξ₂^b



, 0≤a+b≤5.

64 4 Discretization and Numerical Realization

The quadrature pointsξ^k and corresponding wages ωk are defined as follows;

ξ¹ := 1 For the derivation of the presented 7-point rule and less accurate quadratures see [17], Section C.1.

4.4.2 Analytic Integration

Considering a triangulation of ∂Ω, the integration over the boundary reduces to integra-tion over the set of individual elements. For this purpose, it is convenient to introduce a new coordinate system corresponding to a general boundary element τ ⊂R³ with vertices x¹,x²,x³. The vertex x¹ will be considered as the origin of the coordinate system. Or-thogonal coordinate vectors r¹,r²,nintroduced in the following paragraph are visualized in Figure 4.5a.

x¹

(a) Coordinate system corresponding toτ.

x¹

x² x³

s t

(b) Triangle parametrization.

Figure 4.5: Analytic integration over triangles.

We define the vector r² as

r²:= 1 tτ

(x³−x²),

wheretτ denotes the length of the side x²x³, i.e., tτ :=∥x³−x²∥.

Before we define the vectorr¹, we introduce an auxiliary pointx^∗ lying on the intersection of the line defined by the side x²x³ and the altitude coming through x¹. Thus, we can write

x^∗:=x²+t∗r², wheret∗ can be computed as

t∗ :=⟨x¹−x²,r²⟩.

Then we define

r¹ := 1

s_τ(x^∗−x¹) with

sτ :=∥x^∗−x¹∥. (4.23)

The last coordinate vector orthogonal to both r¹ andr² is defined as n:=r¹×r².

Note that forn we have an alternative representation n= (x²−x¹)×(x³−x²)

∥(x²−x¹)×(x³−x²)∥. Any point x∈R³ can now be expressed as

x=x¹+s_xr¹+t_xr²+u_xn, where

sx :=⟨x−x¹,r¹⟩, tx :=⟨x−x¹,r²⟩, ux :=⟨x−x¹,n⟩.

To derive the parametrization of the triangle we introduce two more quantities, namely the tangents of the angles between the altitudex¹x^∗and the sidesx¹x²andx¹x³, respec-tively. These values can be expressed as

α₁ :=−t∗

s_τ, α₂:= t_τ−t∗

s_τ . (4.24)

Using these parameters together with sτ defined in (4.23) we can parametrize the triangle τ as

τ =

y=y(s, t) =x¹+sr¹+tr²: 0< s < s_τ, α₁s < t < α₂s . See Figure 4.5b for the illustration of the parametrization.

Also note that for any point x∈R³ and a pointy∈τ we have

∥x−y∥²=∥(s_x−s)r¹+ (tx−t)r²+uxn∥²= (s−sx)²+ (t−tx)²+u²_x.

66 4 Discretization and Numerical Realization

4.4.3 Single Layer Potential Operator

To compute the matrices arising from the Galerkin equations we need to evaluate the corresponding double surface integrals. In the following three sections we provide analytic evaluations of the singular parts of the inner integrals related to the single layer potential matrixVκ,h, the double layer potential matrixKκ,h and the hypersingular operator matrix D_κ,h. The outer integrals can be computed using the 7-point quadrature rule discussed in Section 4.4.1.

The analytic formulae corresponding to the single layer potential operator and the double layer potential operator are given in [17], Section C.2, and we only provide the final results with some discussion. For the hypersingular operator we provide a more detailed computation.

The entries of the single layer operator matrix are given by V_κ,h[k, ℓ] := 1 For the latter integrand we have

x→ylim

and thus the integral is not singular and can be evaluated numerically. For the remaining singular part of the integral we use an analytic computation based on the local coordinate system introduced in Section 4.4.2 combined with a quadrature scheme. For the inner integral we obtain (see [17], Section C.1.2.)

S_V(τ,x) := 1

with parameters lying on the line defined by the triangle sidesx¹x² and x¹x³, respectively.

The previous formulae provide an analytic computation of the integral for general parameters. However, we need to discuss the singularities that may appear in (4.25).

Although the arguments of the logarithms are non-negative, they may vanish. However, in both cases the situation is similar to

x→0lim+ can set the arctangent part equal to zero due to the jump property of Vκ (3.5) and the boundedness of the arctangent function. Lastly, in the case when s= p, the last part of the function F becomes

2uxarctanαq ux

.

4.4.4 Double Layer Potential Operator

The entries in the double layer potential matrix are given by Kκ,h[k, j] := 1 Again, the latter integrand is bounded and we can use the 7-point scheme to compute the double surface integral. For the first integral it is sufficient to compute the local element contributions

with the affine functionϕ₁ defined onτ by

ϕ1(y) :=

68 4 Discretization and Numerical Realization

Remark 4.2. For integrals with other affine basis functions it is sufficient to rotate the triangle.

Remark 4.3. Forxlying in the plane defined by the triangle τ (and especially in the case of x∈τ) the integrand is equal to zero and we have D(τ,x) = 0.

Using the local coordinates and assuming that the local coordinate vectorncorresponds to the unit outward normal vector toτ, we have

D(τ,x) = 1

For the function FK we have (see [17], Section C.2.2.) F_K(s, α) :=− αu_x

with parametersp and q defined by (4.26) and A1 :=−2αq√

Singularities in the arguments of the logarithms in (4.27) and in the computation of the coefficients A_1/2, B_1/2 and G_1/2 can only occur when u_x = 0. In such case, however, we haveD(τ,x) = 0(see Remark 4.3). To conclude, for the situation whens=p∧ux ̸= 0we have v= 0.

4.4.5 Hypersingular Integral Operator

The matrix corresponding to the hypersingular integral operator is given by Dκ,h[i, j] := 1

Since the unit outward normals are constant on individual elements, we have for the first integral

and thus the matrix entries are some linear combination of corresponding entries in V_κ,h. Recall that

curl_∂Ωϕ(y) :=n(y)× ∇ϕ(y),

where for the extensions of the basis functions we can choose functions constant along the normal defined by the reference functions

ˆ

ϕ₁(ξ) := 1−ξ1−ξ2, ˆ

ϕ₂(ξ) :=ξ1, ˆ

ϕ₃(ξ) :=ξ2 for ξ = [ξ1, ξ2, ξ3]^T∈τˆ×R. Hence, it is sufficient to compute the integral

D²_κ,h[i, j] := κ²

The latter integrand is bounded and we can use the 7-point quadrature scheme. The remaining integral is similar to the single layer potential matrix entries and although the analytical formula is not given in [17], we will follow the computations in Section C.2.1

70 4 Discretization and Numerical Realization

and reuse some results. For the inner integral we have S_D(τ,x) := κ²

withF_V from (4.25) and parametersp, q given by (4.26) and I1 :=

Using per partes with u= ln

With the substitution

s=p+ q

√1 +α² sinhu, ds

du = q

√1 +α² coshu

proposed in [17], page 247, we obtain

With another substitution from [17]

v= tanhu

For the decomposition of the fractions in I3 we get (see [17], pages 247–248) 4q²

The fractions from the remaining integral can be decomposed as 4qu²_x v

72 4 Discretization and Numerical Realization Collecting all terms together we finally obtain

F_D(s, α) =F_V(s, α)

Note that the previous formula can be directly evaluated for u_x ̸= 0 corresponding to the case when x does not lie in the plane defined by τ. In the case of ux = 0 the jump property of the hypersingular operator (3.14) allows us to take the limit u_x → 0. The terms with singularity for v =±1 all vanish for u_x → 0 and it only remains to treat the term u²_xln|A₂v²+A1v+A0|. For the discriminant ofA2v²+A1v+A0 we have

corresponding to the vertex of the parabola. In the case ofαsx=tx we have K := (1 +α²)q²−(αs_x−t_x)²

(1 +α²)q−(αs_x−t_x) =|u_x| and for the limit we get

ulimx→0



u²_xln|u_x|

= 0.

Now we consider the situationαsx̸=tx and α̸= 0. For K we have

Now we consider the situationαsx̸=tx and α̸= 0. For K we have

In document 2011JanZapletal TheBoundaryElementMethodfortheHelmholtzEquationin3DMetodahraničníchprvkůproHelmholtzovurovnicive3D VŠB–TechnicalUniversityofOstravaFacultyofElectricalEngineeringandComputerScienceDepartmentofAppliedMathematics (Stránka 58-102)