Multidimensional interpolation and numerical integration by boundary element method

Engineering Analysis with Boundary Elements 25 (2001) 697±704

www.elsevier.com/locate/enganabound

Multidimensional interpolation and numerical integration by boundary element method Yoshihiro Ochiai* Department of Mechanical Engineering, Kinki University, 3-4-1 Kowakae, 577-8502 Higashi-Osaka City, Japan Received 7 January 1999; revised 20 September 2000; accepted 4 April 2001

Abstract This paper presents a multidimensional interpolation method and a numerical integration for a bounded region using boundary integral equations and a polyharmonic function. In the method using B-spline, points must be assigned in a gridiron layout for two-dimensional cases. In the method presented in this paper, using the polyharmonic function, arbitrary points can be assigned instead of using a gridiron layout, making interpolation an easy process. This method requires the use of boundary geometry of the region and arbitrary internal points. Values at an arbitrary point and the integral value are calculated after solving the discretized boundary integral equations. Numerical integration is performed using this interpolation. In order to investigate the ef®ciency of this method, several examples are given. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Interpolation; Integration; Boundary element method; Computational mechanics; Numerical analysis; Harmonic function; Multiple-reciprocity method

1. Introduction This paper presents a multidimensional interpolation method for distribution in an arbitrary domain and a numerical multiple integration using the interpolation. In general, a multidimensional interpolation is obtained using the B-spline. When using the B-spline, points must be assigned in a gridiron layout for two-dimensional cases as shown in Fig. 1. On the other hand, the radial basis function (radius spline) has been proposed for use in interpolating an arbitrary distribution, and is useful for scattered data [1±3]. Partial differential equations (PDE) are used to generate many kinds of free-form surfaces (smooth curved surface), which are practical examples of smooth interpolation [4]. However, scattered data cannot be used for PDE. Optimization techniques are used to generate a smooth free-form surface, which is a B-spline surface [5]. The author has already proposed another method of using boundary integral equations and polyharmonic functions for two- and threedimensional interpolations [6±10]. In this paper, a numerical multiple integration method is proposed using the above-mentioned multidimensional interpolation method, which is based on boundary integral equations and polyharmonic functions. * Tel.: 181-66-721-2332; fax: 181-66-727-4301. E-mail address: [email protected] (Y. Ochiai).

If a given region can be transformed into several standard regions for which a rule of approximate integration is available, then integration in the given region can be computed numerically. However, it is dif®cult to decompose a bounded region into elementary standard regions [11]. In order to avoid decomposition, multiple integration by sampling is proposed. However, Monte Carlo computation requires excessively long CPU time for convergence of a result [11]. In this paper, multiple integration is numerically calculated by using boundary integral equations, and domain integration is transformed into a boundary integration. One dimension is reduced by this method. This method is based on an improved multiple-reciprocity boundary element method (BEM) for heat conduction analysis with heat generation [6,10]. This integration is useful for boundary integrals in BEM. Speci®cally, this integration is effective for domain integrals in BEM because internal cells are not required. Several examples are given in order to investigate the ef®ciency of this method. 2. Theory 2.1. Interpolation using polyharmonic function First, the one-dimensional case is considered. A linear interpolation, as shown in Fig. 2 by a polygonal line, can

0955-7997/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0955-799 7(01)00056-X

698

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

Fig. 1. Conventional two-dimensional interpolation.

be obtained by 2 7 …1† W1S ˆ 2W2P ;

…1†

2 where 7 …1† ˆ d2 =dx2 : Eq. (1) is the same type of equation as for the steady heat conduction of temperature …W1S † with a point heat source …W2P †: An unknown point heat source (or heat sink) W2P is obtained inversely from Eq. (1) in order to pass through the given points shown in Fig. 2. Many kinds of distributions, including those of sine and exponential function, can be roughly expressed by Eq. (1). Smoother interpolation, as shown in Fig. 3, can be obtained by 2 7 …1† W1S ˆ 2W2S ; 2 7 …1† W2S

ˆ

…2†

2W3P :

…3†

From Eqs. (2) and (3), the following equation can be obtained: 4 7 …1† W1S ˆ W3P :

…4†

This equation is similar to the equation for the deformation W1S of a beam with an unknown point load …W3P †: The natural spline, which is expressed by Eq. (4), results from a deformation of a simply supported thin beam in order to interpolate the one-dimensional distribution. The term W2S corresponds to the moment of the beam, and the term W2S at each end of the beam is assumed to be 0. If the conventional B-spline is used for two-dimensional interpolation, points must be assigned in a gridiron layout, as shown in Fig. 1. Herein, the deformation of the imaginary thin plate is utilized in order to interpolate the distribution of W1S in the two-dimensional case. The following equations

Fig. 2. Linear interpolation.

Fig. 3. Beam for one-dimensional interpolation.

can be used for the two-dimensional interpolation: 7 2 W1S ˆ 2W2S ;

…5†

7 2 W2S ˆ 2W3P :

…6†

The term W2S of Eq. (5) corresponds to the sum of 2 W1S =2x2 and 22 W1S =2y2 : The term W3P is the unknown strength of a Dirac-type function. From Eqs. (5) and (6), the following equation can be obtained: 2

7 4 W1S ˆ W3P :

…7†

This equation corresponds to the following equation for the deformation W1S of a thin plate with a point load P: 7 4 W1S ˆ

P ; D

…8†

where Poisson's ratio n ˆ 0 and the ¯exural rigidity D ˆ 1: In this paper, the deformation of the imaginary thin plate is utilized in order to interpolate the two-dimensional distribution W1S : The deformation wS1 is given, but the force of the point load P is unknown. The force of the point load P is obtained inversely from the deformation wS1 of the thin plate, as shown in Fig. 4. The term W2S corresponds to the moment of the thin plate. The imaginary moment W2S on the boundary is assumed to be 0, which is the same as that of the natural spline. This indicates that the thin plate is simply supported. In order to interpolate W1S by means of the boundary integral equations, the harmonic function T1 …P; Q† in the

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

699

Generally, distributed function W1S …x; y† for the twodimensional case is assumed to satisfy Poisson's equation 7 2 W1S ˆ 2W2S 2 W2L 2 W2P 2 W2D ;

where the function W2S is a distributed one. The function W2P is a Dirac-type function and is de®ned as

Fig. 4. Imaginary thin plate for two-dimensional interpolation.

W2P ˆ lim W2A A:

  1 1 T1 …P; Q† ˆ ln ; 2p r

…9†

where r is the distance between the observation point p and the loading point q. In addition, the biharmonic function T2 …P; Q† in the steady state is used     r2 1 ln 11 r 8p

…10†

The biharmonic function T2 …P; Q† is de®ned by 7 2 T2 …p; q† ˆ T1 …p; q†:

…11†

On the other hand, the quantity of W1S is given by Green's second identity and Eqs. (5), (6) and (11) as follows: cW1S …P† ˆ 2 (

2 X

…21†f

Z

f ˆ1

2

mˆ1

…12†

where the number of W3P is M. Moreover, W2S in Eq. (6) is similarly given by cW2S …P† ˆ

( G

1

T 1 …P; Q† M X

mˆ1

2W2S …Q†

P T 1 W3…m†

2n

)

2

2T1 …P; Q† S W2 …Q† dG 2n

…17†

The term Wf 11 of Eq. (17) corresponds to the sum of curvatures 22 WfS =2x2 and 22 WfS =2y2 : In this method, the distribution is assumed to be that of an free-form surface S [6±9]. The high-order sum of curvatures WF11 ˆ 0 is assumed to interpolate the distribution, and the sum of curvatures WFS may be assumed to approximately satisfy the following equation: L P D 2 WF11 2 WF11 : 7 2 WFS ˆ 2WF11

…18†

7 2 Tf ˆ Tf 21 :

cW 1S …P† ˆ 2

…19†

This function corresponds to the biharmonic function given by Eq. (10).

…21† f

f ˆ1

Z G

) 2WfS …Q† 2Tf …P; Q† S 2 Wf …Q† dG £ T f …P; Q† 2n 2n 2

…14†

F X

(

(13)

Integral Eqs. (12) and (13) are used to interpolate the distribution W1S : The thin plate spline F…p; q†; which is used to make a free-form surface, is de®ned as [1±3]   1 F…p; q† ˆ r 2 ln : r

7 2 WfS ˆ 2WfS11 2 WfL11 2 WfP11 2 WfD11

Let the number of the function WfP11 be M and the shape of the line function WfL11 be G L. Using Green's theorem, the next equation is obtained:

P T 2 W3…m† ;

Z

This function W2P corresponds to the point heat source in heat conduction theory. Similarly, W2L corresponds to the line heat source and W2D to the double-layer source. The next equation is given as

On the other hand, the polyharmonic function Tf is de®ned by

G

) 2WfS …Q† 2Tf …P; Q† S £ T f …P; Q† 2 Wf …Q† dG 2n 2n . M X

…16†

A!0

steady state is used

T2 …P; Q† ˆ

…15†

F X

…21† f

mˆ1

f ˆ1

2

F X f ˆ1

M X

…21† f

Z GD

Tf WfP11…m† 2

F X f ˆ1

2Tf D W dG ; 2n f 11 D

…21† f

Z GL

Tf WfL11 dG L

(20)

where c ˆ 0:5 on the smooth boundary and c ˆ 1 in the domain. The notation G represents the boundary. The notations p and q become P and Q on the boundary. Moreover,

700

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

W fS … f . 1† in Eqs. (17) and (18) is similarly given by cWfS …P† ˆ

F X

…21†e2f

£ T e2f 11 …P; Q† F X

F X

mˆ1

…21†e2f

Z GL

eˆf

1

F X

2WeS …Q†

M X

…21† e2f

eˆf

1

G

eˆf

(

1

Z

…21†e2f

2n

2

2Te2f 11 …P; Q† S We …Q† dG 2n

P Te2f 11 We11…m†

GD

…21† The arbitrary function W1S is interpolated by the above integral equations. The polyharmonic function Tf in n space dimensions is obtained as  Z 1 Z n21 Tf ˆ …22† r Tf 21 dr dr: r n21 The polyharmonic function of order f Tf for the twodimensional case and its normal derivative are generally given by [6,7,12] "   # fX 21 r 2… f 21† 1 1 ; …23† Tf ˆ ln 1 sgn… f 2 1† r e 2p‰2f 2 2†!!Š2 eˆ1

2f 24

2Tf …2f 2 3†r 2r : ˆ 2n 4p…2f 2 2†! 2n

…24†

…25†

…26†

The polyharmonic function Tf in the one-dimensional problem and its normal derivative are generally given by 2r 2f 21 ; Tf ˆ 2…2f 2 1†! 2Tf 2r 2f 22 2r : ˆ 2n 2…2f 2 2†! 2n

) 2T f …P; Q† 2WfS …Q† 22 Tf …P; Q† S Wf …Q† dG 2 £ 2xi 2n 2xi 2n F X

…21† f

2

F X

…21† f

f ˆ1

2

F X f ˆ1

…21† f

M X 2Tf P W 2xi f 11…m† mˆ1

Z GL

Z GD

…29†

2Tf L W dG 2xi f 11 L 22 Tf D W dG ; 2xi 2n f 11 D

where r 2f 23 r;i 2Tf ˆ 2xi 2p‰2f 2 2†!!Š2 " #   fX 21 1 1 £ 2… f 2 1† ln 2 1 1 2… f 2 1† r e eˆ1 22 Tf r 2f 24 ˆ 2xi 2n 2p‰…2f 2 2†!!Š2

("

2r ni 1 2… f 2 2†r;i 2n

…30†

#

"

where sgn( ) is the sign function, …2f †!! ˆ 2f …2f 2 2†¼4 £ 2 and 0!! ˆ 1 [13]. In the three-dimensional case, polyharmonic function Tf and its normal derivative are given by r 2f 23 Tf ˆ 4p…2f 2 2†!;

(

f ˆ1

2Te2f 11 D We11 dG D … f ˆ 2; 3; ¼; F†: 2n

2Tf r 2f 23 2r ˆ 2n 2p‰2f 2 2!!Š2 2n " #   fX 21 1 1 ; £ 2… f 2 1† ln 2 1 1 2… f 2 1† r e eˆ1

F Z X 2W1S …P† ˆ2 …21† f 2xi G f ˆ1

2

L Te2f 11 We11 dG L

Z

eˆf

)

In CAD, derivative values of the surface are important [14]. The gradient of the distribution can be expressed as follows:

…27†

…28†

# ! fX 21 1 1 £ 2… f 2 1† ln 2 1 1 2… f 2 1† r e eˆ1 ) 2r ; 2 2… f 2 1†r;i 2n

…31†

where r;i ˆ 2r=2xi : The functions 2Tf =2n and 22 Tf =2xi 2n in the three-dimensional problem are generally given by 2Tf …2f 2 3†r;i r 2f 24 ˆ 2xi 4p…2f 2 2†!;

…32†

  22 Tf …2f 2 3†r 3f 25 2r ˆ ni 1 …2f 2 5†r;i : 2n 2xi 2n 4p…2f 2 2†!

…33†

The functions in one-dimensional problem are generally given by 2Tf r;i r 2f 22 ; ˆ 2xi 2…2f 2 2†!

…34†

  22 Tf r 2f 23 2r ˆ n 1 …2f 2 1†r;i : 2n 2xi 2n 2…2f 2 2†! i

…35†

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

2.2. Numerical integration Numerical integration for a function de®ned in an arbitrary domain is presented. After interpolation using Eqs. (20) and (21), numerical integration in domain V is obtained using interpolated values. De®ning a new function c 1 ˆ 1; the function c f is introduced. 2

7 c f 11 ˆ c f :

…36†

Using Eq. (36) and Green's theorem, the next equation is obtained because the function W1S satis®es Eqs. (15)±(18). Z V

W 1S …p† dV ˆ

F X

…21† f

Z

f ˆ1

G

(

) 2WfS …Q† 2c f 11 …p; Q† S £ c f 11 …p; Q† 2 Wf …Q† dG…Q† 2n 2n 2

FX 11

…21† f

Z GL

f ˆ2

2

FX 11

…21† f

mˆ1

f ˆ2

2

FX 11

M X

…21† f

Z

f ˆ2

GD

…44†

2c f r 2f 23 2r : ˆ 2n …2f 2 3†! 2n

…45†

Internal cells are used to obtain the domain integral in conventional BEM. Using this integration method, internal cells are not required. The domain integral for steady heat conduction with distributed heat source W1S is transferred to the boundary integral [6] Z V

T 1 …p; q†W1S …q† dV ˆ

c f …p; q†WfP…m†

1

M X mˆ1

f 2f 22 ; ‰2f 2 2†!!Š2

2c f r 2f 23 2r : ˆ 2n …2f 2 2†!!…2f 2 4†!! 2n

(37)

…39†

…40†

Using Eqs. (37) and (41), the mean value of the distribution can be obtained. In the three-dimensional problem, the function c f and its normal derivative are given by

2c f 2… f 2 1†r 2f 23 2r : ˆ 2n …2f 2 1†! 2n

…21† f

f ˆ1

Z G

…42†

P T f 13 …p; q†W3…m† ;

(46)

where W1S is interpolated using Eqs. (12) and (13). 2.3. Numerical procedure of integral equation For practicality, Eqs. (12) and (13), where F ˆ 2; can be used to interpolate a two-dimensional smooth distribution, instead of Eqs. (20) and (21) [6±10]. It is very dif®cult to solve Eqs. (20) and (21), because they have too many unknown terms. Therefore, the next equation is used for the integration instead of Eq. (25) Z V

Using Eq. (37), the area of region V is obtained as follows: Z 2c …p; Q† Z 2 dG: dV ˆ …41† 2n V

f 2f 22 ; …2f 2 1†!

2 X

) 2WfS …Q† 2Tf 11 …P; Q† S £ T f 11 …P; Q† 2 Wf …Q† dG…Q† 2n 2n

In the case of the two-dimensional problem, the function c f and its normal derivative are given by

cf ˆ

r 2f 22 ; …2f 2 2†!

cf ˆ

(

The above integral does not depend on the position of point p. The function c f in n dimensions is obtained as  Z 1 Z n21 cf ˆ c dr dr: …38† r f 21 r n21

cf ˆ

In the case of the one-dimensional problem, the function c f and its normal derivative are given by

c f …p; q†WfL dG L

2c f …p; q† D Wf d G D : 2n

701

W 1S …p† dV ˆ

2 X

…21† f

f ˆ1

Z G

(

) 2WfS …Q† 2c f 11 …P; Q† S £ c f 11 …P; Q† 2 Wf …Q† dG…Q† 2n 2n 1

M X mˆ1

P c 3 …p; q†W3…m† :

In practice, W1S in Eqs. (12) and (13) is given, however, 2WfS =2n; and W3P must be obtained numerically. Constant elements are used for the boundary. The upper index P corresponds to W3P : Replacing WfS and 2WfS =2n by vectors Wf and Vf, respectively, and discretizing Eq. (12), we obtain

W2S ,

H1 W 1 ˆ G1 V1 1 H2 W 2 2 G2 V2 2 GP2 W P3 ; …43†

(1)

…48†

where H1, G1, H2, G2 and GP2 are the matrices with the

702

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

ing elements: H3ij ˆ G3ij ˆ

H4ij ˆ G4ij ˆ

Z Gj

Z Gj

Z Gj

Z Gj

2T1 …pP ; Q† dG j ; 2n

…57†

T1 …pP ; Q† dG j ;

…58†

2T2 …pP ; Q† dG j ; 2n

…59†

T2 …pP ; Q† dG j ;

…60†

GP3ij ˆ T2 …pP ; qP †:

Fig. 5. Interpolation using integral equations (Human face). (a) Given data, (b) Obtained shape.

following elements for a given boundary point i [15]: Z 2T …P; Q† 1 1 dG j ; …49† H1ij ˆ dij 1 2 2n Gj G1ij ˆ

H2ij ˆ

G2ij ˆ

Z Gj

Z Gj

Z Gj

T1 …P; Q† dG j ;

…50†

2T2 …P; Q† dG j ; 2n

…51†

T2 …P; Q† dG j ;

…52†

GP2ij ˆ T2 …P; qP †:

…53†

For F ˆ 2; W2 is obtained by Eq. (13) as follows: H1 W 2 ˆ G1 V2 1 GP1 W P3 ;

…54†

where GP1 is a matrix with the following elements: GP1ij ˆ T1 …P; qP †:

…55†

Moreover, using internal points with the value of function W…pP †; from Eq. (12) we obtain W…pP † ˆ 2H3 W 1 1 G3 V1 1 H4 W 2 2 G4 V2 2 GP3 W P3 ; …56† where H3, G3, H4, G4 and GP3 are the matrices with follow-

…61†

Assuming W 2 ˆ 0 in order to interpolate the distribution, the following equation is obtained by Eqs. (48), (54) and (56): 2 32 3 2 3 V1 H1 W 1 G1 2G2 2GP2 6 76 7 6 7 6 76 7: …62† V 7ˆ6 0 6 0 G1 GP1 7 5 4 54 2 5 4 P P P H3 W 1 1 W…p † W3 G3 2G4 2G3 From Eq. (62), we can obtain V1, V2 and W P3 : In the same manner, V2 ˆ 0 can be assumed on a symmetric axis in order to interpolate the symmetric distribution. On the other hand, contour lines of distribution can be used for W3L [6,7]. When using the constant elements and dividing the boundary into N0 elements and N1 internal points, the simultaneous linear algebraic equations with …2N0 1 N1 † unknowns must be solved. 3. Numerical example Examples for this interpolation using the boundary value and internal points are shown. Fig. 5(a) shows the boundary geometry and internal points of a human face. Fig. 5(b) shows a surface interpolated by using Eq. (62). Using this method, the obtained free-form surface exactly passes through the given boundary geometry and internal points. The integral value obtained by Eq. (47) is 1.688 £ 10 6 mm 3. The next example is the face of a Buddha statue, as shown in Fig. 6(a). Fig. 6(b) shows an interpolated surface. The integral value obtained by Eq. (47), which is a volume, is 4.0822 £ 10 5 mm 3. An example of numerical integration, which has an exact value, is shown. To com®rm the accuracy of the present method, the integral of the distributed function in two dimensions for a square region with a side L ˆ 10 mm; as shown in Fig. 7, is obtained. The distributed function is given by ZL ZL  px   py  Vˆ sin sin dx dy: …63† L L 0 0

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

703

Fig. 6. Interpolation using integral equations (face of Buddha statue). (a) Given data, (b) obtained shape.

Fig. 7(a) and (b) shows the boundary shape and internal points used for interpolation. The length L ˆ 10 mm is assumed. The exact value obtained from Eq. (63) is 40.528, and the numerical value obtained by Eq. (47) is 40.531 in the case of Fig. 7(a) and 40.516 in the case of Fig. 7(b). Constant boundary elements are used for these calculations. Next, the volume of a semisphere is obtained. Fig. 8 shows boundary and internal points. For the radius R ˆ 10 mm; the exact volume is 2094.4 mm 3, and the calculated numerical value is 2092.0 mm 3. To con®rm the accuracy of the present method in a threedimensional case, the integral in a cubic region with a side L is obtained as follows: Vˆ

Fig. 7. Region and internal points. (a) 81 points, (b) 60 points.

ZL ZL ZL 0

0

0

 sin

     px py pz sin sin dx dy dz: …64† L L L

Fig. 8. Volume of semisphere.

704

Y. Ochiai / Engineering Analysis with Boundary Elements 25 (2001) 697±704

method, it has been shown that it is possible to express the distributed function by an integral equation and polyharmonic function. It has been demonstrated that the integration using this interpolation is useful for a distributed function in an arbitrary domain, because this method does not require decomposition of a bounded region into elementary standard regions. Numerical examples have shown that the presented interpolation and integral method will be useful in industrial ®elds.

References

Fig. 9. Triple integral. (a) Boundary discretization, (b) internal points.

F ˆ 3 is used to interpolate a three-dimensional distribution, and the area integral is used for boundary G . Fig. 9(a) and (b) shows the discretized boundary and internal points, respectively. The length L ˆ 10 mm is assumed. The exact value obtained using Eq. (64) is 258.0, and the calculated numerical value is 256.4. Constant boundary elements are used. A Gauss integral is used for area integral. Constant boundary elements are used in order to easily produce a three-dimensional program. A higher accuracy of numerical integration is achieved by using higher order boundary elements. 4. Conclusions By improving the multiple-reciprocity boundary element

[1] Micchelli CA. Interpolation of scattered data. Constr Approx 1986;2:12±22. [2] Dyn N. Interpolation of scattered data by radial functions. In: Chui CK, Schumaker LL, Utreras FI, editors. Topics in multivariate approximation. London: Academic Press, 1987. p. 47±61. [3] Powell MJ. The theory of radial basis function approximation in 1990. Light W, editor. Adv Numer Anal 1992;II:105±210. [4] Bloor MIG, Wilson MJ. Using partial differential equations to generate free-form surfaces. Comput-Aided Des 1990;22(4):202±12. [5] Welch W, Witkin A. Variational surface modeling. Comput Graphics 1992;26(2):157±66. [6] Ochiai Y, Sekiya T. Steady heat conduction analysis by improved multiple-reciprocity boundary element method. Engng Anal Bound Elem 1996;18:111±7. [7] Ochiai Y, Sekiya T. Generation of free-form surface in CAD for dies. Adv Engng Software 1995;22:113±8. [8] Ochiai Y. Generation method of distributed data for FEM analysis. JSME Int J 1996;39(1):93±98. [9] Ochiai Y, Yasutomi Z. Improved method generating a free-form surface using integral equations. Comput Aided Geometric Des 2000;17:233±45. [10] Ochiai Y, Sekiya T. Steady thermal stress analysis by improved multiple-reciprocity boundary element method. J Thermal Stress 1995;18(6):603±20. [11] Davis PJ, Rabinowitz P. Methods of numerical integration. London: Academic Press, 1984. p. 344±417. [12] Nowak AJ, Neves AC. The multiple reciprocity boundary element method. Southampton, Boston: Computational Mechanics Publication, 1994. [13] Abramowitz M, Stegun A. Handbook of mathematical functions. New York: Dover, 1970. p. 258. [14] Roth SD. Ray casting for modeling solids. Comput Graphics Image Process 1982;18:109±44. [15] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniquestheory and applications in engineering. Berlin: Springer, 1984. p. 46±70.