Semi-analytic treatment of nearly-singular Galerkin surface integrals

Applied Numerical Mathematics 60 (2010) 974–993

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Semi-analytic treatment of nearly-singular Galerkin surface integrals ✩ S. Nintcheu Fata Computer Science and Mathematics Division, Oak Ridge National Laboratory, P.O. Box 2008, MS 6367, Oak Ridge, TN 37831-6367, USA

a r t i c l e

i n f o

Article history: Received 14 October 2009 Received in revised form 15 February 2010 Accepted 1 June 2010 Available online 11 June 2010 Keywords: Galerkin approximation Nearly-singular integrals Singular integrals Hyper-singular integrals Boundary element method Potential theory

a b s t r a c t On utilizing piecewise linear test and interpolation functions over flat triangles, three semi-analytic techniques for evaluating nearly-singular surface integrals of potential theory are developed and investigated in the context of a three-dimensional Galerkin approximation. In all procedures, the inner surface integral is calculated exactly via recursive expressions defined over the edges of the integration triangle. In addition, the proposed methods respectively employ a Duffy transformation, sinh transformations, and polynomial transformations to weaken or smooth out near singularities, followed by a standard product Gauss–Legendre quadrature rule to compute the outer boundary integral. Numerical experiments and comparisons of the proposed approaches are included to provide more insight into the performance of these techniques. Computational tests have demonstrated that the algorithm based on a Duffy transformation is the most efficient and effective method for dealing with nearly-singular as well as well-separated integrals. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The numerical solution of three-dimensional elliptic equations using the boundary element method (BEM) usually requires accurate and efficient evaluations of nearly-singular surface integrals. These integrals arise when the source point is close, in some sense, to the integration element or when the distance between the interacting boundary elements is small compared to the element characteristic lengths. Such situations frequently occur when modeling engineering problems involving thin structures [2], or simply when designing a general-purpose boundary integral equation (BIE) algorithm for three-dimensional potential problems [14]. In the context of collocation BEMs, several approaches to effectively deal with nearly-singular integrals have been proposed in the literature. Analytical techniques involving flat boundary elements can be found in [11,13]. Numerical integration methods based on various nonlinear transformations are also available, for example, from [7–9,18,21]. So far, however, the computational treatment of nearly-singular integrals stemming from Galerkin approximations of BIEs has been given very limited attention. Their computation has been impeded, perhaps, by the existence of additional difficulties associated with the calculation of outer surface integrals. Indeed, even with the use of exact formulae [13] to evaluate the inner boundary integrals, direct application of standard cubature rules to estimate the outer surface integrals are often ineffective owing to the rapid variations of the integrands. In view of these challenges, existing techniques developed for nearly-singular integrals arising in collocation BEM have not yet been generalized or extended to the Galerkin discretization of BIEs. As it stands, few exceptions are the variable order composite quadrature method [17] and the polynomial transformation procedure [18] that have been extended to Galerkin BEM, respectively, in [16] and [19].

✩ The U.S. Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes. E-mail address: [email protected].

0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.06.003

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

975

To remedy this situation, three semi-analytic methods for evaluating nearly-singular surface integrals of potential theory are developed and examined within the framework of a three-dimensional Galerkin approximation. By use of piecewise linear shape functions defined over flat triangles and the recursive expressions explicitly provided in [13], the inner surface integral is carried out exactly. Moreover, the proposed approaches are respectively based on a Duffy transformation, sinh transformations, and polynomial transformations to reduce or smooth out near singularities before applying a standard product Gauss–Legendre quadrature rule to approximate the outer boundary integral. Details of the foregoing semi-analytic analysis are given and clarified. In addition, explicit formulae for various surface potentials relevant to the treatment of inner boundary integrals are also provided. Numerical experiments and comparisons of the proposed techniques are included to illustrate their performances. 2. Problem formulation Let E p and E q represent two arbitrarily oriented open and flat triangles in R3 such that E p ∩ E q = ∅. Here, E p stands for the closure of E p . Consider the numerical treatment of the following Galerkin integrals



pq

Gi j =

  ψ i ( x) G (x, y )ψ j ( y ) dΓ y dΓx ,

Ep

pq



Hi j =

Eq

Ep

  ψ i ( x) H (x, y ) · n( y )ψ j ( y ) dΓ y dΓx ,

(1)

Eq

and pq Ti j



 =

ψi (x)n(x) · Ep

 T (x, y ) · n( y ) ψ j ( y ) dΓ y

dΓx ,

(2)

Eq

where ψi is a linear test function and ψ j is a linear shape function; the weakly-singular kernel G, the strongly-singular kernel H and the hyper-singular kernel T are given respectively by

G ( x, y ) =

1

1

4π  x − y 

,

H ( x, y ) =

1

x− y

, 4π x − y 3

T ( x, y ) =

1 4π



I2

x − y 3

−3

(x − y ) ⊗ (x − y ) x − y 5

 (3)

with x = y and x, y ∈ R3 ; n(x) and n( y ) are unit normal vectors on E p and E q respectively. Moreover, i and j are indices associated with the vertices of E p and E q respectively. In the kernel expressions (3), I2 is the symmetric, second-order identity tensor on R3 ;  ·  used throughout this article, unless otherwise specified, is the usual Euclidean norm in R3 ; pq pq x ∈ R3 is the source point and y ∈ R3 is the field or receiver point. The integrals Gi j and Hi j are respectively the local contributions of the flux and potential influence matrices when solving the Laplace equation via a singular Galerkin BIE pq method [12]. On the other hand, Ti j simply represents the local contribution of the potential influence matrix when solving the Laplace equation via a hyper-singular Galerkin BIE approach [14]. With the assumption that the triangles E p and E q do not intersect, the contributions (1) and (2) are not singular. However, experience has demonstrated that the integrands of the Galerkin integrals vary rapidly as the distance between E p and E q decreases. In such situation, straightforward application of standard cubature methods to evaluate these integrals can become ineffective. For this reason, the computational treatment of the non-singular Galerkin integrals (1) and (2) should be decomposed into a nearly-singular (or quasi-singular) case, and a well-separated case. To aid the foregoing analysis, let h p denote the diameter of an arbitrary, open and flat triangle E p in R3 defined as

h p = sup x − y .

(4)

x, y ∈ E p

In addition, let dist( E p , E q ) be the distance between two open flat triangles E p and E q in R3 specified as





dist( E p , E q ) = min x − y : x ∈ E p , y ∈ E q ,

(5)

and introduce the near-singular ratio D pq as

D pq =

dist( E p , E q ) max{h p , hq }

.

(6)

With these settings, a non-singular Galerkin integral is said to be nearly singular if 0 < D pq < 1, and well separated if D pq  1. The well-separated contributions can be accurately calculated to any precision by use of a standard product Gauss–Legendre quadrature rule. However, the accuracy of singular and nearly-singular integrals is essential to obtain a meaningful solution of a typical boundary-value problem via a Galerkin BEM. Since the treatment of singular integrals has been extensively studied in [4,6,14], this article will focus on the computation of nearly-singular Galerkin integrals.

976

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

Fig. 1. Local coordinate systems associated with a triangle E q .

3. Treatment of inner integrals In view of the double integrals (1) and (2), one can conveniently introduce the single-layer potential



q

G j ( x) =

G (x, y )ψ j ( y ) dΓ y ,

x ∈ R3 ,

(7)

Eq

the double-layer potential



q

H j ( x) =

H (x, y ) · n( y )ψ j ( y ) dΓ y ,

x ∈ R3 ,

(8)

x ∈ R3 ,

(9)

Eq

and the “quadruple-layer” potential q



T j ( x) =

T (x, y ) · n( y )ψ j ( y ) dΓ y , Eq

for an arbitrary and fixed source point x ∈ R3 . For a linear shape function ψ j defined on a flat triangle E q , it had been demonstrated in [13] that surface potentials (7)–(9) admit exact representations in terms of analytic formulae over the sides of E q . These explicit expressions will be employed to compute the inner integrals of (1) and (2). 3.1. Analytic formulae for surface potentials To provide exact expressions for surface potentials (7)–(9), let x ∈ R3 represent an arbitrary and fixed source point in space. In addition, assume that the flat triangle E q is also fixed in the space R3 , and let L = ∂ E q . With reference to Fig. 1(a), let y 1 , y 2 , y 3 ∈ R3 be the position vectors of the vertices of E q and introduce the oriented segments L 1 = [ y 1 , y 2 ], L 2 = [ y 2 , y 3 ], L 3 = [ y 3 , y 1 ] so that the boundary of E q is specified as L = L 1 ∪ L 2 ∪ L 3 . As usual, the orientation of a triangle is given by the traversal of its vertices. Further, let {x; e 1 , e 2 , e 3 } be a local orthonormal companion reference of R3 such that (i) the unit vector e 1 is parallel to L 1 and points in the direction of L 1 , (ii) the unit vector e 2 is selected on the plane of E q as dictated by the orientation of E q (see Fig. 1(a)), and (iii) e 3 = e 1 × e 2 . In other words, the orthonormal vectors e 1 and e 2 can be prescribed as

e1 =

y2 − y1

 y2

−

y1

,

e2 =

( y3 − y1) − σ ( y2 − y1) , ( y 3 − y 1 ) − σ ( y 2 − y 1 )

σ=

( y3 − y1) · ( y2 − y1) .  y 2 − y 1 2

Since E q is a flat triangle, the vectors e i (i = 1, 2, 3) are constant unit vectors associated with E q . In particular, e 3 = n( y ), y ∈ E q . Now, let P E q be the plane of E q spanned by {e 1 , e 2 }. Furthermore, let x ∈ P E q be the orthogonal projection of the

source point x on P E q in the direction e 3 . In this setting, {x; e 1 , e 2 } is the induced reference on P E q .

These definitions introduce in R3 a local coordinate system {x; ξ, ζ, η˜ } associated with the triangle E q . In {x; e 1 , e 2 , e 3 }, the position vector r of an arbitrary field point y ∈ P E q can be expressed as

r = y − x = ξ e 1 + ζ e 2 + η˜ e 3 .

(10) j

j

More specifically, the position vectors r of the respective vertices y of E q are decomposed as

r j = y j − x = ξ j e 1 + ζ j e 2 + η˜ e 3 ,

j = 1, 2, 3.

(11)

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

977

By construction, ζ1 = ζ2 . Moreover, the parameter η˜ represents the relative distance between the source point x ∈ R3 and the plane of E q . Indeed, it is readily seen from (10) that x −x = η˜ e 3 . In {x; ξ, ζ, η˜ }, the distance between the field point

y ∈ P E q and the source point x ∈ R3 admits the representation

r =  y − x =



ξ 2 + ζ 2 + η˜ 2 .

(12)

Next, it can be shown that, in {x; ξ, ζ, η˜ }, the linear shape functions ψ j defined on E q take the form

ψ j ( y ) ≡ ψ j (ξ, ζ ) = a j ξ + b j ζ + c j ,

j = 1, 2, 3, y ∈ E q ,

(13)

where the coefficients a j and b j ( j = 1, 2, 3) do not depend on the source point x ∈ R3 , and are expressed as

a1 = −

1

ξ2 − ξ1

a2 =

,

1

ξ2 − ξ1

a 3 = 0;

,

b1 =

ξ3 − ξ2 , Δ

b2 =

ξ1 − ξ3 , Δ

b3 =

1

ζ3 − ζ 1

with Δ = (ξ2 − ξ1 )(ζ3 − ζ1 ); the free terms c j ( j = 1, 2, 3) depend on x, i.e. c j = c j (x), and are given as

c1 =

ξ 2 ζ3 − ξ 3 ζ2 , Δ

c2 =

ξ 3 ζ1 − ξ 1 ζ3 , Δ

c3 = −

ζ1 . ζ3 − ζ 1

(14)

Now, by use of (3), (10), (12) and (13), it can be shown that the single-layer potential (7) is expressible in the local coordinate system {x; ξ, ζ, η˜ } as

1  ξ ζ a j I1 + b j I1 + c j I1 , 4π

q

G j ( x) = where

 I1 =

1 r



ξ

ξ

I1 =

ds,



ζ

Eq

(15)

ζ

I1 =

ds,

r

Eq

x ∈ R3 ,

r

ds.

(16)

Eq

Similarly, the double-layer potential (8) can be written in {x; ξ, ζ, η˜ } as q

H j ( x) = where

 I3 =

−η˜ 

1

ξ



ξ

I3 =

ds,

r3



ζ

x ∈ R3 ,

a j I3 + b j I3 + c j I3 ,

4π

Eq

ξ r3



ζ

ζ

I3 =

ds,

Eq

(17)

ds.

r3

(18)

Eq

Finally, the “quadruple-layer” potential (9) can be decomposed in the local reference {x; e 1 , e 2 , e 3 } as

1 

q

T j ( x) =

4π

q

q



q

x ∈ R3 ,

T j1 (x)e 1 + T j2 (x)e 2 + T j3 (x)e 3 ,

(19)

where



q

ξξ

ξ

ξζ

q

ξ

ξ



q

T j1 (x) = −3η˜ a j I 5 + b j I 5 + c j I 5 ,

ξζ

ζ

ζ



T j3 (x) = a j I 3 − 3η˜ 2 I 5 + b j I 3 − 3η˜ 2 I 5 + c j I 3 − 3η˜ 2 I 5 , ξ

ζζ

ζ

T j2 (x) = −3η˜ a j I 5 + b j I 5 + c j I 5 , x ∈ R3 ,

x ∈ R3 , (20)

ζ

employing the integrals I 3 , I 3 , I 3 given in (18) together with the following integrals



I5 =

1 r5

ds,



ξ

I5 =

Eq

ξξ



I5 =

ξ

r5 Eq

r5 Eq

2

ds,

ξζ

ξ 

I5 =



ζ r5

Eq

ξζ r5

Eq

ζ

I5 =

ds,

ds,

ζζ



I5 =

ds,

ζ2 r5

ds.

(21)

Eq

The dependence of the generic integrals (16), (18) and (21) on x ∈ R3 is omitted here for brevity. According to [13], there exist three local coordinate systems with origin at the source point x ∈ R3 in which the generic integrals (16), (18) and (21) are expressed as analytic formulae over the edges of E q . In order to utilize these exact expressions, let θi ∈ (0, π ) be the inclusion angle at vertex y i (i = 1, 2, 3) of E q such that θ1 + θ2 + θ3 = π (see Fig. 1(a)), and set α1 ≡ 0, α2 = π − θ2 , α3 = π + θ1 . By use of these angles αi (i = 1, 2, 3), one can introduce three local orthonormal companion references of R3 with origin at the source point x, {x; e p i , eqi , e 3 }, such that

978

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

(i) e p i is parallel to L i (i = 1, 2, 3) and points in the direction of L i , (ii) the unit vector eqi ∈ P E q is specified as dictated by the orientation of E q , and (iii) e 3 = e p i × eqi . References {x; e p i , eqi , e 3 } induce respectively in R3 three local coordinate systems {x; p i , qi , η˜ } associated with each side L i (i = 1, 2, 3) of E q such that p i = ξ cos αi + ζ sin αi , and qi = −ξ sin αi + ζ cos αi . In addition, {x; p i , qi } are local coordinate

systems on P E q . j

j

Now, let ( p i , qi , η˜ ) be the coordinate in {x; e p i , eqi , e 3 } of the vector r j = y j − x. With these specifications and in view of (11), the distance between the source point x ∈ R3 and the j-th vertex y j of E q can be effectively computed as



ρj =

j 2

j 2 + qi + η˜ 2 ,

pi

i , j = 1, 2, 3.

(22)

By construction, in the local coordinate systems {x; p i , qi , η˜ } defined for each edge L i (i = 1, 2, 3) of E q , one can write q11 = q21 ≡ q1 , q22 = q32 ≡ q2 , q33 = q13 ≡ q3 . Further, it is useful to introduce the ensuing generic functions for every oriented edge L i = [ y i , y i +1 ] of E q as









χi (x) = ln p ii + ρi − ln p ii+1 + ρi+1 , 

γi (x) = arctan



˜ i −2p ii qi ηρ (qi )2 (ρi )2 − ( p ii )2 η˜ 2

δ i ( x) =  − arctan

p ii

ρi

−

p ii +1

ρ i +1

L i ( x) =

,

˜ i +1 −2p ii +1 qi ηρ

(qi )2 (ρi +1 )2 − ( p ii +1 )2 η˜ 2

1

−

ρi

 ,

1

ρ i +1

,

ρ˜i (x) = ρi − ρi+1 , di = (qi )2 + η˜ 2 , i = 1, 2, 3,

(23) j

j

where the subscript or superscript “4” should be replaced by “1”. Note that p i = p i (x), qi = qi (x), η˜ = η˜ (x), di = di (x). This dependence of Finally, define the function

j p i , qi ,

ρ i = ρ i ( x) , η˜ , ρi , di on the source point x ∈ R is not included here and in the sequel for brevity. 3

 3 1

θ(x) =

2 1 2

γk (x) + θo (x), if η˜ > 0 ˜ 0 k=1 γk (x) − θo (x), if η k =1

3

(24)

when evaluating the fluxes and potentials on the boundary Γ , where

θo (x) =

⎧ 0, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ θi ,

if x ∈ P E q \ E q ,

if x = y i ,

i = 1, 2, 3,

⎪ π , if x ∈ ∂ E q \ { y 1 , y 2 , y 3 }, ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 2π , if x ∈ E q .

On the basis of these definitions, one can write explicit expressions for all integrals featured in (16), (18) and (21). Details on the derivation of these recursive formulae can be found in [13]. With these exact formulae and the relation (15), one can express the single-layer potential (7) as



q G j ( x)



3 3   1  1 = q˜ k ρk (x)ˆakj − dk χk (x)bˆ kj + c j ( x) qk χk (x) − η˜ θ(x) , 8π 4π k =1

x ∈ R3 ,

(25)

k =1

where aˆ kj = a j cos αk + b j sin αk , and bˆ kj = a j sin αk − b j cos αk (k, j = 1, 2, 3). In an analogous manner, it can be shown using (17) and the exact formulae for the generic integrals (18) that the double-layer potential (8) take the form q

H j ( x) = −

1 4π

η˜

3 

χk (x)bˆ kj −

k =1

1 4π

c j (x)θ(x),

x ∈ R3 .

(26) q

Lastly, with the aid of (20) and the analytic expressions for (21), one can verify that the components T jl of the “quadruplelayer” potential T q

T j1 (x) = η˜

q j

given in (19) read

3  k =1

 sin αk Lk (x)ˆakj +

qk dk

 3  δk (x) δk (x)bˆ kj − η˜ c j (x) sin αk − a j θ(x), k =1

dk

x ∈ R3 ,

(27)

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

q T j2 (x)

 3  3   qk δk (x) = η˜ Lk (x) sin αk − δk (x) cos αk bˆ kj + η˜ c j (x) cos αk − b j θ(x), dk

k =1

q

T j3 (x) =

3  

χk (x) −

k =1

η˜ 2 dk

 3  qk δk (x) bˆ kj − c j (x) δk (x), k =1

dk

k =1

dk

979

x ∈ R3 ,

x ∈ R3 .

(28)

(29)

In the formulae (24)–(29), the summation is carried out over the sides L k (k = 1, 2, 3) of E q . 4. Treatment of nearly-singular Galerkin integrals With the analytic expressions for the surface potentials (7)–(9), one can now provide a comprehensive examination of the Galerkin integrals (1) and (2) over a pair of non-intersecting and flat triangles E p and E q . However, evaluations of the outer integrals will be carried out using suitable numerical techniques. To describe the algorithm, it is important to utilize (7) and (8), and to rewrite the Galerkin integrals (1) as pq



Gi j =

q

pq



q

Hi j =

ψi (x)G j (x) dΓx , Ep

ψi (x) H j (x) dΓx ,

i , j = 1, 2, 3,

(30)

Ep

q G j ( x)

q H j (x),

x ∈ E p are given by (25) and (26); i and j are local indices associated with the vertices of the flat where and triangles E p and E q respectively. Likewise, to express (2) as an integral over E p , first assume that the flat triangle E p is fixed in R3 , and let L = ∂ E p . With this assumption, let x1 , x2 , x3 ∈ R3 be the position vectors of the vertices of E p and define the oriented segments L 1 = [x1 , x2 ], L 2 = [x2 , x3 ], L 3 = [x3 , x1 ] so that the boundary of E p is given by L = L 1 ∪ L 2 ∪ L 3 . With these settings, it is useful to construct on E p a local orthonormal companion basis {e 1 , e 2 , e 3 } of R3 such that (i) the unit vector (ii) the unit vector (iii) e 3 = e 1 × e 2 .

e 1 is parallel to L 1 and points in the direction of L 1 , e 2 is prescribed on the plane of E p as dictated by the orientation of E p , and

Owing to the fact that E p is a flat triangle, the basis vectors e i (i = 1, 2, 3) are constant unit vectors associated with E p . Hence, e 3 = n(x), x ∈ E p . On employing the orthonormal basis {e 1 , e 2 , e 3 } associated with E q (see Section 3.1), and the basis {e 1 , e 2 , e 3 } related to E p , one can introduce in R3 an orthogonal matrix C ∈ R3×3 that characterizes the transition from the companion basis {e 1 , e 2 , e 3 } of E q to the companion basis {e 1 , e 2 , e 3 } of E p . The components C i j of C ∈ R3×3 are given as

e j = C 1 j e 1 + C 2 j e 2 + C 3 j e 3 ,

j = 1, 2, 3.

(31)

By virtue of these definitions and taking into account (9), (19) and (31), the Galerkin integral (2) can be rewritten as pq

Ti j = where pq

1  4π



T i jl =

pq

pq

pq



T i j1 C 13 + T i j2 C 23 + T i j3 C 33 ,

q

ψi (x) T jl (x) dΓx ,

(32)

i , j , l = 1, 2, 3;

(33)

Ep q

as before, i and j are local indices associated with the vertices of E p and E q respectively; T jl (x), x ∈ E p are provided by (27)–(29). Thus, the computation of

pq Ti j

has now been reduced to that of (33).

4.1. Preliminaries To design an effective treatment for the Galerkin integrals (30) and (33), it is desirable to incorporate into the algorithm an estimate of the distance between the integration triangles E p and E q . However, as defined in (5), computing the distance x ∈ E p and dist( E p , E q ) between two triangles E p and E q is non-trivial and can be very expensive since the closest points   y ∈ E q for which dist( E p , E q ) =  x− y  may not be unique. Fortunately, for flat triangles, the closest points can always be realized in such a way that one point lies on the boundary of one of the triangles [5]. With this in mind, let for instance  y be located at a vertex y j ( j = 1, 2, 3) of E q : x is an interior point of E p (i.e.  x ∈ E p ), then E p is conveniently partitioned into three sub-triangles given by the (i) if  x, x1 , x2 }, { x, x2 , x3 }, and { x, x3 , x1 } (see Fig. 2(a)); following oriented triplets { (ii) if  x is on an edge of E p (i.e.  x ∈ ∂ E p \ {x1 , x2 , x3 }), then E p is decomposed into two sub-triangles each containing not x as a vertex but also a vertex of E p that is the intersection point of the remaining edges of E p (see Fig. 2(b)). only 

980

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

Fig. 2. Closest points  x ∈ E p and  y ∈ E q between two flat triangles E p and E q : (a)  x ∈ E p and  y = y 1 leading to the decomposition of E p into three sub-triangles; (b)  x is on the edge (x1 , x2 ) of E p and  y = y 1 resulting in the subdivision of E p into two sub-triangles with  x and x3 as vertices.

On the other hand, when  y is situated on an edge of E q (i.e.  y ∈ ∂ E q \ { y 1 , y 2 , y 3 }), the triangle E q is first subdivided into y as a vertex and a vertex of E q residing on the opposite side to the designated edge; next, the two sub-triangles with  x ∈ ∂ E p . Finally, analogous decomposition of E p and E q should be splitting of E p can proceed depending on the position of  carried out when one of the closest point lies on the boundary of E p instead. With these observations, one can assume without loss of generality that the disjoint triangles E p and E q are oriented in R3 such that the closest point  x coincides with vertex x1 of E p (i.e.  x = x1 ), and the closest point  y is positioned at vertex y = y 1 ). With this assumption, one can write y 1 of E q (i.e., 





ε = dist( E p , E q ) = x1 − y 1 ,

(34)

where ε is employed here and in the sequel for brevity. Since E p ∩ E q = ∅, ε > 0. In view of (34), it will be useful to introduce the vector ε = x1 − y 1 . For simplicity in notation, let h = max{h p , hq } and rewrite the near-singular ratio (6) simply as

D=

ε h

(35)

,

where h p and hq represent respectively the diameters of E p and E q as defined by (4). Now, consider the coordinate system { y 1 ; ξ, ζ, η} associated with E q such that the origin is situated at vertex y 1 of E q and the ξ -axis is aligned with L 1 as depicted in Fig. 1(b). In this setting, { y 1 ; e 1 , e 2 , e 3 } is an orthonormal companion reference of R3 and { y 1 ; e 1 , e 2 } is the induced reference on the plane of E q denoted as P E q . In { y 1 ; e 1 , e 2 , e 3 }, the position vector R of an arbitrary point x ∈ R3 can be evaluated as

R = x − y 1 = ξ e1 + ζ e2 + ηe3 .

(36)

In particular, the position vectors R j of the vertices y j ( j = 2, 3) of E q can be specified as

R 2 = y 2 − y 1 = be 1 ,

R 3 = y 3 − y 1 = ce 1 + ae 2 ,

(37)

where R b > 0 is the base of the triangle E q or the length of the edge L 1 ; the parameter R a > 0 is the height of E q in { y 1 ; e 1 , e 2 }, and c ∈ R is the relative position of vertex y 3 in the oriented triplet { y 1 , y 2 , y 3 }. Likewise, introduce a local coordinate system {x1 ; ξ , ζ , η } associated with E p such that the origin is located at vertex x1 of E p and the ξ -axis is aligned with L 1 . In the induced companion reference {x1 ; e 1 , e 2 , e 3 }, the position vector R of an arbitrary point x ∈ R3 can be calculated as

R = x − x1 = ξ e 1 + ζ e 2 + η e 3 .

(38)

More specifically, the position vectors R j of the vertices x j ( j

R 2 = x2 − x1 = b e 1 ,

= 2, 3) of E p can be prescribed as

R 3 = x3 − x1 = c e 1 + a e 2 ,

(39)

where R b > 0 is the base of the triangle E p or the length of the edge L 1 ; the parameter R a > 0 is the height of E p in {x1 ; e 1 , e 2 }, and c ∈ R is the relative position of vertex x3 in the oriented triplet {x1 , x2 , x3 }. By construction, the triangle E lies in the ξ ζ -plane (i.e. the plane of E referred to as P ) that is specified by η = 0. p

p

Ep

With these definitions, it can be shown using (11), (36), (37) and the translation from {x; e 1 , e 2 , e 3 } to the fixed reference { y 1 ; e 1 , e 2 , e 3 } that

p 11 ≡ ξ1 = −ξ,

q1 ≡ ζ1 = −ζ ;

ξ2 = b − ξ,

ζ2 = −ζ ;

ξ3 = c − ξ,

ζ3 = a − ζ ;

η˜ = −η. (40)

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

981

Now, (40) can be employed in (14) to express the functions c j ( j = 1, 2, 3) in { y 1 ; ξ, ζ, η} as

c 1 (ξ, ζ ) = −

ξ b

+

c−b ab

ζ + 1,

c 2 (ξ, ζ ) =

ξ b

−

c ab

ζ,

c 3 (ξ, ζ ) =

ζ a

(41)

.

In addition, the shape functions ψ j ( j = 1, 2, 3) introduced in (13) can be given directly in {x1 ; ξ , ζ , η } as

 ξ c − b ψ1 ξ , ζ = − + ζ + 1, b

Moreover, the vector in { y 1 ; e 1 , e 2 , e 3 } as

 ξ c ψ2 ξ , ζ = − ζ ,

ab

b

ab

 ζ ψ3 ξ , ζ = .

(42)

a

ε = x1 − y 1 , characterizing the distance between E p and E q as expressed in (34), can be decomposed

ε = x1 − y 1 = ε1 e1 + ε2 e2 + ε3 e3 .

(43)

For x ∈ P E p , it can be established using (31), (36), (38), (43), the translation R =

ξ = C 11 ξ + C 12 ζ + ε1 ,

ζ = C 21 ξ + C 22 ζ + ε2 ,

R

+ε

and the condition η

η = C 31 ξ + C 32 ζ + ε3 .

= 0 that (44)

Formula (44) will be used to convert any function defined on the triangle E p in {ξ, ζ, η}-coordinate into a function expressed in {ξ , ζ }-local coordinate on E p . Remark 1. With these preliminary settings, it is now clear that the treatment of nearly-singular integrals within a Galerkin framework can be reduced to the treatment of near singularities concentrated at a pair of closest points between the integration domains. In addition, note that the procedure presented herein automatically degenerates to the singular study for the vertex-adjacent case (see [14]) when the distance between the integration triangles ε → 0. In this situation, (34) and (43) reveal that y 1 = x1 and that ε1 = ε2 = ε3 = 0. 5. Mappings and integrations for nearly-singular Galerkin integrals To effectively deal with quasi-singular Galerkin integrals, the key idea is to (i) introduce various mappings or transformations to weaken or remove concentrated near singularities at vertex x1 of E p and at vertex y 1 of E q , and (ii) evaluate the transformed integrals via a standard product Gauss–Legendre quadrature rule. To this end, it is useful to employ the Duffy transformation [3] as the basis to all other mappings that will be investigated in this study. 5.1. Duffy transformation In the local coordinate system {x1 ; ξ , ζ } defined on the plane of E p , the following Duffy transformation

 

ξ = x b + c − b y ;

ζ = a xy

(45)

with 0  x, y  1 is first considered. The degenerate transformation (45) maps the standard square {(x, y ) ∈ R2 : 0 < x, y < 1} onto E p so that x, y can be used as dimensionless parameters on E p . In addition, (45) together with its Jacobian a b x reduce or remove any singularity at (x = 0, y = 0) which corresponds to vertex x1 of E p . Now, let f = f (x) ≡ f (ξ , ζ ) be an arbitrary function prescribed on E p . With the mapping (45), a typical Galerkin integral can be transformed as

 Fij =



1 1

ψi (x)c j (x) f (x) dΓx = a b Ep

xψi (x, y )c j (x, y ) g (x, y ) dx dy , 0

i , j = 1, 2, 3,

(46)

0

where g (x, y ) = f (x(b + (c − b ) y ), a xy ); c j (x, y ) are given by (41), (44) and (45); the shape functions ψi (x, y ) are provided by (42) and (45). The transformed double integral in (46) can now be accurately calculated by use of a standard product Gauss–Legendre quadrature rule. Remark 2. Since the Duffy transformation is the foundation of the numerical treatment presented herein, below are some additional derivations that are necessary for any subsequent algorithm further relying on a change of variable. Based on the mapping (45), let ξˆ = b + (c − b ) y and ζˆ = a y, and introduce the quantities

pˆ 11 = −C 11 ξˆ − C 12 ζˆ ,

qˆ 1 = −C 21 ξˆ − C 22 ζˆ ,

ηˆ = −C 31 ξˆ − C 32 ζˆ .

(47)

982

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

Moreover, since

e=

ε > 0, one can make use of (34) and (43) to define the unit direction vector

x − y1 1

 x1 − y 1 

= εˆ 1 e 1 + εˆ 2 e 2 + εˆ 3 e 3 .

(48)

Next, by use of (40), (44), (45), (47) and noticing that

p 11 = x pˆ 11 − ε εˆ 1 ,

q1 = xqˆ 1 − ε εˆ 2 ,

εi = εεˆ i (i = 1, 2, 3) given by (43) and (48), one can write

η˜ = xηˆ − εεˆ 3 .

Since near singularities are concentrated at vertex



(49)

y1

E q , it is important to use (22) and rewrite

of

+ (q1 )2 + η˜ 2 , a quantity representing the distance from a source point on E p to vertex y 1 of E q , in terms of the Duffy variables x and y. Indeed, one can utilize (49), the condition that e  = 1 obtained from (48) and the fact that ek · el = δkl (k, l = 1, 2) from (31) to show that

ρ1 =

( p 11 )2



ρ1 = ρˆ1 (x − A )2 + B 2 ,

(50)

ˇ, A = εA

(51)

where

B = ε Bˇ ,

and



2

2

ρˆ1 = (ξˆ ) + (ζˆ ) ,

bˆ 1 = εˆ 1 pˆ 11 + εˆ 2 qˆ 1 + εˆ 3 ηˆ ,

ˇ= A

bˆ 1

(ρˆ1

)2

,

Bˇ =

(ρˆ1 )2 − (bˆ 1 )2 (ρˆ1 )2

.

(52)

Here δkl is the Kronecker delta. Recall that ξˆ = ξˆ ( y ), ζˆ = ζˆ ( y ), pˆ 11 = pˆ 11 ( y ), qˆ 1 = qˆ 1 ( y ), ηˆ = ηˆ ( y ). For this reason, ρˆ1 = ρˆ1 ( y ),

ˇ = Aˇ ( y ), Bˇ = Bˇ ( y ). The dependence of all these quantities on the variable y is omitted in this study for brevity. bˆ 1 = bˆ 1 ( y ), A By construction, ρ1  ε . Furthermore, note that ρˆ1 > 0 and that

ρˆ1  |bˆ 1 | for all y ∈ [0, 1].

(53)

Indeed, (53) can easily be verified by (i) recognizing that bˆ 1 = e ·( pˆ 11 e 1 + qˆ 1 e 2 + ηˆ e 3 ), (ii) applying the Cauchy–Bunyakovsky–

( pˆ 11 )2 + (ˆq1 )2 + ηˆ 2 . It then follows from (53) that B  0 for all y ∈ [0, 1]. The condition (53) equivalently expresses the fact that ρ1  0 for all x, y ∈ [0, 1].

Schwarz inequality, and (iii) noticing that ρˆ1 =

5.2. Sinh transformation To successfully integrate expressions containing arguments of the form (x − A )2 + B 2 , it is suggested in [10] to employ a sinh transformation with the feature that its Jacobian be proportional to the potentially small parameter ε . Moreover, according to [10], the parameter A can be viewed as the position of the nearly-singular point on the integration interval and B can be interpreted as the distance from a source point to the point A. In view of (53), one must therefore consider two cases: (i) ρˆ1 > |bˆ 1 |. This situation corresponds to the case where B > 0 for all y ∈ [0, 1]. Thus, a nonlinear transformation

x = A + B sinh( F 1 z − F 2 ),

z ∈ [0, 1]

(54)

is introduced whereby F 1 and F 2 are determined under the requirement that (54) maps the interval [0, 1] onto [0, 1]. More specifically, by solving the system A − B sinh( F 2 ) = 0 and A + B sinh( F 1 − F 2 ) = 1, and taking into account (51) and (52), one can write

 F 2 = ln

bˆ 1 + ρˆ1

(ρˆ1 )2 − (bˆ 1 )2



 (ρˆ1 )2 + ρˆ1 (ρˆ1 )2 − 2εbˆ 1 + ε 2 − εbˆ 1

F 1 = F 2 + ln . ε (ρˆ1 )2 − (bˆ 1 )2 

,

In view of (51), the Jacobian of transformation (54) is calculated as

dx dz

= ε Bˇ F 1 cosh( F 1 z − F 2 ),

z ∈ [0, 1].

(55)

It is now clear from (55) that the Jacobian of (54) is proportional to ε . In addition, it can be shown using (52), (54) in (50) that

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

ρ1 = ε Thus,

(ρˆ1 )2 − (bˆ 1 )2

ρˆ1

cosh( F 1 z − F 2 ).

983

(56)

ρ1 is also proportional to ε .

(ii) ρˆ1 = |bˆ 1 |. In this situation, B = 0. Consequently, the following nonlinear transformation

x= A+

ε sinh( J 1 z − J 2 ), z ∈ [0, 1] ρˆ1

(57)

is used to carry out the integration. Here, the parameters J 1 and J 2 are chosen such that (57) maps [0, 1] onto itself. Using this argument,

J 2 = ln

ˆ

b1

ρˆ1

+

√









ρˆ1 bˆ 1 ρˆ1 bˆ 1 J 1 = J 2 + ln − + 1+ − ε ρˆ1 ε ρˆ1

2 ,

2  .

The Jacobian of transformation (57) is evaluated as

dx dz

=ε

J1

ρˆ1

cosh( J 1 z − J 2 ),

z ∈ [0, 1].

(58)

It is evident from (58) that the Jacobian of (57) is proportional to ε . Moreover, with the aid of the condition B = 0 and the mapping (57) in (50), it can be established that





ρ1 = ε sinh( J 1 z − J 2 ).

(59)

Owing to (51), it is convenient to rewrite the transformations (54) and (57) simply as

x = ε xˇ ,

(60)

where xˇ is indeed a function of z ∈ [0, 1] and the parameter y ∈ [0, 1] explicitly specified as



xˇ = xˇ ( z, y ) =

if ρˆ1 > |bˆ 1 |, ˇA ( y ) + sinh( J 1 ( y ) z − J 2 ( y ))/ρˆ1 ( y ), if ρˆ1 = |bˆ 1 |.

ˇ ( y ) + Bˇ ( y ) sinh( F 1 ( y ) z − F 2 ( y )), A

Similarly, on inspecting (55), (56) and (58), (59), one can express dx and

dx = ε Kˇ dz,

ρ1 as

ρ1 = ερˇ1 ,

(61)

where Kˇ = Kˇ ( z, y ) and ρˇ1 = ρˇ1 ( z, y ) are defined appropriately for both cases ρˆ1 > |bˆ 1 | and ρˆ1 = |bˆ 1 |. Numerical experiments and the study conducted in [9] have indicated that a sinh transformation should also be used for the Duffy variable y. With these observations, the following change of variable

y = sinh( D 1 t ),

t ∈ [0, 1]

is adopted for the Duffy parameter y ∈ [0, 1]. The coefficient D 1 = ln[1 + maps [0, 1] onto [0, 1]. In view of (62), one can write

√

(62) 2 ] is selected such that the transformation (62)

d y = D 1 cosh( D 1 t ) dt .

(63)

By use of (60) through (63), the Galerkin integral (46) can be converted to

F i j = a b D1ε

1 2

1  cosh( D 1 t ) x¯ ( z, t ) K¯ ( z, t )ψ¯ i ( z, t )¯c j ( z, t ) g¯ ( z, t ) dz dt ,

0

i , j = 1, 2, 3,

(64)

0

where a substitution y = y (t ) provided by (62) has been utilized; x¯ ( z, t ) = xˇ ( z, y (t )), K¯ ( z, t ) = Kˇ ( z, y (t )); the transformed shape functions ψ¯ i ( z, t ) = ψi (ε x¯ ( z, t ), y (t )); similarly c¯ j ( z, t ) = c j (ε x¯ ( z, t ), y (t )) and g¯ ( z, t ) = g (ε x¯ ( z, t ), y (t )). The factor ε 2 in (64) serves to regularize the integral prior to application of a standard product Gauss–Legendre quadrature rule. To illustrate this feature, it is useful to employ (60) and to rewrite (49) as

p 11 = ε pˇ 11 ,

q1 = εqˇ 1 ,

η˜ = εηˇ ,

(65)

984

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 3. Relative errors for Gi j ≡ Gi j on parallel triangles with

ε = 10−4 and D = 7.071 × 10−5 .

where pˇ 11 = xˇ pˆ 11 − εˆ 1 , qˇ 1 = xˇ qˆ 1 − εˆ 2 and ηˇ = xˇ ηˆ − εˆ 3 . In addition, it can also be shown that

p 13 = ε pˇ 13 ,

d1 = ε 2 dˇ 1 ,

q3 = εqˇ 3 ,

d3 = ε 2 dˇ 3 .

(66)

where pˇ 13 = pˇ 11 cos α3 + qˇ 1 sin α3 and qˇ 3 = − pˇ 11 sin α3 + qˇ 1 cos α3 ; dˇ i = (ˇqi )2 + (ηˇ )2 (i = 1, 3) with α3 = π + θ1 as introduced in Section 3.1. By use of (60)–(63), (65) and (66), a Galerkin integral obtained from (29) and (33) with near-singular terms such as q1 /d1 and q3 /d3 can be transformed as



  q1 q3 ψi (x)c j (x) δ1 (x) + δ3 (x) dΓx d1

Ep



1

= a b D1ε 0

d3

1   1     pˇ 13 p 21 qˇ 1 pˇ 1 qˇ 3 p 33 ¯ ¯ ¯ cosh( D 1 t ) x( z, t ) K ( z, t )ψi ( z, t )¯c j ( z, t ) − − dz dt . + ρ2 ρˇ1 dˇ 1 ρˇ1 dˇ 3 ρ3

(67)

0

The transformed double integral in (67) has been regularized, i.e., it no longer contains integrand terms that is explicitly divided by the potentially small parameter ε . A standard product Gauss–Legendre quadrature rule can now be used to efficiently evaluate the right-hand side of (67). Remark 3.Instead of (62), the following transformation y = D sinh( S 1 t ), t ∈ [0, 1] with D specified by (35) and S 1 = ln[1/D + 1 + 1/D 2 ] has also been examined. In the experiments conducted to validate this change of variable, no improvement in the accuracy of the non-singular Galerkin integrals has been observed over the transformation (62). 5.3. Polynomial transformation As proposed in [1,19], one can also use a polynomial transformation to carry out integrals containing expressions of the form (x − A )2 + B 2 . By virtue of (53), there are two cases to consider:

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 4. Relative errors for Gi j ≡ Gi j on parallel triangles with

985

ε = 10−2 and D = 7.071 × 10−3 .

(i) ρˆ1 > |bˆ 1 |. Hence, B > 0. Thus, one can introduce a polynomial transformation

x = A + B (C 1 z − C 2 )mx ,

z ∈ [0, 1]

(68)

where mx is an arbitrary odd natural number; R C 1 > 0 and C 2 ∈ R are numbers to be determined by requiring that (68) maps the interval [0, 1] onto [0, 1]. With the aid of (51) and (52), one can resolve the system A − B (C 2 )mx = 0 and A + B (C 1 − C 2 )mx = 1, and determine that



C2 =

bˆ 1



 m1

(ρˆ1 )2 − (bˆ 1 )2

x

,

C1 = C2 +

(ρˆ1 )2 − εbˆ 1

ε (ρˆ1 )2 − (bˆ 1 )2

 m1

x

.

By use of (68) in (50) and taking into account (51) and (52), one can show that

ρ1 = ε

(ρˆ1 )2 − (bˆ 1 )2 

ρˆ1

1 + (C 1 z − C 2 )2mx .

(69)

Moreover, it directly follows from (51) and (68) that

ˇ 1 (C 1 z − C 2 )mx −1 dz. dx = mx ε BC

(70)

(ii) ρˆ1 = |bˆ 1 |. In this case, B = 0. Thus, the following polynomial transformation

x= A+

ε ( M 1 z − M 2 )mx , z ∈ [0, 1] ρˆ1

(71)

986

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 5. Relative errors for Hi j ≡ Hi j on parallel triangles with

ε = 10−4 and D = 7.071 × 10−5 .

is employed to carry out the integration. Here, the real parameters M 1 > 0 and M 2 are chosen such that (71) maps [0, 1] onto itself. Using this argument,

M2 =

bˆ 1

ρˆ1



,

ρˆ1 bˆ 1 M1 = M2 + − ε ρˆ1

 m1

x

.

Next, by use of (71) in (50) and taking into account the fact that B = 0, one can write

ρ1 = ε| M 1 z − M 2 |mx .

(72)

Further, one can infer from (71) that

dx = mx

ε M 1 ( M 1 z − M 2 )mx −1 dz, z ∈ [0, 1]. ρˆ1

(73)

In view of (51), one can conveniently express the transformations (68) and (71) as

x = ε xˇ ,

(74)

where xˇ depends on z ∈ [0, 1] and the Duffy parameter y ∈ [0, 1] as



xˇ = xˇ ( z, y ) =

ˇ ( y ) + Bˇ ( y )(C 1 ( y ) z − C 2 ( y ))mx , A ˇ ( y ) + ( M 1 ( y ) z − M 2 ( y ))mx /ρˆ1 ( y ), A

if ρˆ1 > |bˆ 1 |, if ρˆ1 = |bˆ 1 |.

In an analogous manner, one can rewrite (69), (70) and (72), (73) as

ρ1 = ερˇ1 , dx = ε Kˇ dz, where ρˇ1 = ρˇ1 ( z, y ) and Kˇ = Kˇ ( z, y ) are suitably defined for both cases ρˆ1 > |bˆ 1 | and ρˆ1 = |bˆ 1 |.

(75)

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 6. Relative errors for Hi j ≡ Hi j on parallel triangles with

987

ε = 10−2 and D = 7.071 × 10−3 .

Finally, the following polynomial transformation

y = tm y ,

t ∈ [0, 1]

(76)

is utilized for the Duffy variable y ∈ [0, 1], where m y is an arbitrary positive integer. By use of (74) through (76), the Galerkin integral (46) can be transformed to

Fij = a b myε

1 2

t 0

m y −1

 1

 x¯ ( z, t ) K¯ ( z, t )ψ¯ i ( z, t )¯c j ( z, t ) g¯ ( z, t ) dz dt ,

i , j = 1, 2, 3,

(77)

0

where a substitution y = y (t ) given by (76) has been used; x¯ ( z, t ) = xˇ ( z, y (t )), K¯ ( z, t ) = Kˇ ( z, y (t )); the transformed shape functions ψ¯ i ( z, t ) = ψi (ε x¯ ( z, t ), y (t )); also c¯ j ( z, t ) = c j (ε x¯ ( z, t ), y (t )) and g¯ ( z, t ) = g (ε x¯ ( z, t ), y (t )). Similar to the sinh transformation case, the factor ε 2 in (77) can be used to regularize the transformed integral before applying the standard product Gauss–Legendre rule. Remark 4. It is important to note that when mx = m y = 1, then the numerical procedure for computing nearly-singular Galerkin integrals employing the polynomial mapping is equivalent to that using Duffy’s transformation. Moreover, since the parameters mx and m y are not known a priori, the task of selecting the optimal exponents for any given configuration of the integration triangles E p and E q is still an open question. Finally, the methodology elaborated in this section for the computation of nearly-singular Galerkin integrals can also be used to effectively evaluate well-separated Galerkin integrals. 6. Results To illustrate the performance and to compare the efficacy of the various mappings described in Section 5, several numerical examples involving parallel and orthogonal triangles, as well as complete boundary-value problems for the Laplace equation have been performed. Recall that to calculate non-singular Galerkin integrals, smoothing transformations are first

988

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 7. Relative errors for Ti j ≡ Ti j on parallel triangles with

ε = 10−4 and D = 7.071 × 10−5 .

applied, followed by an ( N G × N G )-point Gauss–Legendre rule. In the figures, relative errors for all components of the pq pq pq Galerkin integrals Gi j , Hi j and Ti j are plotted against the number of quadrature points N G in [0, 1]. In the numerical experiments, relatively small values of N G have been employed. Namely, N G = 7, 14, 21, 28, 35, 42 have been used for the treatment of nearly-singular integrals. To compute the relative errors presented on the figures, all reference values have been obtained from the algorithm utilizing the Duffy transformation with N G = 1000. In the legends of the figures and in the tables, DT stands for the numerical approach using the Duffy transformation; ST is the procedure employing the sinh transformation; PT is the computational scheme using the polynomial transformation. As to the PT algorithm itself, numerical experiments based on full three-dimensional boundary-value problems for the Laplace equation have indicated that the smoothing exponents mx and m y should not be too high. For this reason and Remark 4, one has set mx = 3 and m y = 3 for all tests reported in this article. 6.1. Parallel triangles With reference to a Cartesian coordinate system {0; x1 , x2 , x3 }, let the outer triangle E p represent the unit twodimensional simplex in R3 with vertices given as x1 = (0, 0, 0)τ , x2 = (0, 1, 0)τ , x3 = (1, 0, 0)τ , where the superscript “τ ” represents the transpose symbol. Next, let the vertices of the inner triangle E q be prescribed as y 1 = (0, 0, ε )τ , √ y 2 = (1, 0, ε )τ , y 3 = (0, 1, ε )τ , where ε is the distance between E p and E q . With this setting, h p = hq = 2, where h p √ and hq denote respectively the diameters of E p and E q . It now follows from (35) that the near-singular ratio D = ε / 2. pq Upon specifying ε = 0.0001, Fig. 3 depicts the relative errors for all components of the Galerkin integrals Gi j involving the weakly-singular kernel G. In this situation, D = 7.071 × 10−5 . It is seen from the figure that the DT algorithm outperforms both the PT and ST procedures. Moreover, one can also infer from Fig. 3 that the PT scheme is more effective than the ST approach. Similar conclusions have been observed for cases when ε = 0.001, and ε = 0.01 as shown on Fig. 4. pq On Fig. 5, the relative errors for all components of the Galerkin integrals Hi j featuring the strongly-singular kernel H are displayed for the case when ε = 0.0001 corresponding again to D = 7.071 × 10−5 . Except for the components H13 and H21 where the PT method is clearly superior, all transformation techniques shown on Fig. 5 give comparable results. Although not presented in this publication, when ε = 0.001, the DT approach is more effective than the remaining methods except

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 8. Relative errors for Ti j ≡ Ti j on parallel triangles with

989

ε = 10−2 and D = 7.071 × 10−3 .

again for H13 and H21 where the superiority of the PT algorithm has significantly decreased. In the situation where pq as depicted on Fig. 6, the DT method is more effective for all components of Hi j .

ε = 0.01

pq

On Fig. 7, the relative errors for all components of the Galerkin integrals Ti j associated with the hyper-singular kernel T are plotted for the case when ε = 0.0001. One can conclude from the figure that the PT algorithm is very efficient on the components T13 and T21 only. However, this performance of the PT method on T13 and T21 is greatly reduced when ε = 0.001, while the DT approach provides better results on the remaining components. In the case when ε = 0.01 as shown pq on Fig. 8, The DT scheme is the most effective method for all components of Ti j . 6.2. Orthogonal triangles With reference to the Cartesian frame {0; x1 , x2 , x3 }, let x1 = (0, 0, 0)τ , x2 = (0, 1, 0)τ , x3 = (1, 0, 0)τ specify the vertices of the outer triangle E p . Now, let y 1 = (0, 0, ε )τ , y 2 = (1, 0, ε )τ , y 3 = (0, 0, 1 + ε )τ be the vertices of the inner triangle E q , where ε is again the distance between E p and E q ; the superscript “τ ” is the transpose symbol. With this specification, E p and E q are perpendicular and each triangle possesses an edge that is parallel to the x1 -axis. This particular configuration will be referred to as x1 -orthogonal triangles. In addition, it is also important to consider a configuration called x2 -orthogonal triangles in which E p is prescribed as before but E q is given by y 1 = (0, 0, ε )τ , y 2 = (0, 0, 1 + ε )τ , y 3 = (0, 1, ε )τ . pq On letting ε = 0.0001, Fig. 9 shows the relative errors for all components of the hyper-singular Galerkin integrals Ti j on x1 -orthogonal triangles. It is evident from the figure that the DT algorithm is more effective than the ST and PT methods. Here, again, the near-singular ratio D = 7.071 × 10−5 . Moreover, the ST scheme performs better than the PT approach on most components. When ε is raised to 0.001 and 0.01, an analogous superiority of the DT and ST techniques is observed over the PT approach. Although, the similarity of the DT and ST is more pronounced. pq On Fig. 10, the relative errors for the hyper-singular Galerkin integrals Ti j on x2 -orthogonal triangles are plotted for the case when ε = 0.0001. For this configuration, the PT method outperforms both the DT and ST procedures. Note that the behavior of the DT and ST approaches remains virtually unchanged at corresponding vertices when compared with the configuration depicted on Fig. 9. More specifically, T11 , T12 , T13 , T21 , T22 , T23 , T31 , T32 , T33 from Fig. 9 are compared respectively with T11 , T13 , T12 , T31 , T33 , T32 , T21 , T23 , T22 from Fig. 10. In other words, nodes {1, 2, 3} of element E q from the x1 -orthogonal configuration correspond respectively to nodes {1, 3, 2} of the same inner triangle in the x2 -orthogonal

990

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 9. Relative errors for Ti j ≡ Ti j on x1 -orthogonal triangles with

ε = 10−4 and D = 7.071 × 10−5 .

setting. Similar conclusions of the PT scheme have been obtained when ε is increased to 0.001 and 0.01. However, the gap on each graph between the PT algorithm and the DT or ST method is greatly reduced. 6.3. Mixed problem on a standard cube To further expose the behavior of the different strategies for computing non-singular Galerkin integrals, consider a mixed boundary-value problem for the Laplace equation in the standard cube {(x1 , x2 , x3 ) ∈ R3 : 0 < x1 , x2 , x3 < 1} based on an oscillatory potential u (x1 , x2 , x3 ) = e x1 sin x3 + e x3 cos x2 . The known harmonic function is prescribed on the face {x1 = 0, 0  x2 , x3  1}, and the flux t = n · ∇ u is specified on the remaining surface components, where ∇ u = (e x1 sin x3 , −e x3 sin x2 , e x1 cos x3 + e x3 cos x2 ) and n is the unit outward normal to the boundary of the cube. This problem will be resolved using a singular Galerkin BIE method involving only the weakly-singular kernel G and the strongly-singular kernel H given by (3). For the treatment of the inherent singular integrals, the numerical codes employ a semi-analytic integration technique for singular Galerkin BIEs that was originally developed in [14] for hyper-singular Galerkin BIEs. In this semi-analytic singular treatment, closed-form expressions are obtained for the coincident double integrals. In addition, in view of Remark 4, the non-singular Galerkin integrals are carried out separately with the DT, ST and PT procedures elucidated in Section 5. In doing so, three computer codes with the same semi-analytic singular treatment are generated. The difference among the codes resides in the way the non-singular Galerkin integrals are calculated. Moreover, the resulting dense linear system in each code is solved iteratively via the BiCGSTAB(3) [20] with a general-purpose sparse preconditioner for BEM detailed in [15]. In the linear solver, the relative tolerance on the BiCGSTAB(3) is set to 10−7 . As commonly done in practice in order to solve such problem, the surface of the cube is partitioned into non-overlapping triangles called boundary elements, each of which is characterized by its vertices or nodes. This discretization process produces a set of boundary nodes and a set of boundary elements. Next, boundary conditions of the problem are assigned at every surface nodes accordingly. On employing various discretizations consisting of uniform triangles, results of this mixed problem are provided on Table 1 in terms of the relative error on the boundary potential (e u ) and relative error on the boundary flux (et ). In the table, N is the total number of boundary nodes; N E is the total number of triangles covering the surface of the cube. The relative errors are expressed as e u = u − u e /u e  and et = t − t e /t e , where u e ∈ R N and t e ∈ R N are vectors representing respectively the exact potential and flux evaluated at all boundary nodes via the known harmonic function;

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

pq

Fig. 10. Relative errors for Ti j ≡ Ti j on x2 -orthogonal triangles with

991

ε = 10−4 and D = 7.071 × 10−5 .

Table 1 Mixed problem on a standard cube; Gauss–Legendre integration points and weights: N G 1 = 20, N G 2 = 14, N G 3 = 4. N

486 2646 5400 7350 10 086 14 406

NE 768 4800 10 092 13 872 19 200 27 648

DT

ST

PT

eu

et

eu

et

eu

et

2.207 × 10−3 3.454 × 10−4 1.629 × 10−4 1.182 × 10−4 8.511 × 10−5 5.887 × 10−5

3.064 × 10−3 6.742 × 10−4 3.764 × 10−4 2.985 × 10−4 2.359 × 10−4 1.789 × 10−4

2.207 × 10−3 3.454 × 10−4 1.629 × 10−4 1.182 × 10−4 8.512 × 10−5 5.886 × 10−5

3.067 × 10−3 6.752 × 10−4 3.760 × 10−4 2.991 × 10−4 2.355 × 10−4 1.795 × 10−4

2.067 × 10−3 3.288 × 10−4 1.721 × 10−4 1.294 × 10−4 1.081 × 10−4 9.322 × 10−5

5.189 × 10−3 4.127 × 10−3 5.412 × 10−3 3.918 × 10−3 4.038 × 10−3 5.274 × 10−3

u ∈ R N and t ∈ R N are vectors whose entries are made up respectively with the computed potential and flux at all surface nodes using the proposed techniques; here  ·  is the usual Euclidean vector norm in R N . In addition, N G 1 is the number Gauss–Legendre quadrature points in [0, 1] for the evaluation of edge-adjacent and vertex-adjacent integrals, N G 2 is the number Gauss–Legendre integration points in [0, 1] for the calculation of nearly-singular integrals, and N G 3 is the number Gauss–Legendre points in [0, 1] for the computation of well-separated integrals. It is clear from Table 1 that both the DT and ST approaches give almost identical accuracies on all mesh refinements. Moreover, the DT and ST methods outperform the PT scheme on full three-dimensional boundary-value problems. The accuracy on the surface flux utilizing the PT algorithm does not seem to converge or exhibits a zigzag pattern. These conclusions are even more accentuated on Table 2, where the results of this mixed problem are given in terms of the total execution time (Time) expressed in seconds, approximate rate of convergence of the boundary potential (βu ), and approximate rate of convergence of the boundary flux (βt ). Recall that an estimate of the rate of convergence can be computed as β = 2 ln(e 2 /e 1 )/ ln( N 1 / N 2 ), where e 1 and e 2 are the errors between two consecutive discretizations, and N 1 and N 2 are their respective number of nodes [14]. As can be seen from Table 2, the convergence rate of the DT or ST method is 2 for the boundary potential and 1.5 for the boundary flux in the L 2 norm. The discrepancy in execution time between the DT and ST procedures at every discretization is the result of functions calls to compute sinh, arcsinh and cosh that are required in the ST method. It is now quite evident from Table 2 that the PT strategy does not converge with the current

992

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

Table 2 Mixed problem on a standard cube; Gauss–Legendre integration points and weights: N G 1 = 20, N G 2 = 14, N G 3 = 4. N 486 2646 5400 7350 10 086 14 406

NE 768 4800 10 092 13 872 19 200 27 648

DT

ST

PT

Time

βu

βt

Time

βu

βt

Time

βu

βt

2.8 83 354 668 1276 2649

– 2.189 2.107 2.081 2.076 2.068

– 1.787 1.634 1.504 1.488 1.552

3.4 104 450 846 1612 3341

– 2.189 2.107 2.081 2.075 2.070

– 1.786 1.641 1.484 1.512 1.523

3 .4 103 442 834 1588 3296

– 2.170 1.815 1.850 1.137 0.831

– 0.270 −0.760 2.096 −0.191 −1.498

set of parameters, i.e., mx = 3, m y = 3, N G 1 = 20, N G 2 = 14, N G 3 = 4. One has to only increase the number of quadrature points for well-separated integrals to N G 3 = 7 in order to achieve the expected rate of convergence as in the DT and ST approaches. As a consequence, the cost of a convergent PT technique becomes very high. For instance, the execution time for the largest mesh on Table 2 with 14 406 nodes has increased from its current value of 3296 seconds to 8688 seconds when N G 3 = 7. Furthermore, setting the exponents mx and m y to higher values requires augmenting N G 2 and N G 3 to even greater numbers. Thus, the PT method proves to be very expensive on full three-dimensional boundary-value problems. 7. Conclusions In this article, three semi-analytic techniques for computing nearly-singular surface integrals of potential theory are developed and investigated within the framework of a three-dimensional Galerkin approximation. The proposed approaches are designed for boundaries discretized using flat triangles. All strategies separately rely on some type of nonlinear transformation to mollify near singularities, followed by a product Gauss–Legendre quadrature rule to calculate the outer boundary integral. In addition, the inner surface integral is evaluated exactly using explicit expressions that are defined over the edges of the integration triangle. Numerical experiments have demonstrated that the algorithm based on a Duffy transformation is the most efficient and effective method for dealing with nearly-singular as well as well-separated integrals. On the other hand, the procedure utilizing sinh transformations, while less effective on individual tests involving parallel and orthogonal pairs of triangles, is equally accurate and robust as compared to the Duffy transformation scheme on full three-dimensional boundary-value problems. Lastly, the technique employing polynomial transformations is the most expensive approach on complete threedimensional boundary-value problems owing to two adjustable parameters (smoothing exponents) in the formulation. The choice of these free parameters depends on the configuration between the interacting triangles. The potential benefit of this strategy will be realized only when an inexpensive method for selecting, at runtime, the optimal smoothing exponents will be developed. Acknowledgements This work was supported by the Office of Advanced Scientific Computing Research, U.S. Department of Energy, under contract No. DE-AC05-00OR22725 with UT-Battelle, LLC. References [1] A. Aimi, M. Diligenti, Numerical integration in 3D Galerkin BEM solution of HBIEs, Comput. Mech. 28 (2002) 233–249. [2] F.C. Araújo, L.J. Gray, Analysis of thin-walled structural elements via 3D standard BEM with generic substructuring, Comput. Mech. 41 (5) (2008) 633–645. [3] M.G. Duffy, Quadrature over a pyramid or cube of integrands with singularity at a vertex, SIAM J. Numer. Anal. 19 (6) (1982) 1260–1262. [4] S. Erichsen, S.A. Sauter, Efficient automatic quadrature in 3-d Galerkin BEM, Comput. Meth. Appl. Mech. Eng. 157 (3-4) (1998) 215–224. [5] C. Ericson, Real-Time Collision Detection, Morgan Kaufmann, San Francisco, CA, 2005. [6] L.J. Gray, A. Salvadori, A.-V. Phan, V. Manti˘c, Direct evaluation of hypersingular Galerkin surface integrals II, Electron. J. Boundary Elem. 4 (3) (2006) 105–130. [7] W. Hackbusch, S.A. Sauter, On numerical cubatures of nearly singular surface integrals arising in BEM collocation, Computing 52 (2) (1994) 139–159. [8] K. Hayami, H. Matsumoto, A numerical quadrature for nearly singular boundary element integrals, Eng. Anal. Boundary Elem. 13 (1994) 143–154. [9] B.M. Johnston, P.R. Johnston, D. Elliott, A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals, Int. J. Numer. Meth. Eng. 69 (7) (2007) 1460–1479. [10] P.R. Johnston, D. Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Int. J. Numer. Meth. Eng. 62 (4) (2005) 564–578. [11] D.E. Medina, J.A. Liggett, Exact integrals for three-dimensional boundary element potential problems, Commun. in Appl. Num. Methods 5 (1989) 555– 561. [12] S. Nintcheu Fata, Fast Galerkin BEM for 3D-potential theory, Comput. Mech. 42 (3) (2008) 417–429. [13] S. Nintcheu Fata, Explicit expressions for 3D boundary integrals in potential theory, Int. J. Numer. Meth. Eng. 78 (1) (2009) 32–47. [14] S. Nintcheu Fata, L.J. Gray, Semi-analytic integration of hypersingular Galerkin BIEs for three-dimensional potential problems, J. Comput. Appl. Math. 231 (2) (2009) 561–576. [15] S. Nintcheu Fata, L.J. Gray, On the implementation of 3D Galerkin boundary integral equations, Eng. Anal. Boundary Elem. 34 (1) (2010) 60–65. [16] S.A. Sauter, C. Schwab, Quadrature for hp-Galerkin BEM in R3 , Numer. Math. 78 (2) (1997) 211–258.

S. Nintcheu Fata / Applied Numerical Mathematics 60 (2010) 974–993

[17] [18] [19] [20] [21]

C. Schwab, Variable order composite quadrature for singular and nearly singular integrals, Computing 53 (2) (1994) 173–194. L. Scuderi, On the computation of nearly singular integrals in 3D BEM collocation, Int. J. Numer. Meth. Eng. 74 (11) (2008) 1733–1770. L. Scuderi, A new smoothing strategy for computing nearly singular integrals in 3D Galerkin BEM, J. Comput. Appl. Math. 225 (2) (2009) 406–427. G.L.G. Sleijpen, D.R. Fokkema, BiCGSTAB(l) for linear equations involving unsymmetric matrices with complex spectrum, ETNA 1 (1993) 11–32. Y. Wenjing, A new transformation technique for evaluating nearly singular integrals, Comput. Mech. 42 (3) (2008) 457–466.

993