Treatment of domain integrals in boundary element methods

Applied Numerical Mathematics 62 (2012) 720–735

Contents lists available at ScienceDirect

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

Treatment of domain integrals in boundary element methods ✩ 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 4 January 2010 Received in revised form 17 June 2010 Accepted 6 July 2010 Available online 13 July 2010 Keywords: Domain integral Newton potential Poisson equation Fundamental theorem of calculus Poincaré lemma Boundary element method

a b s t r a c t A systematic and rigorous technique to calculate domain integrals without a volumefitted mesh has been developed and validated in the context of a boundary element approximation. In the proposed approach, a domain integral involving a continuous or weakly-singular integrand is first converted into a surface integral by means of straightpath integrals that intersect the underlying domain. Then, the resulting surface integral is carried out either via analytic integration over boundary elements or by use of standard quadrature rules. This domain-to-boundary integral transformation is derived from an extension of the fundamental theorem of calculus to higher dimension, and the divergence theorem. In establishing the method, it is shown that the higher-dimensional version of the first fundamental theorem of calculus corresponds to the well-known Poincaré lemma. The proposed technique can be employed to evaluate integrals defined over simplyor multiply-connected domains with Lipschitz boundaries which are embedded in an Euclidean space of arbitrary but finite dimension. Combined with the singular treatment of surface integrals that is widely available in the literature, this approach can also be utilized to effectively deal with boundary-value problems involving non-homogeneous source terms by way of a collocation or a Galerkin boundary integral equation method using only the prescribed surface discretization. Sample problems associated with the three-dimensional Poisson equation and featuring the Newton potential are successfully solved by a constant element collocation method to validate this study. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The mathematical modeling of many scientific and engineering problems often requires the computation of integrals involving continuous or weakly-singular integrands defined over compact domains which are embedded in an Euclidean space of arbitrary but finite dimension. Such is the case when analyzing problems governed by elliptic partial differential equations featuring non-trivial source or body-force terms using the boundary integral equation (BIE) approach [18]. In the spirit of a boundary element method (BEM), an accurate technique which employs a surface-only discretization to evaluate these domain integrals is preferable. Such integration procedure could naturally allow an extension of BEM approaches to non-homogeneous or nonlinear problems [6]. Over the past years, several techniques have been proposed in the literature to deal with domain integrals without discretizing the computational domain into volume elements or internal cells [8,15,21,26,33]. To date, the dual reciprocity method [26] is the most popular approach used to transform domain integrals into boundary integrals. However, the accuracy and efficiency of the dual reciprocity scheme depend largely on the distribution and location of the radial basis functions employed to approximate the source term [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.07.003

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

721

Not long ago, an effective technique for evaluating two- and three-dimensional domain integrals with boundary-only discretization called the radial integration method has been presented in [8]. In this approach, a domain integral is transformed into an equivalent boundary integral by use of straight-path integrals emanating from the source point. In general, these straight-path integrals are carried out numerically [9]. However, the development of the radial integration method still requires an in-depth mathematical scrutiny. A somewhat similar approach can be found in [33] wherein the straight-path integral or radial integral is calculated exactly owing to the simple form of the radial basis functions used to expand the source term (see also [1]). The treatment proposed in this study originates from classical problems of vector analysis which attempt to construct the right inverse of differential operators such as gradient, curl and divergence operators [34]. Existing results for star-shaped domains rely essentially on certain straight-path integrals collectively known as Poincaré maps [5,12,16,17]. These Poincaré maps have recently been utilized in finite element methods to construct conforming finite element subspaces [4,12,16]. Of particular interest to the transformation of domain integrals into boundary integrals is the right inverse of the divergence operator. Namely, the problem is to find a continuous and differentiable vector potential on a compact Lipschitz domain such that its divergence coincides with a prescribed continuous function on the same domain. Once such a vector potential has been constructed, one can use the well-known divergence theorem [7] to convert a domain integral into a boundary counterpart. Proceeding in this manner, the proposed approach employs existing mathematical tools and an extension of the fundamental theorem of calculus to higher dimension to achieve the domain-to-boundary integral transformation. Within the same framework, it is shown that integrals involving weakly-singular integrands can be converted into weakly-singular surface integrals, the treatment of which is well established in the literature [22,23]. In addition, this approach is a systematic and rigorous generalization to any dimension of the radial integration method introduced in [8]. Details of the proposed technique, including the formulation and proof of the higher-dimensional version of the fundamental theorem of calculus, are provided in this article. In addition, mixed boundary-value problems for the threedimensional Poisson equation featuring the Newton potential are successfully dealt with via a constant element collocation BEM to validate this investigation. 2. Problem formulation Let Ω ⊂ Rk , N  k  1 be a bounded domain with Lipschitz-continuous boundary Γ . In addition, let n be the unit normal to Γ directed towards the exterior of Ω . Moreover, let x be a fixed point in Rk and assume that f (x, ·) is a continuous or weakly-singular function defined on Ω = Ω ∪ Γ . In the BEM community, the fixed point x ∈ Rk is usually called source point. Now, consider the task of computing the domain integral



f (x, y ) dΩ y

(1)

Ω

without partitioning Ω into internal cells. In other words, the problem is to calculate (1) by first representing it into a boundary integral. This problem can be conveniently reformulated as a task of constructing or finding a continuously differentiable vector field v (x, ·) such that

div y v (x, y ) = f (x, y ).

(2)

In Eq. (2), div y denotes the divergence operator with respect to the spatial variable y. Once such a vector potential v (x, ·) has been obtained, one can employ the divergence theorem [7] to resolve (1) as desired. Note that the solution to exercise (2) is not unique. Indeed, if v is a solution to (2), then the vector field v + u, where div y u = 0, is also a solution to problem (2). Remark 1. The interest in problem (2) goes well beyond a mere evaluation of domain integrals via surface integrals. It has been considered by many authors investigating fundamental questions of fluid dynamics [11,25,27]. In general, it can be viewed simply as an application of the classical Poincaré lemma involving differential forms [34]. A candidate solution to  (2) for star-shaped domains is provided by the Bogovski˘ı operator [2] with an additional integrability requirement that Ω f dΩ y =  0 (see also [5,10]). Moreover, for a given vector a on the boundary Γ satisfying the necessary compatibility condition Ω f dΩ y = Γ a · n dΓ y , the Bogovski˘ı operator can be used to construct a solution to (2) with the boundary condition v = a on Γ . However, the Bogovski˘ı operator is explicitly expressed as an integral over Ω , a representation which makes the Bogovski˘ı solution unsuitable to deal with problem (1). 3. Domain integrals with continuous integrands To successfully invert the partial differential equation (2) when f is continuous on Ω , the key idea is to construct the unknown vector field v by means of certain path integrals much like those that appear in the proof of Poincaré’s lemma in differential geometry [7,30,32]. To this end, let h(x, ·) be an extension of f (x, ·) into an arbitrary ball B centered at the source point x and containing Ω such that B \Ω h dΩ y is finite. Then, the sought vector field v (x, ·) is provided by the following result:

722

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

Theorem 3.1 (Generalized fundamental theorem of calculus: Part I). Let x be a fixed source point in Rk , k  1, and let Ω ⊂ Rk be a bounded domain with Lipschitz boundary Γ . In addition, assume that f (x, ·) is a continuous function defined on the compact Ω . Further, let h(x, ·) be a continuation of f (x, ·) into any ball B centered at x and containing Ω such that h(x, ·) is bounded and almost everywhere continuous in B \ Ω . Now, consider the function  y −x

R k−1 h(x, x + Re ) dR ,

F ( x, y ) =

y ∈ Ω,

(3)

0

where e ∈ R is a unit direction specified as k

e=

y−x .  y − x

(4)

Then, the vector field

v ( x, y ) =

y−x F (x, y ),  y − xk

y ∈ Ω,

(5)

is continuously differentiable with respect to y on Ω , and

div y v (x, y ) = f (x, y ),

y ∈ Ω.

(6)

Proof. A complete proof for the three-dimensional case, i.e. when k = 3, will be given below. The proof for the k-dimensional case is analogous to that of the three-dimensional situation. Detailed comments to establish the validity of this theorem for an arbitrary natural number k  1 will also be outlined at the end of this proof. Proof for k = 3: Since Ω is a bounded domain in R3 , there exists a ball B ⊂ R3 such that Ω ⊂ B. Indeed, one can set

M x = sup  y − x + 1,

(7)

y ∈Γ

and specify B = { y ∈ R3 :  y − x < M x }. Namely, B is a ball centered at the source point x ∈ R3 and containing Ω . The radius of B is M x . Without loss of generality, the prescribed function f can be extended from Ω into B, for example, as



h(x, y ) =

f (x, y ),

y ∈ Ω,

0,

y ∈ B \ Ω.

With the introduction of a continuation h(x, y ) of f (x, y ), the function  y −x

R 2 h(x, x + Re ) dR ,

F ( x, y ) =

y ∈ B,

(8)

0

and the vector field

v ( x, y ) =

y−x F (x, y ),  y − x3

y ∈ B.

(9)

are well defined on B. Let e be an arbitrary unit direction emanating from the source point x ∈ R3 . Now, consider the spherical coordinate system {x; R , θ 2 , θ 3 } with origin at x such that the induced local orthonormal reference is {x; e , e 2 , e 3 }. Namely, e is the unit vector in the direction of increasing radius R, and e i (i = 2, 3) is the unit tangent vector in the direction of increasing angle θ i . If {x; α 1 , α 2 , α 3 } is a local Cartesian coordinate system with origin at x, then

⎧ 1 2 ⎨ α = R cos θ , 2 2 3 ⎩ α = R sin θ cos θ ,

(10)

α = R sin θ sin θ . 3

2

3

With such definitions, the position vector R of an arbitrary point y ∈ B can be decomposed as

R = y − x = α 1 s1 + α 2 s2 + α 3 s3 ,

(11)

where si is the fixed unit basis vector in the Cartesian

r1 =

∂R ; ∂R

r2 =

∂R ; ∂θ 2

r3 =

∂R ∂θ 3

α i (i = 1, 2, 3) direction. Let

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

723

represent the tangent vectors at an arbitrary point y ∈ B in the R, θ 2 , θ 3 directions respectively. By use of (10) and (11), it can be shown that

⎧ e = rr 1  = cos θ 2 s1 + sin θ 2 cos θ 3 s2 + sin θ 2 sin θ 3 s3 , ⎪ ⎪ 1 ⎨ e 2 = rr 2  = − sin θ 2 s1 + cos θ 2 cos θ 3 s2 + cos θ 2 sin θ 3 s3 , 2 ⎪ ⎪ ⎩ e = r 3 = − sin θ 3 s + cos θ 3 s . 3

2

r 3 

(12)

3

In these settings, scalar fields h, F , and vector field v can be conveniently analyzed in the local orthonormal reference {x; e , e 2 , e 3 }. In fact, h, F , and v can be viewed as scalar and vector fields along the ray characterized by e. In {x; e , e 2 , e 3 },

h(x, y ) ≡ h(x, x + Re ),

R ∈ [0, M x ],

(13)

and, by virtue of (8),

R t 2 h(x, x + te ) dt ,

F (x, y ) ≡ F (x, x + Re ) =

R ∈ [0, M x ].

(14)

0

In view of (13), the continuity of h(x, ·) = f (x, ·) on Ω and the assumption of the theorem that h(x, ·) is bounded and almost everywhere continuous in B \ Ω , it follows that h(x, x + Re ) is bounded and almost everywhere continuous on [0, M x ]. This statement implies that R 2 h(x, x + Re ) is Riemann integrable on [0, M x ]. In other words, the function F (x, x + R e ) expressed via (14) is well defined on [0, M x ]. In fact, F (x, x + Re ) is Lipschitz-continuous on [0, M x ]. This assertion can be seen from the following estimate

R+ R



 2 F x, x + ( R +  R )e − F (x, x + Re ) = t h(x, x + te ) dt  C x  R R 2 + R  R + ( R )2 /3 , R

where R , R +  R ∈ [0, M x ] and C x = sup R ∈[0, M x ] |h(x, x + Re )|. It then follows from the Lipschitz continuity argument that F (x, x + Re ) is almost everywhere differentiable on [0, M x ] with

dF (x, x + R e ) dR

= R 2 h(x, x + R e ) almost everywhere on [0, M x ].

(15)

Since e is an arbitrary unit direction in the ball B, one can use (4) to specify e for an arbitrary point y ∈ B along the ray [0, M x ], and write



∇ y F ( x, y ) =

dF (x, x + Re ) =  y − xh(x, y )( y − x) e dR R = y −x

almost everywhere on B .

(16)

With these results, one can conclude that F (x, y ) is continuous with respect to y on B, and that F (x, y ) is almost everywhere differentiable with respect to y on B. In other words, F (x, y ) is differentiable with respect to y on B except on a subset of B \ Ω with Lebesgue measure zero, where h(x, y ) has a finite jump. From these arguments and letting y ∈ Ω , one can infer that F (x, y ) is continuously differentiable with respect to y on Ω , and one can write from (16) and the definition h(x, y ) = f (x, y ), y ∈ Ω , that

∇ y F (x, y ) =  y − x f (x, y ) ( y − x),

y ∈ Ω.

(17)

Similarly, (9) can be expressed in {x; e , e 2 , e 3 } as

v ( x, y ) ≡ v ( x, x + R e ) =

1 R2

F (x, x + Re )e ,

R ∈ [0, M x ].

(18)

In addition, on the interval [0, M x ], one can write

∇ v (x, x + Re ) = e ⊗

∂ v (x, x + Re ) 1 ∂ v (x, x + Re ) ∂ v (x, x + Re ) 1 + e2 ⊗ + e3 ⊗ . 2 2 ∂R R ∂θ R sin θ ∂θ 3

(19)

On employing (18), and the formulae

∂e = 0; ∂R

∂e = e2 ; ∂θ 2

∂e = sin θ 2 e 3 ∂θ 3

obtained directly from (12), one can verify that (19) takes the form

∇ v (x, x + Re ) =

1 dF (x, x + Re ) R2

dR

e⊗e+

F (x, x + Re ) R3

[I2 − 3e ⊗ e ],

R ∈ [0, M x ],

(20)

724

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

where I2 = e ⊗ e + e 2 ⊗ e 2 + e 3 ⊗ e 3 . I2 represents the symmetric second-order identity tensor on R3 . It was established that F (x, x + Re ) is Lipschitz-continuous on [0, M x ], thus almost everywhere differentiable on [0, M x ] with its derivative given by (15). One can now use (15) to rewrite (20) as

∇ v (x, x + Re ) = h(x, x + Re )e ⊗ e +

F (x, x + Re ) R3

[I2 − 3e ⊗ e ],

almost everywhere on [0, M x ].

(21)

It remains to show that F (x, x + Re )/ R 2 and F (x, x + Re )/ R 3 are continuous on [0, M x ]. Indeed, by use of the ensuing linear transformation t = R z, z ∈ [0, 1], one can express (14) as

F (x, x + Re ) = R 3 g (x, x + Re ),

R ∈ [0, M x ],

(22)

where

1 z2 h(x, x + R ze ) dz,

g (x, x + Re ) =

R ∈ [0, M x ].

0

As an integral of a bounded and almost everywhere continuous function on [0, 1] (i.e. z2 h(x, x + R z e ) is bounded and continuous everywhere on [0, 1] except on a set of Lebesgue measure zero), g (x, x + Re ) is continuous on [0, M x ]. With this latter argument and taking into consideration the relation (22), one can conclude that F (x, x + Re )/ R 2 and F (x, x + Re )/ R 3 are continuous on [0, M x ]. Hence, by virtue of (18), v (x, x + Re ) is continuous on [0, M x ], and the second term of ∇ v (x, x + Re ) in (21) is also continuous on [0, M x ]. As before, since e is an arbitrary unit direction in the ball B, one can again use (4) to prescribe the unit direction e for an arbitrary point y ∈ B along the ray [0, M x ], and utilize (21) to reveal that

∇ y v (x, y ) = ∇ v (x, x + Re )| R = y −x = h(x, y )



( y − x) ⊗ ( y − x) F ( x, y ) ( y − x) ⊗ ( y − x) + I − 3 almost everywhere on B . 2  y − x2  y − x3  y − x2

(23)

With these results, it is now clear that the vector field v (x, y ) given by (9) is continuous with respect to y on B, and that ∇ y v (x, y ) is continuous everywhere on B except on a subset of B \ Ω with Lebesgue measure zero, where h(x, y ) has a finite jump. From these statements and letting y ∈ Ω , it follows that v (x, y ) is continuous and possesses a continuous gradient with respect to y on Ω , and one can write from (23) and the condition h(x, y ) = f (x, y ), y ∈ Ω that



( y − x) ⊗ ( y − x) F ( x, y ) ( y − x) ⊗ ( y − x) , ∇ y v (x, y ) = f (x, y ) + I2 − 3  y − x2  y − x3  y − x2

y ∈ Ω.

(24)

Finally, with the aid of the formulae ( y − x) · ( y − x) =  y − x2 and tr(I2 ) = 3, one can directly obtain from (24) that div y v (x, y ) = f (x, y ), y ∈ Ω . This statement evidently coincides with (6). Here, tr(I2 ) is the trace of I2 . This concludes the proof for the three-dimensional case, i.e. when k = 3. The proof for the k-dimensional case is a straightforward adaptation of the three-dimensional case. Assume that B is any ball centered at the source point x ∈ Rk such that Ω ⊂ B ⊂ Rk . Here, one must consider the spherical coordinate system {x; R , θ 2 , θ 3 , . . . , θ k } with origin at the source point x ∈ Rk such that the induced local orthonormal reference is {x; e , e 2 , e 3 , . . . , ek }. Again, e is the unit vector in the direction of increasing radius R, and e i (i = 2, 3, . . . , k) is the unit tangent vector in the direction of increasing angle θ i . On denoting the local Cartesian coordinate system with origin at x ∈ Rk as {x; α 1 , α 2 , α 3 , . . . , α k }, one must employ the coordinate relations

⎧ 1 α = R cos θ 2 , ⎪ ⎪ ⎪ ⎪ ⎪ α 2 = R sin θ 2 cos θ 3 , ⎪ ⎪ ⎪ ⎪ 3 ⎨ α = R sin θ 2 sin θ 3 cos θ 4 , ⎪ .. ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ ⎪ αk−1 = R sin θ 2 sin θ 3 · · · sin θ k−1 cos θ k , ⎪ ⎩ αk = R sin θ 2 sin θ 3 · · · sin θ k−1 sin θ k .

(25)

In addition, note that the position vector R of an arbitrary point y ∈ B is specified, in this case, as

R = y−x=

k 

α i si ,

(26)

i =1

where si is the fixed unit basis vector in the Cartesian

α i (i = 1, 2, 3, . . . , k) direction. 2

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

725

Remark 2. The validity of Theorem 3.1 is greatly simplified if f (x, ·) is continuously extendible into any ball centered at x ∈ Rk and containing the compact Ω . In addition, note that the unit direction e ∈ Rk , i.e. with e  = 1, expressed via (4) is not defined at the source point x ∈ Rk , but F (x, x) = 0 by construction. Remark 3. It is easy to verify that the vector field (5) also satisfies

∂ v i ( x, y ) ∂ v j ( x, y ) − = 0, ∂yj ∂ yi

y ∈ Ω ⊂ Rk , i , j = 1, 2, 3, . . . , k.

Remark 4. Upon making the change of variable R = t  y − x, t ∈ [0, 1], in the straight-line integral (3), it can be shown using (4) that the vector field (5) admits the representation

1 v ( x, y ) = ( y − x)





t k−1 h x, x + t ( y − x) dt ,

y ∈ Ω.

(27)

0

The representation (27) is obviously void of the apparent singularity featured in (5). Moreover, formula (27) corresponds to the well-known expression in the finite element literature for the right inverse of the divergence operator commonly defined for domains that are star-shaped with respect to x [12]. As can be seen from Theorem 3.1, the star-shaped requirement can be relaxed since Ω can be viewed as a k-dimensional manifold that is embedded in Rk instead of existing on its own. In this regard, the straight-path from the source point x to any field point y ∈ Ω always exists. The integral operator (27) is often called Poincaré map for the divergence operator, and it is usually expressed in terms of differential forms [16,17]. In the context of boundary elements, the path integral (3) was proposed in [8,9] to evaluate domain integrals with boundary-only discretization for two- and three-dimensional problems with the source point x ∈ Ω . Remark 5. In the language of differential calculus, the continuous function f featured in Theorem 3.1 can be identified with a k-form ω . In this correspondence, ω is a closed form, i.e. dω = 0, where d is the exterior derivative operator. In addition, the continuously differentiable vector field v can be identified with a (k − 1)-form ξ . In this consideration, Theorem 3.1 establishes that any closed k-form ω on a compact domain Ω ⊂ Rk with Lipschitz-continuous boundary is exact, i.e. there exists a (k − 1)-form ξ defined on Ω such that ω = dξ . This latter assertion is analogous to the statement of Poincaré’s lemma without the star-shaped condition (see Remark 4). With the result of Theorem 3.1 which claims that any continuous scalar field f (x, ·) defined on a compact domain with Lipschitz boundary is the divergence of some vector potential v (x, ·), one can now transform the domain integral (1) into an equivalent boundary integral. To this end, one can formulate the following: Theorem 3.2 (Generalized fundamental theorem of calculus: Part II). Let x be a fixed source point in Rk , k  1, and let Ω ⊂ Rk be a bounded domain with Lipschitz boundary Γ . In addition, let n be the unit normal to Γ directed towards the exterior of Ω . Moreover, assume that f (x, ·) is a continuous function defined on the compact Ω , and let h(x, ·) be a continuation of f (x, ·) into any ball B centered at x and containing Ω such that h(x, ·) is bounded and almost everywhere continuous in B \ Ω . Now, consider the function  y −x

R k−1 h(x, x + Re ) dR ,

F ( x, y ) =

y ∈ Ω,

0

where e ∈ R is a unit direction specified as k

e= Then

y−x .  y − x



 f (x, y ) dΩ y = Ω

Γ

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Rk .

(28)

Proof. Theorem 3.2 directly follows from Theorem 3.1 and the divergence theorem [7,11]. Indeed, from the assumptions of Theorem 3.2, one can infer from Theorem 3.1 that there exists on Ω a continuously differentiable with respect to y vector field

v ( x, y ) =

y−x

 y − xk

F (x, y ),

y ∈ Ω,

726

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

provided by (5) such that div y v (x, y ) = f (x, y ), y ∈ Ω (see (6)). Next, application of the divergence theorem yields (28), i.e.





f (x, y ) dΩ y = Ω



div y v (x, y ) dΩ y = Ω

v (x, y ) · n( y ) dΓ y . Γ

This completes the proof of the theorem.

2

4. Domain integrals with weakly-singular integrands In an approximation of partial differential equations by the BIE method, it is often required to compute field variables at a set of internal points. These calculations involve domain integrals with weakly-singular kernels. In mathematical terms, let again Ω ⊂ Rk , k  1, be a bounded domain with Lipschitz boundary Γ . A kernel g is said to be weakly singular in Ω if it is defined and continuous for all x, y ∈ Ω , x = y, and there exist positive constants C and β ∈ [0, k) such that

g (x, y )  C x − y −β ,

(29)

for all x, y ∈ Ω , x = y. Now, let the source point x ∈ Ω be selected and fixed. In addition, let the integrand of the domain integral (1) be of the form

f (x, y ) = g (x, y )b( y ),

x, y ∈ Ω,

(30)

where g (x, y ) is a weakly-singular kernel in any ball centered at x and containing Ω and b( y ) is a continuous function prescribed on Ω . Further, let B ε be a ball centered at the source point x with radius ε > 0 specified as B ε = { y ∈ Rk :  y − x < ε }, and consider the domain Ωε = B ε ∩ Ω . With these definitions, the domain integral, that is of interest in this section, can be defined as an improper integral as





g (x, y ) b( y ) dΩ y := lim

ε →0

Ω

g (x, y ) b( y ) dΩ y ,

x ∈ Ω.

(31)

Ω\Ω ε

For a weakly-singular kernel g and a continuous function b defined on Ω , it can be shown that the improper integral on the right-hand side of (31) exists [14,20]. With these settings, one can now establish the following result: Theorem 4.1. Let Ω ⊂ Rk , k  1, be a bounded domain with Lipschitz boundary Γ . In addition, let n be the unit normal to Γ directed towards the exterior of Ω , and let x ∈ Ω be a fixed source point. Moreover, let g (·,·) be a weakly-singular kernel in any ball B centered at x and containing Ω , and assume that b(·) is a continuous function defined on the compact Ω . Further, let h(·) be a continuation of b(·) into B such that h(·) is bounded and almost everywhere continuous in B \ Ω . Now, consider the function  y −x

R k−1 g (x, x + Re )h(x + Re ) dR ,

F ( x, y ) =

y ∈ Ω,

(32)

0

where e ∈ Rk is a unit direction specified as

e= Then

y−x .  y − x

(33)



 g (x, y ) b( y ) dΩ y = Ω

Γ

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Ω.

(34)

Proof. Since Ω is a bounded domain in Rk , one can set

M x = sup  y − x + 1, y ∈Γ

and specify a ball B ⊂ Rk containing Ω as B = { y ∈ Rk :  y − x < M x }. Moreover, in view of the continuity of h = b on Ω and the assumption of the theorem that h is bounded and almost everywhere continuous in B \ Ω , it follows that h is bounded and almost everywhere continuous in the ball B, and one can introduce the quantity

h∞ = sup h( y ) . y∈ B

(35)

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

727

Fig. 1. Weakly-singular integration: (a) the source point x is an interior point of Ω and Γε = S kε−1 ; (b) the source point x is on the boundary Γ .

Further, owing to the fact that g (·,·) is a weakly-singular kernel in B ⊃ Ω , one can employ the definition (29) and the bound (35) to write an estimate for the line integral (32) as k−β F (x, y )  C h∞  y − x , k−β

y ∈ B.

(36)

Now, since the integrand g (x, y )b( y ) is continuous on Ω \ Ωε , one can use Theorem 3.2 to write





g (x, y )b( y ) dΩ y = Ω\Ω ε

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Ω,

(37)

∂(Ω\Ω ε )

where Ωε is defined in the preamble of this section. With reference to Fig. 1(a), let the source point x be an interior point of Ω , i.e. x ∈ Ω . In such situation, it is useful to select the radius ε so small that Ωε = B ε ⊂ Ω , where B ε is a ball centered at x with radius ε . With this choice, ∂(Ω \ Ω ε ) = Γ ∪ S kε−1 , where S kε−1 = ∂ B ε is the (k − 1)-sphere of radius ε in Rk . On account of this surface decomposition, one can express (37) as





( y − x) · n ( y ) F (x, y ) dΓ y +  y − xk

g (x, y )b( y ) dΩ y = Γ

Ω\Ω ε



S kε−1

( y − x) · n( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Ω.

(38)

It now remains to estimate the last integral featured in (38). Indeed, by use of (36), one can write



 C h∞ ( y − x) · n ( y ) F (x, y ) dΓ y  k−β  y − xk

S kε−1

S kε−1

|( y − x) · n( y )| dΓ y ,  y − xβ

x ∈ Ω.

(39)

But, on the sphere S kε−1 , it is useful to make the following change of variable y = x + ε e, where e ∈ S k1−1 and S k1−1 is the (k − 1)-sphere of radius 1 in Rk . Under this change of variable, dΓ y = εk−1 dS, where dS simply denotes the surface element on S k1−1 . In addition, note that the unit outward normal n( y ) = −e on S kε−1 . With these settings, one can verify that



k −1

( y − x) · n( y ) dΓ y = −εk−β ωk−1 ,  y − xβ

x ∈ Ω,

(40)

Sε

where

ωk−1 represents the surface area of the unit sphere S k1−1 . In view of (40), one can rewrite (39) as  C h∞ ωk−1 k−β ( y − x) · n ( y ) F (x, y ) dΓ y  ε , x ∈ Ω. k k−β  y − x

(41)

k −1

Sε

Since 0  β < k, it is clear that the integral appearing in (41) goes to zero when (38), and considering (41) and the definition (31), one can write



ε →0

Ω





g (x, y ) b( y ) dΩ y ≡ lim

Ω\Ω ε

g (x, y )b( y ) dΩ y = Γ

ε → 0. Now, taking the limit as ε → 0 in

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Ω.

(42)

Now, let the source point x be situated on the boundary Γ of the domain Ω , i.e. x ∈ Γ . With reference to Fig. 1(b), let Γε be the boundary of B ε lying in Ω , and let Γε be the portion of Γ residing in the ball B ε . With these definitions,

728

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

∂(Ω \ Ω ε ) = (Γ \ Γε ) ∪ Γε . By virtue of this boundary partition, (37) takes the form





g (x, y ) b( y ) dΩ y = Γ \Γε

Ω\Ω ε

( y − x) · n( y ) F (x, y ) dΓ y +  y − xk



Γε

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

x ∈ Γ.

(43)

Next, since |( y − x) · n( y )|/ y − xβ is a positive function that is obviously defined on the sphere S kε−1 ⊃ Γε , it follows that



Γε

|( y − x) · n( y )| dΓ y   y − xβ



S kε−1

|( y − x) · n( y )| dΓ y = εk−β ωk−1 ,  y − xβ

x ∈ Γ,

(44)

where ωk−1 denotes again the surface area of the unit sphere S k1−1 in Rk . In deriving (44), a change of variable similar to the one used to obtain (40) was employed. With the aid of (36) and (44), one can provide an estimate for the last integral of (43) as

 ( y − x) · n ( y ) C h∞ ωk−1 k−β F (x, y ) dΓ y  ε , k k−β  y − x

x ∈ Γ.

(45)

Γε

With the estimate (45), the fact that 0  β < k and the definition (31), one can now transition to the limit as to reveal that





g (x, y )b( y ) dΩ y ≡ lim

ε →0

Ω



g (x, y )b( y ) dΩ y =

Ω\Ω ε

Γ

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − xk

Lastly, on inspecting (42) and (46), one can retrieve formula (34).

x ∈ Γ.

ε → 0 in (43) (46)

2

5. Illustrations Demonstrative examples are given below to facilitate the interpretation of Theorems 3.1 and 3.2 established in Section 3, and to provide some insights into their application. In addition, three-dimensional examples related to the utilization of Theorem 4.1 are also included. 5.1. One-dimensional case Let the source point x ≡ 0 coincides with the origin of the real axis R, and assume that f is a continuous function defined on the segment [a, b] ⊂ R. As indicated in Theorem 3.1 and also in Remark 2, it is useful to extend f either continuously or simply by zero into the interval [−b, b] containing [a, b]. Let h be such continuation of f into [−b, b], and for all y ∈ [−b, b] such that y = 0, let e = y /| y |. Now, introduce the function

| y | F ( y) =

h( R e ) dR ,

y ∈ [a, b].

(47)

0

Then, Theorem 3.1 states that

v ( y) =

y

| y|

F ( y ),

y ∈ [a, b],

(48)

is continuously differentiable on [a, b], and

dv dy

= f ( y ),

y ∈ [a, b].

(49)

Here, differentiability at the left-end point a is usually understood as the right-hand derivative, and differentiability at the right-end point b means the left-hand derivative. Moreover, since in the one-dimensional situation, the unit normal direction is conventionally prescribed at endpoints of [a, b] as



n( y ) =

−1, if y = a, +1, if y = b,

Theorem 3.2 simply says that

b f ( y ) d y = v (b) − v (a). a

(50)

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

729

Fig. 2. One-dimensional illustration.

Expression (50) corresponds to the well-known Newton–Leibniz formula wherein v is often referred to as an antiderivative of f . With these remarks, it is now clear that Theorems 3.1 and 3.2 can simply be viewed as a generalization to higher dimension of the fundamental theorem of calculus [28,31]. In particular, let f be specified, for example, as the function shown in Fig. 2(a). Then, one can write from (47) and (48) that

a v (a) = F (a) =

b h( R ) dR ;

v (b) = F (b) =

0

h( R ) dR .

(51)

0

By use of (51) in (50), one can show that

b f ( y ) d y = F (b) − F (a). a

On the other hand, if f is given as the function depicted in Fig. 2(b), where the interval [a, b] contains the origin x = 0, then an extension of the integrand f outside [a, b] is not necessary in this case and one can simply set h ≡ f in (47). With these remarks, one can infer from (47) and (48) that

|a| − v (a) = F (a) =

b f (− R ) dR ;

v (b) = F (b) =

0

f ( R ) dR .

(52)

0

Finally, by use of (52), formula (50) becomes

b f ( y ) d y = F (b) + F (a). a

5.2. Three-dimensional case Let the source point x ≡ 0 coincides with the origin of some reference Cartesian coordinate system {0; y 1 , y 2 , y 3 } in R3 . Now, assume that f ( y ) is prescribed as n f ( y ) = yl1 ym 2 y3

(53)

on a certain compact domain Ω ⊂ R3 with Lipschitz boundary Γ , where l, m, n are any non-negative integer. In this case, it is useful to extend f continuously outside Ω via the same formula given by the right-hand side of (53). Then, it can be shown via direct integration of (3) that

F ( y) =

n yl1 ym 2 y3

l+m+n+3

 y 3 .

(54)

Thus, with the aid of (54), the vector potential (5), using k = 3, is expressible as

v ( y) =

n yl1 ym 2 y3

l+m+n+3

y.

Lastly, by use of (54) in (28) with k = 3, one can write



n yl1 ym 2 y 3 dΩ y =

Ω

1 l+m+n+3



n yl1 ym 2 y 3 y · n( y ) dΓ y .

Γ

(55)

730

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

In particular, if f ( y ) = 1, i.e. l = m = n = 0, one can obtain from (55) that



dΩ y =

1



y · n( y ) dΓ y .

3

Ω

Γ

Furthermore, when f ( y ) = c · y with c describing a constant vector in R3 , it can be established, via sequential substitution in (55) the parameters (l = 1, m = n = 0), (m = 1, n = l = 0), (n = 1, l = m = 0) and summing up the results, that



c · y dΩ y =

1



(c · y ) y · n( y ) dΓ y .

4

Ω

Γ

As to the case with weakly-singular integrands, let the source point x be fixed in R3 and let

f i ( x, y ) =

yli

x − y 

i = 1, 2, 3,

,

(56)

define three functions on a compact domain Ω with Lipschitz boundary Γ , where l is a non-negative integer and y i (i = 1, 2, 3) are the components of the field point y ∈ Ω . Here again, f i (x, ·), i = 1, 2, 3, are simply redefined everywhere in R3 via the same formulae expressed by the right-hand side of (56). Then, one can introduce, using brute force integration of (3) or (32), the functions

F i ( x, y ) =

l  y − x2  (l + 1 − k) yli−k xki , 2 + 3l + l2

(57)

k =0

where xi (i = 1, 2, 3) are the components of the source point x ∈ R3 in the reference Cartesian frame. By virtue of (57) in expression (28) when the source point x ∈ R3 \ Ω , and by use of (57) in (34) when x ∈ Ω , one can compute the integral of f i (x, y ) over the domain Ω for all x ∈ R3 . More specifically, by use of (57) in (28) and (34) with k = 3, and the condition that l = 0 from (56), one can write



Ω

1

x − y 

1

dΩ y =



2 Γ

( y − x) · n ( y ) dΓ y ,  y − x

x ∈ R3 .

(58)

Moreover, by setting, for example, l = 1 and i = 1 in (56), one can again employ (28), (34), and (57) to reveal that



Ω

y1

x − y 

1

dΩ y =



6 Γ

( y − x) · n ( y ) (2 y 1 + x1 ) dΓ y ,  y − x

x ∈ R3 .

(59)

In general, let the weakly-singular kernel g (x, y ) appearing in Theorem 4.1 be expressed as

g (x, y ) =

1

x − y 

,

x, y ∈ R3 , x = y .

(60)

In addition, assume that the function h(·), also from Theorem 4.1, is an extension of b(·) into any ball centered at x and containing Ω . Then, by use of (60) in the straight-line integral (32) with k = 3, the change of variable R = z y − x, z ∈ [0, 1] and expression (33) for the unit direction e, one can verify that

F (x, y ) =  y − x2 F (x, y ),

y ∈ Ω,

(61)

where

1 F (x, y ) =





zh x + z( y − x) dz,

y ∈ Ω.

(62)

0

One can now utilize (61) in (28) when x ∈ R3 \ Ω , and (61) in (34) when x ∈ Ω , again with k = 3, to establish that

 Ω

b( y )

x − y 

 dΩ y = Γ

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − x

x ∈ R3 .

(63)

Formula (63) illustrates the transformation of a weakly-singular domain integral into a weakly-singular boundary integral.

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

731

Fig. 3. Left: Discretization of a cylindrical tube using 1600 triangles and 896 nodes. Right: Triangulation of an L-shaped domain with 578 triangles and 401 nodes.

6. Applications Consider the computational treatment of the Poisson equation ∇ 2 u + b = 0 in a bounded domain Ω ⊂ R3 with Lipschitzcontinuous boundary Γ using a collocation BIE method, where b is a continuous source function prescribed on Ω . To successfully deal with this problem, let t = n · ∇ u be the flux associated with the solution u, where n is the unit normal to Γ directed towards the exterior of Ω . In addition, for some R  ε > 0, let xε ∈ R3 \ Ω . On employing the Green’s representation formula [3,13,18], the unknown functions u and t on the surface Γ can be obtained by solving the singular BIE



 G (xε , y ) t ( y ) dΓ y −

lim

xε →x∈Γ

Γ



 H (xε , y ) · n( y )u ( y ) dΓ y +

Γ

G (xε , y )b( y ) dΩ y

= 0,

(64)

Ω

where G (x, y ) is the fundamental solution of the Laplace operator, and H (x, y ) = ∇ y G (x, y ). These kernels are given respectively by

1

G ( x, y ) =

1

4π  x − y 

,

H ( x, y ) =

1

x− y

4π x − y 3

,

x, y ∈ R3 , x = y .

(65)

Toapproximate (64), it is useful to assume that Γ is the surface of a polyhedral domain and consider a triangulation of Γ j into closed and non-overlapping boundary elements such that Γ j is an open flat triangle (see Fig. 3). In addition, assume that on each triangle Γ j , the boundary unknowns u and t are constants. Now, let N be the total number of boundary elements on the surface Γ . With these assumptions, a collocation method for solving (64) requires that the singular BIE be satisfied exactly at a set of collocation points {xi }iN=1 resting on Γ . In practice, a collocation point xi is often specified as the centroid of the triangle Γi . The above requirement leads to a dense linear system of algebraic equations for boundary unknowns u and t as

Γ =

G{t } − H{u } = −{ B },

(66)

where {u } and {t } are vectors whose components respectively represent unknown functions u and t evaluated at the centroid of each boundary element Γ j ( j = 1, 2, . . . , N). In (66), the components of influence matrices G and H can be expressed as j

Gi j = lim g j (xε ), xε → x i

Hi j = lim h j (xε ), xε → x i

j

xi ∈ Γ,

(67)

where the single-layer potential g j and double-layer potential h j are specified as





g j ( x) =

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

h j ( x) =

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

x ∈ R3 .

(68)

Γj

Moreover { B } is a vector characterizing the contributions of the source function b with component given by

B i = lim V(xε ), xε →xi

xi ∈ Γ,

(69)

where V is the Newton potential specified as



V( x) =

G (x, y )b( y ) dΩ y , Ω

x ∈ R3 .

(70)

732

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

Upon prescribing the boundary conditions (e.g. Dirichlet, Neumann or mixed type) for the specific boundary-value problem associated with the Poisson equation, the linear system (66) can be rearranged as

A{ z} = { f },

(71)

where { z} ∈ R is a vector containing unknown quantities u or t on the surface Γ , and { f } ∈ R is a vector whose entries are obtained from known boundary data and the contributions of the source term { B }. For an arbitrary source point x ∈ R3 , it was shown in [22] that surface potentials g j (x) and h j (x) admit exact representations in terms of recursive formulae expressed over the edges of the flat triangle Γ j . These analytic expressions will be utilized to compute the entries of the influence matrices G and H implicitly featured in (71). It now remains to evaluate the volume potential (70) at every collocation point xi ∈ Γ (i = 1, 2, . . . , N) using the methodology elucidated in Sections 3 and 4 to completely generate the right-hand side of the linear system (71). To this end, one can (i) substitute in (70) the kernel G (x, y ) by its corresponding expression provided by (65), and (ii) make use of (63) to reveal that N

V( x) =

N

1



4π Γ

( y − x) · n ( y ) F (x, y ) dΓ y ,  y − x

x ∈ R3 ,

(72)

where F (x, y ) is given by (62). With this boundary representation of the Newton potential V(x), and in view of the decomposition of the surface Γ into non-overlapping flat triangles Γ j ( j = 1, 2, . . . , N), one can also adopt a constant element approximation for F (x, ·) on each triangle Γ j so that (72) is estimated as

V( x) =

N 1 

4π



j

A j (x)F x, yˆ ,

x ∈ R3 ,

(73)

j =1

j

where yˆ denotes the centroid of the flat triangle Γ j and



A j ( x) = Γj

( y − x) · n ( y ) dΓ y ,  y − x

x ∈ R3 .

(74)

Following the treatment of surface potentials thoroughly investigated in [22], it can be established, using the fact that Γ j is a flat triangle, that A j (x) = η j I 1 j (x) resulting in a representation of V(x) as

V( x) =

N 1 

4π





η j I 1 j (x)F x, yˆ j , x ∈ R3 .

(75)

j =1

In expression (75), the parameter η j describes the relative distance between the source point x ∈ R3 and the plane of the triangle Γ j , and it is constructed using a local orthonormal companion reference of R3 associated with Γ j and having its origin at x such that ( y − x) · n( y ) = η j , y ∈ Γ j ;



I 1 j ( x) = Γj

1

 y − x

dΓ y ,

x ∈ R3 .

Exact formula for the generic integral I 1 j (x), x ∈ R3 , is also provided explicitly in [22]. In view of the continuity of the Newton potential V across the surface Γ (see, e.g. (72)), the passage to the limit in (69) is straightforward and one can simply write B i = V(xi ), xi ∈ Γ (i = 1, 2, . . . , N). In this section, formula (75) will be is utilized to approximately calculate V on Γ . With an estimate of the Newton potential at hand, the linear system (71) can now be resolved to obtain the solution u j and its normal derivative t j at the centroid of each boundary element Γ j ( j = 1, 2, . . . , N). Finally, the remaining task of computing the solution u (x) at an arbitrary interior point x in Ω , in this constant element framework, can be accomplished by use of the Green’s representation formula [18], (68) and (75) as

u ( x) =

N  j =1

t j g j ( x) −

N 

u j h j (x) + V(x),

x ∈ Ω.

(76)

j =1

All examples given below are stated and computed in a reference Cartesian coordinate system {0; y 1 , y 2 , y 3 }. For comparison purposes, test problems are formulated using exact solutions for the Poisson equation in the domain of interest Ω . In the tables, the accuracy of these numerical examples are characterized by the relative errors e u = u − u e /u e  and et = t − t e /t e , where u ∈ R N and t ∈ R N are vectors representing respectively the computed solution u and its flux t at all collocation points (i.e. centroids of boundary elements) on Γ ; u e ∈ R N and t e ∈ R N are vectors describing respectively

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

733

Table 1 Mixed problem on a cylindrical tube: Boundary solution. N

eu

βu

et

βt

400 1600 6400 25 600

3.401 × 10−2 1.220 × 10−2 4.501 × 10−3 1.659 × 10−3

– 1.479 1.439 1.440

1.616 × 10−2 9.283 × 10−3 5.601 × 10−3 3.657 × 10−3

– 0.800 0.729 0.615

Table 2 Mixed problem on a cylindrical tube: Evaluation at interior points. Coordinates

(3.99, 0, 0)

(−1.5, 2.598, 2)

(−1.05, −1.819, −3)

400 1600 6400 25 600

−0.6274742 −57.227771 −58.003332 −58.173962

7.890677 7.828704 7.762403 7.729389

5.501348 5.572604 5.640466 5.643616

Exact values

−58.214499

7.708201

5.628045

the exact solution and its normal derivative also evaluated at all collocation points on Γ . With these errors, one can provide an empirical rate of convergence of the proposed algorithm as β = 2 ln(e 2 /e 1 )/ ln( N 1 / N 2 ), where e 1 and e 2 are the relative errors between two consecutive discretizations N 1 and N 2 of the surface Γ . In addition, the dense linear system (71) is solved iteratively using the BiCGSTAB(3) [29] with a general-purpose sparse preconditioner for BEM detailed in [24]. In the linear solver, the relative tolerance on the BiCGSTAB(3) is set to 10−7 . 6.1. Mixed problem on a cylindrical tube Consider a mixed boundary-value problem for the Poissonequation ∇ 2 u + b = 0 in a cylindrical tube Ω centered at

the origin 0 = (0, 0, 0) such that Ω = {( y 1 , y 2 , y 3 ) ∈ R3 : 2 <

y 21 + y 22 < 4, −5 < y 3 < 5}. In this test problem, the exact

solution is specified as

u ( y ) = − y 31 +

y 21 + y 22 + y 23 3

y ∈ Ω,

,

(77)

leading to a source function given by

b ( y ) = 6 y 1 − 2,

y ∈ Ω.

On the top and bottom faces {2 



(78) y 21 + y 22  4, y 3 = ±5} of the boundary Γ , the known solution u given via (77) is

employed to specify the Dirichlet boundary conditions. On the remaining surface components of Γ , Neumann boundary conditions are prescribed as t = n · ∇ u, where ∇ u = (−3 y 21 + 23 y 1 , 23 y 2 , 23 y 3 ) and n denotes the unit outward normal on the surface Γ . A typical triangulation of Γ is depicted in Fig. 3 (left). In view of the source term (78) and the boundary representation for the Newton potential (72), one can use (58) and (59) to show that, in this situation, the path integral (62) can be carried out exactly as F (x, y ) = 2 y 1 + x1 − 1, where x1 is the component of the source point x ∈ R3 in the first coordinate direction. On employing several mesh refinements, results of this mixed problem can be seen in Tables 1 and 2. In Table 1, the relative errors on the boundary solutions (i.e. e u and et ) indicate that the proposed method is convergent. Moreover, the table shows the approximate rate of convergence for the function u (βu ) and that for the flux t (βt ) on the boundary Γ . The computation of the solution u at selected interior points (3.99, 0, 0), (−1.5, 2.598, 2) and (−1.05, −1.819, −3), is provided in Table 2, where the first column contains the number of surface elements employed in the calculation; the last row of the table includes the exact values of u evaluated at interior points using (77). Note that the first interior point (3.99, 0, 0) is very close the discretized surface. It is therefore evident from Tables 1 and 2 that the computation of the Newton potential is successful and accurate. 6.2. Mixed problem on an L-shaped domain With reference to Fig. 3 (right), consider a mixed boundary-value problem for the Poisson equation ∇ 2 u + b = 0 in an L-shaped domain Ω = Ω1 ∪ Ω2 , such that Ω1 = {( y 1 , y 2 , y 3 ) ∈ R3 : 0 < y 1 , y 2 < 2, 0 < y 3 < 6} and Ω2 = {( y 1 , y 2 , y 3 ) ∈ R3 : 0 < y 1 , y 3 < 2, 0 < y 2 < 4}. In this example, the exact solution is prescribed as

u ( y ) = y 32 e y 1 + y 3 ,

y ∈ Ω,

resulting in a source function given by

(79)

734

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

Table 3 Mixed problem on an L-shaped domain using exact solution (79) and source term (80): Boundary solution. N

eu

βu

et

βt

578 2312 9248 36 992

4.979 × 10−2 1.719 × 10−2 5.909 × 10−3 2.038 × 10−3

– 1.534 1.541 1.536

1.279 × 10−1 6.185 × 10−2 3.249 × 10−2 1.852 × 10−2

– 1.048 0.923 0.811

Table 4 Mixed problem on an L-shaped domain using exact solution (79) and source term (80): Evaluation at interior points. Coordinates

(1.5, 0.5, 0.01)

(0.7, 3, 1.3)

(1, 1, 5)

578 2312 9248 36 992

90.183668 31.682404 10.839389 3.870195

242.011190 214.482206 204.476363 201.110028

445.976760 421.044699 409.801924 405.603457

0.565841

199.504515

403.428793

Exact values

Table 5 Mixed problem on an L-shaped domain using exact solution (77) and source term (78): Boundary solution. N

eu

βu

et

βt

578 2312 9248 36 992

3.468 × 10−2 1.216 × 10−2 4.172 × 10−3 1.423 × 10−3

– 1.512 1.543 1.552

1.040 × 10−2 8.170 × 10−3 5.921 × 10−3 4.182 × 10−3

– 0.348 0.464 0.502

Table 6 Mixed problem on an L-shaped domain using exact solution (77) and source term (78): Evaluation at interior points. Coordinates

(1.5, 0.5, 0.01)

(0.7, 3, 1.3)

(1, 1, 5)

578 2312 9248 36 992

−2.231350 −2.417487 −2.472667 −2.519760

3.479950 3.415469 3.393762 3.386818

8.108386 8.035923 8.011438 8.003582

Exact values

−2.541633

3.383667

8.000000





b( y ) = − 2 y 32 + 6 y 2 e y 1 + y 3 ,

y ∈ Ω.

(80)

On the surfaces {( y 1 , y 2 , y 3 ) ∈ R : 0  y 1 , y 2  2, y 3 = 6} and {( y 1 , y 2 , y 3 ) ∈ R : 0  y 1 , y 3  2, y 2 = 4}, the exact solution (79) is utilized to specify the Dirichlet boundary conditions. On the remaining faces of Γ , Neumann boundary conditions are given as t = n · ∇ u, where ∇ u = ( y 32 e y 1 + y 3 , 3 y 22 e y 1 + y 3 , y 32 e y 1 + y 3 ) and n is the unit normal to Γ directed towards the exterior of the domain Ω . Since the path integral (62) does not always admit an analytic representation, one will employ a Gauss–Legendre quadrature rule to compute the integral F (x, y ) entering the Newton potential (75). In addition, quadrature formulae can be directly applied in this test problem owing to the fact that an explicit expression for the source function b( y ) is available. This exact expression is given by the right-hand side of (80). Indeed, a continuous extension, h( y ), of the source function b( y ) outside Ω is accomplished by the same formula (80). This extension guarantees that the path integral (62) is well defined (see e.g. Theorem 4.1 and also Remark 2). In this example, 7 Gauss–Legendre quadrature points and weights are employed in the calculation of the straight-line integral F (x, y ). Results of this test problem are illustrated in Tables 3 and 4 for different discretizations. One can conclude from both tables that the proposed approach for solving the Poisson equation employing only the underlying surface triangulation is convergent. However, one can also notice from Table 4 that the accuracy of the solution u at interior points is still deficient even for the largest mesh (containing 36992 elements) used in the table. This is particularly evident at (1.5, 0.5, 0.01), a point that is very close to the surface Γ of the domain Ω . The reason for this lack of accuracy is related to the fact that the numerical procedure is trying to recover a highly-varying function (79) given (80) using only constant element approximations. Indeed, if one simply replaces, in this mixed problem on an L-shaped domain, the exact solution (79) with (77) and similarly swaps the source function (80) with (78), then one obtains the results shown in Tables 5 and 6 after application of the proposed numerical algorithm. Comparison of Table 3 and Table 5 reveals that the rate of convergence on the solution (βu ) has settled around 1.5, whereas the convergence rate on the normal derivative (βt ) is approaching its limiting value from above as in Table 3 and from below as in Table 5. Furthermore, it is seen from Table 6 that the computation of a mildly-varying function u at interior points is now sufficiently accurate. 3

3

S. Nintcheu Fata / Applied Numerical Mathematics 62 (2012) 720–735

735

7. Conclusions In this study, an effective treatment of domain integrals with surface-only discretization has been exposed and its usefulness demonstrated in the context of a collocation boundary element technique. In the proposed method, an integral, defined over a simply- or multiply-connected domain and involving a continuous or weakly-singular integrand, is first transformed into a boundary integral by means of straight-line integrals that intersect the underlying domain. Then, the resulting surface integral is calculated via standard methods. This transformation of an arbitrary domain integral into an equivalent boundary integral is based on a rigorous mathematical foundation and is derived from an extension of the fundamental theorem of calculus to higher dimension, and the divergence theorem. Combined with the singular treatment of surface integrals, this technique can be employed to successfully deal with boundary-value problems featuring internal sources using a collocation or a Galerkin boundary integral equation method without the need to partition the underlying domain into volume cells. Indeed, given only the surface triangulation, the boundary data of mixed-type and the source function, a complete treatment of the three-dimensional Poisson equation via a constant element collocation method has indicated that the proposed approach is accurate and robust. 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] E.L. Albuquerque, P. Sollero, W. Portilho de Paiva, The radial integration method applied to dynamic problems of anisotropic plates, Commun. Numer. Meth. Engng. 23 (9) (2007) 805–818. [2] M.E. Bogovski˘ı, Solution of the first boundary value problem for an equation of continuity of an incompressible medium, Dokl. Akad. Nauk SSSR 248 (5) (1979) 1037–1040. [3] M. Bonnet, Boundary Integral Equation Methods for Solids and Fluids, Wiley & Sons, New York, 1995. [4] M. Costabel, M. Dauge, L. Demkowicz, Polynomial extension operators for H 1 , H (curl) and H (div)-spaces on a cube, Math. Comp. 77 (264) (2008) 1967–1999. [5] M. Costabel, A. McIntosh, On Bogovski˘ı and regularized Poincaré operators for de Rham complexes on Lipschitz domains, Math. Z. 265 (2) (2010) 297–320. [6] J. Ding, W. Ye, L.J. Gray, An accelerated surface discretization-based BEM approach for non-homogeneous linear problems in 3-D complex domains, Internat. J. Numer. Methods Engrg. 63 (12) (2005) 1775–1795. [7] W.H. Fleming, Functions of Several Variables, Addison-Wesley, Massachusetts, 1965. [8] X.-W. Gao, The radial integration method for evaluation of domain integrals with boundary-only discretization, Eng. Anal. Boundary Elem. 26 (10) (2002) 905–916. [9] X.-W. Gao, Evaluation of regular and singular domain integrals with boundary-only discretization-theory and fortran code, J. Comput. Appl. Math. 175 (2) (2005) 265–290. [10] M. Geissert, H. Heck, M. Hieber, On the equation div u = g and Bogovski˘ı’s operator in Sobolev spaces of negative order, in: Partial Differential Equations and Functional Analysis, in: Operator Theory: Advances and Applications, vol. 168, Birkhäuser, Basel, 2006, pp. 113–121. [11] G.P. Gladi, An Introduction to the Mathematical Theory of the Navier–Stokes Equations. Vol. I: Linearized Steady Problems, Springer Tracts in Natural Philosophy, vol. 38, Springer-Verlag, New York, 1994. [12] J. Gopalakrishnan, L.F. Demkowicz, Quasioptimality of some spectral mixed methods, J. Comput. Appl. Math. 167 (1) (2004) 163–182. [13] W. Hackbusch, Elliptic Differential Equations: Theory and Numerical Treatment, Springer, New York, 1992. [14] W. Hackbusch, Integral Equations: Theory and Numerical Treatment, ISNM, Birkhäuser, Basel, 1995. [15] M.R. Hematiyan, A general method for evaluation of 2D and 3D domain integrals without domain discretization and its application in BEM, Comput. Mech. 39 (4) (2007) 509–520. [16] R. Hipmair, Canonical construction of finite elements, Math. Comp. 68 (228) (1999) 1325–1346. [17] R. Hipmair, Finite elements in computational electromagnetism, Acta Numer. 11 (2002) 237–339. [18] G.C. Hsiao, W.L. Wendland, Boundary Integral Equations, Applied Mathematical Sciences, vol. 164, Springer, New York, 2008. [19] S.-C. Hsiao, A.A. Mammoli, M.S. Ingber, The evaluation of domain integrals in complex multiply-connected three-dimensional geometries for boundary element methods, Comput. Mech. 32 (4–6) (2003) 226–233. [20] R. Kress, Linear Integral Equations, 2nd ed., Applied Mathematical Sciences, vol. 82, Springer, New York, 1999. [21] A.C. Neves, C.A. Brebbia, Multiple reciprocity boundary element method in elasticity. A new approach for transforming domain integrals to the boundary, Internat. J. Numer. Methods Engrg. 31 (4) (1991) 709–727. [22] S. Nintcheu Fata, Explicit expressions for 3D boundary integrals in potential theory, Internat. J. Numer. Methods Engrg. 78 (1) (2009) 32–47. [23] 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. [24] S. Nintcheu Fata, L.J. Gray, On the implementation of 3D Galerkin boundary integral equations, Eng. Anal. Boundary Elem. 34 (1) (2010) 60–65. [25] A. Novotný, I. Straškraba, Introduction to the Mathematical Theory of Compressible Flow, Oxford Lecture Series in Mathematics and its Applications, vol. 27, Oxford University Press, Oxford, 2004. [26] P.W. Partridge, C.A. Brebbia, L.C. Wrobel, The Dual Reciprocity Boundary Element Method, Computational Mechanics Publications, Southampton, 1992. [27] L. Poul, On dynamics of fluids in meteorology, Cent. Eur. J. Math. 6 (3) (2008) 422–438. [28] D.B. Scott, S.R. Tims, Mathematical Analysis: An Introduction, Cambridge University Press, London, 1966. [29] G.L.G. Sleijpen, D.R. Fokkema, BiCGSTAB(l) for linear equations involving unsymmetric matrices with complex spectrum, ETNA 1 (1993) 11–32. [30] M. Spivak, Calculus on Manifolds, Westview Press, Boulder, 1965. [31] J. Stewart, Calculus Early Transcendentals, 5th ed., Brooks/Cole, Belmont, 2003. [32] M.E. Taylor, Partial Differential Equations: Basic Theory, Springer, New York, 1996. [33] P.H. Wen, M.H. Aliabadi, D.P. Rooke, New method for transformation of domain integrals to boundary integrals in boundary element method, Commun. Numer. Meth. Engng. 14 (11) (1998) 1055–1065. [34] S.L. Yap, The Poincaré lemma and an elementary construction of vector potentials, Amer. Math. Monthly 116 (3) (2009) 261–267.