r3 using Gauss quadrature

Efficient evaluation of integrals of order l/r l/r 2, l/r 3 using Gauss quadrature LIU JUN,* G. B E E R a n d J. L. M E E K Department of Civil Engineering, University of Queensland, Australia An efficient method is presented for the numerical evaluation of integrals as they are encountered in the two- and three-dimensional implementation of the Boundary Element Method where the integrands are of 0(1/r), 0(1/r 2) and 0(1/ra). Key Words: Boundary Element Method, ~ntegration

INTRODUCTION The two- and three-dimensional implementation of the Boundary Element Method (BEM) requires the evaluation of integrals where the integrand is 0(l/r), 0(1/r 2) and 0(1/ra), where r is a distance from a source point. For a general implementation using isoparametric boundary elements these integrands have to be evaluated numerically. The Gauss-Quadrature (G-Q) formula ~ has been used almost exclusively in the past) -s For the efficient implementation of the BEM it is very important to select the minimum number of Gaussian points required to give a sufficient precision of integration. This may include the subdivision of the element into a number of sub-elements. This selection of the integration order has been discussed in previous papers) 'a In this paper, we examine these methods in more detail and present a new method for one- and two-dimensional integration. Some examples are presented and comparison is made with the methods given in reference 2.

ONE-DIMENSIONAL GAUSS-QUADRATURE The one-dimensional formula for G-Q can be written as: 1

/(x) dx= y.A./(x.)

(1)

For the simple one-dimensional integration shown in Fig. 1, the value of H can be determined in the following way (here the Iacobian is constant):

1 f=-r

(2N)! L 2N

H > 22Nr2N+l He

e =--

/

4

~

(4,)2N/L2~v

r=R

gives maximum error

If we choose a value of r = R the error will be overestimated and the estimate for the number of integration points will probably be conservative. This is shown in Fig. 2 where the predicted (curve 1) and actual (curve 2) number of integration points required for a given precision as a function of R]L is shown. Obviously, we must work with a value ofr/L = r*/L where R/L
R

. . . .

L

a ~ -

L

(5)

L

In equation (5) all is a translation of the curve towards the N-axis as shown in Fig. 2 (curve 3). Substituting (5) into equation (4) and taking eo = 4(L/4R) 2N, the equation (4) for the function f l = 1/r can be written as:

1

e(ax) = eo

f(x) dx -- ~ Anf(Xn) < e H

(4)

r = R + L gives minimum error

--1

In the above N is the total number of integration (Gaussian) points, x n their co-ordinates and A n the weight coefficients. The number of points needed depends on the required precision of the integration and the nature of the integrand. Stroud and Secrest I give an estimate of the error as:

(3)

T w o extreme values can be substituted for r:

n=l

--1

L r=--(1--~)+R 2

(6)

(2)

n=l

where H~ld2Nf[dxZN[ and e ~ 4[(22N(2N)!) *Visiting scholar, Institute of Geology, Academia Simica, Bering, China. Paper presented at the 7th Boundary Element Forum in Japan sponsored by JASCOME. Discussion closes November 1985.

Figure 1.

One-dimensional element 0264-682Xl85/030118-06

118

Engineering Analysis, 1985, Vol. 2, No. 3

© 1985 CMLPubUcations

$02.00

Efficient evaluation o f integrals using Guass quadrature: Liu Jun, G. Beer and J. L. Meek N 10

I

6

" . CURvz3 2

1.0

20

30

"

Figure 2. Predicted error curve 1 (r*/L = R/L)~ curve 3 (r*]L = R [L + O.3 5) - and actual error (curve 2)

N

Using expressions (6)-(8) we plot diagrams in Figs. 3, 4 and 5 which give the Gaussian points needed for the integration of these functions versus the closeness o f the source point (R/L) for an allowable error e~n) o f 10-a and l0 -4 respectively. In all of the expressions, a l l is taken to be 0.35. For a given ratio R / L the number, N, o f Gaussian points depends on both the chosen allowable error and the nature o f the integrand. To check the equations (6)-(8), as an example, Fig. 4 shows a comparison curve o f the predicted and actual number o f integration points needed to achieve a certain accuracy (Ca(2)= 10-a) for function f2 = 1/ta. We can see that both curves are very close to each other. Therefore using the above equation in the numerical integration, we can obtain a result which meets the accuracy required and is economical in computer resources. Tables 1-3 show the limiting values o f R[L for particular numbers of integration points and the accuracy of inte-

N

10 6

6 E~U a-"

2

10-

10-4

10-2

lO-Z

lO-,a

ZO

" R/L

Figure 3. Predicted number o f integration points, N for function 1/r

F ~ r e S. tion 1/ra

Predicted number o f integration points for func-

Table 1. Number of integration points needed for precision,, e = 10-3

N

f = 1/r

f = l/r 2

f = 1/r3

1.1180 (0.7680) 0.6786 (0.3286) 0.5287 (0.1787) 0.4551 (0.1051) 0.4119 (0.0619) 0.3835 (0.0335) 0.3636 (0.0136) 0.3487 (0.0000) 0.3374

1.6719 (1.3219) 0.9386 (0.5886) 0.6958 (0.3458) 0.5785 (0.2285) 0.5100 (0.1600) 0.4654 (0.1154) 0.4340 (0.0840) 0.4107 (0.0607) 0.3928 (0.0428)

2.2003 (1.8503) 1.1825 (0.8325) 0.8508 (0.5008) 0.6920 (0.3420) 0.5998 (0.2498) 0.5399 (0.1899) 0.4979 (0.1479) 0.4668 (0.1168) 0.4428 (0.0928)

f = 1/r

jr= 1#2

f = 1#3

1.9882 (1.6382) 0.9961 (0.6461) 0.7050 (0.3550) 0.5730 (0.2230) 0.4990 (0.1490) 0.4521 (0.1021) 0.4198 (0.0698) 0.3963 (0.0463) 0.3785 (0.0285)

2.9730 (2.6230) 1.3776 (1.0276) 0.9279 (0.5779) 0.7283 (0.3783) 0.6179 (0.2679) 0.5486 (0.1986) 0.5012 (0.1512) 0.4668 (0.1168) 0.4407 (0.0907)

3.9127 (3.5627) 1.7357 (1.3857) 1.1346 (0.7846) 0.8712 (0.5212) 0.7267 (0.3767) 0.6364 (0.2864) 0.5749 (0.2249) 0.5305 (0;1805) 0.4968 (0.1468)

10 2 3 4 5 6 7 8 9 10

6 10-4

lO-Z_ 10-z :

Table 2.

e = lO-3

Figure 4. Predicted and actual number o f integration points for function 1/r 2

Similarly, the expression for the function 1"2 = 1/r 2 can be written as: e (2) =

(2N +

1) eo

(7)

and for the functionfa = 1/ra: e (3) = ( N + I ) ( 2 N + 1) eo

Where e0), e(2) and e0 ) are allowable errors.

(8)

2 3 4 5 6 7 8 9 10

Efficient evaluation of integrals using Guass quadrature: Liu Jun, G. Beer and J. L. Meek Table 3.

N~

2 3 4 5 6 7 8 9 10

e = 10-"

R/L

f = 1/r

f = 1/r2

f = 1# 3

3.5355(3.1855) 1.4620(1.1120) 0.9402 (0.5902) 0.7213(0.3713) 0.6046(0.2546) 0.5329 (0.1829) 0.4848(0.1348) 0.4504 (0.1004) 0.4247(0.0747)

5.2869 (4.9369) 2.0221(1.6721) 1.2373(0.8873) 0.9168 (0.5668) 0.7486(0.3986) 0.6466 (0.2966) 0.5787 (0.2287) 0.5305(0.1805) 0.4945(0.1445)

6.9479 (6.6079) 2.5477(2.1977) 1.5130(1.1630) 1.0967(0.7467) 0.8804 (0.5304) 0.7502 (0.4002) 0.6639 (0.3139) 0.6028 (0.2528) 0.5575(0.2075)

r--. L~

_Ro L~ ~i" L;~L L

Lo F ~ r e 6.

- I

the limiting value listed in Tables 1-3. Using the expression (11) and Tables 1-3, the number M is obtained. Then the number Nq) for each sub-element is found using (10), (11) and (6), (7), (8). The numerical integration is carried out using equation (9). Two examples are analysed below. The first example is: 1

I=

(1.1 - - x ) a --1

with the allowable error ea(a) = 10-4. The different results are listed in Table 4: We can see that if the integration region is not subdivided (M = 0), the results is obtained according to the limiting value (= 10) in the program, so that the computation error is greater by far than e~a)(= 10-a). On the contrary, i f M = 3 the result is very good. It should be pointed out that the computer time depends on the limitation of integration points used in the program as the integration region has to be subdivided. Let us consider a second example: 1

An element and its sub-elements I--

(1.1 -- x) 2 -1

gration. The values of the R/L shown in the brackets are obtained according to curve 3, r*/L = R[L + 0.35. Otherwise they are according to curve 1, (r*[L = R[L), i.e. the value used in reference 2. For the case where the source points is very near to the integration region and the number of the integration points which can be used in a program is limited, the region has to be divided into some sub-elements. This has been suggested by Watson. 2,a However, the scheme presented here is different from Watson's, as we subdivide this element unequally instead of equally. This is explained for the simple example in Fig. 1. Figure 6 shows an element which is subdivided into four sub-elements. Each sub-element has its own local co-ordinate which varies from -- 1 to + 1. According to Fig. 6, the general formulae can be obtained as follows: +1

M+I N(/) E Hnq)f(~On D)

I= f / ( ~ ) d ~ = E d

--1

(9)

]=1 n = l

The result is listed in Table 5: In Table 5, M, N and e* are the same as the above example. N L is a limitation of the integration points which can be used in the program. We can find that if the largest NL(= 10) is used, the M, N are smaller when e* < e[2). That means computer time is shortened.

Table 4 M

N

e*

0

10

10 -2

3

21

1.5 X10-6

Herein, M + 1 is total number of sub-elements, N is total number of integration points, and e* is computation error

Table 5.

e~~)'~ = 10-2

NL~

M+1

N

e* (%)

6 5 4 4 4 3 3 3 2

18 16 14 14 12 12 12 12 12

0.04 0.05 0.I0 0.10 0.I0 0.40 0.40 0.40 0.40

l

where M + 1 is total number of the sub-elements. For all sub-elements except the last, we have: Lo

3

~ = 2-/~(D =-~(1-- 2/) R~/, _ 1 + 2 / ( ~ o )

(10)

2 3 4 5 6 7 8 9 10

for the last sub-element: 1

TWO-DIMENSIONAL GAUSS-QUADRATURE (11)

LM+I

The two-dimensional G-Q formula can be written as: 1

\Lo]

If the source point is very near to the element, the total number of the integration points needed is more than

f

1

f

-1 -1

f(~(1)' ~(2)) d~(l) d~(2)

Efficient evaluation o f integrals using Guass quadrature: Liu Jun, G. Beer and J. L. Meek N"t N 2

= E E A0)A(2)Cte0)g2)~

(12)

1

l=1]=1

In the above, N.r (7 = 1, 2) is the number of integration points in ~('t) (7 = 1, 2) direction. The number N~ depends on the required precision of integration. The order (N1, N2) of the formula to be used is determined from the upper bound for the error given by Stroud and Secrest:

I

•

"

[ -

"

"

•

•

•

J

•

•

•

•

•

•

•

J

•

•

•

•

•

I

•

•

D

•

•

]

•

•

•

•

•

2

2

IA~I < 2 ~ e,r/-Lr

(13)

3,=1

where


•

•

and

It

4 e v - 22N,r (2N,r)!

,, ,O

The norrnalised error is ea = A~ff. The function 1/r2 is taken to be representative of actual integrands and the case considered is a rectangular region of integration of dimensions L1 xL2 (Fig. 7). Proceeding as before, it is found that N~ should be chosen such that: e(2) = (2N,r + 1)eo

ff

1

for last one: 2\

2 ~/

(17)

R(M'v+1) = 2M7 __R~ L~.y+l) L.v Some examples are now given for the two-dimensional problem. The following examples can be solved exactly in order to check whether the method is correct or not. Firstly, the integral: 1

1

-1

--1

f(/j(1),/j(2)) d/j0) d/j(2)

--1 --1 M,+I =

.,

Distribution o f Gaussian points - example 1

(14)

For the case x a qt Sb, where is a quadrilateral element, N~ is chosen by R/L~ in terms of (14) as for the onedimensional problem. Being similar to the one.dimensional problem, when N~ exceeds the limiting value of the number of Gaussian points which can be used in the program, the region has to be divided into some sub-regions (Fig. 7). Then the integration formula must be written as: 1

2 Figure 8.

M2+I N z N 2 a(1)~(2)¢{~(1)~(2h

E E EE-,-j,,.,,',i-

(15)

!=1 k=l i=1/=1

for all of the sub-regions except the last one, we have: Sj~= 2-(/' k)~r + ~ (1 -- 2 g ~ ) (16)

k)( R'r ) R g ~) L (ct,k) = 1 + 2 (£ k-~-~,/

is computed. The integration region and the source point x a are shown in Fig. 8. The allowable error e~3) is 10-2. The exact solution is 2.327344790. In the process of numerical integration, the region is divided into four sub-regions. Fifty-eight integration points are used which are distributed in Fig. 8. The result is 2.32729505, and computation error e* is 0.21 x 10-3. Obviously, the result is sufficiently accurate. In the second example, the integral is: +1

1

- a t dr/ r -i

I [.2

_i

t !

L1 Figure 7. An element and sub-regions

×a

-i

the integration region is shown in Fig. 9. The exact solution is 3.343673397. We compute this integral according to two different cases: (1) the allowable error is 10-2, the result is 3.34523876 without any subdivision. The number of the integration points is 81 (Fig. 9a), the computation error e* is equal to 0.49 x 10- 3. (2) The allowable error is taken to be 10- 4. In the computation process, the region is divided into ten sub-regions, the total number of the integration points is 187 distribution of which is shown in Fig. 9b the result is 3.34369565 e* = 0.68 x 10- 6. The above results show that this method given by authors is very efficient for the functions 0(I/r), 0(I/r 2)

Efficient evaluation of integrals using Guass quadrature." Liu Jun, G. Beer and J. L. Meek We compute the integral: +1

+1

dx

ff(1.4-x) a

- - 1 --1

with the method given by the authors and in reference 2. The exact solution is 0.9447190. The results are listed in Table 6. The distribution of integration points are shown in

1]t

Fig. 11.

N _

L

10

2 (a)

+ ~

I

9

I I I 19

I

16

~1~1)=

2

io-~

IO-L

16

26 : R/L.

1.o Figure 10.

1116

Table 6

I - - - -I-

Methods's~'-...~

M

N

I

e*

Author's Reference 2

4 9

84 188

0.9441790 0.9441792

0.64 X 10-s 0.20 X10-6

%

I 16 19q16-

16

(b) Figure 9.

~

16

_

16

16

116 I

Distribution o f integration point -- example 2

and 0(1/r a) with different singularity and for different allowable error.

COMPARISON Lachat and Watson 2'a give the following formula for the integral of (1/r2):

(2Np + 1)Otp2NP< K

(18)

where ap = (ds/dr/p) = Lp[4R and K is a constant which is chosen in the light of experience. It should be pointed out that the physical meaning of K is explicit. The value is believed to be dependent on the required integration error and a value of K = 10-s is suggested in reference 5. The plot of expression (18) computed with this value is shown in Fig. 10. Comparing this with the curves if Fig. 5 it is found that the humber of integration points given by equation (18) is always larger than by equation (7) when (R[L) is specified. This means that expression (18) is conservative.

\

.4

Xo

Xa

(a)

(b)

Figure 11. Distribution of Gaussian points: (a) author's method; (b) method o f ref 2

SUMMARY AND CONCLUSIONS The efficient evaluation of integrals where the integrand varies rapidly and is of order (l/r), (1/r2), (1/r a) has been investigated in detail. Integrands of this type occur in the numerical implementation of the BEM method in two- and three-dimensional elastostatics, where r is the distance to a source point.

Efficient evaluation o f integrals using Guass quadrature: Liu Jun, G. Beer and J. L. M e e k

It has been found that in particular for three-dimensional Boundary Element analysis the computational effort and therefore computing time, needed to integrate the KernelShape function products is significant and in some cases may be more significant than the other stages of computation including the solution o f the system o f equations. The method presented in this paper allows to select, from a table, the number o f integration points needed for the numerical integration to achieve a prescribed precision. This number is dependent on the order of the integrand and the proximity o f the source point to the integration region. The values o f the tables are checked by simple examples where the integral can be evaluated analytically and good agreement was found between predicted and actual precision o f the integration. When the source point is very close to the integration region a subdivision into subregions has been suggested as the number o f integration points needed may exceed the number available in the computer program. In the present paper the authors are presenting a new and more efficient method of less integration points for a more accurate result. The methods of integration presented here are expected to make the evaluation o f Boundary Element coefficient matrices more efficient and save a significant amount o f computational effort especially for three-dimensional analysis. They have been implemented in the computer program BEFE which was developed by Beer. s Recent examples on three-dimensional stress analysis problems have given excellent results with regard to accuracy and efficiency. There is a possibility that the integration can be made even more efficient by using coordinate transformations near the singular point. This has been proposed recently by Higashimachi et al. 9 and Mustoe. I° These methods have been investigated in more detail by the authors and are discussed in a sequel to this paper.

ACKNOWLEDGEMENTS The authors would like to thank the Australia/China Council for providing support for Mr Liu Jun.

REFERENCES 1 Stroud, A. H. and Secrest, D. Gaussian Quadrature Formulae, Prentice-Hall, Engelwood Cliffs, NJ, 1966 2 Lachat, J. C. and Watson, J. O. Efffective numerical treatment of boundary integral equation, J. N u n Meth. Eng. 1976, 10, 981 3 Lachat, J. C. and Watson, J. O. In Boundary Integral Equation Method: Computational Application in Applied Mechanics, eds T. A. Cruse and F. J. Rizzo, AMD 11, ASME, New York, 1975 4 Chadouet, A. and Afzal, M. CA, ST, or 3D: Three-dimensional boundary element analysis computer code, in Boundary Elements, Prec. 5th Int. Conf., Hiroshima, eds. C. A. Brebbia, T. Futagami and M. Tanaka, Springer-Vedag, Japan, 1983, pp. 787-809 5 Higashimachi, T., Ezawa, Y., Okamoto, Y. and Kamimiya, Y. Some investigations concerning precision o f three-dimensional analysis by the Boundary Element Method, 60th General Meeting o f Japan Soc. Mech. Engrs, No. 830-1, April, 1983, pp. 7-9 (in Japanese) 6 Cristessu, M. and Loubignac, G. Gaussian quadrature formulae functions with singularities in 1/R over triangles and quadranges, in Recent Advances in Boundary Element Methods, ed. C. A. Brebbia, Pentech Press, London, 1978, pp. 375-390 7 Pina, H. L. G. and Fernandes, J. L. M. Some numerical integration formulae over triangles and squares with a 1]R singularity, Appl. Math. Modelling 1981, 5, 209 8 Beer, G. BEFE- a combined boundary element-f'mite element computer program. Adv. in Eng. Software 1983, 6, 103 9 Higashimachi, T., Okamoto, N., Ezawa Y.I Aizawa, T. and Ire, A. Interactive structural analysis system using the advanced boundary element method, Prec. 5th Int. Conf., Hiroshima, eds C. A. Brebbia, T. Futagami and M. Tanaka, SpringerVerlag, Japan, 1983, pp. 847-856 10 Mustoe, G. G. W. Advanced integration schemes over boundary elements and volume cells for two- and three-dimensional non-linear analysis, in Developments in Boundary Element Methods- 3, eds I. K. Banerjee and S. Mukherjee, Applied Science, London, 1984, pp. 213--269