An h-hierarchical Galerkin BEM using Haar wavelets

Engineering Analysis with Boundary Elements 25 (2001) 581±591

www.elsevier.com/locate/enganabound

An h-hierarchical Galerkin BEM using Haar wavelets Kazuhisa Abe*, Kazuhiro Koro, Kayoko Itami Department of Civil Engineering and Architecture, Niigata University, 8050 Igarashi 2-Nocho, Niigata 950-2181, Japan Received 24 May 2000; revised 2 November 2000; accepted 15 November 2000

Abstract Application of Haar wavelet to h-hierarchical adaptive method is attempted to improve the performance. The wavelet boundary element method (BEM) is developed by means of the Galerkin method. In order to save memory requirement and computation time, sparse matrices are obtained by truncation of matrix coef®cients. Error indicator and estimator are constructed based on the difference between the true solution and its projection onto the mesh. These values are effectively evaluated using the boundary element solution. The adaptive method is developed for 2D Laplace problems. Through numerical examples, performance of the present method is investigated concerning the sparseness of matrices and the computation time devoted to calculation of matrix coef®cients and to equation solving process. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: H-hierarchical method; Haar wavelet; Galerkin BEM; Sparse matrix

1. Introduction In recent years, the boundary element methods (BEM) have been applied to a wide range of engineering problems. Since accuracy of the boundary element solution depends on the element distribution, error estimation and adaptive re®nement will be very useful strategies in practical analysis. The adaptive boundary element methods based on h-, p-, r-, hp-, or hr-versions have been proposed by many researchers [1]. Although these methods allow us to obtain a good solution with reducing the degrees of freedom (DOF), at least a part of the previously constructed matrices must be updated at each re®nement step. It will be a drawback of those methods. Hierarchical method will be a remedy for this problem. Adaptive boundary element methods based on h-hierarchical version have been proposed by Parreira and Dong [2], Guiggiani and Lombardi [3], Chara® et al. [4] and Chara® and Wrobel [5]. They have proved the ef®ciency of the hierarchical method. However even if the h-hierarchical method is employed, coef®cient matrices of the boundary element equations become fully populated ones. It implies that the application of the conventional BEM is dif®cult to large problems. To overcome this dif®culty, in this paper, the Haar wavelet is employed as h-hierarchical basis. Since the wavelets are bases having the property of multiresolution analysis, * Corresponding author. Tel.: 181-25-262-7028; fax: 181-25-262-7021. E-mail address: [email protected] (K. Abe).

they are also available as the h-hierarchical bases. The Haar wavelet is the simplest orthonormal wavelet having vanishing moments of ®rst order. Owing to the properties of vanishing moments and local support of the wavelet, matrices of the boundary element equations have many small matrix coef®cients and few large ones. Hence a large number of matrix coef®cients can be truncated without affecting the accuracy of the solution. Koro and Abe [6] have constructed an h-hierarchical BEM using Haar wavelets in which the boundary element equations are derived based on collocation method. In the collocation method, the wavelet basis is adopted only to the boundary elements but not to the weighting functions. It is dif®cult for this method to save memory requirement and computation time. In order to improve the performance, in this study, employment of the Galerkin method is attempted, so that the wavelet can be applied not only to boundary elements but to weighting functions. Ef®ciency of the adaptive method is investigated through numerical examples.

2. Haar wavelet 2.1. Wavelet Expansion of a function with wavelet basis is achieved based on a scaling function f (x) and a wavelet c (x). The scaling function is de®ned as a function such that integer translates of the scaling function {f (x 2 k)}k[Z form basis

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

582

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

in a subspace V0 , L 2(R). Then the functions

f m;k …x† ˆ 2m=2 f…2m x 2 k†

…m; k [ Z†

…1†

are basis in Vm. Here Vm is the subspace of resolution level m such that ¼ , V21 , V0 , V1 , ¼ , L2 …R†:

…2†

Let us consider an orthogonal decomposition of Vm Vm ˆ Vm21 % Wm21 ;

…3†

where Wm21 is the orthogonal complementary subspace of Vm21. Basis in the subspace Wm is de®ned by the wavelet

c m;k …x† ˆ 2

m=2

m

c…2 x 2 k†

…m; k [ Z†:

…4†

Since V0 and W0 are subspaces of V1, f (x) and c (x) can be expanded using f 1;k …x†…ˆ f…2x 2 k†† 9 1 X > f…x† ˆ ak f…2x 2 k†; > > > = kˆ2 1 …5† 1 > X > > > c…x† ˆ b k f…2x 2 k†: ; kˆ2 1

From Eq. (3) subspaces Wm are orthogonal. The subspace Vm with resolution level m can also be decomposed into a hierarchical form Vm ˆ V0 % W0 % ¼ % Wm21 :

…6†

Fig. 1. The Haar scaling function and wavelet. (a) Scaling function. (b) Wavelet.

Wavelet coef®cients c0,1 and dm,k can be obtained by c0;1 ˆ

Z1 0

f f0;1 dx;

dm;k ˆ

Z1 0

f c m;k dx:

…11†

In Eq. (10) c m,k is generated by dilation and translation of the Haar wavelet c , and the function f is expanded in a hierarchical form. 2.3. Vanishing moment of wavelets A wavelet has a property of vanishing moments such that Z1 x p c…x† dx ˆ 0 …p ˆ 0; 1; ¼; n 2 1†: …12† 21

The Haar wavelet has vanishing moments of ®rst order …n ˆ 1†:

2.2. Haar wavelet

3. Galerkin BEM with Haar wavelet

The Haar wavelet is the simplest orthonormal wavelet. The scaling function is de®ned as ( 1; 0 # x , 1 f…x† ˆ …7† 0; otherwise:

Although very few adaptive BEMs have been developed based on the Galerkin formulation [7,8], the Galerkin method is suitable for wavelet BEM. In the Galerkin version, the wavelet basis is used not only for the boundary elements but the weighting functions. Hence the property of the vanishing moments in wavelets exhibits the validity in double integrations, and the effect of truncation of matrix coef®cients becomes more remarkable than for the collocation method. Therefore the Galerkin method is employed popularly in the wavelet BEM [9±14]. In this paper, we consider 2D Laplace problem. The direct formulation of the boundary integral equation is expressed as Z Z c…x† u…x† 1 q p …x; y† u…y† dG y ˆ u p …x; y† q…y† dG y ;

Note that the function of Eq. (7) is the basis of piecewise constant discretization. The Haar wavelet is de®ned as 8 1; 0 # x , 1=2 > > < c…x† ˆ 21; 1=2 # x , 1 …8† > > : 0; otherwise:

c (x) can also be expressed using the scaling function c…x† ˆ f…2x† 2 f…2x 2 1†

…9†

The Haar scaling function and wavelet are depicted in Fig. 1. A function f(x) [ L 2 [0,1] can be expressed by the Haar wavelet basis, m

f …x† ˆ c 0;1 f0;1 …x† 1

1 X 2 X mˆ0 kˆ1

dm;k c m;k …x†:

…10†

G

G

…x [ G† …13† p

p

where u is the potential and q is the ¯ux. u and q are the fundamental solutions corresponding to potential and ¯ux, respectively, and c(x) is the free-term. G is the boundary, and x and y are points on G .

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

On G , u and q are discretized using the Haar wavelets ~ j† ˆ u…

nb X jˆ1

~ j† ˆ q…

nb X jˆ1

u0; j f0; j …j† 1

q0; j f0; j …j† 1

nm M X X mˆ0 kˆ1 nm M X X mˆ0 kˆ1

Z G

Z G

(14)

~ dG y : u p …x; y† q…y†

…15† In the Galerkin method, boundary element equations are derived from the following equation Z Z r f 0; j dG ˆ 0; r c m;k dG ˆ 0: …16† G

G

Eq. (16) leads to " #( ) H ff Hfc u~ f Hcf

Hcc

u~ c

" ˆ

Gff

Gfc

Gcf

Gcc

#(

q~ f q~ c

hcc

gfc ˆ gcc ˆ

4.1. Error indicator and error estimator For the sake of brevity, we express the bases f 0,i and c m,k as wi's, i.e. w 1 ˆ f0;1 ; ¼; wnb ˆ f0;nb ; wnb 11 ˆ c 0;1 ; ¼:

;

…17†

…20†

Hence the boundary element solutions u~ and q~ are expanded as u~ ˆ

N X i

Z

q~ ˆ

u~i wi ;

N X i

q~ i wi :

…21†

Eq. (16) is expressed as G

rw i dG ˆ 0

…i ˆ 1; ¼; N†:

…22†

Since the true solutions u and q satisfy Eqs. (13) and (22) leads to ZZ Z ~ i dG 1 ~ i dG 2 c…u 2 u†w q p …u 2 u†w G

G

ZZ

~ i dG 2 ˆ 0 up …q 2 q†w

…i ˆ 1; ¼; N†:

…23†

Here we introduce functions u^ and q^ de®ned as u^ ˆ

N X i

q^ ˆ

u^i wi ;

N X i

where Z uw i dG; u^i ˆ G

q^ i wi ;

q^ i ˆ

Z G

…24†

qw i dG:

…25†

^ we can rewrite u 2 u~ and q 2 q~ as Using u^ and q;

ZZ

p

G

2

f 0;i c m;k q dG ;

ZZ 1 ˆ dml dkj 1 c m;k c l; j qp dG 2 ; 2 G

gff ˆ

4. Error estimation and adaptive re®nement

2

)

where u~ f ; u~ c and q~ f ; q~ c are subvectors consisting of the coef®cients associated with the scaling functions f 0,j and wavelets c m,k. The elements of submatrices Hff , Gff ,¼ are de®ned as ZZ 1 h ff ˆ dij 1 f 0;i f0; j qp dG 2 ; …18† 2 G hfc ˆ

…19†

where A is a matrix obtained from H and G, z is an unknown vector and b is a known vector.

q~m;k c m;k …j†;

~ dG y 2 q p …x; y† u…y†

Arranging Eq. (17), we obtain a matrix equation Az ˆ b;

u~m;k c m;k …j†;

where u~ and q~ are approximations of u and q, respectively, and u0, j, u~ m;k and q0, j, q~ m;k are wavelet coef®cients of u~ ~ j is an intrinsic coordinate de®ned on the boundary. and q. M is the number of hierarchical resolution levels and n m ˆ 2m nb : Here nb is the number of subboundaries divided by the scaling functions. The degrees of freedom of the problem is N ˆ 2M11 nb : For the sake of simplicity, the intrinsic coordinate j and x are not distinguished in the following description. Substitution of the approximate solution Eq. (14) into Eq. (13) yields a residual r ~ 1 r…x† U c…x† u…x†

583

ZZ G

ZZ G

ZZ G

^ 1 u^ e ; u 2 u~ ˆ …u 2 u†

u^e U u^ 2 u~ ˆ q^e U q^ 2 q~ ˆ

c m;k c l; j up dG 2 ;

where d ij is the Kronecker delta.

…26†

where

f 0;i f0; j up dG 2 ; f 0;i c m;k up dG 2 ;

^ 1 q^ e ; q 2 q~ ˆ …q 2 q†

N X i N X i

…u^ i 2 u~ i †wi ; …27† …q^ i 2 q~ i †wi :

Since the Haar wavelet is orthonormal basis, u^i 2 u~i and q^ i 2 q~ i in Eq. (27) can be expressed as follows Z Z ~ i d G; ~ i d G: …u 2 u†w q^ i 2 q~ i ˆ …q 2 q†w u^i 2 u~ i ˆ G

G

…28†

584

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

Eq. (28) implies that u^ i 2 u~ i and q^i 2 q~ i correspond to the wavelet coef®cients of error. Therefore u^ e and q^ e can be regarded as the projections of error onto the subspace VM. Substituting Eq. (26) into Eq. (23), we have Z ZZ ZZ cu^ e wi dG 1 q p u^ e wi dG 2 2 u p q^ e wi dG 2 G

G

ˆ2

ZZ G

G

^ i dG 2 1 q p …u 2 u†w

ZZ G

^ i dG 2 : u p …q 2 q†w

Eq. (29) corresponds to Eq. (9) of [15] in which the collocation method is employed. In this paper, the error indicator and estimator are de®ned based on the right-hand side of Eq. (29). u 2 u^ and q 2 q^ are approximated by X X u 2 u^ < q 2 q^ < …30† u j c j ; q j c j ; j

where c j denotes the wavelet at the next hierarchical resolution level. u j and q j are wavelet coef®cients corresponding to c j, Z Z ^ c j dG ˆ u j ˆ …u 2 u† u c j d G; G

q j ˆ

Z G

G

^ c j dG ˆ …q 2 q†

Z G

…31† q c j d G:

The de®nitions of a solution f, f^ and f j c j are depicted in Fig. 2. Substitution of Eq. (30) into Eq. (29) provides an equation Z ZZ cu^ e wi dG 1 q p u^ e wi dG 2 G

2

G

ZZ G

u p q^ e wi dG 2 <

N X j

lij ;

…32†

where

lij U 2

ZZ G

q p wi c j dG 2 u j 1

ZZ G

u p wi c j dG 2 q j :

…33†

The error indicator re¯ecting the quality of approximation on each interval G j …ˆ Supp c j † is de®ned as

lj U Max{ulij u; i ˆ 1; ¼; N}:

Fig. 2. De®nitions of f ; f^ and fj c j :

hU

N X j

…29†

j

Notice that the number of piecewise constant intervals G j's is equal to the degrees of freedom N. The error estimator evaluating global accuracy of the boundary element solution is de®ned as

…34†

lj:

…35†

For the right-hand side of Eq. (32), the following inequality holds X X N N lij # lj : …36† j j Hence, substitution of Eq. (27) into Eq. (32) leads to X X N N aij ej # lj ; …37† j j where aij denotes a coef®cient of the matrix A in Eq. (19), and ej stands for wavelet coef®cient of the error u^ j 2 u~ j or q^ j 2 q~ j : Eq. (37) implies that thePerror of boundary element solution can be estimated by lj : Therefore, in this paper, h de®ned in Eq. (35) is employed as the error estimator, while in many papers the error estimator is de®ned based on the residual r. 4.2. Estimation of error indicator Since uj and q j are wavelet coef®cients of the true solution u and q, the value of l ij cannot be evaluated exactly. Thus these values are approximated using the boundary ~ In order to achieve this, let us element solution u~ and q. consider a wavelet element wi at the highest resolution level with element length h…ˆ Supp wi †: On this element, we approximate the solution u using a linear function u given by interpolation of the boundary element solution u~ (Fig. 3). The linear function u passing the middle points of each piecewise constant interval on the element is expressed as 4u~  ˆ u0 2 pi x; u…x† h h

…38†

Fig. 3. Boundary element solution u, ~ linear interpolation u and correction u~j c j :

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

where u0 is a constant and u~i is the boundary element solution corresponding to the wavelet coef®cient of wi. As depicted in Fig. 2, we introduce a wavelet c j having the next hierarchical resolution level. Wavelet coef®cient u j of the function u associated with the wavelet c j is given by u j ˆ

Z

4u~ Z u~ xc j dG ˆ pi : uc j dG ˆ 2 pi G 2 2 h h G

…39†

It should be pointed out that the term concerning the constant u0 is disappeared owing to the property of vanishing moments of c j. u j de®ned in Eq. (31) is approximated by uj u j < u j ˆ

u~i p : 2 2

…40†

Similarly, approximation of qj can be obtained, q~ q j < pi  : 2 2

4.3. Adaptive strategy The optimum mesh distribution will be a discretization in which every error indicator l j has the same value. Consequently, in the adaptive meshing, elements having larger error indicator will be re®ned. In the h-hierarchical method, the re®nement is carried out by addition of elements. In this study, the following criterion is employed

lj $ b Max{lj };

the accuracy. This property allows us to save the memory and computation time devoted to solving the matrix equation. Besides, if the value of each matrix coef®cient can be estimated and truncated a priori, we can also shorten the computation time consumed for the calculation of matrix coef®cients. In this paper the estimation of the matrix coef®cient is carried out considering the vanishing moments of the Haar wavelet, similarly to the method employed in Ref. [13]. The matrix coef®cients de®ned in Eq. (18) are estimated as h h fc < h fc U 20 22…3k11†=2 ; r hcc < h cc U

h0 2…3…ki 1kj †12†=2 2 ; r3

gfc < gfc U

h20 2…3k11†=2 2 ; r

gcc < g cc U

h20 2…3…ki 1kj †12†=2 2 ; r2

…41†

From Eqs. (40) and (41) the error indicator and estimator are evaluated using the boundary element solution.

0 # b # 1;

…42†

where b is a speci®ed parameter. Element c j satisfying the condition of Eq. (42) is added to the present discretization. The value of b in¯uences the performance of the adaptation. As for larger values of b , the number of elements added at each adaptive stage decreases. Consequently the ®nal element distribution will be almost optimum one, while the number of adaptive stages increases. The accuracy of the boundary element solution is evaluated by the error estimator h . When h becomes smaller than a tolerance e , the adaptive process is terminated

h # e:

…44†

where h0 is the length of element corresponding to the scaling function, and r ˆ r=h0 : Here r is the distance between two elements. k, ki and kj are the hierarchical resolution level of the wavelets. None of the matrix elements hf f and gff which do not include any wavelet bases c m,k in their integrations are truncated. Since the number of these elements is very small, they don't affect the performance of wavelet BEM. The value of matrix coef®cient will be small for large r; and in such cases it can be considered that h ij , g ij ; where hij and gij denote h fc ; h cc and g fc ; g cc : Then the truncation is accomplished based on the value of g ij : When g ij satis®es the condition gij # tgmax ;

…45†

then hij and gij are truncated, that is to say, these values are

…43†

5. Truncation of matrix coef®cients Although the wavelet BEM has fully populated matrices, they have many small coef®cients and few large ones. Because of this, truncation of the matrix coef®cients for a certain threshold provides sparse matrices without affecting

585

Fig. 4. Boundary conditions and initial discretization of Ex.1.

586

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

Fig. 5. Boundary conditions and initial discretization of Ex.2.

set to zero. Here t is a threshold parameter and gmax is the largest value of ugff u 0 s. The matrix elements which don't satisfy Eq. (45) are calculated based on Eq. (18). Furthermore, when these

Fig. 7. Relation between the error estimator b and computation time. Results of the adaptive re®nement are shown for b ˆ 0.1, 0.3, 0.5 and 0.7. Uniform re®nement is represented with broken curve. (a) Ex.1. (b) Ex.2.

values satisfy the condition gij # tgmax ;

Fig. 6. Convergence of solution. Results of the adaptive re®nement are shown for b ˆ 0:1; 0:3; 0:5 and 0:7: Convergence for uniform re®nement is represented with broken curve. (a) Ex.1. (b) Ex.2.

hij # thmax ;

…46†

then they are truncated after the calculation. The above strategy is also employed in the calculation of l ij which is used for construction of the error indicator. Since l ij is composed of hij and gij associated with the wavelet basis at the next hierarchical resolution level, when it is decided to add the element to the existing mesh, hij and gij are used again as the matrix coef®cients in the re®ned boundary element equations. To save the computational effort, a list of elements which are not truncated is made when they are calculated in the evaluation of the error indicator. The list allows omission of a part of the truncation decision for the added elements. Matrix entries having larger values than the threshold are stored in the memory. The boundary element equations with a sparse matrix are solved by the GMRES method [16].

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

587

6. Numerical examples

6.2. In¯uence of b on the performance

6.1. Analytical conditions

Fig. 6 is showing the relation between accuracy (h ) and DOF (N) in Ex.1 and Ex.2 for b ˆ 0:1; 0:3; 0:5 and 0:7: In the ®gures results for uniform re®nements are also shown. Here the truncation of small matrix coef®cients is not achieved. It can be seen that the rate of convergence is increased using the adaptive method, especially for Ex.2. From the ®gures it is found that the value of b does not in¯uence the rate of convergence. Fig. 7 depicts the in¯uence of b on computation time. The performance is not changed for b # 0:3; while the computation time increased with increasing value of b for b . 0.3. As for larger values of b , the number of adaptive steps increases. Then the computation time increases. It should be noted that, Ex.1 with b ˆ 0:7 has longer computation time than uniform meshing. In the following analyses b ˆ 0:1 is employed. The ®nal meshes for tolerance

In order to investigate the validity of the developed method, two potential problems are considered. Analytical conditions are shown in Figs.4 (Ex.1) and 5 (Ex.2). In these ®gures initial meshes are also shown. The initial mesh is corresponding to the distribution of scaling functions. Ex.1 has no singularity, while Ex.2 has singularity in ¯ux at points (0,1) and (0.5,1). In the GMRES(m) method, diagonal preconditioning is applied to the right side. In this study GMRES algorithm is restarted after m ˆ 10 steps. Convergence of iterative solution is evaluated based on the residual vector r k rk ˆ b 2 Azk ;

…47†

where z k is k-th iterative solution. The convergent condition is de®ned as uurk uu # m; uubuu

…48†

where m is a tolerance. In the present analyses m ˆ 1 £ 10210 is used.

Fig. 8. Final discretization for e ˆ 2 £ 1025 : (a) Ex.1. (b) Ex.2.

Fig. 9. Accuracy at each re®nement step with threshold parameters of t ˆ 1 £ 1025 (Case 1), 1 £ 10 26 (Case 2), 1 £ 10 27 (Case 3), 1 £ 10 28 (Case 4) and 1 £ 10 29 (Case 5). Results for full matrix are represented with broken curves. (a) Ex.1. (b) Ex.2.

588

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

e ˆ 2 £ 1025 are shown in Fig. 8. For this accuracy about 30% and 70% of DOF are saved in Ex.1 and Ex.2, respectively, comparing with the uniform re®nement. Mesh distribution is concentrating around a region where the boundary solution has large slope. Since Ex.2 has singularity, concentration of the mesh distribution is obtained at these singular points.

Truncation of matrix coef®cients is achieved based on Eq. (46). In this section, ®ve threshold parameters (t ) of 1 £ 10 25 (Case 1), 1 £ 10 26 (Case 2), 1 £ 10 27 (Case 3), 1 £ 10 28 (Case 4) and 1 £ 10 29 (Case 5) are considered. In¯uence of truncation on the accuracy of boundary element solution is shown in Fig. 9. Although the matrices of wavelet boundary element equations have many small elements, they are fully populated. Therefore the truncation of matrix coef®cients, even if for small ones, in¯uences the boundary

element solution. From Fig. 9, it can be seen that DOF at which reduction of accuracy becomes remarkable is increasing with decreasing value of the threshold parameter. Relation between accuracy and memory requirement for the matrix elements is shown in Fig. 10. The matrix entries are stored in double precision. Since the truncated matrix is a sparse one, to store the elements we need at least two vectors: One is for storing the matrix elements and the other is for address (of row or column) of the elements. Hence larger memory is required for compressed matrix than for full matrix, at small DOF. Although sparseness of the matrix increases as the DOF increases, the accuracy decreases for excessive truncation. In the analyses the memory requirement can be saved at a certain accuracy about 55 and 70% for Ex.1 and Ex.2, respectively. Effect of truncation on computation time for the calculation of matrix coef®cients is shown in Fig. 11. In this study, the values of matrix coef®cients are estimated and truncated a priori. Therefore the computation time devoted

Fig. 10. Relation between accuracy and memory requirement for Case 1±5 and full matrix. (a) Ex.1. (b) Ex.2.

Fig. 11. Relation between accuracy and computation time for calculation of matrix coef®cients for Case 1±5 and full matrix. (a) Ex.1. (b) Ex.2.

6.3. Ef®ciency of truncation

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

Fig. 12. Relation between accuracy and computation time for solution of matrix equation for Case 1±5 and full matrix. (a) Ex.1. (b) Ex.2.

to calculation of matrix coef®cients can be shortened. Fig. 11 shows that the computation time for matrix coef®cients can be shortened by a factor of 1/7 to 1/10. Effect of truncation on computation time for solution of matrix equation is shown in Fig. 12. In this ®gure, results for LU factorization is also shown. Although many conventional hierarchical methods employ LU factorization algorithm as a solver, GMRES becomes faster than the LU factorization for larger DOF. The GMRES algorithm includes a multiplication of matrix and vector. The computation time devoted to each iteration can then be quickened by increasing sparseness of the matrix. As shown in Fig. 12, the truncation can quicken the equation solving process. Fig. 13 depicts the relation between accuracy and total computation time. The computation time of the present method is 1/2 to 1/3 that of conventional method. 6.4. Comparison with collocation method Adaptive BEM with wavelets have been developed by

589

Fig. 13. Relation between accuracy and total computation time devoted to the adaptive re®nement for Case 1±5 and full matrix. (a) Ex.1. (b) Ex.2.

Hall and McKenzie [17] and Koro and Abe [6]. They have employed collocation BEM. In the collocation method, however, the wavelets are applied to boundary element but not to weighting functions. Hence the compression of matrices due to the property of vanishing moments of wavelets is exhibited for boundary elements and not for weighting functions. On the other hand, the Galerkin method adopts the wavelet both to boundary elements and weighting functions. Therefore sparseness of the matrices becomes higher for Galerkin method than for collocation method. Besides of it, as a preconditioner incomplete LU factorization (ILU) is necessary for convergence of GMRES in collocation method [6], while diagonal preconditioning is effective for Galerkin method. Those differences will be re¯ected in the computational cost. Fig. 14 is showing relations between memory requirement and DOF for collocation method and Galerkin method. From the ®gure effect of the truncation can be seen. Computational results are plotted until the accuracy is

590

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

Fig. 14. Memory requirement at each re®nement step (Ex.2). The memory required for storing the address of matrix entries and for preconditioning are included. (a) Collocation method. (b) Galerkin method.

decreased due to the truncation of matrix elements. Hence each data shown in the ®gure has an accuracy comparable to that of full matrix. In the ®gure the memory for preconditioner is included. As mentioned above, the collocation method must use ILU as preconditioner, while the Galerkin method can use diagonal preconditioning. For this reason, memory saving is dif®cult for the collocation method. On the other hand, the Galerkin method can save the memory about 60±70%. Fig. 15 is showing the effect of truncation on the relation between computation time for solution of equation and DOF. In the collocation method, computation time consumed in ILU factorization is longer than that in GMRES. As a result, the solution needs long computation time. From the ®gure, it can be seen that the GMRES with ILU preconditioning has computation time comparable to that of LU factorization method used for the full matrix. On the other hand, since the Galerkin method has matrices with

Fig. 15. Computation time for solution of matrix equation at each re®nement step (Ex.2). (a) Collocation method. LU factorization method is used for the case of full matrix. (b) Galerkin method. LU factorization method and GMRES are used for the case of full matrix.

higher sparseness and does not need ILU, the computation time can be shortened remarkably. 7. Conclusions An h-hierarchical adaptive strategy for boundary element method, using Haar wavelet as a basis, has been developed by means of the Galerkin method. Since the wavelets are hierarchical basis they can be adopted readily to h-hierarchical basis. Estimating the matrix coef®cients a priori and truncating them with small values, we can save the memory requirement and computation time for calculation of matrix elements and for equation solving process. Although these effects have been investigated so far for uniform re®nement with wavelets exclusively, it was found that the wavelet is effective to improve the performance of h-hierarchical adaptive re®nement.

K. Abe et al. / Engineering Analysis with Boundary Elements 25 (2001) 581±591

Ef®ciency of Galerkin method in h-hierarchical wavelet BEM was investigated through comparison with collocation method. In the collocation method, it is dif®cult to obtain a sparse matrix comparable to that of Galerkin method. Besides, the collocation method needs ILU as a preconditioner of GMRES. As a result, the collocation method requires larger memory and longer computation time than the case of Galerkin method.

[8] [9] [10] [11]

References [1] Kita E, Kamiya N. Recent studies on adaptive boundary element methods. Adv Engng Soft 1994;19:21±32. [2] Parreira P, Dong YF. Adaptive hierarchical boundary elements. Adv Engng Soft 1992;15:249±59. [3] Guiggiani M, Lombardi F. Self-adaptive boundary elements with h-hierarchical shape functions. Adv Engng Soft 1992;15:269±77. [4] Chara® A, Neves AC, Wrobel LC. h-Hierarchical adaptive boundary element method using local reanalysis. Int J Numer Meth Engng 1995;38:2185±207. [5] Chara® A, Wrobel LC. h-Hierarchical functions for 2D and 3D BEM. Engng Anal Boundary Elem 1995;16:341±9. [6] Koro K, Abe K. H-hierarchical adaptive BEM with Haar wavelet functions for two-dimensional Laplace problems. In: Brebbia CA, Power H, editors. Boundary elements XXI, 1999. [7] Rank E. Adaptivity and accuracy estimation for ®nite element and boundary integral element methods. In: Babuska I, Zienkiewicz OC, Gago J, Oliveira ER de A, editors. Accuracy estimates and adaptive

[12] [13] [14] [15] [16] [17]

591

re®nements in ®nite element computations. New York: John Wiley, 1986 (Chapter 4). Parreira P. Self-adaptive p-hierarchical boundary elements in electrostatics. In: Brebbia CA, Wendland WL, Kuhn G, editors. Boundary elements IX, 1987. Beylkin G, Coifman R, Rokhlin V. Fast wavelet transforms and numerical algorithms I. Comm Pure Appl Math 1991;44:141±83. Sabetfakhri K, Katehi LPB. Analysis of integrated millimeter-wave and submillimeter-wave waveguides using orthonormal wavelet expansions. IEEE Trans Microwave Theor Techniq 1994;42(12):2412±22. Goswami JC, Chan AK, Chui CK. On solving ®rst-kind integral equations using wavelets on a bounded interval. IEEE Trans Antennas Prop 1995;43(6):614±22. Wang G. Application of wavelets on the interval to numerical analysis of integral equations in electromagnetic scattering problems. Int J Numer Meth Engng 1997;40:1±13. Lage Ch, Schwab Ch. On the implementation of a fully discrete multiscale Galerkin BEM. In: Marchetti M, Brebbia CA, Aliabadi MH, editors. Boundary elements XIX, 1997. Rathfeld A. A wavelet algorithm for the boundary element solution of a geodetic boundary value problem. Comp Meth Appl Mech Engng 1998;157:267±87. Abe K. A new residue and nodal error evaluation in h-adaptive boundary element method. Adv Engng Soft 1992;15:231±9. Saad Y, Schultz MH. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986;7(3):856±69. Hall WS, McKenzie RA. A sparse h-adaptive boundary integral equation solution for the 2D Laplace equation using multi-wavelets. In: Marchetti M, Brebbia CA, Aliabadi MH, editors. Boundary elements XIX, 1997.