A Domain Decomposition Method Based on BEM and FEM for Linear Exterior Boundary Value Problems

Journal of Mathematical Analysis and Applications 262, 70᎐86 Ž2001. doi:10.1006rjmaa.2001.7537, available online at http:rrwww.idealibrary.com on

A Domain Decomposition Method Based on BEM and FEM for Linear Exterior Boundary Value Problems1 Gabriel N. Gatica Departamento de Ingenierıa Uni¨ ersidad de Concepcion, ´ Matematica, ´ ´ Casilla, 160-C, Concepcion, ´ Chile E-mail: [email protected]

George C. Hsiao Department of Mathematical Sciences, Uni¨ ersity of Delaware, Newark, Delaware 19716 E-mail: [email protected]

and Mario E. Mellado Departamento de Matematica, Uni¨ ersidad del Bıo-Bıo, ´ ´ ´ Casilla 5-C, Concepcion, ´ Chile E-mail: [email protected] Submitted by Hans-Gorg ¨ Roos Received April 3, 2000

We develop the finite dimensional analysis of a new domain decomposition method for linear exterior boundary value problems arising in potential theory and heat conductivity. Our approach uses a Dirichlet-to-Neumann mapping to transform the exterior problem into an equivalent boundary value problem on a bounded domain. Then the domain is decomposed into a finite number of annular subregions and the local Steklov᎐Poincare ´ operators are expressed conveniently either by BEM or FEM in order to obtain a symmetric interface problem. The global Steklov᎐Poincare ´ problem is solved by using both a Richardson-type scheme and the preconditioned conjugate gradient method, which yield iteration-by-subdomain algorithms well suited for parallel processing. Finally, contractivity results and finite dimensional approximations are provided. 䊚 2001 Academic Press 1 This research was partially supported by CONICYT-Chile through the Program A on Numerical Analysis of the FONDAP in Applied Mathematica and by the Universidad de Concepcion ´ through the Direccion ´ de Investigacion, ´ Escuela de Graduados and Facultad de Ciencias Fisicas y Matematicas. ´

70 0022-247Xr01 $35.00 Copyright 䊚 2001 by Academic Press All rights of reproduction in any form reserved.

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

71

1. INTRODUCTION In this work we develop and analyze a new domain decomposition method for linear exterior boundary value problems appearing in physics and engineering sciences. In order to solve this kind of problem, the well known coupling of boundary element methods ŽBEM. and finite element methods ŽFEM. has been successfully applied during the last two decades Žsee, for instance, w4, 7, 9, 11, 13, 18, 27x and the references therein .. In addition, the so-called uncoupling procedure, which makes use of a Dirichlet-to-Neumann mapping expressed in terms of the hypersingular integral operator, has been shown to be a very effective simplification of the coupling technique Žsee, e.g., w8, 17x.. We apply here the combination of Dirichlet-to-Neumann mappings Žobtained by the uncoupling procedure. and domain decomposition methods Žsee w21, 24x. to solve the exterior Dirichlet problem for the Laplacian in the plane. Although this is a simple model problem in potential theory, it is of sufficient interest and generality to illustrate the main aspects of our approach. Indeed, as shown in w20x, our method is also applicable to a larger class of exterior problems, including the Lame ´ system in elastostatics. However, instead of the uncoupling procedure, one needs to use a Dirichlet-to-Neumann mapping based on Fourier series developments Žsee w5, 14x.. The use of domain decomposition methods to solve linear and nonlinear boundary value problems has recently gained considerable attention due principally to the fast development of the parallel processing field. In fact, its inherent parallelism and its ability to handle complex geometries and discontinuous coefficients equations are very attractive features to be taken into account. A rather complete discussion on these methods, based principally on FEM, is given in w19, 21, 22, 24x. Also, in w6, 10x, we have successfully applied FEM-based domain decomposition techniques to linear exterior problems in the framework of Steklov᎐Poincare ´ operators Žsee w1, 21x.. We propose here the alternative use of boundary integral operators and solution operators of domain based variational problems in order to express the corresponding local Steklov᎐Poincare ´ operators. At the discrete level, this corresponds to a combination of BEM and FEM. We recall that the use of BEM in Steklov᎐Poincare ´ operators has already been applied in several papers Žsee, e.g., w3, 15, 16, 23, 25, 26x.. As usual, our final goal is to obtain domain decomposition algorithms that can be naturally implemented on a parallel computer. Our emphasis is on the iterative solvers of the Steklov᎐Poincare ´ problem and the corresponding finite dimensional approximations. The rest of this paper is organized as follows. In Section 2 we state some previous results and notations. Next, in

72

GATICA, HSIAO, AND MELLADO

Section 3 we discuss the iterative solvers, present the algorithms, and state their contractivity properties. Finally, the finite dimensional approximations are discussed in Section 4.

2. PRELIMINARIES We first recall our model problem. Let D be a bounded and simply connected domain in R 2 with smooth boundary ⌫D . Then, given f g L2 ŽR 2 y D . with compact support, g 0 g H 1r2 Ž ⌫D ., and b g R, we seek u such b that y⌬u s f in R 2 y D, u s g 0 on ⌫D , and uŽ x . s ␲ log < x < q O Ž1. as < x < ª q⬁. We also introduce a circle ⌫N , of radius r, whose interior region contains D, such that the support of f lies inside the annular region ⍀ bounded by ⌫D and ⌫N . Then we apply the boundary integral equation method as in w8, 17x to derive a Dirichlet-to-Neumann mapping on ⌫N . In this way, the variational formulation of the exterior problem is set on the bounded region ⍀ and becomes: find u g H 1 Ž ⍀ . such that u s g 0 on ⌫D and b

H⍀ ⵜu ⭈ ⵜ¨ dx q 2² ¨˙, Vu˙: s H⍀ f¨ dx q ␲ r H⌫ ¨ ds

Ž 1.

N

for all ¨ g H ⌫1DŽ ⍀ ., where H ⌫1DŽ ⍀ . [  ¨ g H 1 Ž ⍀ . : ¨ s 0 on ⌫D 4 . Here, ² ⭈ , ⭈ : denotes the duality pairing between Hy1 r2 Ž ⌫N . and H 1r2 Ž ⌫N . with respect to the L2 Ž ⌫N .-inner product, the dot indicates tangential derivative along ⌫N , and the mapping V : Hy1 r2 Ž ⌫N . ª H 1r2 Ž ⌫N . is the boundary integral operator of the single layer potential, that is, V␴ Ž x . [

H⌫ E Ž x, y . ␴ Ž y . ds

y

N

1 for all x g ⌫N , for all ␴ g Hy1 r2 Ž ⌫N ., where E Ž x, y . [ y 2␲ log < x y y < is the fundamental solution of the two-dimensional Laplacian. We now let p ) 2 be a given integer. Then, in order to define the interface problem, we let ⌫j , j s 1, p y 1 be polygonal closed curves contained in ⍀, such that the interior region of ⌫j contains both D and the interior region of ⌫jy1 , and such that they split ⍀ into p subdomains ⍀ j , j s 1, p. In other words, ⍀ j is the annular region bounded by ⌫jy1 and ⌫j . Here, we adopt the notation ⌫0 s ⌫D and ⌫p s ⌫N . In addition, let ␯ i , i s 1, p, be the unit outward normal to ⭸ ⍀ i .

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

73

Then, we define linear extension operators R i , i s 1, p, as follows: 䢇

R1 : H 1r2 Ž ⌫1 . 2 ␭1 ª R1Ž ␭1 . g H 1 Ž ⍀ 1 . is the solution of y⌬ R1 Ž ␭1 . s 0 R1 Ž ␭1 . s 0 on ⌫D

in ⍀ 1 , R1 Ž ␭1 . s ␭1

and

on ⌫1 .

For all i s 2, p y 1, R i : H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i . 2 Ž ␭ iy1 , ␭ i . ª Ž R i ␭ iy1 , ␭ i . g H 1 Ž ⍀ i . is the solution of 䢇

y⌬ R i Ž ␭ iy1 , ␭ i . s 0

in ⍀ i ,

R i Ž ␭ iy1 , ␭ i . s ␭ iy1 on ⌫iy1 , R i Ž ␭ iy1 , ␭ i . s ␭ i on ⌫i . 䢇

R p : H 1r2 Ž ⌫py1 . 2 ␭ py1 ª R p Ž ␭ py1 . g H 1 Ž ⍀ p . is the solution of y⌬ R p Ž ␭ py1 . s 0

in ⍀ p ,

R p Ž ␭ py1 . s ␭ py1

⭸ ⭸␯

on ⌫py1 ,

R p Ž ␭ py1 . s y2W Ž R p Ž ␭ py1 . .

on ⌫N .

Here, ␯ denotes the unit outward normal to ⌫N and W : H 1r2 Ž ⌫N . ª H ⌫N . is the hypersingular boundary integral operator defined by y1 r2 Ž

W␭Ž x . [ y

⭸ ⭸␯x

H⌫

N

⭸ ⭸␯ y

E Ž x, y . ␭ Ž y . ds y

for all x g ⌫N , for all ␭ g H 1r2 Ž ⌫N ., where ␯z stands for the unit outward normal at z g ⌫N . We note that the above boundary value problems define harmonic extensions to ⍀ i , i s 1, p, of the unknown Dirichlet data on the interfaces. Now, we introduce the bilinear forms A i : H 1 Ž ⍀ i . = H 1 Ž ⍀ i . ª R, i s 1, p, and A : H 1 Ž ⍀ . = H 1 Ž ⍀ . ª R, defined, respectively by A i Ž z i , ¨ i . s H⍀ i ⵜz i . ⵜ ¨ i dx, i s 1, p y 1, A p Ž z p , ¨ p . s H⍀ p ⵜz p ⭈ ⵜ ¨ p dx q p 2² ¨˙p , Vz˙p :, and AŽ z, ¨ . s Ý ks1 A k Ž z k , ¨ k . for all z, ¨ g H 1 Ž ⍀ ., where z k [ z < ⍀ k and ¨ k [ ¨ < ⍀ k . With this notation, the unique solution u g H 1 Ž ⍀ . of Ž1. satisfies u s g 0 on ⌫D and AŽ u, ¨ . s F Ž ¨ . for all ¨ g H ⌫1DŽ ⍀ ., b where F Ž ¨ . [ H⍀ f¨ dx q ␲ r H⌫N ¨ ds. Also, we note that, for all ␭ i g H 1r2 Ž ⌫i ., the harmonic extensions satisfy A1Ž R1Ž ␭1 ., ␸ 1 . s 0 for all ␸ 1 g H01 Ž ⍀ 1 ., A i Ž R i Ž ␭ iy1 , ␭ i ., ␸ i . s 0 for all ␸ i g H01 Ž ⍀ i ., i s 2, p y 1, and A p Ž R p Ž ␭ py1 ., ␸p . s 0 for all ␸p g H ⌫1py 1Ž ⍀ p ..

74

GATICA, HSIAO, AND MELLADO

Hereafter, w⭈, ⭈ x i denotes the usual duality pairing between Hy1 r2 Ž ⌫i . and H 1r2 Ž ⌫i . with respect to the L2 Ž ⌫i .-inner product, and w⭈, ⭈ x iy1, i stands for the corresponding duality pairing between Hy1 r2 Ž ⌫iy1 . = Hy1 r2 Ž ⌫i . and H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i ., i s 2, p y 1, that is,

Ž ␰ iy1 , ␰ i . , Ž ␮ iy1 , ␮ i .

iy1, i

[ w ␰ iy1 , ␮ iy1 x iy1 q w ␰ i , ␮ i x i

for all Ž ␰ iy1 , ␰ i . g Hy1 r2 Ž ⌫iy1 . = Hy1 r2 Ž ⌫i ., Ž ␮ iy1 , ␮ i . g H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i .. Then, we let S1 : H 1r2 Ž ⌫1 . ª Hy1 r2 Ž ⌫1 ., S i : H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i . ª Hy1r2 Ž ⌫iy1 . = Hy1r2 Ž ⌫i . ,

i s 2, p y 1

and S p : H 1r2 Ž ⌫py1 . ª Hy1 r2 Ž ⌫py1 . be the local Steklov᎐Poincare ´ operators defined, respectively, by S1 Ž ␭1 . , ␮ 1 S i Ž ␭ iy1 , ␭ i . , Ž ␮ iy1 , ␮ i . [

ž

⭸ ⭸␯ i

1

⭸

[

⭸␯ 1

R1 Ž ␭1 . , ␮ 1 , 1

iy1, i

R i Ž ␭ iy1 , ␭ i .

⌫ iy 1

,

⭸ ⭸␯ i

R i Ž ␭ iy1 , ␭ i . ⌫i

/

, Ž ␮ iy1 , ␮ i . iy1, i

and S p Ž ␭ py1 . , ␮ py1

py1

[

⭸ ⭸␯p

R p Ž ␭ py1 . , ␮ py1

. py1

Also, from now on we assume, without loss of generality, that p is odd. Then, we express all the local Steklov᎐Poincare ´ operators S i with odd index as S1 Ž ␭1 . , ␮ 1

1

s A1 Ž R1 Ž ␭1 . , R1 Ž ␮ 1 . .

Ž 2.

for all ␭1 , ␮ 1 g H 1r2 Ž ⌫1 ., S i Ž ␭ iy1 , ␭ i . , Ž ␮ iy1 , ␮ i .

iy1, i

s A i Ž R i Ž ␭ iy1 , ␭ i . , R i Ž ␮ iy1 , ␮ i . . Ž 3 .

for all Ž ␭ iy1 , ␭ i ., Ž ␮ iy1 , ␮ i . g H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i ., and S p Ž ␭ py1 . , ␮ py1

py1

for all ␭ py 1 , ␮ py1 g H 1r2 Ž ⌫py1 ..

s A p Ž R p Ž ␭ py1 . , R p Ž ␮ py1 . .

Ž 4.

75

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

On the other hand, for even i we write the local Steklov᎐Poincare ´ operators by using the boundary integral equation method. Then, following w16x, we use the hypersingular representation of S i given by Si [

ž

1 2

I q KX⭸ ⍀ i V⭸y1 ⍀i

/ ž

1 2

I q K ⭸ ⍀ i q W⭸ ⍀ i ,

/

Ž 5.

where the boundary integral operators are given by V⭸ ⍀ i ␴ Ž x . [

H⭸ ⍀ E Ž x, y . ␴ Ž y . ds , y

i

K ⭸ ⍀ i ␭Ž x . [

⭸

H⭸ ⍀ ⭸␯ i

KX⭸ ⍀ i ␴ Ž x . [

E Ž x, y . ␭Ž y . ds y , y

⭸

H⭸ ⍀ ⭸␯ i

E Ž x, y . ␴ Ž y . ds y , i

and W⭸ ⍀ i ␭Ž x . [ y

⭸ ⭸␯x

⭸

H⭸ ⍀ ⭸␯ i

E Ž x, y . ␭Ž y . ds y . y

Now, for even i we define the bilinear form A˜i : H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i . = H 1r2 Ž ⌫iy1 . = H 1r2 Ž ⌫i . ª R as A˜i Ž Ž ␭ iy1 , ␭ i . , Ž ␮ iy1 , ␮ i . . [ S i Ž ␭ iy1 , ␭ i . , Ž ␮ iy1 , ␮ i . s

½ž

1 2

I q KX⭸ ⍀ i V⭸y1 ⍀i

/ ž

1 2

iy1, i

5

q K ⭸ ⍀ i q W⭸ ⍀ i Ž ␭ iy1 , ␭ i . ,

/

Ž ␮ iy1 , ␮ i .

. Ž 6. iy1, i

py 1 Then, denoting ww⭈, ⭈ xx as the duality pairing between Ł is1 Hy1r2 Ž ⌫i . py 1 1r2 Ž . and Ł is1 H ⌫i , our Steklov᎐Poincare ´ problem reads: find ␭ [ py 1 Ž ␭1 , . . . , ␭ py1 . g Ł is1 H 1r2 Ž ⌫i . such that

SŽ ␭. s ␹ ,

Ž 7.

76

GATICA, HSIAO, AND MELLADO

py 1 py 1 where S : Ł is1 H 1r2 Ž ⌫i . ª Ł is1 Hy1r2 Ž ⌫i . is the global Steklov᎐ Poincare ´ operator defined by

SŽ ␭. , ␮

[ A1 Ž R1 Ž ␭1 . , R1 Ž ␮ 1 . . q A p Ž R p Ž ␭ py1 . , R p Ž ␮ py1 . . q

A i Ž R i Ž ␭ iy1 , ␭ i . , R i Ž ␮ iy1 , ␮ i . .

Ý 3FiFpy2 i odd

q

A˜i Ž Ž ␭ iy1 , ␭ i . , Ž ␮ iy1 , ␮ i . .

Ý 2FiFpy1 i even

py 1 for all ␭ [ Ž ␭1 , . . . , ␭ py1 ., ␮ [ Ž ␮ 1 , . . . , ␮ py1 . g Ł is1 H 1r2 Ž ⌫i ., and ␹ py 1 y1r2 Ž . g Ł is1 H ⌫i is given by

w ␹ , ␮ x [ yA1 Ž w˜1 , R1 Ž ␮ 1 . . q H fR1 Ž ␮ 1 . dx y A p Ž wp , R p Ž ␮ py1 . . ⍀1

q

H⍀ f R

p

Ž ␮ py1 . dx

p

q

Ý

V⭸y1 ⍀i

2FiFpy1 i even

q

Ý 3FiFpy2 i odd

H⍀ f Ž y . E Ž x, y . dy, Ž ␮

iy1 ,

␮i .

i

½H

⍀i

iy1, i

f R i Ž ␮ iy1 , ␮ i . dx y A i Ž wi , R i Ž ␮ iy1 , ␮ i . .

5

py 1 for all ␮ g Ł is1 H 1r2 Ž ⌫i .. We end this section by remarking that the interface problem Ž7., which is equivalent with Ž1., is uniquely solvable Žsee w20x..

3. ITERATIVE INTERFACE SOLVERS It is well known that a crucial part in any domain decomposition algorithm is how to efficiently solve the Steklov᎐Poincare ´ problem. The choice of an iterative method must take into account the inherent characteristics of the problem such as symmetry and positiveness. Since in general, at the discrete level, the spectral condition number of the Steklov᎐Poincare ´ operator increases dramatically as the step size of the mesh becomes smaller, an adequate choice of the preconditioner is of central importance.

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

77

3.1. The Richardson Scheme The first iterative algorithm that we use to solve the interface problem py1 Ž7. is the Richardson iteration method. Given ␭0 g Ł is1 H 1r2 Ž ⌫i . and a preconditioner operator P, the Richardson method for solving SŽ ␭. s ␹ constructs the sequence ␭ nq 1 s ␭ n y ␣ Py1w ␹ y SŽ ␭ n .x for all n g N j  04 , where ␣ ) 0 is a relaxation parameter. There are several possibilities for choosing P. In most applications, it reduces to a convex linear combination of S i , i s 1, p. However, in our present case, once we introduce a finite number of annular subregions partitioning ⍀, the above choice leads to an iteration-by-subdomain algorithm in which non-uniquely solvable problems appear. Hence, in order to overcome this difficulty, we use the preconditioner P [ Ýeven S i q T where py 1 py 1 T : Ł is1 H 1r2 Ž ⌫i . ª Ł is1 Hy1r2 Ž ⌫i . is defined by py1

T Ž ␭. , ␮

[

Ý

c i ² ␭ i , ␮ i :L2 Ž⌫i . ,

is1

where c i ) 0, i s 1, p, are arbitrary constants. Alternative choices for T can be seen in w10x. Here we are considering, without loss of generality, that p ) 3 is odd. Then, in order to obtain an iteration-by-subdomain algorithm for the interface problem Ž7., we substitute P into the Richardson sequence obtaining ␭ nq 1 s Ž1 y ␣ . ␭ n q ␣␤ n, where ␤ n satisfies

␤ns

y1

½ Ý S q T5 ½ ␹ y i

even

Ý S i y T Ž ␭n . odd

5

.

The latter equation suggests the following iterative algorithm to upgrade ␭n : Ž1. Solve in parallel the boundary value problems: y⌬u1n s f in ⍀ 1 ,

u1n s ␭1n on ⌫1 ;

u1n s g 0 on ⌫D ,

for i s 3, 5, . . . , p y 2, y⌬u in s f in ⍀ i ,

u in s ␭ niy1 on ⌫iy1 ,

u in s ␭ ni on ⌫i ;

and y⌬u np s f in ⍀ p ,

⭸ u np ⭸␯

u np s ␭ npy1 on ⌫py1 ,

s y2W Ž u np . q

b

␲r

on ⌫N .

78

GATICA, HSIAO, AND MELLADO

Set d1 s yŽ ⭸ u1nr⭸␯ 1 .< ⌫1 , Ž d iy1 , d i . s ŽyŽ ⭸ u inr⭸␯ i .< ⌫iy 1 , yŽ ⭸ u inr⭸␯ i .< ⌫i ., i s 3, 5, . . . , p y 2, and d py1 s yŽ ⭸ u np r⭸␯p .< ⌫py 1. n Ž2. Solve in parallel for Ž ␤iy1 , ␤in ., i s 2, 4, . . . , p y 1, the boundary integral equations,

ž

1 2

I q KX⭸ ⍀ i V⭸y1 ⍀i

/ ž

1 2

s Ž d iy1 , d i . q y

⭸ ⭸␯x

I q K ⭸ ⍀ i q W⭸ ⍀ i

/

1 2

I q KX⭸ ⍀ i

V⭸y1 ⍀i

n n , ␤in . q Ž c iy1 ␤iy1 , c i ␤ in . Ž ␤iy1

H⍀ f Ž y . E Ž x, y . dy i

H⍀ f Ž y . E Ž x, y . dy q Ž c

n iy1 ␭ iy1 ,

c i ␭ ni . ,

i

and compute ␭ nq 1 s Ž1 y ␣ . ␭ n q ␣␤ n. Ž3. If a given stop criterion is not satisfied, set n s n q 1 and go Ž to 1.. Finally, we observe that the convergence of this procedure is guaranteed by the fact that, for a sufficiently small parameter ␣ , the mapping py1

Ł

is1

py1

H 1r2 Ž ⌫i . 2 ␭ ª ␭ q ␣ Py1 ␹ y S Ž ␭ . g

Ł H 1r2 Ž ⌫i .

is1

becomes a contraction. We note that this can be easily proved by using Theorem 10 in w10x. Then, in this case we obtain that the solution u n [ Ž u1Ž ␭1n ., . . . , u p Ž ␭ npy1 .. converges to u in H 1 Ž ⍀ .. Note that, just for the sake of simplicity, we have utilized here the strong formulations of the boundary value problems. In order to deal with the discrete schemes, we must use the corresponding bilinear forms representing the variational formulations of those problems. 3.2. The Preconditioned Conjugate Gradient Method It is well known that fixed point formulations can be accelerated by using conjugate gradient techniques, which only applies to symmetric and py1 positive definite systems. We remark that, given ␭0 g Ł is1 H 1r2 Ž ⌫i ., the preconditioned conjugate gradient algorithm reads Žsee w12x. r 0 s ␹ y S Ž ␭0 . ,

␣0 s

² z 0 , P Ž z 0 .: , ² z 0 , S Ž z 0 .:

z 0 s Py1 Ž r 0 . ,

␭1 s ␭0 q ␣ 0 z 0 , w 1 s 1.

79

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

do n s 1 until convergence PŽ z n . s r n s ␹ y S Ž ␭n . ,

␣n s 1 wnq 1

² z n , P Ž z n .: , ² z n , S Ž z n .:

␣n

² z n , P Ž z n .: 1 s1y , ny 1 ny1 ␣ ny1 ² z , PŽ z .: wn

␭ nq 1 s ␭ ny 1 q wnq 1 Ž ␣ n z n q ␭ n y ␭ ny 1 . . Test on convergence end do. We note that the steps Ž1. and Ž2. of the BEM᎐FEM algorithm described in Section 3.1 will be used to describe the present procedure. Now, if we substitute P s Ýeven S i q T into the preconditioned conjugate py 1 gradient scheme and denote as before z i s z < ⌫i for all z g Ł is1 H 1r2 Ž ⌫i ., we obtain the following BEM᎐FEM iteration-by-subdomain algorithm: Ž1. Calculate ␤ 0 by using Ž1. and Ž2. and set z 0 s ␤ 0 y ␭0 . Ž2. Calculate Žtotally parallel step. ␣0 s

0 0 py1 Ýeven w Ž z iy1 , z i0 . , S i Ž z iy1 , z i0 . x iy1 , i q Ý is1 c i w z i0 , z i0 x i

py2 0 0 0 0 w z10 , S1 Ž z10 . x 1 q Ý is2 w Ž z iy1 , z i0 . , S i Ž z iy1 , z i0 . x iy1 , i q w z py1 , S p Ž z py1 . x py1

.

Ž3. Update ␭0 by means of ␭1 s ␭0 q ␣ 0 z 0 and set w 1 s 1. Ž4. Calculate ␤ n by using Ž1. and Ž2. and set z n s ␤ n y ␭ n. Ž5. Calculate Žtotally parallel step. py1

␥n s

n n , z in . , S i Ž z iy1 , z in . Ý Ž z iy1

even

1

␣n

s

1

␥n

iy1, i q

Ý

c i z in , z in i ,

is1 py2

½

z1n , S1 Ž z1n .

1

q

n n , z in . , S i Ž z iy1 , z in . Ý Ž z iy1

iy1, i

is2 n n q z py 1 , S p Ž z py1 .

py1

5

.

Ž6. Update ␭ n and wn through 1 ␣ n ␥n 1 s1y , wnq 1 ␣ ny1 ␥ny1 wn

␭ nq 1 s ␭ ny 1 q wnq 1 Ž ␣ n z n q ␭ n y ␭ ny 1 . . Ž7. If a given stop criterion is not satisfied, set n s n q 1 and go to Ž4..

80

GATICA, HSIAO, AND MELLADO

4. FINITE DIMENSIONAL APPROXIMATIONS Let Th be a quasi-uniform triangulation of ⍀ made up of triangles ␶ not crossing the interfaces ⌫i , i s 1, p y 1, and set ⍀ h [ D␶ : ␶ g Th 4 . Also, for i s 1, p, we put ⍀ i, h [ D␶ : ␶ g ⍀ i / ⭋4 and denote ⌫D, h and ⌫N, h , the polygonal curves Ždetermined by Th . approximating ⌫D and ⌫N , respectively. Now, we introduce the discrete analogue of A i ; that is, for odd i s 1, p y 2, A i, h : H 1 Ž ⍀ i, h . = H 1 Ž ⍀ i, h . ª R is defined by Ai , h Ž zi , ¨ i . [

H⍀

ⵜz i ⭈ ⵜ¨ i dx i, h

and A p, h : H 1 Ž ⍀ p, h . = H 1 Ž ⍀ p, h . ª R by Ap, hŽ zp , ¨ p . [

H⍀

ⵜz p ⭈ ⵜ¨ p dx q 2 ² ¨˙p , Vh Ž ˙ z p .: h p, h

for all z i , ¨ i g H 1 Ž ⍀ i, h ., where Vh : Hy1 r2 Ž ⌫N, h . ª H 1r2 Ž ⌫N, h . is the single layer operator on ⌫N, h and ² ⭈ , ⭈ : h stands for the duality pairing between Hy1 r2 Ž ⌫N, h . and H 1r2 Ž ⌫N, h . with respect to the L2 Ž ⌫N, h .-inner product. Also, we define the finite element subspaces Hi , h [  ¨ h g C Ž ⍀ i , h . : ¨ h < ␶ g P1 Ž ␶ . ᭙␶ : ⍀ i , h 4 and ⌳ i , h [  ¨ h < ⌫ i : ¨ h g Hi , h 4 , where P1Ž␶ . is the space of polynomials of degree F 1 defined on ␶ . In addition, we set H1,0 h [  ¨ h g H1, h : ¨ h s 0 on ⌫D , h j ⌫1 4 , Hi0, h [  ¨ h g Hi , h : ¨ h s 0 on ⌫iy1 j ⌫i 4 ,

i s 2, p y 1

and H p0, h [  ¨ h g H p , h : ¨ h s 0 on ⌫py1 4 . Then, in analogy with the continuous operators, we define for odd i the discrete harmonic extensions: 䢇

R1, h : ⌳ 1, h 2 ␭1, h ª R1, hŽ ␭1, h . g H1, h solution of A1, h Ž R1, h Ž ␭1, h . , ¨ h . s 0

᭙ ¨ h g H1,0 h ,

R1, h Ž ␭1, h . s ␭1, h

on ⌫1 ,

R1, h Ž ␭1, h . s 0

on ⌫D , h .

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

81

For odd i s 3, p y 2, R i, h : ⌳ iy1, h = ⌳ i, h 2 Ž ␭ iy1, h , ␭ i, h . ª R i, hŽ ␭ iy1, h , ␭ i, h . g Hi, h solution of 䢇

A i , h Ž R i , h Ž ␭ iy1, h , ␭ i , h . , ¨ h . s 0

᭙ ¨ h g Hi0, h ,

R i , h Ž ␭ iy1, h , ␭ i , h . s ␭ iy1, h R i , h Ž ␭ iy1, h , ␭ i , h . s ␭ i , h 䢇

on ⌫iy1 , on ⌫i .

R p, h : ⌳ py1, h 2 ␭ py1, h ª R p, hŽ ␭ py1, h . g Hp, h solution of A p Ž R p , h Ž ␭ py1, h . , ¨ h . s 0

᭙ ¨ h g H p0, h ,

R p , h Ž ␭ py1, h . s ␭ py1 , h

on ⌫py1 .

In what follows we denote ⌳ h [ ⌳ i, h and by ⌳Uh its dual. Moreover, w⭈, ⭈ x i, h Žresp. w⭈, ⭈ x iy1, i, h . denotes the duality pairing between ⌳ i, h Žresp. ⌳ iy1, h = ⌳ i, h . and its dual space ⌳Ui, h Žresp. ⌳Uiy1, h = ⌳Ui, h .. With the above definitions, we can introduce the discrete local Steklov᎐Poincare ´ operators expressed by FEM Žodd index. as py 1 Ł is1

S1, h Ž ␭1, h . , ␮ 1, h

1, h

[ A1, h Ž R1, h Ž ␭1, h . , R1, h Ž ␮ 1, h . .

for all ␭1, h , ␮ 1, h g ⌳ 1, h , S i , h Ž ␭ iy1, h , ␭ i , h . , Ž ␮ iy1, h , ␮ i , h .

iy1, i , h

[ A i , h Ž R i , h Ž ␭ iy1, h , ␭ i , h . , R i , h Ž ␮ iy1, h , ␮ i , h . . for all Ž ␭ iy1, h , ␭ i, h ., Ž ␮ iy1, h , ␮ i, h . g ⌳ iy1, h = ⌳ i, h and S p , h Ž ␭ py1, h . , ␮ py1 , h

py1, h

[ A p , h Ž R p , h Ž ␭ py1, h . , R p , h Ž ␮ py1, h . .

for all ␭ py 1, h , ␮ py1, h g ⌳ py1, h . In addition, for even i the discrete Steklov᎐Poincare ´ operators S i, h are expressed by BEM as in Ž5., that is, S i, hŽ ␭ iy1, h , ␭ i, h . [ S i Ž ␭ iy1, h , ␭ i, h . for all Ž ␭ iy1, h , ␭ i, h . g ⌳ iy1, h = ⌳ i, h , which implies that they preserve the continuity and coerciveness properties. At this point it is very important to remark that, since the interfaces ⌫i , i s 1, p y 1 are polygonals, the bilinear forms A i and A i, h coincide for all odd i s 2, p y 1. Also, the same is true for A˜i and the corresponding discrete bilinear form A˜i, h . However, it is also worth remarking that for practical computations one must replace V⭸y1 ⍀ i in the definition of S i, h by an approximated operator of finite range. More precisely, given finite dimenr2 sional subspaces Hy1 of Hy1r2 Ž ⌫i . with mesh size ˜ h, we define the ˜ i, h operator 1 y1r2 V⭸y1 I q K ⭸ ⍀ i : ⌳ iy1, h = ⌳ i , h ª Hy1r2 ˜ ˜ = Hi , h˜ ⍀ i, h iy1, h 2

ž

/

Ž ␭ iy1, h , ␭ i , h . ª Ž ␰ iy1, h˜ , ␰ i , h˜ . ,

82

GATICA, HSIAO, AND MELLADO

r2 y1r2 where Ž ␰ iy1, h˜ , ␰ i, h˜ . g Hy1 is the unique solution of ˜ = Hi, h˜ iy1, h

Ž ␩iy1, h˜ , ␩i , h˜ . , V⭸ ⍀ iŽ ␰ iy1, h˜ , ␰ i , h˜ . s Ž ␩iy1, h˜ , ␩i , h˜ . ,

ž

1 2

iy1, i

I q K ⭸ ⍀ i Ž ␭ iy1, h , ␭ i , h .

/

iy1, i

r2 y1r2 for all Ž␩iy1, h˜ , ␩i, h˜ . g Hy1 ˜ = Hi, h˜ . iy1, h For further details on this technical analysis we refer to w20x. Now, the discrete Steklov᎐Poincare ´ operator S h : ⌳ h ª ⌳Uh becomes p S h [ Ý is1 S i, h . According to our previous choice, we take constants c i ) 0, i s 1, p y 1, and consider Ph : ⌳ h ª ⌳Uh as follows: Ph [ Ýeven S i, h q Th , where for all ␭ h s Ž ␭1, h , . . . , ␭ py 1, h ., ␮ h s Ž ␮ 1, h , . . . , ␮ py 1, h . g ⌳ h

py1

Th Ž ␭ h . , ␮ h

h

[

Ý ci w ␭i , h , ␮i , h x i , h .

is1

Here, ww⭈, ⭈ xx h stands for the duality pairing between ⌳ h and ⌳Uh . We note that the invertibility of Ph and S h guarantees the convergence of the Richardson algorithm for a sufficiently small parameter ␣ . However, from the numerical point of view, one needs to prove that this convergence does not depend on h. In fact, we first observe that Phy1 S h is symmetric and positive definite with respect to the inner product defined by w ␭ h , ␮ h x Ph [ wwPhŽ ␭ h ., ␮ h xx h for all ␭ h , ␮ h g ⌳ h . It follows that the eigenvalues of Phy1 S h are positive real numbers 0 - ␴ 1 F ␴ 2 F ⭈⭈⭈ F ␴m , where m is the dimension of ⌳ h . Also, we note that ␬ , the spectral condition number of Phy1 S h , is given by ␬ [ ␴m r␴ 1 where ␴m s sup␭h g ⌳ hŽwPhy1 S hŽ ␭ h ., ␭ h x Phrw ␭ h , ␭ h x Ph . and ␴y1 s sup␭ h g ⌳ hŽw ␭ h , ␭ h x Phr 1 wPhy1 S hŽ ␭ h ., ␭ h x P .. h The following two lemmata will be needed to estimate ␬ . LEMMA 4.1. Let ⌫ be a Lipschitz continuous boundary and consider ˜ g H 1r2 Ž ⌫ .rR, and ˜r g R such that ␭ s ␭˜ q ˜r. Then, there ␭ g H 1r2 Ž ⌫ ., ␭ exists C ) 0 such that

˜5 2H 1r 2 Ž⌫ . q 5 ␭ 5 2L2 Ž⌫ . G C 5 ␭ 5 2H 1r2 Ž⌫ . 5␭

᭙␭ g H 1r2 Ž ⌫ . .

Proof. First, by using the triangle inequality, it follows that

˜ q ˜r 5 2H 1r2 Ž⌫ . 5 ␭ 5 2H 1r 2 Ž⌫ . s 5 ␭ ˜5 2H 1r 2 Ž⌫ . q 5 ˜r 5 2H 1r2 Ž⌫ . F 2 5␭

½

5

83

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

and then, by using the fact that 5 ˜ r 5 H 1r 2 Ž⌫ . s 5 ˜ r 5 L2 Ž⌫ . , we can write

˜5 2H 1r2 Ž⌫ . q 5 ˜r 5 2L2 Ž⌫ . 5 ␭ 5 2H 1r 2 Ž⌫ . F 2 5 ␭

½ ˜5 s 2½5 ␭ ˜5 F 2½5 ␭

5

2 H 1r 2 Ž⌫ .

˜5 2L2 Ž⌫ . q 5␭ y ␭

2 H 1r 2 Ž⌫ .

˜5 2L2 Ž⌫ . . q 2 5 ␭ 5 2L2 Ž⌫ . q 2 5 ␭

5 5

Finally, since 5 ␰ 5 L2 Ž⌫ . F 5 ␰ 5 H 1r 2 Ž⌫ . for all ␰ g H 1r2 Ž ⌫ ., we conclude that

˜5 2H 1r2 Ž⌫ . q 5 ␭ 5 2L2 Ž⌫ . 4 5 ␭ 5 2H 1r 2 Ž⌫ . F 6  5 ␭

᭙␭ g H 1r2 Ž ⌫ . ,

which completes the proof. LEMMA 4.2. There exist two positi¨ e constants, ␣ S , MS , independent of h, such that py1

␣S

Ý

5 ␭ i , h 5 2H 1r 2 Ž⌫i . F w S h ␭ h , ␭ h x

is1

py1 h F MS

Ý

5 ␭ i , h 5 2H 1r2 Ž⌫i .

is1

for all ␭ h s Ž ␭1, h , . . . , ␭ py1, h . g ⌳ h . Proof. It follows from the well known extension theorem for FEM Žsee w2, 28x., the properties of the discrete single layer potential Vh , the generalized Poincare ´ inequality, and the trace theorem. We omit further details and refer the reader to w20, 25, 26x. The following result establishes that the convergence of our algorithms is independent of h. THEOREM 4.1. There exists a positi¨ e constant C, independent of h, such that ␬ F C. Proof. By using the definitions given through this section, we can write the discrete preconditioner Ph as Ph Ž ␭ h . , ␭ h

h

s

S i , h Ž ␭ iy1, h , ␭ i , h . , Ž ␭ iy1, h , ␭ i , h .

Ý 2FiFpy1 i even py1

q

Ý

c i 5 ␭ i , h 5 2L2 Ž⌫i .

᭙␭ h g ⌳ h

is1

and then, using the fact that 5 ␭ i , h 5 L2 Ž⌫i . F 5 ␭ i , h 5 H 1r 2 Ž⌫i .

᭙␭ i , h g H 1r2 Ž ⌫i . ,

iy1, i , h

84

GATICA, HSIAO, AND MELLADO

we find that there exists MP ) 0, independent of h, such that py1

Ph Ž ␭ h . , ␭ h

F MP h

5 ␭ i 5 2H 1r 2 Ž⌫i .

Ý

᭙␭ h g ⌳ h .

is1

On the other hand, in order to obtain a lower bound for wwPhŽ ␭ h ., ␭ h xx h , it follows that for all ␭ h g ⌳ h Ph Ž ␭ h . , ␭ h

h

G

˜iy1, h , ␭˜i , h ␣i ␭

Ý 2FiFpy1 1 even

ž

2

/

H 1r 2 Ž ⌫ iy 1 .=H 1r2 Ž ⌫ i .

py1

q

Ý

c i 5 ␭ i , h 5 2L2 Ž⌫i .

is1 py1

G

Ý ␥ i  5 ␭˜i , h 5 2H

1r 2

Ž⌫ i .

q 5 ␭ i , h 5 2L2 Ž⌫i . 4 ,

is1

˜i, h g H 1r2 Ž ⌫i .rR is defined as in Lemma 4.1 and ␥ i [ min ␣ i , c i 4. where ␭ By setting ␥ [ min␥ i 4 and using Lemma 4.1, we conclude that py1

Ph Ž ␭ h . , ␭ h

h

G␥

Ý  5 ␭˜i , h 5 2H

1r 2

Ž⌫ i .

q 5 ␭ i , h 5 2L2 Ž⌫i . 4

is1 py1

G ␣P

Ý

5 ␭ i , h 5 2H 1r 2 Ž⌫i .

is1

for all ␭ h g ⌳ h . Finally, taking into account the above estimates concerning the discrete preconditioner and Lemma 4.2 for the discrete global Steklov᎐Poincare are given by ´ operator, the bounds for ␴m and ␴y1 1 1 Žww ␴m s sup␭ h g ⌳ h S hŽ ␭ h ., ␭ h xx h rwwPhŽ ␭ h ., ␭ h xx h . F MS r␣ P and ␴y s 1 Žww Ž . xx ww Ž . xx . sup␭h g ⌳ h Ph ␭ h , ␭ h hr S h ␭ h , ␭ h h F MPr␣ S . This completes the proof since ␬ s ␴m ␴y1 1 . Now, concerning numerical experiments, we mention that the exterior boundary value problem studied here has already been successfully solved by using FEM-based domain decomposition methods in w6x. Also, several numerical experiments performed on a parallel computer, using the preconditioners of the present paper together with Richardson’s and conjugate gradient methods, are reported in w10x. Hence, the suitable combination of FEM and BEM explained here must produce at least as good results as in w6, 10x, with less computational effort. As a final remark, we note that the extension of these algorithms to nonlinear exterior boundary value problems can be done similarly as shown here by using linear preconditioners on H 1 Ž ⍀ . Žsee w10, 20x..

DOMAIN DECOMPOSITION BASED ON BEM AND FEM

85

ACKNOWLEDGMENTS This research was started while M. E. Mellado visited the Department of Mathematical Sciences of the University of Delaware during the period March 21᎐May 6, 1998. This author expresses his gratitude to Professor G. C. Hsiao and his colleagues for their hospitality and stimulating scientific discussions.

REFERENCES 1. V. I. Agoshkov, Poincare᎐Steklov operators and domain decomposition methods in finite ´ dimensional spaces, in ‘‘Domain Decomposition Methods for Partial Differential Equations’’ ŽR. Glowinski et al., Eds.., pp. 73᎐112, SIAM, Philadelphia. 2. J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring, I, Math. Comp. 47 Ž1986., 103᎐134. 3. C. Carstensen, M. Kuhn, and U. Langer, Fast parallel solvers for symmetric boundary element domain decomposition equations, Numer. Math. 79 Ž1998., 321᎐347. 4. M. Costabel, Symmetric methods for the coupling of finite elements and boundary elements, in ‘‘Boundary Elements IX’’ ŽC. A. Brebbia et al., Eds.., Vol. 1, pp. 411᎐420, Springer-Verlag, Berlin, 1987. 5. G. N. Gatica, Combination of mixed finite element and Dirichlet-to-Neumann methods in nonlinear plane elasticity, Appl. Math. Lett. 10 Ž1997., 29᎐35. 6. G. N. Gatica, E. C. Hernandez, and M. E. Mellado, A domain decomposition method for ´ linear exterior boundary value problems, Appl. Math. Lett. 11 Ž1998., 1᎐9. 7. G. N. Gatica and G. C. Hsiao, On the coupled BEM and FEM for a nonlinear exterior Dirichlet problem in R 2 , Numer. Math. 61 Ž1992., 171᎐214. 8. G. N. Gatica and G. C. Hsiao, The uncoupling of boundary integral and finite element methods for nonlinear boundary value problems, J. Math. Anal. Appl. 189 Ž1995., 442᎐461. 9. G. N. Gatica and G. C. Hsiao, ‘‘Boundary-Field Equation Methods for a Class of Nonlinear Problems,’’ Pitman Research Notes in Mathematics Series, Vol. 331, Longman, HarlowrNew York, 1995. 10. G. N. Gatica and M. E. Mellado, Nonoverlapping domain decomposition methods for linear and nonlinear exterior boundary value problems, in ‘‘Proceedings of the Fourth World Congress on Computational Mechanics, Buenos Aires, Argentina, June 29᎐July 2, 1998’’ ŽS. Idelsohn and E. Onate, Eds.., CIMNE, Barcelona, 1998. 11. G. N. Gatica and W. L. Wendland, Coupling of mixed finite elements and boundary elements for a hyperelastic interface problem, SIAM J. Numer. Anal. 34 Ž1997., 2335᎐2356. 12. G. H. Golub and C. F. Van Loan, ‘‘Matrix Computations,’’ Johns Hopkins Press, Baltimore, 1983. 13. H. Han, A new class of variational formulations for the coupling of finite and boundary element methods, J. Comput. Math. 8 Ž1990., 223᎐232. 14. H. Han and X. Wu, The approximation of the exact boundary conditions at an artificial boundary for linear elastic equations and its applications, Math. Comp. 59 Ž1992., 21᎐37. 15. G. C. Hsiao, B. N. Khoromskij, and W. J. Wendland, Boundary integral operators and domain decomposition, preprint 94-11, Mathematisches Institut A, Universitat ¨ Stuttgart, 1994.

86

GATICA, HSIAO, AND MELLADO

16. G. C. Hsiao and W. L. Wendland, ‘‘Domain Decomposition via Boundary Element Methods,’’ Technical Report 92-10, Department of Mathematical Sciences, University of Delaware, 1992. 17. G. C. Hsiao and S. Zhang, Optimal order multigrid methods for solving exterior boundary value problems, SIAM J. Numer. Anal. 31 Ž1994., 680᎐694. 18. C. Johnson and J. C. Nedelec, On the coupling of boundary integral and finite element methods, Math. Comput. 35 Ž1980., 1063᎐1079. 19. P. Le Tallec, Domain decomposition methods in computational mechanics, Comput. Mech. Ad¨ . 1 Ž1994., 121᎐220. 20. M. E. Mellado, ‘‘Numerical Solution of Exterior Problems in Potential Theory and Elasticity,’’ Ph.D. Thesis, Universidad de Concepcion, ´ 2000. 21. A. Quarteroni and A. Valli, Theory and application of Steklov᎐Poincare ´ operators for boundary value problems, in ‘‘Applied and Industrial Mathematics’’ ŽR. Spigler, Ed.., pp. 179᎐203, Kluwer Academic, Dordrecht, 1991. 22. A. Quarteroni and A. Valli, ‘‘Domain Decomposition Methods for Partial Differential Equations,’’ Oxford Univ. Press, London, 1999. 23. G. Schmidt, Boundary element discretization of Poincare᎐Steklov operators, Numer. ´ Math. 69 Ž1994., 83᎐101. 24. B. Smith, P. Bjorstad, and W. Gropp, ‘‘Domain Decomposition,’’ Cambridge Univ. Press, Cambridge, UK, 1996. 25. O. Steinbach, Boundary elements in domain decomposition methods, in ‘‘Domain Decomposition Methods in Scientific and Engineering Computing, University Park, PA, 1993,’’ Contemporary Mathematics, Vol. 180, pp. 343᎐348, Amer. Math. Soc., Providence, 1994. 26. O. Steinbach and W. L. Wendland, Domain decomposition and preconditioning techniques in boundary element methods, in ‘‘Boundary Element Topics, Stuttgart, 1995,’’ pp. 471᎐490, Springer-Verlag, Berlin, 1997. 27. W. L. Wendland, On asymptotic error estimates for the combined BEM and FEM, in ‘‘Finite Element and Boundary Element Techniques from Mathematical and Engineering Point of View’’ ŽE. Stein and W. L. Wendland, Eds.., CISM Lecture Notes, Vol. 301, pp. 273᎐333, Udine, Springer-Verlag, New York, 1988. 28. O. B. Widlund, An extension theorem for finite element spaces with three applications, in ‘‘Numerical Techniques in Continuum Mechaniques’’ ŽW. Hackbush and K. Witsch, Eds.., BraunschweigrWiesbaden, 1987.