Terminada Catherine (Catalina) 2014 | Michael Alexander - Excel-VBA Power-Programmierung für Dummies | twin peaks 265

The h-p boundary element method for solving 2- and 3-dimensional problems

The h-p boundary element method for solving 2- and 3-dimensional problems

Computer methods in a p p l l o d mochamics and engineering EI~SEVIER Comput. Methods Appl. Mech. Engrg. 133 (1996) 183-208 The h-p boundary elemen...

1MB Sizes 0 Downloads 4 Views

Computer methods in a p p l l o d

mochamics and engineering EI~SEVIER

Comput. Methods Appl. Mech. Engrg. 133 (1996) 183-208

The h-p boundary element method for solving 2- and 3-dimensional problems E.P. Stephan lnstina f iir Angewandte Mathematik. Universitiit Hannover, ttatmover, Germany Received 2I July 1995

Abstract

In this survey paper we report on recent developments of the hp-version of the boundary element method. We consider weakly singular and hypersingular integral equations on polygons and polyhedral surfaces. We show that the Galerkin solutions (computed with the hp-version on geometric meshes) converge exponentially fast towards tile exact solutions of the integral equations, "lhe implementation of the hp-version of the boundary element method is discussed and adaptive algorithms are given. Model problems for the coupling of finite elements and boundary elements are also considered. Numerical results are presented which underline the theoretical results,

1. I n t r o d u c t i o n

The h-p version was first analyzed theoretically for the finite element method ( F E M ) by Babugka and others (for an overview see [2]). Meanwhile, the technique of the h-p version was applied to the b o u n d a r y element method ( B E M ) for solving Dirichlet and N e u m a n n problems for the Laplacian in plane polygonal domains [ 3 , 4 , 9 ] , for corresponding mixed boundary value problems [13] and transmission problems [1I]. Recently, three-dimensional boundary value problems have been treated with the h-p version of the B E M [16, 19]. All boundary value/transmission problems have in c o m m o n that they can be reduced to strongly elliptic boundary integral equations (pseudodifferential equations) on the b o u n d a r y / i n t e r f a c e . The basic idea of the convergence proofs is the observation that for strongly elliptic pseudodifferential equations one obtains quasioptimai convergence in the energy norm for any Galerkin scheme with conforming elements [32]. This result is used to analyze tht~ h-p version which turns out to be the most suitable numerical approach to achieve accuracy and reliable solutions: the h-p version with a geometric mesh refinement yields exponentially fast convergence, This is not restricted to two-dimensional problems. As recently shown, the h-p version of the B E M can also be applied to three-dimensional boundary value problems [16, 19] and transmission problems [12], it is a very robust m e t h o d since a geometric mesh refinement towards edges and corners of the domain together with an appropriate p-distribution suffices to obtain exponentially fast converging numerical solutions, Sometimes, one uses a combination of finite elements and boundary elements, e.g. in transmission problems, where the interior problem is described by a (possibly non-linear) differential o p e r a t o r with variable coefficients, and the exterior problem is given via a linear differential o p e r a t o r with constant coefficients. In this case the symmetric coupling method of F E M and B E M can also be performed with an h-p version, and again exponentially fast convergence can be obtained [14, 15]. 0045-7825/96/$15.00 i~) 1996 Elsevier Science S.A. All rights reserved SSD! 0045-7825(95)00940-X

E.I', Stephan / Cmnpttt. Methods" Appl. Mech. Engrg. 133 (I096) I83-208

184

2. Boundary integral equations We consider integral equation m e t h o d s for solving Dirichlet and N e u m a n n b o u n d a r y value p r o b l e m s for the Laptacian in polygonial domains .(5 c-_ ~-" or polyhedral d o m a i n s .(? c-_ ~ .

Dirichlet: For given f E H ' 2 ( F ) on the b o u n d a r y F of II find u E H~(~I) such that Au=0

in.O,

u=f

on F .

(1)

Neumann" For given g E H - I"-'(F) on F = a~l with J'r g ds = 0 find u E H~(I'~) such that At,=(I

~ll

On- g

in.f~.

onF.

(2)

H e r e . . Q is a polygon or p o l y h e d r o n with piecewise C ~ b o u n d a r y F, cusps are excluded. The n o r m a l vector always points into the exterior of J2. H e r e , the Sobolev spaces It'(F) for s > 0 are defined to be the restriction of H " t " - ( R '') to F, n = 2 c 3"

{ul, "u~H"' "(~:)}

H'(I')'=

for s = 0 to be L : ( F ) and for s < 0 by duality (see [17]). It is well k n o w n [5] that p r o b l e m s (1) and (2) can be reduced to integral e q u a t i o n s of the first kind on the b o u n d a r y F for the C a u c h y data v = ulr and ~,O= iJu/On]r. N a m e l y , we have, respectively

V0=(l+K)f

onl" Wv=(l-K')g o n F

(3) (4)

with the integral o p e r a t o r s Vw(x) "= 2

fi G(x, y ) w ( y ) ds, .

K'w(x) "= 9_ . ~G(X,gm,

K;v(x) "= 2 f~ an,a---mG(x, y)w(y) ds,.

.v)w( vI ds,

Wn,(x) : = 2 0--~i,

.-~., G(x, y)w(y) ds~ "

with the f u n d a m e n t a l solution of the Laplacian

- ~ G(x, y ) ' =

1

1

In Ix - , v l .

n

= 2.

(5)

1

4--~- Ix - y l

"

n = 3.

W h e r e a s for n = 3. -Q can be completely arbitrary, let us assume that for the case of a twodimensional d o m a i n .O its b o u n d a r y F has conformal radius less than one; this can always be achieved -|.~ by an a p p r o p r i a t e scaling of the d o m a i n . T h e n there exist unique solutions qJ ~ H - ( F ) of (3) and o ~ H t '-'(F) of (4), where the solution of (4) is only unique up to constants. Eq. (3) is equivalent to p r o b l e m ( I ) : i.e. let u ~,tli(.(2) solve (1). then ~b = au/Onl,-solves (3). conversely, let t ! , ~ . H l ' " ( F ) solve (3) thcn u defined by u(x)=

.

G(.r, y)4,(y)--b-~,.G(x, y)v(y)

ds.,., x E 1 2

(6)

with v = f solves ( i ) . O n the o t h e r hand, (4) is equivalent to (2) in the sense that let u @ H i ( g 2 ) s o l v e (2). then t, = U[r solves (4), conversely, let v E H t ' 2 ( F ) solve (4), then u defined by (6) with d , = g

E.I" Stephan

/

Comput. Method~ AppI. Meclt. Engrg. 133 (1996) i83-208

185

solves (2). Of course, boundary integral equations can also be used to solve boundary value problems in unbounded regions. For brevity, we consider here only screen problems in l~ ~, i.e. boundary value problems exterior to an open surface I" (for corresponding two-dimensional problems see [25]). !

3

7.

Dirichlet screen" For given f ~ H ~ -(1") find u E HI,,,(II~ '4 ) such that Au =

()

in ~k/-~ .

u = f

onl'.

N e u m a n n screen: For given g C H ,xu=0

in R ~

""

u=O

(T I-V)

aslxl--,~

(7)

' -'(I') find u E HI,,,.(I]~~\/~) such that o n 1".

an au - g

u = o (I--~t)

~s Ix I---" ~

(8)

As shown in [29]. these screen problems can be converted into the integral equations V~,=Zf

onF

PcL' = 2g

on I"

(9) (10)

with the jumps O = l a u / a n l l , ~ f l ' -'(1") and t " = [ U I ] r E H ' 2(F) across F. Again, there holds existence and uniqueness of solutions O E / t i -(i.) and v ElYl ~ "-(F) of (9) and (10). Here, the Sobolev space for an open surface piece F is introduced with the help of a closed surface /~, where F C F, namely H ' ( F ) = {u E H'(/~): supp u C/-'}. Note that n " = ( F ) = H ~ ' " ( F ) (for H'~"-(F) see [17]) and 1 4 - ' " - ( F ) = ( H " Z ( F ) ) ' by duality. In electrostatics qJ describes the ~.[m[g¢ density of the .~uiia,_~: ,~ith given putcnti~! f Cleally, the application of boundary integral equations is not restricted to the above benchmark problems, e.g. they can be applied to the Lamd system of linear clasticity in a similar manner when the corresponding fundamental solution is taken instead of (5). The boundary integral equations (3) and (4) (or (9) and (10)) can be solved approximately by the Galerkin method using co',forming subspaces {XN} of H -I'_(i.) ( o r / 4 - ' / z ( F ) ) and {YN} of H I / 2 ( F ) (or / i ~ ' z ( F ) ) . respectively. In the b o u n d a r y element m e t h o d one uses finite element spaces as such subspaces. The majority of results and boundary element software are based on the classical h-version. Recently, the p-version [311 and the h-p version were developed [30] and here we report on the progress of the latter version. Since the operators of the single layer potential V and the normal derivative of the double layer potential W satisfy a G&rding's inequality in H - J ' z ( F ) and_ H ~ " - ( F ) , respectively, in case of a closed manifold F (and in H - J ' : ( F ) a n d / 4 J " ( F ) . respectively, if F is open), any conforming Galerkin scheme converges quasioptimally [32], i.e. we have with a constant C independent of N IIq, - q,,.ll x ~< C l l ~ , - ~ , t l x

(1I)

I1~ - v,~ll,, ~ cflv - wiiy

(12)

and

for any ~/, ~ XN C H - 1,2([.) =. X and w E Yx C H ~ ' : ( I ") =" Y. Here, 0 solves (3) and v solves (4) and qJN, V,,, are the solutions of the corresponding Galerkin equations with (-, .) denoting the L Z-inner product on F:

(v,/,.,,. 4,) = ( L ,/,)

v,/, ~ x~

(~3)

and

(WvN, W ) - - ( g , w

)

VwEY

N.

(14)

186

E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 133 (1996) 183-208

For the h-version optimal convergence results for the above Galerkin schemes together with mesh grading towards the edges of the p o l y h e d r o n / " have been given in 124]. For the p-version applied to the weakly singular integral equation (3) and the hypersingular equation (4), see [28] and [26], respectively. The key of the error analysis for both h- and p-version are sharp regularity results for the solutions of (3) and (4); these results are derived in [23, 24] and are based on the work in [71. In order to obtain these results we analyse the original boundary value problems (1) and (2). Since the solutions of (3) and (4) are the corresponding unknown Cauchy data, we get their regularity results by taking traces of the solutions of the boundary value problems. Let us shortly review this approach by considering the three-dimensional Dirichlet problem which was analysed in [23]; for the Neumann problem see [24 !. Let us consider the case of a polyhedron 1'1 with F = 0.0 = U,= /~ in ~~ with plane faces F ' . We describe Dirichlet data f E H~(F) where H " ( F ) is defined above. Then. the Neumann data 4, of the solution of (1) has regularity H ' - ' away from the edges and corners. Near an edge with opening w there are singularities of the form c( y)p "" ~"-'p - t Here v = +r/to. m > 0 and p ~>0 are integers and p denotes the distance to the edge while the stress intensity factor c(y) is a function defined on the edge. Near the corners we get additional corner singularities of the form ra'-'wk(~" ) ,

~E~,.

where r denotes the distance to the vertex and w k is a sphere centered in the vertex. ~" = (0. ,#) are function w~ are obtained as follows: Consider operator A,.¢ on F,,. let/~k be the kth eigenvalue

is a function on the spherical polygon F,, - F rl S_~. S_, polar coordinates on S_,. The exponent Ak and the the eigenvalue problem for the L a p l a c e - B e l t r a m i with corresponding eigenfunction v~:

n~.,Vk(O, ~) = ttkvk(O, 'P) V~(O, 'P)I,~, = O. Then Ak =

-- ~" +

k +

~

and

w,~ ' =

~

%(0, ~)

,;,

With this notation there holds the following decomposition result for the solution of (3) into edge and corner singularities where the edge terms are of tensor product form. T H E O R E M I. Let f be smooth on F j, j = 1. . . . . (3) has the decomposition on 1"i On , . , = ~ ' + ,

¢or

~'~ x , r ^ ' - ' w ~ ( O ) + rlvr*;

~

CdgC~

J, F = U i =~

X~c/(y)p~'-'

pJ

Then, the solution d~ = Ou / On on F o f

(15)

where the sums o f edge and corner singularity terms are finite. Here, Xc(Xe) are suitable C ~ cut-offfunctions near corners (edges) and ~h~ consists o f a global remainder term belonging to s o m e H~(F j) plus higher-order.edge and corner singularity terms, c i is a function on the ]th edge belonging to specific weighted Sobolev spaces and t~ := ~ / w / where wp is the opening angle at the jth edge. ( F o r a detailed formulation o f the theorem see [23]). R E M A R K . Since p~, E H " : + ' , - ' ( F ) sufficiently smooth f is ¢,lr, ~ H-"-~+"(F i )

and r ~ ~ H t + ~ ' - ' ( F ) the regularity of the solution ~b of (3) for (i.6)

E.P. Stephan f Comput. Methods Appl. Me,'h. Engrg. 133 (1996) 183-2(18

187

with

(17)

a=min{~.a~+l/2}.

In case of screen p r o b l e m s , the p o l y h e d r o n d c g c n c r a t e s to a surface piece F a n d ,(,~ b e c o m e s the region e x t e r i o r to the screen F. T h e regularity t h e o r y in [23] shows that e v e n for a s m o o t h r i g h t - h a n d .~ide ./~ the solution ~k of the integral e q u a t i o n (9) b e c o m e s u n b o u n d e d n e a r t h e e d g e s of r , it e x h i b i t s so-called edge and c o r n e r singularities. F r o m [23] we cite the regularity result for the s q u a r e plate F = (0, 1 ~, × (0, 1 ). T h e r e f o r e , a c o r r e s p o n d i n g regularity analysis leads to t h e s i t u a t i o n in E x a m p l e 1.

T H E O R E M 2. Let f ~_H~-(F). Then. the unique solution q, E t71 t : ( / , ) o f (9) has near the origin a d e c o m p o s i t i o n into edge a n d corner singularities o f the f o r m q' = ~, + x~(r)a lr~ " 'w(O ) + X~.(O )e l(x)y "'- + X , . ( ' r r / 2 - O ) e , ( y ) x ' ~''a-tj2 +c,, with ~ , , G H I - ' ( I ' ) , atE~, w E H I , [0.~T/21, e ~ ( x ) = b l x A I : + c ~. e . _ ( y ) = b z y Hql~(~'/, b, E ~ , i = I. 2. A corresponding d e c o m p o s i t i o n holds at the other vertices o f 11.

c,E

H e r e , (r, 0 ) are p l a n e p o l a r c o o r d i n a t e s c o n c e n t r a t e d at the origin • X~ a n d ,~~ a r e C ~- c u t - o f f f u n c t i o n s with X~ - 1 for r < 1 / 4 . x~ - I for 0 < ~ / 4 . We have A ~ [).2966 (see [21]) a n d ~ = ~ since % = 2"tr.

3.

Standard

boundary

element

methods

3. I. The h - m e t h o d with graded mesh F o r simplicity we o n l y c o n s i d e r the case w h e r e all the faces m e e t i n g at t h e v e r t e x a r e c o n v e x . T h e n e a c h of t h e s e s e c t o r s can be m a p p e d linearly on the q u a d r a n t ~ , × ~ +t c t ~ - ' . O n t h e s q u a r e [0, 1] 2 we i n t r o d u c e a g r a d e d m e s h given by the lines x t = (ih) t', x , = ( ] h ) ' , h = l / N , i, ] = 0 . . . . . N. F o r t h e c o n s t r u c t i o n o f a g r a d e d m e s h on a general p o l y h e d r a l d o m a i n see [22]. Let X N = S a0 . h be the space o f piccewise c o n s t a n t f u n c t i o n s on such a g r a d e d m e s h . W e c a n c o m p e n s a t e t h e effect o f t h e singularities by an a p v r o p r i a t e l y g r a d e d m e s h a n d get t h e c o n v e r g e n c e r a t e h?.,2.

THEOREM

3 [24]. Let h be sufficiently small a n d ot as in (17). Then we have I

116 - 6 , , 1 1 , , ,

h "t~-"

.-,,, ~ C

[h""

3 i f 13 < 2,~ 3 ;f t 3 > 2---d

where the constant C = C( fl ) is independent o f h. T o illustrate this resuJt we c o n s i d e r t h e screen F - - (0, 1): (cf. T h e o r e m 2). Let g2 be t h e e x t e r i o r o f t h e s q u a r e . T h e n , vt = 1 / 2 , A~ = 0.297 thus, ~ =- 1 / 2 a n d T h e o r e m 3 we have t h e c o n v e r g e n c e r a t e (~(h ° ' : - ' )

if / 3 < 3

and

(3(h 3':)

if / 3 > 3 .

F o r b r e v i t y we have c o n s i d e r e d h e r e o n l y the Dirichlet p r o b l e m ; for t h e N e u m a n n p r o b l e m see [24].

3.2. The h-p version with q u a s i u n i f o r m m e s h C o n s t r u c t a q u a s i u n i f o r m m e s h with width h a n d define the space X N = St,.t,(F ) o | d i s c o n t i n u o u s piecewi.~e p o l y n o m i a l s o f d e g r e e p o n this m e s h . T h e z again, t h e rate o f c o n v e r g e n c e o f t h e G a l e r k i n s o l u t i o n is d e t e r m i n e d by the regularity o f the s o l u t i o n , i.e. the p a r a m e t e r a in (17).

E.P. Stephan t ('omput. Methods Appl. Mecla. Engrg. 133 (1996) lH.]-2Obl

18~

T H E O R E M 4 [19. 2b,']. Let p be sufficiently large and h be small enough. Then, the Galerkin equations (13) are uniquely soh, able in Sp.h(F ). Let 0 ~ H ~t"-(F) be the solution o f tire integral equation (3) with sufficiently smooth rigttt-hand side f, and let ~/'r,.,,E X x ' = Sp.,,(F) be the Galerkin solution of (13), then we have with a in (17) and arbitrary • > i ) -

,,11,,

,

Ch"'p

= " * =" .

(18)

R E M A RK. If N = dim X v denotes the degrees of freedom, then T h e o r e m 4 gives for the h - m e t h o d the rate L(N -") and for the p - m e t h o d the rate eu(N -'"). Corresponding results are given in [8] for the finite element method. In the screen problem (9) with F = (0. 1)-" we have t~ = min{~, i 0.297 + J~ } = i hence T h e o r e m 4 yields the convergence rate ~(h ~-" - ' ) for the h - m e t h o d and b ( p ~ i. ,., ) for the p - m e t h o d . C o r r e s p o n d i n g results can be derived for the N e u m a n n problem, too. For the two-dimensional problems see Stephan and Suri [301.

4. T h e

h-p

version of the B E M for 2D p r o b l e m s

4.1. Exponential convergence Let us consider in more detail the problems (1) and (2) when .(~ is a b o u n d e d , plane domain with polygonal boundary F = U~:~ J /-", F ' being open straight line segments. The interior angle at the vertex tj is denoted by (o~. It is well known that the solutions of (1) and (2) have special singular forms at the corners of ,(/ [10]. (Note for 2D problems a = min ~ in (17) with z,~= ~/to,.) W h e n these problems are converted into the boundary integral equations (3) and (4), then the solutions have corresponding corner singularities which reduce the convergence of the Galerkin schemes (13) and (14) when the h-version or the p-version are performed (cf. T h e o r e m s 3 and 4). Nevertheless the h-p version with geometrically spaced mesh (near the vertices of .O) converges exponentially fast. T o describe the h-p version of the boundary element m e t h o d let us first introduce a geometric mesh F~', on the b o u n d a r y I" of .Q. As in [9, 13], we refine geometrically the mesh towards each vertex and introduce on this geometric mesh FI', the spaces S"-I(F~',) of boundary elements as follows: s"-'trT',).=

{,t,!,!,.,=,

=

1.....

,r + l )

where Fk(F~ff) denotes the space of polynomials of degree <~k on the subinterval F'ff. T h e geometric mesh !'~' is obtained by bisecting each side F ' with length d~ into two pieces F ~_ (containing the vertex t,.~) and i " . (containing the vertex t,). T h e n , d i s t ( t j _ t , F ' f f ) = d j / 2 t r 0-~*~ for k ~ n + 1 w h e r e tr ~ (0, l) and n is an integer. T h e n . with the choice X,,. = S "'-t(l . ,,). .where. N .= dim . S " ~(F,,), we obtain the exponential convergence of the Galerkin method for the interior problem ( l ) as shown in [3, 41.

THEOREM

5. Provided the given data f in (1) is piecewise analytic, then there holds the estimate

lid, - +~,,.ll,, . -'~,,, ~ C e-h,,.v

(19)

between the Galerkin solution ~.,, E S"-~(F',',) o f (13) and tlre exact solution • = ou / onlr of (3) where the positive constants C and b > 0 depend on the mesh parameter tr but not on N = dim S"- ~(,,). F" Next, we con~ider the interior Neumann problem (2) and since W is only positive semidefinite, we consider as in [14] the system

E.I'. Stephan ,' ('omput. Method.s Appl. Meets. i-,gr£. !~3 r1996)183-208

189

(

Wv - - ( 1 - K ' ) g - - a

on I ' ,

j~ v d s = 0

(20)

which has a u n i q u e solution ( v . a ) E H ~ 2 ( ! ' ) × C. T h e c o r r e s p o n d i n g t e r m i n e (v,,,, a,. ) (Z Y.~ x C with

(Wv, , w} + (a,. w) = ((/(v,. t ) = o

Galerkin

scheme

reads:

De-

K')g. w) (21)

for all w ~ Y,. C t t I 2(!'). F o r the h-p version o f (21) we l a k e c o n t i n u o u s , piecewisc p o l y n o m i a l s on t h e g e o m e t r i c mesh 1"i',:

.

(I ,)A("(I')CI-I'

,,)=S

.

:(I')

(22)

wheFt2

S"(r::)=

{vlvi, ,:, e P ~ ( l ' " ~ ) . k

= 1. . . . .

n + I}.

T h e n . with the choice Yx = S"(I'I',) we o b t a i n ' h e e x p o n e n t i a l c o n v e r g e n c e o f the G a l e r k i n m e t h o d (21) fo~ !he i n t e r i o r N e u m a n n p r o b l e m (2) ( c o m p a r e B a b u ~ k a ct ~!. [ 3 . 4 ] ) .

THEOREM

6. P r o v i d e d the given data g in (2) is piece,vise analytic, there h o l d s the e s tima te

Iio - o , i l , , , - . , , ,

+ t- i - , I

~ c~

"' "

(23)

b e t w e e n the G a l e r k i n s o l u t i o n (v x, a . ~ ) E "( FI') x C o f ( 2 I ) a n d tire exact s o l u t i o n (v. a ) ~ H ' ' ( F ) × o f (20). A g a i n , tire p o s i t i v e c o n s t , n t s C, b d e p e n d on o hut trot o n N dim ( ,,).

C

~

4.2. I m p l e m e n t a t i o n F o r X~. C H - t ' : ( F ) the basis f u n c t i o n s can bc d i s c o n t i n u o u s at the grid p o i n t s o f t h e g e o m e t r i c m e s h .A : = F~',. W e can use linearly t r a n s f o r m e d m o n o m i a i s . F o r an clerr, cnt K C F~', with p o l y n o m i a l d e g r e e p~,o. the e n d p o i n t s ( a , . a , ) . ( b , . b~) a n d the t r a n s f o r m a t i o n

T h. " =

~ x

"

-_-

' 2[(hi

-

a,

al

)x+

+b,l

- ~ [ ( b 2 ~ a , _ ) x -l- a , _ +

,,7?

b2l]

we t a k e t h e f u n c t i o n s

(i = () . . . . .

A function f~.X,,, with m e s h F " , , = { K , ' i (Pl ..... P.,,ta~) has the representatiocJ ('n,( J )

K,

-, P:

aK,

w h e r e ajK t (5_ ~ , j = O . . . . .

K

p,. i = I . . . . .

m(J).

PK ) -

= 1

.

.

.

.

.

re(A)}

and

polynomial

degree

vector

p=

190

E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 133 (1996) 1613-20~

For Yx C H t ; - ' ( F ) the basis functions must be c o n t i n u o u s on F. We construct a basis using the affine i m a g e s of antiderivatives of L e g e n d r e polynomials, then these basis functions vanish at the grid points of F',,. T h e missing d e g r e e s of f r e e d o m at those points are delivered by the s t a n d a r d piecewise linear r o o f functions. T h e s t a n d a r d basis functions on the reference interval [ - 1 , 11 arc defined by t -

(-l)'x

O,(x) =

2

(i = O, i )

q,,(x) =

V S

I,_~(t) dt

(i = 2 . 3 . . . . )

-I

with the L e g e n d r e p o l y n o m i a l s I,. T h e r e holds tO,(-+ 1 ) = 0 (i >~ 2) and I

~' = k/2(2i - ~, j

(l, - I,_,.)

(i >~ 2 ) .

H e n c e , for 1"• YN there holds

f-"

I

m{~

I~

)

K,

-,

N p,

with a,h,' = a ~ ' - ' (i-- 1,2 . . . . . rn(A), K , , = K,,,jj) by continuity 1 . . . . . m(ZI)). H e r c . the L,~"s are the t r a n s f o r m e d basis functions K

(a~'~.

j=0 .....

p,.

i=

----, (i = 0 . . . . .

Pk ) •

In o u r c o m p u t a t i o n s there a p p e a r four integral o p e r a t o r s , the single layer p o t e n t i a l V, the h y p e r s i n g u l a r o p e r a t o r W, the d o u b l e layer potential K and its adjoint K'. All test a n d trial functions in the G a l e r k i n s c h e m e s can be r e p r e s e n t e d as c o m b i n a t i o n s of t r a n s f o r m e d m o n o m i a l s e,~ = QkA', ( d e g r e e k. s u b i n t e r v a l K, C FI°). H e n c e , it suffices to consider those functions. T h e single layer p o t e n t i a l t e r m s

are c o m p u t e d as follows: T h e o u t e r integral by 32 point G a u s s i a n q u a d r a t u r e , analytically m a k i n g use of the expression I

1 xt loglx - z I dx - / + 1 [x'+' - z t " l Ioglx - z l -

"

1 ~ / +-----i-

the inner integral

x ~- ~+2z ~-t t - k + 2

k-:l

F o r details sec I9, 251. For the t e r m s involving the h y p e r s i n g u l a r o p e r a t o r there holds d u c to integration by parts

,.e, > =c
.e,

>

with c o n s t a n t s c = c ( K i, Ki). H e r e , t;jt d e n o t e s the derivative with respect to the arc length of e~. N o t e in the t e r m s with the h y p e r s i n g u l a r o p e r a t o r the test and trial functions are c o n t i n u o u s t h r o u g h o u t F.

E.P Stephan

/

('otnpla. Method.~ Appl. Mech. Engrg. 133 (1996) 183-208

191

For the terms with the double layer potential o p e r a t o r we observe the relation

( Kg.

=

h-,,4>

A g a i n . the o u t e r integral is c o m p u t e d by Gauss(an q u a d r a t u r e whereas the inner integral is c o m p u t e d analytically (compare Ervin et al. 191), The linear systems (13) and (21) for the various versions are solved direcllv.

4.3. N u m e r i c a l results EXAMPLE 1. (Dirichlct Probicm) u = Im z-" ~ and 12 be the L-shaped region with vertices ( 0 , 0 1 , (0, 1). ( - 1 , 1). ( - 1 , - - 1 ) , ( I , - 1 ) , (1,(1). In Table 1, N d e n o t e s thc n u m b e r of u n k n o w n s of the G a l e r k i n a p p r o x i m a t i o n @N from (131. n is the n u m b e r of subintervals on the edges joining at the rcentrant corner of .Q. e~.: d e n o t e s the relative error in the energy n o r m

I1 e,: =

- 6.,,1(,:li il,:

with I1 , -

•=

(v(6

-.

0,.).

(6

-

',/,,.))

'

2

EXAMPLE 2. ( N e u m a n n Problem) ,', = Im z ' 7 and ~2 ils in E x a m p l e 1. % = II" - -NIl,.--,, ,/ii-II,.--,,, where u.~. solves (21). Note" In Tables 1 and 2 the experimental convergence rates or~, agree with the theoretical c o n v e r g e n c e rates a * for the h and p versions, w h e r e a s the error of the h-p version d e c r e a s e s exponenti~!!y underlining (19) and (23).

Table l R e l a t i v e e r r o r in e n e r g y n o r m f o r h - . p- a n d //-version. p = 0 a* = 0.66

h-p v e r s i o n o f ( 3 ) h-p v e r s i o n , o- = 0 . 2 5

p-version, 8 elem. ~x * = i .33

N

ej

a~

N

e,

8 16 32 64 128

I1.18811 11. ! 126 0.0710 (i.(t147 I).I1282

0.74 O.t~7 I).67 11.66

8 I6 24 32 40

II. 188(I I).|)874 0.0543 11.(1375 I).I1273

Table 2 R e l a t i v e e r r o r in L " - n o r m f o r h - . p- a n d

~ 1.11 1.I7 1,29 ! .43

N

I 2 3 4

8 18 311 44

I[I.%' ik 188(I O.0699 I).112,12 O.(li(Xl

1.22

2.08 2.31

h-p v e r s i o , o f ( 2 0 ) h-p v e r s i o n , rr = 0 . 1 5

p-version. 8 clem.

h - v e r s i o n , p = 0. " = (I.64

n

N

t°,,

O[,,

C,,

8 16 24 32 4(1 48

0.1)241 0.11149 {I.0115 (I.1x195 (1.0(183 0.(X174

O.70 0.64 0.65 (1.64 0.64

0.(1241 11.(1136 0.01194 o.(x)71 0.1X156 0.01147

O[., (t.83 liAR) 11.98 1.03 1 ,o7

11

N

t

8

2 3 4 5

18 3(1 44 60

e.

11.(12411 II.011(-192 0.00217 O.(RH~7

0.(X1023

0¢%. 1.54

2.27 3.1)7 3.45

E.P. Stephan 1 Comput. Methods Appi. Mech. Engrg. 133 (1996) lg3-20bl

192

5. The

h-p

version of the BEM for 3D problems

5. I. Exponential convergence Let I" bc the boundary of a simply connected bounded polyhedron or a screen 1"~ in ~". We first introduce appropriate countably normed spaces B~(F), (0
B~(F) =

{,,I,,

~ HI-t(F): ::lc. d > 0 i n d e p e n d e n t of k such that

I!rf +".-'(o,(~,,- o,))" +~. •

.!}

.(alar,)~'.(alaoi)".ullt.:,v,,

<~Cd"'+""-t(a,+ao-l)!

BL(r) =

{-I,, ~

tr + a o = k = l . l + i

.....

i=1 .....

m}.

(24)

H ' - ' t r ) , ,, ~ Bt~(F) for all faces F of F }

with (a)+ := max(a, 0). If we would like to emphasize the dependence on the constants C, d we will ! write B ~ ( F ) = B , . c . , ( F ) , etc. Now we define the geometric mesh on the faces of the polyhedron. We assume that the face F is a triangle. There is no loss of generality because every polygonal domain can be d e c o m p o s e d into triangles. We divide this triangle into three parallelograms and three triangles where each parallelogram lies in a corner of the face F and each triangle lies at the edge of F apart of the corners. By linear t r a n s f o r m a t i o n s ,¢, we can t r a n s f o r m the p a r a l l e l o g r a m s on to the reference square Q = [0, 112 such that the vertices of the face F are t r a n s f o r m e d t o ( 0 , 0 ) . T h e triangles can be t r a n s f o r m e d by a linear t r a n s f o l m a t i o n ~i on to the reference triangle Q = {(x, y ) E Qly-<-x} such that the c o m e r point of the triangle in the interior of the face F is t r a n s f o r m e d to ( I , 1) of the reference triangle. By the following definition the g e o m e t r i c mesh and a p p r o p r i a t e spline spaces are defined on the r e f e r e n c e e l e m e n t Q. A n a l o g o u s l y , the g e o m e t r i c mesh can be defined on the reference 'triangle Q (see Fig. 1). Via the t r a n s f o r m a t i o n s ~ ¢ - t , ~ 1 the g e o m e t r i c mesh F~ is also defined on the faces of the p o l y h e d r o n . T h e a p p r o x i m a t i o n on the reference square is the m o r e interesting case because it handles the c o r n e r - e d g e singularities. N o w let us introduce the spaces of spline-functions on g e o m e t r i c meshes. For a given 0 < or < 1 we use the partition I~', of ! = [0, 1] into n subintervals (xj_ i, xj), j = 1 . . . . . n, w h e r e

x. = 0,

x, = ~r"-',

] = 1. . . . .

n.

(25)

L,,, Fig. 1. Geometric mesh for tr = !/2 and n = 5.

E.P. Stephan / Compm. Method,: Appl. Mech. E,grg.

133

(1996) 11q3-208

193

W i t h ii', we associate a d e g r e e - v e c t o r p = (p,, . . . . . p,, ~) and define SV~(l:',) as t h e v e c t o r s p a c e of all piecewisc p o l y n o m i a l s v on / having d e g r e e p, on ( x , . x , , , } , / = 0 . . . . . n - l . i.e. 0 ] ~ , . ~ , . , ~ P,,((x,, x,. t)). For a given {I < er < i we use the p a r t i t i o n 0~', of Q = [0, 1] 2 into n-" r e c t a n g l e s R ~ n t . l = Ixl.

. i X.L ] X. IX.! , . X l ] .

.~ k . l :.

.I .

,~)

O -- I~j /,.t

!¢~, .

(26)

I

W i t h Q,,~l we associate a d e g r e e vector p = (p,, . . . . . P,, t) and define $, p. I ( Q , I1, ) as the v e c t o r s p a c e of all piecewise p o l y n o m i a l s v on O having d e g r e e p~ , in x a n d Pt t in y on [x~_~.xkl x [ x l - , . x t l . k . l = 1 . . . . . n. i.e. vii,, ,.,,I-l,, , . , , i E l ~ ,, ,.p, ,(R tel. F o r the d i f f e r c n c e s h k = x k - x ~ _ t we have

h~ = x ~ - x ~ ,

(5)(5)

"-x~_,

- !

~x

- !

=x7.

VxC[.r~

~.x~],

2~k<~n

(27)

with 3 , = 1 / o ' - 1 . T h e index I.~ {0. t} in St'¢(t~',) (rcsp. St'J(Qi',)) d e t e r m i n e s the r e g u l a r i t y of the piecewise p o l y n o m i a l s , i.e. d i s c o n t i n u i t y in case of I = () and c o n t i n u i t y in case o f I = 1. T h e n we have by c o n s t r u c t i o n ' ~ . , . ,~ s , . , ( Q ; ; ) ~ s ~ . , ( t , ,,. ) x ~~,..,,,,,

(28)

a n d via t h e t r a n s f o r m a t i o n s ,¢,. ~, the c.~rresponding b o u n d a r y c l e m e n t spaces S t'. I ( / ' .n, ) a r e o b t a i n e d o n t h e g e o m e t r i c m e s h F~',. F o r the c o u n t a b l y n o r m e d spaces t h e r e holds the following a p p r o x i m a t i o n results.

T H E O R E M 7 [20. 18]. (i) Let u(x. y ) @ B ~ ( Q ) with [ ) < / 3 < I. /3 # 1 / 2 . Then there is a spline f i m c t i o n u N E S P " ( Q ~ ' , ) and constants C t . b ~ > 0 independent o f N. but deFendent on ~r. Iz, ~ such that

Ilu(x. y) - u , ( x .

y ) l l , -',~,, -<- c , e " " '

'

with p = (p,, . . . . . p,,). p~ = [ # k l f o r a il > 0 and N = dim S'"J(Q",). (ii) Let o ( x . y ) E B ~ ( Q ) f q C ~ I ( Q ) with t ) < / 3 < 1 , /3#1/2. Tiwn there is a spline f u n c t i o n o N E Sp'I(Q~',) and constants C2. be > 0 independent o f N. but dependent on or. IX. [3 such that

llv(x, y) - o,.(x, y ) [ I , , , , , , with p = (p,, . . . . .

_ 1~ .;%, i

4

<- c , e

p,,). p , , = 1. p , = m a x ( 2 . [t~kl + 1) ( k ~

I) for a t~ > 0 and N = d i m S ' " ( Q : : ) .

In t h e following we w a n t to a p p l y these results to the G a l c r k i n s o l u t i o n s of t h e b o u n d a r y integral e q u a t i o n s (3) a n d (4). First. we c o n s i d e r the s c r e e n p r o b l e m (9) a n d take for simplicity I ' = [ 0 . 2 ] " a n d divlde into 4 pieces b y s y m m e t r y (for t h e g e n e r a l case see [ 1 8 . 1 9 ] ) . N o w we c o n s i d e r t h e G a l e r k i n e q u a t i o n s to (9) w h e n X.,; = Sz'"(Q',',). N = d i m S t ' " ' ( Q ~ ' , ) . E s t i m a t e (I I) causes us to look for a g o o d a p p r o x i m a t i o n in X x to o b t a i n an u p p e r b o u n d for t h e G a l e r k i n e r r o r . D u e to t h e t e n s o r p r o d u c t c o n s t r u c t i o n o f St":~(Q~',) a n d t h e special d e c o m p o s i t i o n o f ~/, in T h e o r e m 2 t h e a p p r o x i m a t i o n p r o p e r t i e s of St"°(Q~',) can be studies terra by t e r m . w h a t will be i n d i c a t e d in t h e following. T h u s . (due to T h e o r e m 2) we n e e d only to look -~t the b e h a v i o r n e a r t h e origin. H e n c e . let us c o n s i d e r the m o s t singular part o f t h e exact solution ;:ear the origin, n a m e l y the c o r n e r - e d g e s i n g u l a r i t y in (0, 0) a n d the e d g e - s i n g u l a r i t y for the edge ( 0 . 0 ~ . ( 1 . 0 ) . T h a t m e a n s we have to a p p r o x i m a t e t h e function t/re : = r* - 10 t.,_, + x.~ - t, 23,- t~2 n by a s p l i n e - f u n c t i o n ~ E S p "o (Q,,). The term

(29) r x -- t o t : 2

has a singular b e h a v i o u r for r - ~ 0 a n d 0---* 0 a n d

E.P. Stephan / Comput. Methods Appl, Mech. Engrg. 133 (I906) lbt3-2Obl

194

y-I,~, corresponds to the edge-singularity. The function x ~-I'-" has the same behaviour as the stressintensity distribution. We use a similar decomposition of ~. like (29) to approximate the terms separately, (3o)

gtu "= dp~(r, 0 ) + a ~ ( x ) ~ ( y ) .

Note that ¢b~ is a two-dimensional spline, while a~, and ~e are splincs on a one-dimensional cut of the two-dimensional mesh parallel to the axis, i.e. ~., C S"('(Q~',) and ax, qJ~,E S'"(I~',). From [18, 20 i we know that r"O" C B ~ ( Q ) for a > - B ,

v>I/2-B,

(31) i.e. with , 1 , = A - 1 = - 0 . 7

and v = l / 2

we have

r*-'o"zGB~(Q)

(32)

with B > max( 1 - A, 0) = I --A ~ 0 . 7 . Thus, due to T h e o r e m 7 there exists a $~, ~ S""'(QI',) such that for amax(I/2, l-A)=l-A<13 0

with constants C I , b, independent of N = dim X.,,.. (ii) Let g in (2) be piecewise analytic attd let v be defined by (4) and ux, its Galerkin a p p r o x i m a t i o n , t)e defined by" (14) with Y,,. = S t'. (i ,,) on geometric meshes F,, on F. Then there holds ]'or arty [

iiv - v , I1., ..,,,

c._ e

- --~

Pl

st

+ e4N-").

with constants C,, b: independent o f N = dim Yx. 5.2. Implementation Here, we consider only the case that the polyhedron 1" consists of surfaces which can be completely divided into rectangles (no m a n g l e s are needed). Since we use tensor products of Legendre polynomials as test and trial functions in the Galerkin schemes all arising integrals can be computed analytically.

Appl.

F.P. St,7dmll / ( ' o m p m . Mt,th,Ms

Mech. l-ngrg. 133 (19963 l g 3 - 2 0 g

195

N o w the rectangles which are the s u p p o r t s of trial functions may have 3 different o r i e n t a t i o n s to each o t h e r . N a m c l y , (i) both rectangles arc in the same plane. (ii) the rcctangles are in parallel planes, (iii) the rectangles are p e r p e n d i c u l a r to each other. This leads to 2 different antiderivatives for the e l e m e n t s of the G a l e r k i n matrix: (i)

(ii)

I"~ ~..... ( x ~ . x 2 .

j

Yr. Y 2 " x ~ -

.v ~ ) " =

Fh~,,,,(xi,x Z. y:..v3;x ~. yt )--

fff dx ,

th" :

dy,

IfJf tl.r~

dx 2

-rlx2~ ~ Y2

d.v.

dy2

¢,

","'"

-r~J:z.} 2

dy~

;

*'"'"'

)" s

--. )2

Vz(,, v,

in all cases the integrals are of convoluti(,n form D,,(x. vl =rj__d_x__.[dv x ~ _ v ' D ( x - y). w h e r e D ( x - y) can be a general kcrnel, e.g. here we have D ( . r - 3 ' ) = i / V ( ( x - - 3 ; ) - " + a-)." H e n c e , the n u m b e r of integrations car, be reduced. D,,(.r.y)= -

; fdx d y x ,,y D ( x - 3 " ) i + l! y ,.,f d.r , x

D(x

-" ~. ).

if

. . /. .~- I

d v.

v"'. .f h,

dv . - 1 +I1 { .v ; ' l f .d x . r ~ D. ( x - - v ).+ . r 'I,'

&r

.

.,-'D ( x - v )

~D(x- v)--kD~

i.t t ( x , y ) . . . .

}

In case (i) we apply the reduction twice and obtain integrals of the form "

GZ/(y |,

f dv "J d r .

v.;x..a)"= ""

-

''

J

/

YIY2

__

"- \~y, -.r,)Z +(yz-x._)'+a-"

in case (ii) we can only apply this reduction once and obtain

j

-

6£,,,,(y,.y._.3.,:.r,.x._..~-,)=fdy, f

d)'2

y~

dY3 / - 2 -

t

,,,

Y2Y.~ (X I

.... .Vt ) 2

With the antiderivativc g1(Y" x. a ) ' =

I dy ~ ( y

Y

, x) z + a:

wc can write G , ; and G£~,,, as

G£t()',. y.'x,.x?.,a)=

f

d y : y 2 Ig , ( y , o x ~ . ~ i v , - x ; ) ' - - t - a ' - )

and

f G£I,,,(Y,. Y2. Y3" x , . x.,. 3'3) = j dy3

y';'G;t(y

1,

Y2" "lt°l " "t"2" 3'&

--A'-~)

.

T h u s . with the r e c u r r e n c e relation V - - .It"

g,,(y: x. a) = arsinh "

g l ( y . x . a ) = - [ { y]

i-1

t-I V~-) , - x ) Z + a - + (" 2 1

- l)xgt i ( y ; x . a ) - ( l -

we get also r c c u r r e n c e relations for G£; and G~;,,,.

l)(x2+a')gt.._z(y;x.a)}

E.I: 3tcphan ! Comput. Method.s Appi. Mech. lfngrg. !.1.1 (19961 lSJ-2ob¢

19(,

5.3. Crack problems Crack p r o b l e m s in linear elasticity lead to systems of integral equations of the first kind which can also be solved by the h-p version of G a l e r k i n ' s m e t h o d , e.g. the inclusion (Dirichlet) p r o b l e m on a plane crack surface /" leads for given ff~' to 0

ex-yl CI

o

1 (x 2 --y,_)-" Ix - y[ + c, ix _ y l '

0

c,

(x2 - Yz)(x~ - y~)

-

Ix - yt

(x_. - y_.)(x 3 - Y3) Ix - yl ~ 1

-Ix- - +- cyl, _

~"(Y)) ~:((y)Y)

ds,

(x: - y , ) ' lx._yl~

[ g'l'(x)~ xEF

(34)

\g'?i(x)/ where the constants cz, c, d e p e n d on the L a m d constants. A gencralisation of the analysis in [16] will yield again exponentially fast convergence for the h-p version G a l e r k i n scheme with g e o m e t r i c m e s h as in T h e o r e m 8. Note that the necessary explicit decomposition of the solution into edge and c o r n e r singularities can be found in [22].

5..1. A dapative h-p algorithm for screen problems (a) For the weakly singular integral e q u a t i o n we proceed as follows: Starting with s o m e G a l e r k i n solution ~,,, one c o m p u t e s the local residual

R~ "= fl

(f(x)

- VO:~,(x)) ~

For given 0-< 0 ~ 1, if

dr. ?

R~ <~OR'-,.... the G a l e r k i n solutior~ is already accurate e n o u g h , if R, > OR'-,,.... for

s o m e e l e m e n t s F, we call this e l e m e n t critical. In the case of rectangles (triangles) we can also p e r f o r m a direction control on the critical elements. In the case of rectangles t e m p o r a r i l y we divide every critical e l e m e n t I~ uniformly into four e l e m e n t s (see Fig. 2) and take a closer look at the four local residues

Fi.3

Fi.4

r,. L

F,,2

Fig,. 2. D e c o m p o s i t i o n of an element In determine the refinement direction.

( ,,.~lmt. Meth.d', .,|pI d

I'r/ '~. .(~l~rl'J/¢lfl.

R~.~ =

f

[ ([(x)at;~

Vqrx(x)) 2 d x ,

k = I.....

3 I t c h l-m;t~.

4.

1.13 (H)96) lb13-2tt,'¢

197

(35)

Let /.t > 0 be a given n u m b e r (/.t = 1.5 by default), called the direction p a r a m e t e r , and define the two indicators

ex,:=max

,

R,. I + R 2

RT.,- + R .~

f'~ R,.2 + R~.4

R.t +

R

)

,2 -~



(36)

T h e n , p e r f o r m a refinement in x-direction on I: (x-refinement) if e_r, ~ ey, or ex, > p. and a y - r e f i n e m e n t on 1-~ if ey, ~> ex, or ey, >/~. T h e critical e l e m e n t s F, have to bc refined or the polynomial degrees have to be increased. This is decided as follows. Firstly. one c o m p u t e s a new G a l e r k i n .~olution tkx with the ~ame mesh but with p o l y n o m i a l degree increased by one on the critical e l e m e n t s I~. Now one c o m p a r e s the old with the new residual

/~, ' = f,. (f¢,:) - v ~ , ( x ) ) z ct,-. t

T h e n , i f / ~ ~< y R 7 on a critical e l e m e n t one performs a p - r e f i n e m e n t whereas if/~ ~ > y R ~ on a critical e l e m e n t one p e r f o r m s a h - r e f i n e m e n t , where y is a preset p a r a m e t e r . T h e h - r e f i n e m e n t is d o n e by interval halving when T, is an interval and by subdividing into rectangles (triangles) when C is a rectangle (triangle). If the rectangles are o b t a i n e d by just dividing the critical c l e m e n t I~ into 4 s u b s q u a r e s we s p e a k of an adaptive routine without xy-control (otherwise with xy-con~rol). T h e p - r e f i n e m e n t is d o n e by increasing the polynomial d e g r e e by 1. Dividing the e l e m e n t s into more than 4 subsquares o r increasing the p o l y n o m i a l degree by m o r e than 1 usually gives a lower convergence rate because we introduce in this way m o r e degrees of f r e e d o m away from the singularities of the solution than necessary. If we use a direction (or x y - ) control, then we have to do the refinement in the indicated direction, i.e. a p p r o p r i a t e halving of the rectangle (triangle) in the case of a h - r e f i n e m e n t and increasing the p o l y n o m i a l d e g r e e by I in the case of a p - r e f i n e m e n t . (b) For the hypersingular integral e q u a t i o n wc proceed as follows: Starting with s o m e G a l e r k i n solution v~. we c o m p u t e an a p p r o x i m a t e value for thc local residual n~ "= f,. (g(x)

- Wv,(x))'-

dx .

Since v vanishes on the edges of F = [ - 1 , 1] 2 thcrc holds with the surface Laplacian A r for x E/-" Wv(x;

'-

,~n., - ~

. I. . x .. - y l

v(y)dy

=

--A,

. Ix-y] dy

T h e r e f o r e we a p p r o x i m a t e R~ by / ~ "= f , . ( g ( x ) + A ! r V v x ( x ) ) 2 d x using the 5 point difference star ~1" With this trick we avoid to c o m p u t e directly the action of the hypersingular integral o p e r a t o r W. N o w the adaptive algorithm can be p e r f o r m e d as above (see 1161).

E.P Stephan / Comput. Methods AppL Mech. Engrg. 133 (1996) 183-20~

198

I t • ¢, = 0 . 2 5 °o=0.2

-,m

Iiu,llo.

• ~= e.25 • = 0.2 0.171

10-~.

9

10-'~ i .q

10-3~

I0-~

1O-." ! 1

. . . . . . . . . Fig.

3.

L-'-errors e

"

"

"

~

"

"

"

"

'

1'0

~'o =

Ilu. -...tl,

~... ~ = It~ -~'11,

i, .-

0.04001 0.03001

0.01!

O ¢:: t~

¢o

Q. c¢1

(.3 ¢-

uJ

0.00001

I 445

Number of Unknowns

Z~

T7

<> [3

graded mesh. beta = 4.0 h-vers~..o-O I>-versk>n,I~1 I~VemJon. h=0.5

A

mesh,

~

= 0.17, mU=1.0

II

geomemc mesh. r,~g,-na = 0.5, mu=l.O ~pUveht. h - ~ (.~t~xN xy-comro0 ~ t l v e t y . ~ v e r s ~ (w~h W-oomoO

o

geometric mesh. sigma = o.17. mu=0.5

V

41.

~oUv~.

hp-vorraon (with xy-c~roO

Fig. 4, Dirichlet problem on tile L-shaped pIate.

885

17~5P_205

E.P Stephan i C'otnput. Methods Appl. Mech. Etagrg. 133 (1996) lg3-201g

5.5.

199

N u m e r i c a l results

As described in Sections 3. l, 3.2, 5. i and 5.2, we implement boundary element Galerkin schemes to obtain approximations ~/,,,. to the solution ~ of (9) and approximations u~v to the solution u of (10). We consider the screens given by the L-shaped plate in Fig. 6 and the square plate with sidelength 2. In all cases we choose the f u n c t i o n s f = g = I on the right-hand side in (9) and (10). For the Dirichlet problem we look for the capacitance C of F which is given by C = l / 4 r r fr ~/~ds where the charge density ~, solves the integral equation (9). As approximation we compute Cx = 1/4,rr J'r @N ds with the Galerkin approximation ~.~, of ~. Due to f = 1 there holds with constants c=, c 2 > 0 c,114, -

tc-

c.,.l

-

(37)

,-,,,

0.06001 o.05001 0.04001 -0.03001 -0.02001 -

o.o1~1-

¢D ¢,3 ¢-

o a. rO c lu

0.0o001

! 4

NAI/4, N: Number of Unknowns

A v I t

georn~L,'tc mesh, r,lgma = 0.17, re,u=1.0 g,eometrlc mesh, t ~ = 0 5, mu=l.0 ~qcztm.my. h--ve,-sk~ ( w t U ~ xy-,cx~m0 ,m:tat~ve~. h-ve,,',~,k~(wtth W ~ zm,zq~veq, hD-verr,k ~ ( w ~ xy-oomrm) geomem¢ r n e ~ , r,z,grrm = 0.17, m u ~ . 5 Fig. 5, D i r i c h l e t

problem

~n t h e L - s h a p e d p l a t e .

6

E.P. Swphan / Comput. Methods Appl. Mech, Engrg. 133 (1996) 183-208

200

yielding with (18) the error estimates for the pure h- and p-versions

l

h

| --~

(38)

] C - C N I ~ c ~ p -2÷2e

with c 3 > 0 and arbitrary • > 0 . For the h-p version with geometric mesh, Figs. 4 and 5 show the predicted exponentially fast convergence (cf. Theorem 8), whereas they show only algebraic convergence for the pure h- and p-versions as predicted by (38). In Fig. 6 we present typical mesh sequences and degree distributions for the adaptive hp-version with direction control. For the Neumann problem we consider the relative error in the energy norm i l v for the solution v of (10) and its Galerkin approximation ON. There hold the error estimates [24.26]

llv

fh

-

t/2-e

(39)

VNlln, r--','", ~< C~[p--]+ ; . .

Capacitance=0.6241132

Capacitance=0.6303089

P o4

oo

oO

O0

O0

O0

O0

O0

oo

GO

o0

oo

oo

O0

O0

O0

oo

Capacitance'-0.6411131 |o

oo

oO

O0

0o

[email protected]

Capacitanc¢-~0.6456763

oo

O0 i

!o

oo

1o

oo

i

|0

I1

10

[email protected]

O0 O 0

-

,

Oo 10

=F

Q1

01

Ol

110

I0

20

O0

411

O|

Capacitance=0.6506758

OO

1

10

OJ

CapsrAtance=0.6530956 2lO

'11

O~

"at---

30

O0

10 2O

12

1o 20 02,

:;13

Fig. 6. The adaptive hp-vcrsion with direction control for the Dirichlet problem ( C = (}.659265).

E.P. Stephan

/

201

Comput. Methods Appl. Mech. Engrg. 133 (1996) 183-208

r" UJ

c: L_

o t...

UJ Q

>._ (1)

rr

0.003

|

5O4

::

I 1004

I 1504

Number of Unknowns

A V

gracm~ n'Nmh, l:~ta = 4.o ~ p,,,o

,& v

ge0fTlef]l¢ mesh. ~ mE0,17, t o t a l , 0 geometric mesh, eroma ,: 0-% r n u : l . 0

D

p-vernon, h=1



O e o r ~ mesa, m ~am~,.

,: o.17, mu-o_~ h-~,,mk: (win xy-txxW~

Fig. 7. The N e u m a n n problem on the _~luarc plate.

For the hp-version with geometric mesh, Figs. 7 and 8 show the predicted exponentially fast convergence (cf. Theorem 8), whereas they show only algebraic convergence for the pure h- and p-versions as predicted by (39). We have also plotted in those figures the convergence rate for the graded mesh (cf. Theorem 3). Next, we present some numerical experiments for the weakly singular Eq. (3) on the L-block given in Fig. 10. Again, Fig. 9 shows only algebraic convergence for the pure hand p-versions and exponentially fast convergence for the h-p version on the geometric mesh. Fig. 10 shows the mesh of the adaptive h-version with xy-control. Finally, Fig. 11 shows the numerical experiments for the system of integral equations (34) which corresponds to the inclusion problem for the Lam6 system.

E,P. Stephan / Comput, Methods Appl, Mech. Engrg. 133 (1~u36) !#3-20bi

202 O.S03

"

0 . 4 0 3 -O.303-

0 . 2 0 3 --

E

0.103-

o

c

II) ruJ rL_

o

UJ ¢D

._> G)

tr-

0~0~



1

2

6

NAI/4,

N: Number

of Unknowns

geocrNm~ mesh. sk~ma z 0.17, mu=1.0 g e o c n e v ¢ n~r.h, s~3ma = 0.5. m u = l . 0

¢~omemc me~. r~0ma = o.~ 7. n ~ O . s II

adamh~, h-re.on (webW-tonal Fig. 8. The Neumann problem on the square plate.

T o solve the linear systems (13) and (14) for the different versions of the b o u n d a r y element m e t h o d we used the conjugate gradient ( C G ) algorithm. In case of the p-version we investigate various kinds of preconditioning: (i) diagonal preconditioning, (ii) preconditioning by the diagonal plus the subblock of piecewise constant or piecewise linear functions in (a) one variable or (b) both variables. The numerical experiments show that preconditioning reduces significantly the n u m b e r of iterations which are necessary to reach a given precision. !11 Fig. 12 the n u m b e r of iterations are plotted versus the n u m b e r of u n k n o w n s for the p-version to solve the Dirichlet problem on the square plate. For the h-p version we have always used CG with diagonal preconditioning. For further details and further numerical results, see [16, 18, 19].

2(13

E+F. Stephan / Comput. Methods Appl. Mech. Engrg. 133 (1996) 183-2(18 0.1430.123,

0,083

--

0.063-

0.043o O~ UJ

o.o23 -

6

o UJ

>o q}

ft.

0.003

I 418

22

i 814

m 121o

I

l

I

I

1606 2 0 0 2 2 3 9 ~ 9 4

N u m b e r of U n k n o w n s

LI V []

h-versk)n, p~0 p-vsf~on, h=:l p-Version, rr-0.5 graded mesh. beta = 3.0



geometric mosh. f ~ t a = 0.17

v

geometric mesh, s~mo = o.5

t

~ .

m

~m~.,~. ~ . a o .



~th~,



geometric ~ ,

n-~

(wm~ovtxy-~ntr~) (w,h ~-,x,.t,,~

hp-ve.don(wah xy-oon,~) skarnm =0.17, mue = 0.5

Fig. 9. Dirichlel problem on the L-block.

6. Coupling of F E M and BEM----&-p version

H e r e , we consider the h-p version of the F E M / B E M coupling for transmission p r o b l e m s in ~2 (see [14]). For simplicity, we consider the interface p r o b l e m -Au=f

in 12,

Au=0

in ~ 2 ~

with given j u m p s [-ll,-= -,,

and

r u!l

t 3rt - r = t,,

and growth condition u ( x ) = 13 Ioglxl 4-~(1) where 13 ~ C is not known in advance.

on F = c3.(2 (40)

as I x l - '

E.P. Stephan 1 Comput. Methods Appl. Mech. Engrg. 133 (1996) lbt3-2Olq

21)4

i I 1 f

i

I

!.

_/

J

I

I

Hq

r--7---7--~

,_L _'L

H

j

t r

/

---J--I



J'



3B

Fig. Ilk T h e a d a p t i v e h - v e r s i o n with d i r e c t i o n c o n t r o l for t h e Dirichlct p r o b l e m o n t h e L-block.

It is well known that (40) has the following variational formulation: Find u ~E H~(O). ,b E H - ~ ' - ' ( F ) such that for all o EEH ' ( [ / ) . ~/,E H-~"(F) there holds with a(u, v)"= .[. VuVv dx

B((~,), (~,))= L((:))

(41)

where

,, ((:;). (;,))--~°~,,.

o~ + ~,,'-

,~,~ + ,,,,I,. o> + <~,,-,~.t,

- v,/,, ¢,)

For X x C H ' ( O ) and Y.~tC H-~'2(F) the Galerkin scheme reads: Find uN E X N and @M • YM satisfying

t~ 4,,,, ( ~ ) ) = L ( ( : ) )

V eX,..,t, eY~,

(42)

Next, we consider the h-p version on a geometric mesh. For simplicity we consider only the case where F is an L-shaped polygon. Let .f-t"., be a geometri,- ~nesh on l-/, refined at the reentrant corner. Then, O~', induces a geometric mesh F~', on F = O/-/. O n the finite element mesh .f/,'~ we approximate u by tensor products of antiderivatives of Legendre polynomials together with continuous piecewise linear functions, this gives the space ~,,.l(O~). O n the boundary element mesh F~, we approximate Ou/On as

E.P. Stephon / Comput. Melhod~ Appl. Mech. Engrg. 133 11906) IH3-20H

2115

0.41

E O ¢--

¢-

I.U

i,._

LU G)

.>_ G)

n-

0.005

|

i

I

5t5

4

1026

I

1"

I

!

1537 204~559070

Number of Unknowns

A V ~3

graded mesh. beta = 4.0 ~vemaon. p~O

,k me

p-Vemk~n.~ i

9)

mesh, ~

= 0 . 1 7 , tlnU~=1.0

adapttv~. ~ v e ~ o n (v ~t~ xy-oomra) geometric m e ~ , .s~gnw~.17. mu=0.5

Fig. I1. Lame problem (inclusion) on the square plate.

in Sections 4.1 and 4.2 by elements in S"-2(F~). T h e n , the above symmetric coupling procedure converges exponentially (see [151). THEOREM Ilu -

9. L e t f. u o, t o p i e c e w i s e analytic, then there h o l d s the e s t i m a t e

u~ll.,,,., ÷ 114, - 4,MII,,-,

,-,,, <~ C(e

-~'~ + e -~'~)

b e t w e e n the G a l e r k i n s o l u t i o n u N ~ XN = ~,,. t (-Q~r), ,, 4,,u ~ Y.~ = S " - ~( F ~ ) a n d the exact s o l u t i o n u, ¢k = ou/ anlr o f ( 4 1 ) , w h e r e the p o s i t i v e c o n s t a n t s C, b~, b 2 are i n d e p e n d e n t o f N = dim X~ a n d M = dim YM. REMARK. For transmission problems in ~s, G u o and Stephan [12] show exponential convergence for the h - p version of the F E M / B E M coupling, Ilu - - N i l . , , . , ,

+ 114, - 4,MIl..-. ~,,-, ~< f i e - ~ ' ~

+ e-b~"~)



with corresponding finite element and boundary element spaces X~ and YM. T h e coupling a p p r o a c h can

E.P. Stephan / Comput. Methods Appl. Mech. Engrg. 133 (1996) IKI-20B

2116

(/1

£: o

a)

.,...

fD

~_~--:~ -~_---~-~,~i~-~---~~ ~___:~~--:--~ ~_~L-:-2~ ~ f • 1~ 16

!

l

645

11~

l'

|

!

t

!

I



Number of Unknowns

Fig. 12. Dirichlet problem on the square plate in R" 3. Number o1" iterations for the CG-algorithm (p-version. 4 elements).

also be applied to interface problems with a non-linear differential o p e r a t o r in the b o u n d e d domain g~ (see [27]), e.g. the Hencky problem in non-linear elasticity. For brevity we consider here the problem

-div(p[~Tutt)Vu ,) + ut = f

in O ,

-Au2 = 0

in ~z~.0

with U, = U, + u , ~ .

~3u t

Ot~ z

P ( I V U ' I ) ~ n -- ~n- +t'~

on F,

a n d . : - - Ioglx! as ixl - " = .

This leads to (41) with a(u,

w):=

,

2+

1 + Iv,,I

VuVw + u . w

dx

p ( t ) ' = 2 + (1 + t)-'

E. P, .~'tephan ! Comput. Methods Appl. Mech, Engrg. 133 (1996) 183-208

207

2 Let , ¢ 2 = [ - 1 , 1 ] " and tt 1 =(2--X21--X2)' 4'3` U 2 = ~ Iog(x "~ 7 + x 2 ) , and the mesh ,O,~ be geometrically graded towards the corners of 12. T h e n we obtain exponentially fast c o n v e r g e n c e in the energy norm for the G a l e r k i n solution of (42). Fig. 3 shows L-'-errors for meshes with a different a'. EXAMPLE.

Acknowledgmenls T h e a u t h o r would like to thank M. Maischak for performing the p r e s e n t e d numerical results. T h e a u t h o r was partly s u p p o r t e d by the G e r m a n Research F o u n d a t i o n u n d e r G r a n t Nr. Ste 573/1-1.

Reference~ I!1 B. Andersson+ I. Babugka and T+ yon Petersdo,'ff. Numerical treatment of vertex singularities and intensity factors for mixed boundary value problems for the Laplace equatio', in ~ ' . SIAM J. Numcr. Anal+ 31(5) (1994) 1265-1288.

121 I. Babu~ka and B. Guo. The h. p and h-p version of the finite clement method: basis theory and applications. Adv. Engrg. Software I5 (1992) 159-174.

131 I. Babugka. B.Q. Guo and E.P. Stephan. On the exponential conucrgence of the h-p version for boundary element Galerkin methods on polygons. Math. Methods Appl. Sci. 12 (1990) 413-427.

141 !. Babugka. B.Q. Guo and E.P. Stephan. The h-p version of the boundary element method with geometric mesh on polygonal domains. Comput. Methods Appl. Mech. Engrg. 80 (I99t}J 319-325.

I51 M. Costabel+ Boundary integral operators on Lipschitz domains. Elementary results. SIAM J. Math. Anal. 19 (1988) 613-626.

161 M. Costabel and E.P. Stephan. Coupling of finite and boundary element methods for an elastoplastic interface problem, SIAM J. Numer. Anal. 27 (19911) 1212-1226.

171 M. Dauge, Elliptic Boundary Value Problems on Corner Domains (Springer-Verlag, Lecture Notes in Mathematics, 1341, 1988).

181 M.R. Dorr. The approximation of solutions of elliptic botmdary-value problems via the p-Version of the finite element method, SlAM .I. Numer. Anal. 23 (1986) 58-77.

191 V.J. Ervin. N. Heuer and E.P. Stephan, On the h-p version of the boundary element method for Symm's integral equation on polygons, Comput. Methods Appl. Mech. Engrg. 110 (I993) 25-38.

[to] P. Grisvard. Elliptic Problems in Nonsmooth Domains. Pitman Publishing Program (L985). !!I! B. Guo. N. bhzuer and E.P. Stephan. The h-p version of the b,~undary element method for transmission problems with piecewlse analytic data. to appear.

1121 B. Guo and E.P. Stephan, The h-p version of the coupling of finite element and boundary element methods for transmission problems in polyhedral domains, to appear.

It31 B.Q. Guo, T. yon Petersdorff and E.P. Stephan, An h-p version of the boundary, element method for plane mixed boundary value problems, C.A. Brebbia and J.J. Connor, eds.. Proc. Conf. BEM-I !. 1989 Cambridge, USA. Advances in Boundary Elements, %/ol. 1 (1989) 95-11)3. 114] N. Heuer, hp-Versionen der Randelementmethods. Ph.D. Thesis, Hannover, 1992. 1151 N. Heuer and E.P. Stephan, Coupling of the finite element and boundary element method, The h-p version, Z A M M 71(6) ( 1991 ) T584-Tr586. 1161 N. Heuer+ M. Maischak and E.P. Stephan, The h-p version of the boundary element method for screen problems, Preprint Institut f/it Angewandte Mathematik Universit/it Hannover. 1994. [t71 J.L. Lions and E. Magenes, Non-homogeneous Boundary Value Problems and Applications, %/ol. 1 (Springer, BerlinHeidellrerg-New York. 1972). [I81 M. Maischak, hp-Methoden fiir Randintegralgleichungen bet 3D-Problemen+ Theorie und Impleraentierung, Ph.D. "lWaesis, Hannover, 1995. 1191 M. Maischak and E.P. Stephan, The h-p version of the boundary element method for first kind integral equaiions on polyhedral surfaces, Preprint Institut fiJr Angewandte Mathematik Universitat Hannover, 1994. I2Ol M. Maisehak and E.P. Stephan, The hp-Version of the boundary element method in R'. The basic approximation results, in press. I211 J.A. Morrison and J.A. Lewis. Singularity at the corner o ~ a flat plate. SIAM J. Appi. Math. 31 (1976) 233-250. 1221 T. yon Petersdhorff, Randwertprobleme der Elastizit~itstheorie fiir Polyeder~Singularitfiten und Approximation mit Randelemethoden, Ph.D. Thesis, Darmstadt, 1989. I231 T. yon Petersdorff and E.P. Stephan, Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron, Math. Nachrichten 149 (I990) 7I-104+ 1241 T. yon Petersdorff and E.P. Stephan, Regularity of mixed boundary value problems in R + and boundary, element methods on graded meshes, Math. Methods Appl. Sci. 12 (1990) 229-249.

208

E,P. 5;tephan / Comput. Methods Appi. Mech. Engrg. 133 (1996) 183-208

[25] F.V. Postell and E.i'. Stephan, On the h-, p- and h-p versions of the boundary integral element method--Numerical results, Comput. Methods Appi. Mech. Engrg. 83 (i99t!) 69-89. 126j Ch. Schwab and M. Suri. The optimal p-version approximation of singularities on polyhedra in the boundary element method. Uni. Maryland, Baltimore County, r~!:~.'n.',r~ .MD. USA, 1993. 1271 E.P. Stephan, Coupling of finite elements an0 boundary elemenen for ~olm: ~,onlinear interface problems. Comput. Methods Appl. Mech. Engrg. 101 (1992)61-72. [281 E.P. Stephan, Improved Galerkin methods for integral equations on polygons and on polyhedral surfaces. First loint Japan/US Syrup. on Boundary Element Methods. Tokyo (1988) 73-8(I. 129] E.P. Stephan, Boundary integral equations for screen problems in [~. J. Integral Eqns Operator Theory I(I (1987) 467-5(14. [3ill E.P. Stephan and M. Suri, The h-p version of the boundat3' element method on polygonal domains with quasiuniform meshes, Math. Modeling Numer. Anal. (RAIRO) 25(6) (1991) 783-8(17. [311 E.P. Stephan and M. Suri, On the convergence of the p-version for some boundary element Galerkin methods, Math. Comput. 52 (I989) 31-48. [321 E.P. Stephan and W.L. Wendland, Remarks to Galerkin and least squares methods with finite elements for general elliptic problems, Manuscripta Geodaetiea ! (1976} 93-123.