Computer methods in a p p l l o d
mochamics and engineering EI~SEVIER
Comput. Methods Appl. Mech. Engrg. 133 (1996) 183208
The hp boundary element method for solving 2 and 3dimensional 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 hpversion 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 hpversion on geometric meshes) converge exponentially fast towards tile exact solutions of the integral equations, "lhe implementation of the hpversion 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 hp 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 hp 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, threedimensional boundary value problems have been treated with the hp 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~ hp version which turns out to be the most suitable numerical approach to achieve accuracy and reliable solutions: the hp version with a geometric mesh refinement yields exponentially fast convergence, This is not restricted to twodimensional problems. As recently shown, the hp version of the B E M can also be applied to threedimensional 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 pdistribution 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 nonlinear) 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 hp version, and again exponentially fast convergence can be obtained [14, 15]. 00457825/96/$15.00 i~) 1996 Elsevier Science S.A. All rights reserved SSD! 00457825(95)00940X
E.I', Stephan / Cmnpttt. Methods" Appl. Mech. Engrg. 133 (I096) I83208
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=(lK')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,amG(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) i83208
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 twodimensional 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 IV)
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 hversion. Recently, the pversion [311 and the hp 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 Zinner 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) 183208
For the hversion 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 pversion 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 pversion 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 threedimensional 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 ~ cutofffunctions near corners (edges) and ~h~ consists o f a global remainder term belonging to s o m e H~(F j) plus higherorder.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) 1832(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 socalled 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 ' ~''atj2 +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 > 2d
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 hp 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 rigttthand 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 twodimensional problems see Stephan and Suri [301.
4. T h e
hp
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 hversion or the pversion are performed (cf. T h e o r e m s 3 and 4). Nevertheless the hp version with geometrically spaced mesh (near the vertices of .O) converges exponentially fast. T o describe the hp 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 eh,,.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)183208
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 hp 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')CII'
,,)=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) 161320~
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) 183208
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 Lshaped 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,.,, ,/iiII,.,,, 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 hp 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
hp v e r s i o n o f ( 3 ) hp v e r s i o n , o = 0 . 2 5
pversion, 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
hp v e r s i o , o f ( 2 0 ) hp v e r s i o n , rr = 0 . 1 5
pversion. 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) lg320bl
192
5. The
hp
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,,
~ HIt(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,+aol)!
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 splinefunctions 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) 11q3208
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,, ,.,,Il,, , . , , 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) lbt32Obl
194
yI,~, corresponds to the edgesingularity. 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 twodimensional spline, while a~, and ~e are splincs on a onedimensional cut of the twodimensional 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/2B,
(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, lA)=lA<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. lngrg. 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 +(yzx._)'+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]
i1
tI 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 lSJ2ob¢
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 hp 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
exyl 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 hp 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 hp 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 lm;t~.
4.
1.13 (H)96) lb132tt,'¢
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 xdirection on I: (xrefinement) 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 xycontrol (otherwise with xycon~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,
. Ixy] 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) 18320~
198
I t • ¢, = 0 . 2 5 °o=0.2
,m
Iiu,llo.
• ~= e.25 • = 0.2 0.171
10~.
9
10'~ i .q
103~
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 hvers~..oO 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 xycomro0 ~ t l v e t y . ~ v e r s ~ (w~h WoomoO
o
geometric mesh. sigma = o.17. mu=0.5
V
41.
~oUv~.
hpvorraon (with xyc~roO
Fig. 4, Dirichlet problem on tile Lshaped pIate.
885
17~5P_205
E.P Stephan i C'otnput. Methods Appl. Mech. Etagrg. 133 (1996) lg3201g
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 Lshaped 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 righthand 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. hve,sk~ ( w t U ~ xy,cx~m0 ,m:tat~ve~. hve,,',~,k~(wtth W ~ zm,zq~veq, hDverr,k ~ ( w ~ xyoomrm) 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) 183208
200
yielding with (18) the error estimates for the pure h and pversions
l
h
 ~
(38)
] C  C N I ~ c ~ p 2÷2e
with c 3 > 0 and arbitrary • > 0 . For the hp 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 pversions as predicted by (38). In Fig. 6 we present typical mesh sequences and degree distributions for the adaptive hpversion 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/2e
(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 hpvcrsion with direction control for the Dirichlet problem ( C = (}.659265).
E.P. Stephan
/
201
Comput. Methods Appl. Mech. Engrg. 133 (1996) 183208
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
pvernon, h=1
•
O e o r ~ mesa, m ~am~,.
,: o.17, muo_~ h~,,mk: (win xytxxW~
Fig. 7. The N e u m a n n problem on the _~luarc plate.
For the hpversion 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 pversions 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 Lblock given in Fig. 10. Again, Fig. 9 shows only algebraic convergence for the pure hand pversions and exponentially fast convergence for the hp version on the geometric mesh. Fig. 10 shows the mesh of the adaptive hversion with xycontrol. 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) !#320bi
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~, hre.on (webWtonal 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 pversion 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 pversion to solve the Dirichlet problem on the square plate. For the hp 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) 1832(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 []
hversk)n, p~0 pvsf~on, h=:l pVersion, rr0.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,,~
hpve.don(wah xyoon,~) skarnm =0.17, mue = 0.5
Fig. 9. Dirichlel problem on the Lblock.
6. Coupling of F E M and BEM&p version
H e r e , we consider the hp 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) lbt32Olq
21)4
i I 1 f
i
I
!.
_/
J
I
I
Hq
r77~
,_L _'L
H
j
t r
/
JI
•
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 Lblock.
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 hp version on a geometric mesh. For simplicity we consider only the case where F is an Lshaped polygon. Let .ft"., 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) IH320H
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
pVemk~n.~ i
9)
mesh, ~
= 0 . 1 7 , tlnU~=1.0
adapttv~. ~ v e ~ o n (v ~t~ xyoomra) 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  ~ ' ~
+ eb~"~)
•
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) IKI20B
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 CGalgorithm (pversion. 4 elements).
also be applied to interface problems with a nonlinear 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 nonlinear 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) 183208
207
2 Let , ¢ 2 = [  1 , 1 ] " and tt 1 =(2X21X2)' 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/11.
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) 12651288.
121 I. Babu~ka and B. Guo. The h. p and hp version of the finite clement method: basis theory and applications. Adv. Engrg. Software I5 (1992) 159174.
131 I. Babugka. B.Q. Guo and E.P. Stephan. On the exponential conucrgence of the hp version for boundary element Galerkin methods on polygons. Math. Methods Appl. Sci. 12 (1990) 413427.
141 !. Babugka. B.Q. Guo and E.P. Stephan. The hp version of the boundary element method with geometric mesh on polygonal domains. Comput. Methods Appl. Mech. Engrg. 80 (I99t}J 319325.
I51 M. Costabel+ Boundary integral operators on Lipschitz domains. Elementary results. SIAM J. Math. Anal. 19 (1988) 613626.
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) 12121226.
171 M. Dauge, Elliptic Boundary Value Problems on Corner Domains (SpringerVerlag, Lecture Notes in Mathematics, 1341, 1988).
181 M.R. Dorr. The approximation of solutions of elliptic botmdaryvalue problems via the pVersion of the finite element method, SlAM .I. Numer. Anal. 23 (1986) 5877.
191 V.J. Ervin. N. Heuer and E.P. Stephan, On the hp version of the boundary element method for Symm's integral equation on polygons, Comput. Methods Appl. Mech. Engrg. 110 (I993) 2538.
[to] P. Grisvard. Elliptic Problems in Nonsmooth Domains. Pitman Publishing Program (L985). !!I! B. Guo. N. bhzuer and E.P. Stephan. The hp version of the b,~undary element method for transmission problems with piecewlse analytic data. to appear.
1121 B. Guo and E.P. Stephan, The hp 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 hp version of the boundary, element method for plane mixed boundary value problems, C.A. Brebbia and J.J. Connor, eds.. Proc. Conf. BEMI !. 1989 Cambridge, USA. Advances in Boundary Elements, %/ol. 1 (1989) 9511)3. 114] N. Heuer, hpVersionen der Randelementmethods. Ph.D. Thesis, Hannover, 1992. 1151 N. Heuer and E.P. Stephan, Coupling of the finite element and boundary element method, The hp version, Z A M M 71(6) ( 1991 ) T584Tr586. 1161 N. Heuer+ M. Maischak and E.P. Stephan, The hp 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, Nonhomogeneous Boundary Value Problems and Applications, %/ol. 1 (Springer, BerlinHeidellrergNew York. 1972). [I81 M. Maischak, hpMethoden fiir Randintegralgleichungen bet 3DProblemen+ Theorie und Impleraentierung, Ph.D. "lWaesis, Hannover, 1995. 1191 M. Maischak and E.P. Stephan, The hp 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 hpVersion 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) 233250. 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) 7I104+ 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) 229249.
208
E,P. 5;tephan / Comput. Methods Appi. Mech. Engrg. 133 (1996) 183208
[25] F.V. Postell and E.i'. Stephan, On the h, p and hp versions of the boundary integral element methodNumerical results, Comput. Methods Appi. Mech. Engrg. 83 (i99t!) 6989. 126j Ch. Schwab and M. Suri. The optimal pversion 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)6172. [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) 738(I. 129] E.P. Stephan, Boundary integral equations for screen problems in [~. J. Integral Eqns Operator Theory I(I (1987) 4675(14. [3ill E.P. Stephan and M. Suri, The hp version of the boundat3' element method on polygonal domains with quasiuniform meshes, Math. Modeling Numer. Anal. (RAIRO) 25(6) (1991) 7838(17. [311 E.P. Stephan and M. Suri, On the convergence of the pversion for some boundary element Galerkin methods, Math. Comput. 52 (I989) 3148. [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} 93123.