- Email: [email protected]
Optimal extensions on tensor-product meshes
Applied Numerical Mathematics 54 (2005) 391–405 www.elsevier.com/locate/apnum
Optimal extensions on tensor-product meshes Sven Beuchler ∗,1 , Joachim Schöberl 2 Institute for Numerical Mathematics, Kepler-University Linz, Altenberger Strasse 69, A-4040 Linz, Austria
Abstract In this paper, a uniformly elliptic second order boundary value problem in 2D is discretized by the p-version of the finite element method. An inexact Dirichlet–Dirichlet domain decomposition preconditioner for the system of linear algebraic equations is investigated. The ingredients of such a preconditioner are a preconditioner for the Schur complement, a preconditioner for the subdomains and an extension operator operating from the edges of the elements into their interior. Using methods of multi-resolution analysis, we propose a new method in order to compute the extension efficiently. Numerical experiments show the optimal performance of the described extension. 2004 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Domain decomposition; Wavelets; Solution of discretized equations; Finite elements
1. Introduction We investigate the following boundary value problem. Let Ω ⊂ R2 be a domain which can be decom01 (Ω) = {u ∈ H 1 (Ω), u = 0 on Γ1 }, Γ1 ∩ Γ2 = ∅, Γ1 ∪ Γ2 ⊂ ∂Ω such posed into elements Rs . Find u ∈ H
* Corresponding author.
E-mail addresses: [email protected] (S. Beuchler), [email protected] (J. Schöberl). URLs: http://www.ricam.oeaw.ac.at/people/page/beuchler (S. Beuchler), http://www.hpfem.jku.at/joachim/index.html (J. Schöberl). 1 This work is supported by the Austrian Academy of Sciences. 2 This work is supported by the Sonderforschungsbereich F013 of the Austrian FWF with grant Y192 Start project hp-FEM. 0168-9274/$30.00 2004 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2004.09.010
392
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
that
a (u, v) :=
∇u · ∇v =
Ω
fv +
Ω
f1 v := f, v + f1 , v Γ2
(1.1)
Γ2
01 (Ω). The set Γ1 denotes the part of ∂Ω having a Dirichlet boundary condition u = 0, holds for all v ∈ H = f1 . Problem (1.1) and the set Γ2 denotes the part of ∂Ω having a Neumann boundary condition ∂u ∂n will be discretized by means of the p-version of the finite element method (FEM) using quadrilateral elements, see [25,34]. Let R2 = (−1, 1)2 be the reference element and Φs : R2 → Rs be the mapping to 01 (Ω), u|Rs = u(Φ ˜ s−1 (x, y)), u˜ ∈ Qp }, the element Rs . We define the finite element space M := {u ∈ H where Qp is the space of all polynomials of maximal degree p in each variable. Let (ζ1 , . . . , ζN ) be a basis for M. The Galerkin projection of (1.1) onto M leads to the linear system of algebraic finite element equations N N (1.2) Au = f , where A = a∆ (ζj , ζi ) i,j =1 , f p = f, ζi + f1 , ζi Γ2 i=1 . Using the vector u, an approximation of the exact solution u of (1.1) can be built by the usual finite element isomorphism. For the h-version of the FEM, the polynomial degree p of the shape functions on the elements is kept constant and the mesh-size h is decreased. This is in contrast to the p-version of the FEM in which the polynomial degree p is increased and the mesh-size h is kept constant. The advantage of the p-version in comparison to the h-version is that the approximate solution up converges faster to the exact solution u with respect to the number of unknowns N , if u is sufficiently smooth. Both ideas, mesh refinement and increasing the polynomial degree, can be combined. This is called the hp-version of the FEM. In both cases, h- and p-version of the FEM, the stiffness matrix is usually ill-conditioned. Therefore, efficient preconditioning techniques are required in order to solve (1.2). Preconditioners for (1.2) can be built by substructuring techniques (DD-methods). We refer to [9–12, 28–31], in the case of the h-version of the FEM and to [2,3,16,22–24,27] for the p-version of the FEM. From the algebraic point of view, the symmetric positive definite (spd) stiffness matrix A is split into I 0 S 0 A11 A12 I A12 A−1 22 (1.3) = A= 0 A22 A21 A22 A−1 0 I I 22 A21 with the Schur-complement S = A11 − A12 A−1 22 A21 . Then, the preconditioner CS 0 I −E T I 0 C= 0 I −E I 0 C22
(1.4)
is introduced. The following result has been proved in [18,28,29,32] and is the key in order to analyze the preconditioner C. Lemma 1.1. Let CS and C22 be spd preconditioners for S and A22 , i.e., cS (CS g, g) (Sg, g) CS (CS g, g) ∀g ∈ Rm , cI (C22 v, v) (A22 v, v) CI (C22 v, v) ∀v ∈ Rn . Moreover, let I I g cE2 (Sg, g) A g, E E
∀g ∈ Rm .
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
393
Then, the inequalities c(Cx, x) (Ax, x) C(Cx, x) ∀x ∈ Rn+m hold with c = 12 cE−2 min{cs , cI } and C = 2cE2 max{CS , CI }. In the p-version of the FEM, the splitting of the matrix A is naturally given by a splitting of the ansatz functions. Let nv , ne , and ni be the number of vertices not belonging to Γ 1 , number of edges not belonging to Γ 1 , and number of elements, respectively. One basis function corresponds to each vertex, p − 1 basis functions correspond to each edge, (p − 1)2 basis functions correspond to each element. Thus, the dimension of the ansatz space is N = nv + (p − 1)ne + (p − 1)2 ni . We define the functions ζ1 , . . . , ζnv as the usual piecewise bilinear hat functions. The functions ζnv +(j −1)(p−1)+1 , . . . , ζnv +j (p−1) correspond to the edge ej of the mesh, and vanish on all other edges, i.e., satisfy the condition ζnv +(j −1)(p−1)+k−1 |el = δj,l pk , where pk is a polynomial of degree p, k = 2, . . . , p. The support of an edge function is formed by those two elements, which have this edge ej in common. The remaining basis functions are interior bubble functions consisting of a support containing one element only. Now, the basis functions ζi are divided into three groups, • the vertex functions, • the edge bubble functions, • the interior bubbles. Corresponding to the division of the shape functions, the matrix A is split into the blocks Av Av,e Av,i A = Ae,v Ae Ae,i , Ai,v Ai,e Ai where the indices v, e, and i denote the blocks corresponding to the vertex, edge bubble, and interior bubble functions, respectively. In a first step, the simpler matrix Av 0 0 Ae Ae,i , C= 0 (1.5) 0 Ai,e Ai is investigated. It has been proved in [3] that the condition number κ(C −1 A) grows as 1 + log p. Therefore, the vertex unknowns can be determined separately. Preconditioners for Av are multi-grid methods, [20], or BPX-preconditioners, [13,35]. Computing the other unknowns, we factorize the remaining 2 by 2 block I 0 S 0 Ae Ae,i I Ae,i A−1 i (1.6) = 0 Ai Ai,e Ai 0 I I A−1 i Ai,e with the Schur-complement S = Ae − Ae,i A−1 i Ai,e . The matrix Ai is a block diagonal matrix, each block ARs corresponds to one element Rs , i.e., i . Ai = blockdiag[ARs ]ns=1
(1.7)
394
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
Therefore, in order to compute the interior unknowns, we have to solve a Dirichlet problem on each quadrilateral. The edge unknowns are computed via the Schur-complement S. Due to Lemma 1.1, an inexact DD-preconditioner for (1.1) includes a preconditioner for Ai , a preconditioner for the Schurcomplement S and an extension operator E operating from the edges of the quadrilateral into its interior in order to replace the matrix A−1 i Ai,e by a matrix −E.
which is the result of assembling the element In [24], Jensen/Korneev have introduced a matrix A stiffness matrices on the reference element R2 = [−1, 1]2 instead of the element stiffness matrix corre −1 A) = O(1). This result allows sponding to the element Rs . Then, they proved in Lemma 4.2., that κ(A us to restrict ourselves to the case of the reference element in order to derive preconditioners for Ai , S and to find the extension operator E. In [22–24], two preconditioners for the Schur-complement are proposed. The preconditioners use basis transformations from the integrated Legendre polynomials to the Chebyshev polynomials or a Lagrange basis. Another preconditioning techniques can be found in [1,17]. Moreover, an approach using multiresolution bases for the Schur-complement in 3D, [27], can be applied to the 2D case as well. For the interior solver, several preconditioners Ci are developed in [5,6,8,24,26]. The condition number of Ci−1 Ai is bounded independent of the polynomial degree, [7]. The operation Ci−1 w requires O(p2 ) operations. All preconditioners use interpretations of the matrix AR2 as an h-version FEM discretization matrix of a degenerated elliptic problem, published in [5]. In order to replace A−1 i Ai,e in (1.6), an efficient extension operator E ↔ E is required, cf. Lemma 1.1. Given a polynomial g ∈ Pp on one edge of R2 , find a polynomial u ∈ Qp with uH 1 (R2 ) cE gH 0.5 (∂ R2 )
∀g ∈ Pp
(1.8)
under the constraint u = g on ∂R2 , where the constant cE is independent of the polynomial degree p. Usually, the system Ai u = −Ai,e g
or
Ai u = −Ai,f g
(1.9)
can be treated by a preconditioned Chebyshev iteration with the preconditioner Ci for Ai , see, e.g., [27]. Since this operation is required in each preconditioning step, cf. (1.3), the system solve (1.9) is too expensive. In the case of the h-version of the FEM, i.e., the function g is a piecewise linear function, Nepomnyaschikh, [19,32], derived several techniques in order to develop such an extension operator without the solution of (1.9). For the p-version of the FEM, a pioneering work is done by [3]. The practical and efficient implementation of this extension operator in certain polynomial bases, i.e., integrated Legendre polynomials, has been an open question. This paper is dedicated to the development of fast solvers for the system (1.2) arising from the discretization by the p-version of the FEM. Efficient tools are domain decomposition methods. We mainly focus on one ingredient of the inexact domain decomposition preconditioner, the extension operator. On the one hand, we will prove the estimate (1.8). On the other hand, we will show that this extension can be computed in O(p2 ) operations. Our technique uses a basis for the given function g in (1.8) which is stable in L2 and H 1 . With methods of multi-resolution analysis, cf. [14,15,33] for the h-version of the FEM, and [8] for the p-version of the FEM, such bases can be determined. The paper is organized as follows. In Section 2, the model problem in order to derive the extension operator is introduced. Then, the algorithm for the general definition of extension operators is defined. In Section 3, some examples of extension operators for the h- and the p-version of the FEM in 2D are given. In Section 4, we show the numerical performance of the proposed extension operators.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
395
Throughout this paper, the notation g is used for vectors of the Euclidean space, where the corre sponding function g(x) = nj=1 gj φj (x) is denoted by g. We write g ↔ g. Moreover, let I = (−1, 1). The parameter p denotes the polynomial degree and Pp and Qp denote the space of all polynomials of degree p in one and two variables. Pp,00 denotes the subspace of Pp satisfying q(±1) = 0. 2. Construction of extension operators 2.1. Extensions from one face or edge We consider the following problem: −u = 0 in Ω = I × ω, u(1, y) = 0, u(−1, y) = g(y) u(x, y) = 0 on I × ∂ω,
(2.1) (2.2) (2.3)
on ω,
where I = (−1, 1). Problem (2.1)–(2.3) is the typical model problem for the harmonic extension of a function g on ω into the interior of the domain Ω in the tensor product case. The weak formulation of (2.1)–(2.3) is: find u ∈ Hg (Ω) = {u = uˆ + g, uˆ ∈ H01 (Ω)} such that (2.4) a(u, v) = ∇u · ∇v = 0, ∀v ∈ H01 (Ω) Ω
holds. Problem (2.1)–(2.3) is solved approximately by a tensor product version of the finite element method (FEM). Let Ξ = {ξj (x)}m j =1 be a basis for the discretization of the interval I = (−1, 1) with ξj (±1) = 0 for j = 2, . . . , m, ξ1 (1) = 0 and ξ1 (−1) = 1. Moreover, let X = {χj (y)}nj=1 be a basis for the discretization of ω satisfying χj (∂ω) = 0. On Ω, we use the approximation space m,n V = span φij (x, y) i,j =1 , where φij (x, y) = ξi (x)χj (y). In order to satisfy the boundary condition (2.2), the functions φ1j , 1 j n, are defined, i.e., we assume that g(y) = nj=1 gj φ1 (−1, y) with the given coefficients gj . of problem (2.4) onto the space V leads to the following problem: find u = nThe Galerkinprojection m n g φ + u j 1j j =1 i=2 j =1 ij φij such that a(u, φij ) = 0,
∀i = 2, . . . , m, j = 1, . . . , n.
Introducing the matrices (m,n) K22 = a(φij , φlk ) (i,j ),(k,l)=(2,1) ,
n K11 = a(φ1j , φ1k ) k;j =1 ,
(2.5) m,n;n K21 = a(φij , φ1k ) i=2,j =1;k=1 ,
problem (2.5) is equivalent to solve the system of linear algebraic equations K22 u = −K21 g
(2.6)
with u = [u21 , . . . , umn ]T and g = [g1 , . . . , gn ]T . We are not interested in solving the problem (2.5), (2.6) exactly, we only try to find a vector u satisfying the inequality T K11 K21 (2.7) x, x = (Kx, x) cE2 (Sg, g) ∀g K21 K22
396
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
−1 T under the constraint x = gu . The matrix S = K11 − K21 K22 K21 denotes the Schur complement. The constant cE should not depend on the discretization parameter. In the case of the h-Version of the FEM, several techniques in order to construct a vector u ↔ u satisfying (2.7) are given in [32], where the system (2.6) has not to be solved. By the usual FEM-isomorphism, one has g 2 (Kx, x) =| u |H 1 (Ω) with x = ↔u (2.8) u between the vector x = gu and the corresponding function u. Due to (2.2) and the Poincare–Friedrichs −2 inequality, the H 1 -seminorm is equivalent to the H 1 -norm. Thus, one obtains CΩ u2H 1 (Ω) (Kx, x) 2 uH 1 (Ω) for all u, where the constant CΩ depends on the domain Ω only. Concerning a relation between the function g and the vector g, the following assumption is required.
Assumption 2.1. Let us assume that for every g ∈ span X there exists a function u such that uH 1 (Ω) cH gH 0.5 (∂Ω) ,
(2.9)
under the constraints u(−1, y) = g(y), u(1, y) = 0, u(x, ∂ω) = 0, where the constant c does not depend on the discretization parameter. Then, the following result can be proved. Lemma 2.2. Let us assume that Assumption 2.1 is satisfied and that ∂Ω is Lipschitz continuous. Let g = nj=1 gj φj and g = [gj ]nj=1 . Then, we have the norm equivalence relation 2 cT2 g2H 0.5 (∂Ω) (Sg, g) cH g2H 0.5 (∂Ω) ,
(2.10)
where the constant cH is the constant in (2.9) and cT is the constant from the trace theorem. Proof. The proof is similar to the arguments in order to derive formula (3.12) in [24]. By the properties of the Schur-complement, we have (Sg, g) = minu (Ku, u), where x = gu . Let u∗ ↔ x ∗ = ug∗ be the optimal solution of this minimization problem. Thus, by the trace theorem, cT2 g2H 0.5 (∂Ω) u∗ 2H 1 (Ω) = Kx ∗ , x ∗ = (Sg, g). Moreover, let u0 ↔ x0 = ug0 be the extension of Assumption 2.1. Then, 2 g2H 0.5 (∂Ω) , (Sg, g) (Kx0 , x0 ) = u0 2H 1 (Ω) cH
which proves the lemma.
2
This lemma is important in order to derive a preconditioner for the Schur-complement embedded in the domain decomposition preconditioner (1.4). Using this norm equivalence, most preconditioners for 0.5 S can be constructed by finding a stable basis in H (∂Ω). Moreover, we are able to define (Sg, g) as an equivalent norm on H 0.5 (∂Ω). The next lemma gives a relation for the extension and the corresponding norms.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
Lemma 2.3. Let us assume that the extension satisfies (2.7). Then, x =
g u
397
↔ u satisfies
uH 1 (Ω) CΩ cE cH gH 0.5 (∂Ω) .
(2.11)
Proof. Use relations (2.10), (2.7), (2.8) and the Poincare–Friedrichs inequality.
2
In order to define the extension operator, we use a technique which is similar to the multi-level decomposition technique derived in [19]. In the space Y = span X, a so-called stable basis is used. Assumption 2.4. Let us assume that there exists a basis Ψ = {ψj }nj=1 in the space Y which is stable in L2 (ω) and H 1 (ω), i.e., for any function g = nj=1 g˜ j ψj (y) the norm equivalence relations c1,1
n
dj(1) g˜ j2 g2L2 (ω) c2,1
j =1
c1,2
n
n
dj(1) g˜ j2 ,
(2.12)
j =1
dj(2) g˜ j2
g2H 1 (ω)
c2,2
j =1
n
dj(2) g˜ j2
(2.13)
j =1
hold with some numbers dj(i) > 0, j = 1, . . . , n, i = 1, 2. The constants ci,j , i, j = 1, 2, do not depend on the discretization parameter n. The computation of the nearly discrete harmonic extension of the function g(y) = nj=1 gj χj (y) consists of three steps: Algorithm 2.5. 1. Transform the function g into the basis Ψ = {ψj }nj=1 , i.e., g(y) =
n
g˜ j ψj (y).
j =1
2. For j = 1, . . . , n, let u∗j be the solution of the problem |uj|H1 (Ω) → min,
where uj (x, y) = ψj (y) ξ1 (x) +
fi
m
ξi (x)fi .
(2.14)
i=2
Put u(x, y) =
n
g˜ j u∗j (x, y).
j =1
3. Transform u(x, y) =
n
˜ j uj (x, y) j =1 g
(2.15) into the basis {φij (x, y)}i,j .
We introduce the following real numbers, vectors and matrices: let m m a 1 = ξ1 , ξj I j =2 , A1,0 = ξi , ξj I i,j =2 , βj,1 = ψj , ψj ω , m m βj,2 = Dψj , Dψj ω , a 2 = Dξ1 , Dξj I j =2 , A2,0 = Dξi , Dξj I i,j =2 ,
(2.16)
398
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
where ·, · Ω denotes the L2 (Ω)-scalar product and Du is the gradient of u. In order to compute the extension efficiently, the following lemma and remark are useful. Lemma 2.6. The computation of uj in (2.14) requires the system solve (βj,1 A1,0 + βj,2 A2,0 )f = −βj,1 a 1 − βj,2 a 2 ,
(2.17)
where βj,1 , βj,2 > 0. Proof. Inserting (2.14) into uj leads to the optimization problem 2 m ξi (x)fi → min . ψj (y)ξ1 (x) + fi 1 i=2
H (Ω)
Using (2.16) and tensor product arguments, this is equivalent to 2 2 m m βj,1 ξ1 (x) + ξi (x)fi + βj,2 ξ1 (x) + ξi (x)fi 1 i=2
i=2
H (I )
L2 (I )
→ min . fi
The solution of this quadratic optimization problem with respect to f = [fk ]m k=2 leads to the system βj,1 Dξi , Dξk I + βj,2 ξi , ξk I fk = −βj,1 Dξi , Dξ1 I − βj,2 ξi , ξ1 I , i = 2, . . . , m. Inserting now (2.16) proves the lemma.
2
Remark 2.7. In (1), one basis change from the basis X to the basis Ψ has to be done, in (3), n basis changes from the Ψ to the basis X have to be done. (2) means from the algebraic point of view n system solves with positive linear combinations of the matrices A1,0 and A2,0 , cf. Lemma 2.6. The matrices A1,0 and A2,0 are the 1D mass and 1D stiffness matrix with respect to the basis Ψ . Remark 2.8. Moreover, we investigate the following modification of Algorithm 2.5. In this second approach, instead of the systems (2.17) in order to compute uj (2.14), the systems (1) dj A1,0 + dj(2) A2,0 f = −dj(1) a 1 − dj(2) a 2 , with the constants dj(i) in (2.12) and (2.13) have to be solved. This is a modification of Algorithm 2.5. 2.2. The general case In the last subsection, only an extension of a function living on one edge (in 2D) or face (in 3D) of the domain Ω has been investigated. This subsection deals with the general case of an extension from ∂Ω into Ω. Usually, the boundary of Ω is split into several faces (edges in 2D) ωi , i = 1, . . . , r, r i=1 ωi = ∂Ω, ωi ∩ ωj = ∅ for i = j . On each ωi , a function gi is given which satisfies the boundary condition gi (∂ωi ) = 0. Then, we consider the problem: find u ∈ H 1 (Ω) such that u2H 1 (Ω) cg2H 0.5 (∂Ω) ,
u|ωi = gi = g|ωi , i = 1, . . . , r.
(2.18)
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
399
This problem can be solved by the method described in the previous subsection (Algorithm 2.5). Let ui be the extension of ui |ωj = gi δij
on ωj and ui 2H 1 (Ω) ci gi 2H 0.5 (ωi ) .
(2.19)
Then, we define u = ri=1 ui as the extension of the function g into Ω. Obviously, u|ωj = ri=1 ui |ωj = r i=1 gi δij = gj , i.e., condition (2.18) is fulfilled. Moreover, by the Cauchy–Schwarz inequality and (2.19), we obtain u2H 1 (Ω)
r 2 = uj 1 j =1
r
H (Ω)
r
r
uj 2H 1 (Ω) r
j =1
r j =1
cj gj 2H 0.5 (ωj ) r max {cj } j =1,...,r
r j =1
gj 2H 0.5 (ωj )
max {cj }g2H 0.5 (∂Ω) . j =1,...,r
In the last estimate, a relation between the norms · 2H 0.5 (∂Ω) and rj =1 · 2H 0.5 (ωj ) is required. In the second norm, the coupling between two neighbouring faces or edges ωi and ωj has no influence in the norm, whereas an additional connectivity term is involved in the definition of the first norm. So, with an optimal extension for problem (2.19), cf. Algorithm 2.5, an optimal extension for (2.18) can be computed.
3. Applications In this section, several examples for extension operators of the type proposed in Algorithm 2.5 are given. So, Assumption 2.4 has to be verified, i.e., a polynomial basis Ψ has to be derived which is stable in H 1 (I ) and L2 (I ). Using methods of multi-resolution analysis, [8,14,15,33], such a basis Ψ can be constructed in the case of h- and p-version of the FEM. We start with the p-version of the FEM. 3.1. Extension operators for the p-version of the FEM In this subsection, we describe a method deriving an extension operator E, i.e., u = Eg in matrix representation, in the case of the p-version of the FEM. Usually, such an extension operator is embedded in a domain decomposition preconditioner of Dirichlet–Dirichlet-type, cf. Section 1, or [5]. The operator u = Eg will require O(p2 ) arithmetical operations in 2D, i.e., it is arithmetically optimal. Moreover, numerical experiments indicate the estimate I I g, g c(Sg, g) ∀g ∈ Rm , (Kx, x) = K E E where c is independent of the polynomial degree p. A system solve K22 u = −K21 g is not required. In Section 3.1.1 a method in order to construct stable polynomial bases in the spaces H 1 (I ) and L2 (I ) is shown. In Section 3.1.2, we derive the extension operators for the 2D case.
400
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
3.1.1. A stable polynomial basis of Pp,00 in H 1 (I ) and L2 (I ) In this subsubsection, we propose a basis in Pp,00 , which is stable in the spaces L2 (I ) and H 1 (I ). For i ∈ N, let Li (x) = ( dxd )i (x 2 − 1)i be the ith Legendre polynomial and let √ x (2i − 3)(2i − 1)(2i + 1) i (x) = Li−1 (s) ds (3.1) L 2 −1
i }pi=2 is a basis in i (±1) = 0. Thus, {L be the ith integrated Legendre polynomial (i 2). Note that L the space Pp,00 . Moreover, we introduce the 1D mass and 1D stiffness matrix with respect to the basis i }pi=2 , i.e., {L j , L i I p j , D L i I p . and H1,I = D L (3.2) L2,I = L i,j =2 i,j =2 In [7], we have proved that 1 0 −γ2 0 1 0 −γ2 0 1 L2,I = .. . 0
...
0
0 −γ3 0 .. .
−γ4
−γp−2
0
2 0 −1 ∼ .. . 1
0
0 2 0
−1 0 0 −1 2 0 −1 .. .
...
−1
Moreover, an easy calculation shows, [24], (2j − 3)(2j + 1) p . H1,I = diag 2 j =2
0
.
(3.3)
2
(3.4)
Here and in the following, the relation A ∼ B means that c1 (Av, v) (Bv, v) c2 (Av, v) ∀v, where the constants c1 and c2 do not depend on the polynomial degree p. Using a permutation P of rows and columns, one obtains T L1 T H1 P and H1,I = P P, (3.5) L2,I ∼ P L1 H1 where L1 is a tridiagonal matrix with 2’s on the main diagonal and −1’s on the first subdiagonal. Obi }pi=2 is stable in L2 (I ). In i }pi=2 is stable in H1 (I ). However, we are not able to prove that {L viously, {L order to find a stable basis, we use interpretations of the matrices L1 and H1 as discretization matrices on the unit interval (0, 1) using piecewise linear ansatz functions. ], i = 0, . . . , n − 1, be an equidistant Let k ∈ N, n = 2k and p = 2n − 1. Moreover, let τi = [ ni , i+1 n be the basis of the usual hat functions mesh on the interval [0, 1] and Θ = {θi }n−1 i=1 nx − (i − 1) on τi−1 , θi (x) = (i + 1) − nx on τi , i = 1, . . . , n − 1. (3.6) 0 else, We introduce the matrices 1 n−1 1 n−1 θi (x)θj (x)x 2 dx and T Θ = θi (x)θj (x) dx . MΘ = −1
i,j =1
In [4], we have proved the following result.
−1
i,j =1
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
401
Lemma 3.1. The spectral equivalence relations H1 ∼ n3 M Θ and L1 = n1 T Θ are valid. Concerning the weighted mass matrix M Θ and the unweighted stiffness matrix T Θ , it is known from the wavelet theory, [8], that it can be constructed a multi-level basis in which M Θ and T Θ are spectrally equivalent to diagonal matrices. Let J = (j, l), where 0 l k and j = 1, 3, 5, . . . , 2l − 1, card(J ) = n − 1. The diagonal matrices (3.7) D1,0 = diag 22l J =(j,l) and D2,0 = diag j 2 2−2l J =(j,l) are introduced. Theorem 3.2. There are multi-level bases Ψ˜ = {ψ˜ J }J with Ψ˜ = ΘQ such that the matrices T Ψ = QT T Θ Q and M Ψ = QT M Θ Q satisfy the spectral equivalence relations T Ψ ∼ D1,0 and M Ψ ∼ D2,0 . The transformation v = QT w can be performed in O(n) operations. Examples for wavelets bases are given in [8]. Now, let Q 0 , W = PT 0 Q
(3.8)
where P denotes the permutation matrix in (3.5) and Q denotes the fast wavelet transform in Theorem 3.2. Theorem 3.3. The matrices L2,I and H1,I satisfy the relations 1 −T D1,0 0 −1 3 −T D2,0 W and H1,I ∼ n W L2,I ∼ W 0 D1,0 0 n
0 W −1 . D2,0
(3.9)
Proof. Using (3.5), Lemma 3.1 and Theorem 3.2, we have PT TΘ 1 T Q−T T Ψ Q−1 0 0 0 T L1 P= P= P P L2,I ∼ P 0 L1 0 TΘ 0 Q−T T Ψ Q−1 n n −1 P T Q−T 0 0 0 D1,0 Q ∼ P. 0 Q−T 0 D1,0 0 Q−1 n Since P is a orthogonal matrix, the first assertion follows from (3.8). The second assertion can be proved in the same way. 2 Via the matrix W (3.8), the basis 2 (x), . . . , L p (x) W Ψ = ψ1 (x), . . . , ψp−1 (x) = L
(3.10)
is introduced. Corollary 3.4. The basis Ψ is stable in L2 (I ) and H 1 (I ). Moreover, the operation v = W w can be performed in O(p) arithmetical operations. p j (x) ∈ Pp,00 from the So, we have found an algorithm which transforms a function g(x) = j =2 gj L basis of the integrated Legendre polynomials into the basis Ψ ⊂ Pp,00 which is stable in L2 (I ) and H 1 (I ) and requires O(p) operations.
402
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
3.1.2. The 2D case We consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the p-version of the FEM with ω = (−1, 1), Ω = (−1, 1)2 . We use the integrated Legendre polynomials (3.1) for the discretization in x- and y-direction, i.e., j (x), ξj (x) = L j +1 (y), χj (y) = L
j = 2, . . . , p, ξ1 (x) =
1−x 2
and
j = 1, . . . , p − 1.
Using the notation (2.16), an easy calculation shows, cf. [24], T 5 7 , , 0, . . . , 0 , A1,0 = L2,I , a1 = − 12 60 a 2 = [0, . . . , 0]T ,
A2,0 = H1,I .
(3.11)
We use Algorithm 2.5 in order to compute a nearly optimal discrete harmonic extension. Since Corollary 3.4, the basis Ψ (3.10) is stable in L2 (I ) and H1 (I ). Thus, Assumptions 2.4 are satisfied. Lemma 3.5. Let x = gu ↔ u be defined via Algorithm 2.5 or the generalization of Remark 2.8. Then, the computation of u requires O(p2 ) operations. Proof. Since the multiplication w = Qv requires O(p) operations, the cost of p − 1 operations w = W v is of order p2 . Due to (3.3) and (3.4), the matrices dj(1) L2,I + dj(2) H1,I are block diagonal matrices of two tridiagonal matrices after reordering the unknowns, dj(i) ∈ R, dj(i) > 0, i = 1, 2, j = 1, . . . , p − 1. Using Cholesky decomposition, the cost of (dj(1) L2,I + dj(2) H1,I )−1 v is O(p). Solving p − 1 systems, cf. Lemma 2.6, the total cost is of step 2 in Algorithm 2.5 of order p2 . For the operation W −1 g, an inverse wavelet transform is required, where W ∈ Rp−1×p−1 . Since W is a product of [log2 p] matrices with a bandwidth O(1), the cost for W −1 g is of order [log2 p]p using Gaussian elimination. 2 Remark 3.6. This technique in order to construct the extension can be applied to the 3D-case. Then, we consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the p-version of the FEM with ω = (−1, 1)2 , Ω = (−1, 1)3 . With tensor product arguments, the basis Ψ × Ψ is stable in the tensor product spaces L2 (I × I ) = L2 (I ) × L2 (I ) and H 1 (I × I ) = H 1 (I ) × L2 (I ) + L2 (I ) × H 1 (I ). 3.2. h-version of the FEM in 2D and 3D We consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the h-version of the FEM ], i = with ω = (−1, 1), Ω = (−1, 1)2 . Let k ∈ N, n = 2k+1 − 1, s = 2k . Moreover, let τi = [ si , i+1 s s−1 −s, . . . , s − 1, be an equidistant mesh on the interval I and Θ = {θi }i=−s+1 be the basis of the usual hat functions (3.6). Then, we choose χj (y) = θj +s (y) for j = −s + 1, . . . , s − 1. Moreover, let τ˜i , i = 1, . . . , m be an arbitrary 1D finite element mesh of (−1, 1). Then, Ξ is the basis of the usual hat functions according to the mesh τ˜i . It is known from the theory of multi-resolution bases, [14,15,33], that several wavelet bases on the interval are stable in H 1 and L2 . Hence, Assumption 2.4 is satisfied. Moreover, it has been shown
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
403
in [19] that there exists a piecewise linear harmonic extension u of the piecewise linear function g = s−1 j =−s+1 gj θj into the interior of the domain Ω. Thus, we have verified the second Assumption 2.1, too. So, the optimal extension can be computed via Algorithm 2.5. For the computation of u, n systems of linear algebraic equations with a linear combination Mj , j = 1, . . . , n, of the one-dimensional mass matrix and one-dimensional stiffness matrix with respect to the basis Ξ have to be solved, cf. Remark 2.7. Since Ξ is the basis of the usual hat-functions, the system matrices Mj are tridiagonal matrices. Thus, one system solve requires O(m) operations and the total cost is O(mn). The approach can be extended to the 3D-case, where ω is a surface in R3 , [21]. 4. Numerical examples In this section, we compute the smallest possible constant cE in (2.7) for the proposed extension operator in the case of the p-version of the finite element method in two dimensions. In all experiments, the wavelet basis described in [7] is used with two types of scalings, an explicit one given by Theorem 3.3 and a implicit one given by the numbers dj(i) = βj,i , i = 1, 2, j = 1, . . . , n. In the first experiment, we use the numbers dj(i) , i = 1, 2, j = 1, . . . , n, defined via (3.7) and (3.9) instead of the numbers βj,i (2.16). This is a slight modification of Algorithm 2.5, cf. Remark 2.8. In the second example, we take dj(i) = βj,i , i = 1, 2, j = 1, . . . , n. Using (3.2) and (3.8), these entries can explicitly be computed in O(p2 ) operations. Table 1 displays the quality of the extension of Algorithm 2.5, i.e., Remark 2.8, in the first example. Moreover, the best constants c1,1 , c1,2 , c2,1 and c2,2 in (2.12) and (2.13) are computed. For all p, the constant cE2 lies between 1.4 and 1.9 and is bounded by a constant independent of p. Table 2 displays the constants cE , ci,j , i, j = 1, 2, for the second example. Now, the constants cE and c1,1 , c2,1 and c2,2 are closer to one. The only difference between the two proposed versions of extension is the type of diagonal scaling of the basis Ψ in L2 (ω) and H 1 (ω), i.e., the choice of the parameters dj(i) . So, the scaling factors dj(i) > 0 in Table 1 Smallest possible constant of (2.7) for the p-version extension operator p 3 7 15 31 2 in (2.7) cE 1.45 1.89 1.85 1.79 c1,1 0.22 0.20 0.18 0.17 0.94 0.37 0.25 0.21 c1,2 2.23 2.35 2.39 2.41 c2,1 3.94 3.94 3.42 3.23 c2,2
63 1.76 0.17 0.19 2.42 3.15
Table 2 Smallest possible constant of (2.7) for the p-version extension operator with implicitly given βj,i p 3 7 15 31 63 2 in (2.7) cE 1.00 1.13 1.26 1.34 1.37 c1,1 0.36 0.36 0.36 0.36 0.36 1.00 0.33 0.25 0.21 0.19 c1,2 1.63 1.68 1.77 1.83 1.88 c2,1 1.00 1.96 2.16 2.27 2.35 c2,2
127 1.75 0.16 0.18 2.43 3.11
127 1.39 0.36 0.18 1.92 2.40
404
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
order to solve dj(1) L2,I + dj(i) H1,I , cf. the proof of Lemma 3.5, are different between the first and second experiment. Since these systems are solved by a direct method and dj(i) > 0, the arithmetical cost for the computation of the extension itself is the same. For the second experiment, an additional initialization has to be done, i.e., the entries dj(i) = βj,i have to computed which is not required for the first experiment. Usually, the extension is embedded in a domain decomposition preconditioner C for A (1.2). Since the quality of the extension is not so good for the first experiment, the number of the pcg iterations may be higher than using the extension of the second experiment. The additional cost for the second experiment, i.e., the computation of the entries dj(i) = βj,i , can be done before starting the pcg method. Moreover, the arithmetical cost in order to compute βj,i is comparable with the cost for the computation of the extension one times. Summarizing, the method which uses dj(i) = βj,i , i.e., the second experiment should be preferred.
References [1] M. Ainsworth, A preconditioner based on domain decomposition for h–p finite element approximation on quasi-uniform meshes, SIAM J. Numer. Anal. 33 (4) (1996) 1358–1376. [2] M. Ainsworth, B. Gao, An additive Schwarz preconditioner for p-version boundary element approximation of the hypersingular operator in three dimensions, Numer. Math. 85 (3) (2000) 343–366. [3] I. Babuška, A. Craig, J. Mandel, J. Pitkäranta, Efficient preconditioning for the p-version finite element method in two dimensions, SIAM J. Numer. Anal. 28 (3) (1991) 624–661. [4] S. Beuchler, Preconditioning for the p-version of the FEM by bilinear elements, Technical Report SFB393 01-17, Technische Universität Chemnitz, May 2001. [5] S. Beuchler, Multi-grid solver for the inner problem in domain decomposition methods for p-FEM, SIAM J. Numer. Anal. 40 (3) (2002) 928–944. [6] S. Beuchler, AMLI preconditioner for the p-version of the FEM, Numer. Linear Algebra Appl. 10 (8) (2003) 721–732. [7] S. Beuchler, Optimal preconditioners for the p-version of the FEM, Technical Report SFB393 03-03, Technische Universität Chemnitz, March 2003. [8] S. Beuchler, R. Schneider, C. Schwab, Multiresolution weighted norm equivalences and applications, Numer. Math. 98 (1) (2004) 67–97. [9] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, I, Math. Comp. 47 (175) (1986) 103–134. [10] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, II, Math. Comp. 49 (179) (1987) 1–16. [11] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, III, Math. Comp. 51 (184) (1988) 415–430. [12] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, IV, Math. Comp. 53 (187) (1989) 1–24. [13] J. Bramble, J. Pasciak, J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (191) (1991) 1–22. [14] A. Cohen, I. Daubechies, P. Vial, Biorthogonal spline wavelets on the interval—stability and moment conditions, Appl. Comput. Harm. Anal. 6 (2) (1998) 132–196. [15] W. Dahmen, Wavelet and multiscale methods for operator equations, Acta Numerica 6 (1997) 55–228. [16] B. Guo, W. Cao, A preconditioner for the hp-version of the finite element method in two dimensions, Numer. Math. 75 (1996) 59–77. [17] B. Guo, W. Cao, An iterative and parallel solver based on domain decomposition for the hp-version of the finite element method, J. Comput. Appl. Math. 83 (1997) 71–85. [18] G. Haase, U. Langer, A. Meyer, The approximate Dirichlet domain decomposition method. Part I: An algebraic approach, Computing 47 (1991) 137–151.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
405
[19] G. Haase, S.V. Nepomnyaschikh, Explicit extension operators on hierarchical grids, East–West J. Numer. Math. 5 (4) (1998) 231–284. [20] W. Hackbusch, Multigrid Methods and Applications, Springer, Heidelberg, 1985. [21] H. Harbrecht, R. Schneider, Biorthogonal wavelet bases for the boundary element method, Math. Nachr. 269–270 (1) (2004) 167–188. [22] S.A. Ivanov, V.G. Korneev, On the preconditioning in the domain decomposition technique for the p-version finite element method, Part I, Technical Report SPC 95-35, Technische Universität Chemnitz-Zwickau, December 1995. [23] S.A. Ivanov, V.G. Korneev, On the preconditioning in the domain decomposition technique for the p-version finite element method, Part II, Technical Report SPC 95-36, Technische Universität Chemnitz-Zwickau, December 1995. [24] S. Jensen, V.G. Korneev, On domain decomposition preconditioning in the hierarchical p-version of the finite element method, Comput. Methods Appl. Mech. Engrg. 150 (1–4) (1997) 215–238. [25] G.M. Karniadakis, S.J. Sherwin, Spectral/HP Element Methods for CFD, Oxford University Press, Oxford, 1999. [26] V.G. Korneev, Poqti optimalny metod rexeni zadaq Dirihle na podoblasth dekompozicii ierarhiqesko hp-versii, Differencialnye Uravneni 37 (7) (2001) 1–15. [27] V. Korneev, U. Langer, L. Xanthis, On fast domain decomposition methods solving procedures for hp-discretizations of 3d elliptic problems, Technical Report 03-07, Universität Linz, March 2003. [28] P.L. Lions, On the Schwarz alternating method I, in: R. Glowinski, G.H. Golub, G.A. Meurant, J. Périaux (Eds.), Proc. 1st Int. Symp. on Domain Decomposition Methods, SIAM, Philadelphia, PA, 1988, pp. 1–42. [29] P.L. Lions, On the Schwarz alternating method, II, in: Proc. 1st Internat. Symp. on Domain Decomposition Methods, Los Angeles, 1989, pp. 47–70. [30] A.M. Matsokin, S.V. Nepomnyaschikh, The Schwarz alternation method in a subspace, Izv. Vyssh. Uchebn. Zaved. Mat. 29 (10) (1985) 61–66. [31] A.M. Matsokin, S.V. Nepomnyaschikh, Norms in the space of traces of mesh functions, Sov. J. Numer. Anal. Math. Modelling 3 (3) (1988) 199–216. [32] S.V. Nepomnyaschikh, Method of splitting into subspaces for solving elliptic boundary value problems in Complex-form domains, Sov. J. Numer. Anal. Math. Modelling 6 (1991) 151–168. [33] R. Schneider, Multiskalen- und Wavelet Matrixkompression, Teubner, Leipzig, 1998. [34] C. Schwab, p- and hp-Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics, Clarendon, Oxford, 1998. [35] X. Zhang, Multilevel Schwarz methods, Numer. Math. 63 (1992) 521–539.
Optimal extensions on tensor-product meshes Sven Beuchler ∗,1 , Joachim Schöberl 2 Institute for Numerical Mathematics, Kepler-University Linz, Altenberger Strasse 69, A-4040 Linz, Austria
Abstract In this paper, a uniformly elliptic second order boundary value problem in 2D is discretized by the p-version of the finite element method. An inexact Dirichlet–Dirichlet domain decomposition preconditioner for the system of linear algebraic equations is investigated. The ingredients of such a preconditioner are a preconditioner for the Schur complement, a preconditioner for the subdomains and an extension operator operating from the edges of the elements into their interior. Using methods of multi-resolution analysis, we propose a new method in order to compute the extension efficiently. Numerical experiments show the optimal performance of the described extension. 2004 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Domain decomposition; Wavelets; Solution of discretized equations; Finite elements
1. Introduction We investigate the following boundary value problem. Let Ω ⊂ R2 be a domain which can be decom01 (Ω) = {u ∈ H 1 (Ω), u = 0 on Γ1 }, Γ1 ∩ Γ2 = ∅, Γ1 ∪ Γ2 ⊂ ∂Ω such posed into elements Rs . Find u ∈ H
* Corresponding author.
E-mail addresses: [email protected] (S. Beuchler), [email protected] (J. Schöberl). URLs: http://www.ricam.oeaw.ac.at/people/page/beuchler (S. Beuchler), http://www.hpfem.jku.at/joachim/index.html (J. Schöberl). 1 This work is supported by the Austrian Academy of Sciences. 2 This work is supported by the Sonderforschungsbereich F013 of the Austrian FWF with grant Y192 Start project hp-FEM. 0168-9274/$30.00 2004 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2004.09.010
392
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
that
a (u, v) :=
∇u · ∇v =
Ω
fv +
Ω
f1 v := f, v + f1 , v Γ2
(1.1)
Γ2
01 (Ω). The set Γ1 denotes the part of ∂Ω having a Dirichlet boundary condition u = 0, holds for all v ∈ H = f1 . Problem (1.1) and the set Γ2 denotes the part of ∂Ω having a Neumann boundary condition ∂u ∂n will be discretized by means of the p-version of the finite element method (FEM) using quadrilateral elements, see [25,34]. Let R2 = (−1, 1)2 be the reference element and Φs : R2 → Rs be the mapping to 01 (Ω), u|Rs = u(Φ ˜ s−1 (x, y)), u˜ ∈ Qp }, the element Rs . We define the finite element space M := {u ∈ H where Qp is the space of all polynomials of maximal degree p in each variable. Let (ζ1 , . . . , ζN ) be a basis for M. The Galerkin projection of (1.1) onto M leads to the linear system of algebraic finite element equations N N (1.2) Au = f , where A = a∆ (ζj , ζi ) i,j =1 , f p = f, ζi + f1 , ζi Γ2 i=1 . Using the vector u, an approximation of the exact solution u of (1.1) can be built by the usual finite element isomorphism. For the h-version of the FEM, the polynomial degree p of the shape functions on the elements is kept constant and the mesh-size h is decreased. This is in contrast to the p-version of the FEM in which the polynomial degree p is increased and the mesh-size h is kept constant. The advantage of the p-version in comparison to the h-version is that the approximate solution up converges faster to the exact solution u with respect to the number of unknowns N , if u is sufficiently smooth. Both ideas, mesh refinement and increasing the polynomial degree, can be combined. This is called the hp-version of the FEM. In both cases, h- and p-version of the FEM, the stiffness matrix is usually ill-conditioned. Therefore, efficient preconditioning techniques are required in order to solve (1.2). Preconditioners for (1.2) can be built by substructuring techniques (DD-methods). We refer to [9–12, 28–31], in the case of the h-version of the FEM and to [2,3,16,22–24,27] for the p-version of the FEM. From the algebraic point of view, the symmetric positive definite (spd) stiffness matrix A is split into I 0 S 0 A11 A12 I A12 A−1 22 (1.3) = A= 0 A22 A21 A22 A−1 0 I I 22 A21 with the Schur-complement S = A11 − A12 A−1 22 A21 . Then, the preconditioner CS 0 I −E T I 0 C= 0 I −E I 0 C22
(1.4)
is introduced. The following result has been proved in [18,28,29,32] and is the key in order to analyze the preconditioner C. Lemma 1.1. Let CS and C22 be spd preconditioners for S and A22 , i.e., cS (CS g, g) (Sg, g) CS (CS g, g) ∀g ∈ Rm , cI (C22 v, v) (A22 v, v) CI (C22 v, v) ∀v ∈ Rn . Moreover, let I I g cE2 (Sg, g) A g, E E
∀g ∈ Rm .
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
393
Then, the inequalities c(Cx, x) (Ax, x) C(Cx, x) ∀x ∈ Rn+m hold with c = 12 cE−2 min{cs , cI } and C = 2cE2 max{CS , CI }. In the p-version of the FEM, the splitting of the matrix A is naturally given by a splitting of the ansatz functions. Let nv , ne , and ni be the number of vertices not belonging to Γ 1 , number of edges not belonging to Γ 1 , and number of elements, respectively. One basis function corresponds to each vertex, p − 1 basis functions correspond to each edge, (p − 1)2 basis functions correspond to each element. Thus, the dimension of the ansatz space is N = nv + (p − 1)ne + (p − 1)2 ni . We define the functions ζ1 , . . . , ζnv as the usual piecewise bilinear hat functions. The functions ζnv +(j −1)(p−1)+1 , . . . , ζnv +j (p−1) correspond to the edge ej of the mesh, and vanish on all other edges, i.e., satisfy the condition ζnv +(j −1)(p−1)+k−1 |el = δj,l pk , where pk is a polynomial of degree p, k = 2, . . . , p. The support of an edge function is formed by those two elements, which have this edge ej in common. The remaining basis functions are interior bubble functions consisting of a support containing one element only. Now, the basis functions ζi are divided into three groups, • the vertex functions, • the edge bubble functions, • the interior bubbles. Corresponding to the division of the shape functions, the matrix A is split into the blocks Av Av,e Av,i A = Ae,v Ae Ae,i , Ai,v Ai,e Ai where the indices v, e, and i denote the blocks corresponding to the vertex, edge bubble, and interior bubble functions, respectively. In a first step, the simpler matrix Av 0 0 Ae Ae,i , C= 0 (1.5) 0 Ai,e Ai is investigated. It has been proved in [3] that the condition number κ(C −1 A) grows as 1 + log p. Therefore, the vertex unknowns can be determined separately. Preconditioners for Av are multi-grid methods, [20], or BPX-preconditioners, [13,35]. Computing the other unknowns, we factorize the remaining 2 by 2 block I 0 S 0 Ae Ae,i I Ae,i A−1 i (1.6) = 0 Ai Ai,e Ai 0 I I A−1 i Ai,e with the Schur-complement S = Ae − Ae,i A−1 i Ai,e . The matrix Ai is a block diagonal matrix, each block ARs corresponds to one element Rs , i.e., i . Ai = blockdiag[ARs ]ns=1
(1.7)
394
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
Therefore, in order to compute the interior unknowns, we have to solve a Dirichlet problem on each quadrilateral. The edge unknowns are computed via the Schur-complement S. Due to Lemma 1.1, an inexact DD-preconditioner for (1.1) includes a preconditioner for Ai , a preconditioner for the Schurcomplement S and an extension operator E operating from the edges of the quadrilateral into its interior in order to replace the matrix A−1 i Ai,e by a matrix −E.
which is the result of assembling the element In [24], Jensen/Korneev have introduced a matrix A stiffness matrices on the reference element R2 = [−1, 1]2 instead of the element stiffness matrix corre −1 A) = O(1). This result allows sponding to the element Rs . Then, they proved in Lemma 4.2., that κ(A us to restrict ourselves to the case of the reference element in order to derive preconditioners for Ai , S and to find the extension operator E. In [22–24], two preconditioners for the Schur-complement are proposed. The preconditioners use basis transformations from the integrated Legendre polynomials to the Chebyshev polynomials or a Lagrange basis. Another preconditioning techniques can be found in [1,17]. Moreover, an approach using multiresolution bases for the Schur-complement in 3D, [27], can be applied to the 2D case as well. For the interior solver, several preconditioners Ci are developed in [5,6,8,24,26]. The condition number of Ci−1 Ai is bounded independent of the polynomial degree, [7]. The operation Ci−1 w requires O(p2 ) operations. All preconditioners use interpretations of the matrix AR2 as an h-version FEM discretization matrix of a degenerated elliptic problem, published in [5]. In order to replace A−1 i Ai,e in (1.6), an efficient extension operator E ↔ E is required, cf. Lemma 1.1. Given a polynomial g ∈ Pp on one edge of R2 , find a polynomial u ∈ Qp with uH 1 (R2 ) cE gH 0.5 (∂ R2 )
∀g ∈ Pp
(1.8)
under the constraint u = g on ∂R2 , where the constant cE is independent of the polynomial degree p. Usually, the system Ai u = −Ai,e g
or
Ai u = −Ai,f g
(1.9)
can be treated by a preconditioned Chebyshev iteration with the preconditioner Ci for Ai , see, e.g., [27]. Since this operation is required in each preconditioning step, cf. (1.3), the system solve (1.9) is too expensive. In the case of the h-version of the FEM, i.e., the function g is a piecewise linear function, Nepomnyaschikh, [19,32], derived several techniques in order to develop such an extension operator without the solution of (1.9). For the p-version of the FEM, a pioneering work is done by [3]. The practical and efficient implementation of this extension operator in certain polynomial bases, i.e., integrated Legendre polynomials, has been an open question. This paper is dedicated to the development of fast solvers for the system (1.2) arising from the discretization by the p-version of the FEM. Efficient tools are domain decomposition methods. We mainly focus on one ingredient of the inexact domain decomposition preconditioner, the extension operator. On the one hand, we will prove the estimate (1.8). On the other hand, we will show that this extension can be computed in O(p2 ) operations. Our technique uses a basis for the given function g in (1.8) which is stable in L2 and H 1 . With methods of multi-resolution analysis, cf. [14,15,33] for the h-version of the FEM, and [8] for the p-version of the FEM, such bases can be determined. The paper is organized as follows. In Section 2, the model problem in order to derive the extension operator is introduced. Then, the algorithm for the general definition of extension operators is defined. In Section 3, some examples of extension operators for the h- and the p-version of the FEM in 2D are given. In Section 4, we show the numerical performance of the proposed extension operators.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
395
Throughout this paper, the notation g is used for vectors of the Euclidean space, where the corre sponding function g(x) = nj=1 gj φj (x) is denoted by g. We write g ↔ g. Moreover, let I = (−1, 1). The parameter p denotes the polynomial degree and Pp and Qp denote the space of all polynomials of degree p in one and two variables. Pp,00 denotes the subspace of Pp satisfying q(±1) = 0. 2. Construction of extension operators 2.1. Extensions from one face or edge We consider the following problem: −u = 0 in Ω = I × ω, u(1, y) = 0, u(−1, y) = g(y) u(x, y) = 0 on I × ∂ω,
(2.1) (2.2) (2.3)
on ω,
where I = (−1, 1). Problem (2.1)–(2.3) is the typical model problem for the harmonic extension of a function g on ω into the interior of the domain Ω in the tensor product case. The weak formulation of (2.1)–(2.3) is: find u ∈ Hg (Ω) = {u = uˆ + g, uˆ ∈ H01 (Ω)} such that (2.4) a(u, v) = ∇u · ∇v = 0, ∀v ∈ H01 (Ω) Ω
holds. Problem (2.1)–(2.3) is solved approximately by a tensor product version of the finite element method (FEM). Let Ξ = {ξj (x)}m j =1 be a basis for the discretization of the interval I = (−1, 1) with ξj (±1) = 0 for j = 2, . . . , m, ξ1 (1) = 0 and ξ1 (−1) = 1. Moreover, let X = {χj (y)}nj=1 be a basis for the discretization of ω satisfying χj (∂ω) = 0. On Ω, we use the approximation space m,n V = span φij (x, y) i,j =1 , where φij (x, y) = ξi (x)χj (y). In order to satisfy the boundary condition (2.2), the functions φ1j , 1 j n, are defined, i.e., we assume that g(y) = nj=1 gj φ1 (−1, y) with the given coefficients gj . of problem (2.4) onto the space V leads to the following problem: find u = nThe Galerkinprojection m n g φ + u j 1j j =1 i=2 j =1 ij φij such that a(u, φij ) = 0,
∀i = 2, . . . , m, j = 1, . . . , n.
Introducing the matrices (m,n) K22 = a(φij , φlk ) (i,j ),(k,l)=(2,1) ,
n K11 = a(φ1j , φ1k ) k;j =1 ,
(2.5) m,n;n K21 = a(φij , φ1k ) i=2,j =1;k=1 ,
problem (2.5) is equivalent to solve the system of linear algebraic equations K22 u = −K21 g
(2.6)
with u = [u21 , . . . , umn ]T and g = [g1 , . . . , gn ]T . We are not interested in solving the problem (2.5), (2.6) exactly, we only try to find a vector u satisfying the inequality T K11 K21 (2.7) x, x = (Kx, x) cE2 (Sg, g) ∀g K21 K22
396
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
−1 T under the constraint x = gu . The matrix S = K11 − K21 K22 K21 denotes the Schur complement. The constant cE should not depend on the discretization parameter. In the case of the h-Version of the FEM, several techniques in order to construct a vector u ↔ u satisfying (2.7) are given in [32], where the system (2.6) has not to be solved. By the usual FEM-isomorphism, one has g 2 (Kx, x) =| u |H 1 (Ω) with x = ↔u (2.8) u between the vector x = gu and the corresponding function u. Due to (2.2) and the Poincare–Friedrichs −2 inequality, the H 1 -seminorm is equivalent to the H 1 -norm. Thus, one obtains CΩ u2H 1 (Ω) (Kx, x) 2 uH 1 (Ω) for all u, where the constant CΩ depends on the domain Ω only. Concerning a relation between the function g and the vector g, the following assumption is required.
Assumption 2.1. Let us assume that for every g ∈ span X there exists a function u such that uH 1 (Ω) cH gH 0.5 (∂Ω) ,
(2.9)
under the constraints u(−1, y) = g(y), u(1, y) = 0, u(x, ∂ω) = 0, where the constant c does not depend on the discretization parameter. Then, the following result can be proved. Lemma 2.2. Let us assume that Assumption 2.1 is satisfied and that ∂Ω is Lipschitz continuous. Let g = nj=1 gj φj and g = [gj ]nj=1 . Then, we have the norm equivalence relation 2 cT2 g2H 0.5 (∂Ω) (Sg, g) cH g2H 0.5 (∂Ω) ,
(2.10)
where the constant cH is the constant in (2.9) and cT is the constant from the trace theorem. Proof. The proof is similar to the arguments in order to derive formula (3.12) in [24]. By the properties of the Schur-complement, we have (Sg, g) = minu (Ku, u), where x = gu . Let u∗ ↔ x ∗ = ug∗ be the optimal solution of this minimization problem. Thus, by the trace theorem, cT2 g2H 0.5 (∂Ω) u∗ 2H 1 (Ω) = Kx ∗ , x ∗ = (Sg, g). Moreover, let u0 ↔ x0 = ug0 be the extension of Assumption 2.1. Then, 2 g2H 0.5 (∂Ω) , (Sg, g) (Kx0 , x0 ) = u0 2H 1 (Ω) cH
which proves the lemma.
2
This lemma is important in order to derive a preconditioner for the Schur-complement embedded in the domain decomposition preconditioner (1.4). Using this norm equivalence, most preconditioners for 0.5 S can be constructed by finding a stable basis in H (∂Ω). Moreover, we are able to define (Sg, g) as an equivalent norm on H 0.5 (∂Ω). The next lemma gives a relation for the extension and the corresponding norms.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
Lemma 2.3. Let us assume that the extension satisfies (2.7). Then, x =
g u
397
↔ u satisfies
uH 1 (Ω) CΩ cE cH gH 0.5 (∂Ω) .
(2.11)
Proof. Use relations (2.10), (2.7), (2.8) and the Poincare–Friedrichs inequality.
2
In order to define the extension operator, we use a technique which is similar to the multi-level decomposition technique derived in [19]. In the space Y = span X, a so-called stable basis is used. Assumption 2.4. Let us assume that there exists a basis Ψ = {ψj }nj=1 in the space Y which is stable in L2 (ω) and H 1 (ω), i.e., for any function g = nj=1 g˜ j ψj (y) the norm equivalence relations c1,1
n
dj(1) g˜ j2 g2L2 (ω) c2,1
j =1
c1,2
n
n
dj(1) g˜ j2 ,
(2.12)
j =1
dj(2) g˜ j2
g2H 1 (ω)
c2,2
j =1
n
dj(2) g˜ j2
(2.13)
j =1
hold with some numbers dj(i) > 0, j = 1, . . . , n, i = 1, 2. The constants ci,j , i, j = 1, 2, do not depend on the discretization parameter n. The computation of the nearly discrete harmonic extension of the function g(y) = nj=1 gj χj (y) consists of three steps: Algorithm 2.5. 1. Transform the function g into the basis Ψ = {ψj }nj=1 , i.e., g(y) =
n
g˜ j ψj (y).
j =1
2. For j = 1, . . . , n, let u∗j be the solution of the problem |uj|H1 (Ω) → min,
where uj (x, y) = ψj (y) ξ1 (x) +
fi
m
ξi (x)fi .
(2.14)
i=2
Put u(x, y) =
n
g˜ j u∗j (x, y).
j =1
3. Transform u(x, y) =
n
˜ j uj (x, y) j =1 g
(2.15) into the basis {φij (x, y)}i,j .
We introduce the following real numbers, vectors and matrices: let m m a 1 = ξ1 , ξj I j =2 , A1,0 = ξi , ξj I i,j =2 , βj,1 = ψj , ψj ω , m m βj,2 = Dψj , Dψj ω , a 2 = Dξ1 , Dξj I j =2 , A2,0 = Dξi , Dξj I i,j =2 ,
(2.16)
398
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
where ·, · Ω denotes the L2 (Ω)-scalar product and Du is the gradient of u. In order to compute the extension efficiently, the following lemma and remark are useful. Lemma 2.6. The computation of uj in (2.14) requires the system solve (βj,1 A1,0 + βj,2 A2,0 )f = −βj,1 a 1 − βj,2 a 2 ,
(2.17)
where βj,1 , βj,2 > 0. Proof. Inserting (2.14) into uj leads to the optimization problem 2 m ξi (x)fi → min . ψj (y)ξ1 (x) + fi 1 i=2
H (Ω)
Using (2.16) and tensor product arguments, this is equivalent to 2 2 m m βj,1 ξ1 (x) + ξi (x)fi + βj,2 ξ1 (x) + ξi (x)fi 1 i=2
i=2
H (I )
L2 (I )
→ min . fi
The solution of this quadratic optimization problem with respect to f = [fk ]m k=2 leads to the system βj,1 Dξi , Dξk I + βj,2 ξi , ξk I fk = −βj,1 Dξi , Dξ1 I − βj,2 ξi , ξ1 I , i = 2, . . . , m. Inserting now (2.16) proves the lemma.
2
Remark 2.7. In (1), one basis change from the basis X to the basis Ψ has to be done, in (3), n basis changes from the Ψ to the basis X have to be done. (2) means from the algebraic point of view n system solves with positive linear combinations of the matrices A1,0 and A2,0 , cf. Lemma 2.6. The matrices A1,0 and A2,0 are the 1D mass and 1D stiffness matrix with respect to the basis Ψ . Remark 2.8. Moreover, we investigate the following modification of Algorithm 2.5. In this second approach, instead of the systems (2.17) in order to compute uj (2.14), the systems (1) dj A1,0 + dj(2) A2,0 f = −dj(1) a 1 − dj(2) a 2 , with the constants dj(i) in (2.12) and (2.13) have to be solved. This is a modification of Algorithm 2.5. 2.2. The general case In the last subsection, only an extension of a function living on one edge (in 2D) or face (in 3D) of the domain Ω has been investigated. This subsection deals with the general case of an extension from ∂Ω into Ω. Usually, the boundary of Ω is split into several faces (edges in 2D) ωi , i = 1, . . . , r, r i=1 ωi = ∂Ω, ωi ∩ ωj = ∅ for i = j . On each ωi , a function gi is given which satisfies the boundary condition gi (∂ωi ) = 0. Then, we consider the problem: find u ∈ H 1 (Ω) such that u2H 1 (Ω) cg2H 0.5 (∂Ω) ,
u|ωi = gi = g|ωi , i = 1, . . . , r.
(2.18)
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
399
This problem can be solved by the method described in the previous subsection (Algorithm 2.5). Let ui be the extension of ui |ωj = gi δij
on ωj and ui 2H 1 (Ω) ci gi 2H 0.5 (ωi ) .
(2.19)
Then, we define u = ri=1 ui as the extension of the function g into Ω. Obviously, u|ωj = ri=1 ui |ωj = r i=1 gi δij = gj , i.e., condition (2.18) is fulfilled. Moreover, by the Cauchy–Schwarz inequality and (2.19), we obtain u2H 1 (Ω)
r 2 = uj 1 j =1
r
H (Ω)
r
r
uj 2H 1 (Ω) r
j =1
r j =1
cj gj 2H 0.5 (ωj ) r max {cj } j =1,...,r
r j =1
gj 2H 0.5 (ωj )
max {cj }g2H 0.5 (∂Ω) . j =1,...,r
In the last estimate, a relation between the norms · 2H 0.5 (∂Ω) and rj =1 · 2H 0.5 (ωj ) is required. In the second norm, the coupling between two neighbouring faces or edges ωi and ωj has no influence in the norm, whereas an additional connectivity term is involved in the definition of the first norm. So, with an optimal extension for problem (2.19), cf. Algorithm 2.5, an optimal extension for (2.18) can be computed.
3. Applications In this section, several examples for extension operators of the type proposed in Algorithm 2.5 are given. So, Assumption 2.4 has to be verified, i.e., a polynomial basis Ψ has to be derived which is stable in H 1 (I ) and L2 (I ). Using methods of multi-resolution analysis, [8,14,15,33], such a basis Ψ can be constructed in the case of h- and p-version of the FEM. We start with the p-version of the FEM. 3.1. Extension operators for the p-version of the FEM In this subsection, we describe a method deriving an extension operator E, i.e., u = Eg in matrix representation, in the case of the p-version of the FEM. Usually, such an extension operator is embedded in a domain decomposition preconditioner of Dirichlet–Dirichlet-type, cf. Section 1, or [5]. The operator u = Eg will require O(p2 ) arithmetical operations in 2D, i.e., it is arithmetically optimal. Moreover, numerical experiments indicate the estimate I I g, g c(Sg, g) ∀g ∈ Rm , (Kx, x) = K E E where c is independent of the polynomial degree p. A system solve K22 u = −K21 g is not required. In Section 3.1.1 a method in order to construct stable polynomial bases in the spaces H 1 (I ) and L2 (I ) is shown. In Section 3.1.2, we derive the extension operators for the 2D case.
400
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
3.1.1. A stable polynomial basis of Pp,00 in H 1 (I ) and L2 (I ) In this subsubsection, we propose a basis in Pp,00 , which is stable in the spaces L2 (I ) and H 1 (I ). For i ∈ N, let Li (x) = ( dxd )i (x 2 − 1)i be the ith Legendre polynomial and let √ x (2i − 3)(2i − 1)(2i + 1) i (x) = Li−1 (s) ds (3.1) L 2 −1
i }pi=2 is a basis in i (±1) = 0. Thus, {L be the ith integrated Legendre polynomial (i 2). Note that L the space Pp,00 . Moreover, we introduce the 1D mass and 1D stiffness matrix with respect to the basis i }pi=2 , i.e., {L j , L i I p j , D L i I p . and H1,I = D L (3.2) L2,I = L i,j =2 i,j =2 In [7], we have proved that 1 0 −γ2 0 1 0 −γ2 0 1 L2,I = .. . 0
...
0
0 −γ3 0 .. .
−γ4
−γp−2
0
2 0 −1 ∼ .. . 1
0
0 2 0
−1 0 0 −1 2 0 −1 .. .
...
−1
Moreover, an easy calculation shows, [24], (2j − 3)(2j + 1) p . H1,I = diag 2 j =2
0
.
(3.3)
2
(3.4)
Here and in the following, the relation A ∼ B means that c1 (Av, v) (Bv, v) c2 (Av, v) ∀v, where the constants c1 and c2 do not depend on the polynomial degree p. Using a permutation P of rows and columns, one obtains T L1 T H1 P and H1,I = P P, (3.5) L2,I ∼ P L1 H1 where L1 is a tridiagonal matrix with 2’s on the main diagonal and −1’s on the first subdiagonal. Obi }pi=2 is stable in L2 (I ). In i }pi=2 is stable in H1 (I ). However, we are not able to prove that {L viously, {L order to find a stable basis, we use interpretations of the matrices L1 and H1 as discretization matrices on the unit interval (0, 1) using piecewise linear ansatz functions. ], i = 0, . . . , n − 1, be an equidistant Let k ∈ N, n = 2k and p = 2n − 1. Moreover, let τi = [ ni , i+1 n be the basis of the usual hat functions mesh on the interval [0, 1] and Θ = {θi }n−1 i=1 nx − (i − 1) on τi−1 , θi (x) = (i + 1) − nx on τi , i = 1, . . . , n − 1. (3.6) 0 else, We introduce the matrices 1 n−1 1 n−1 θi (x)θj (x)x 2 dx and T Θ = θi (x)θj (x) dx . MΘ = −1
i,j =1
In [4], we have proved the following result.
−1
i,j =1
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
401
Lemma 3.1. The spectral equivalence relations H1 ∼ n3 M Θ and L1 = n1 T Θ are valid. Concerning the weighted mass matrix M Θ and the unweighted stiffness matrix T Θ , it is known from the wavelet theory, [8], that it can be constructed a multi-level basis in which M Θ and T Θ are spectrally equivalent to diagonal matrices. Let J = (j, l), where 0 l k and j = 1, 3, 5, . . . , 2l − 1, card(J ) = n − 1. The diagonal matrices (3.7) D1,0 = diag 22l J =(j,l) and D2,0 = diag j 2 2−2l J =(j,l) are introduced. Theorem 3.2. There are multi-level bases Ψ˜ = {ψ˜ J }J with Ψ˜ = ΘQ such that the matrices T Ψ = QT T Θ Q and M Ψ = QT M Θ Q satisfy the spectral equivalence relations T Ψ ∼ D1,0 and M Ψ ∼ D2,0 . The transformation v = QT w can be performed in O(n) operations. Examples for wavelets bases are given in [8]. Now, let Q 0 , W = PT 0 Q
(3.8)
where P denotes the permutation matrix in (3.5) and Q denotes the fast wavelet transform in Theorem 3.2. Theorem 3.3. The matrices L2,I and H1,I satisfy the relations 1 −T D1,0 0 −1 3 −T D2,0 W and H1,I ∼ n W L2,I ∼ W 0 D1,0 0 n
0 W −1 . D2,0
(3.9)
Proof. Using (3.5), Lemma 3.1 and Theorem 3.2, we have PT TΘ 1 T Q−T T Ψ Q−1 0 0 0 T L1 P= P= P P L2,I ∼ P 0 L1 0 TΘ 0 Q−T T Ψ Q−1 n n −1 P T Q−T 0 0 0 D1,0 Q ∼ P. 0 Q−T 0 D1,0 0 Q−1 n Since P is a orthogonal matrix, the first assertion follows from (3.8). The second assertion can be proved in the same way. 2 Via the matrix W (3.8), the basis 2 (x), . . . , L p (x) W Ψ = ψ1 (x), . . . , ψp−1 (x) = L
(3.10)
is introduced. Corollary 3.4. The basis Ψ is stable in L2 (I ) and H 1 (I ). Moreover, the operation v = W w can be performed in O(p) arithmetical operations. p j (x) ∈ Pp,00 from the So, we have found an algorithm which transforms a function g(x) = j =2 gj L basis of the integrated Legendre polynomials into the basis Ψ ⊂ Pp,00 which is stable in L2 (I ) and H 1 (I ) and requires O(p) operations.
402
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
3.1.2. The 2D case We consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the p-version of the FEM with ω = (−1, 1), Ω = (−1, 1)2 . We use the integrated Legendre polynomials (3.1) for the discretization in x- and y-direction, i.e., j (x), ξj (x) = L j +1 (y), χj (y) = L
j = 2, . . . , p, ξ1 (x) =
1−x 2
and
j = 1, . . . , p − 1.
Using the notation (2.16), an easy calculation shows, cf. [24], T 5 7 , , 0, . . . , 0 , A1,0 = L2,I , a1 = − 12 60 a 2 = [0, . . . , 0]T ,
A2,0 = H1,I .
(3.11)
We use Algorithm 2.5 in order to compute a nearly optimal discrete harmonic extension. Since Corollary 3.4, the basis Ψ (3.10) is stable in L2 (I ) and H1 (I ). Thus, Assumptions 2.4 are satisfied. Lemma 3.5. Let x = gu ↔ u be defined via Algorithm 2.5 or the generalization of Remark 2.8. Then, the computation of u requires O(p2 ) operations. Proof. Since the multiplication w = Qv requires O(p) operations, the cost of p − 1 operations w = W v is of order p2 . Due to (3.3) and (3.4), the matrices dj(1) L2,I + dj(2) H1,I are block diagonal matrices of two tridiagonal matrices after reordering the unknowns, dj(i) ∈ R, dj(i) > 0, i = 1, 2, j = 1, . . . , p − 1. Using Cholesky decomposition, the cost of (dj(1) L2,I + dj(2) H1,I )−1 v is O(p). Solving p − 1 systems, cf. Lemma 2.6, the total cost is of step 2 in Algorithm 2.5 of order p2 . For the operation W −1 g, an inverse wavelet transform is required, where W ∈ Rp−1×p−1 . Since W is a product of [log2 p] matrices with a bandwidth O(1), the cost for W −1 g is of order [log2 p]p using Gaussian elimination. 2 Remark 3.6. This technique in order to construct the extension can be applied to the 3D-case. Then, we consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the p-version of the FEM with ω = (−1, 1)2 , Ω = (−1, 1)3 . With tensor product arguments, the basis Ψ × Ψ is stable in the tensor product spaces L2 (I × I ) = L2 (I ) × L2 (I ) and H 1 (I × I ) = H 1 (I ) × L2 (I ) + L2 (I ) × H 1 (I ). 3.2. h-version of the FEM in 2D and 3D We consider the discretization (2.5) of problem (2.1)–(2.3) in the case of the h-version of the FEM ], i = with ω = (−1, 1), Ω = (−1, 1)2 . Let k ∈ N, n = 2k+1 − 1, s = 2k . Moreover, let τi = [ si , i+1 s s−1 −s, . . . , s − 1, be an equidistant mesh on the interval I and Θ = {θi }i=−s+1 be the basis of the usual hat functions (3.6). Then, we choose χj (y) = θj +s (y) for j = −s + 1, . . . , s − 1. Moreover, let τ˜i , i = 1, . . . , m be an arbitrary 1D finite element mesh of (−1, 1). Then, Ξ is the basis of the usual hat functions according to the mesh τ˜i . It is known from the theory of multi-resolution bases, [14,15,33], that several wavelet bases on the interval are stable in H 1 and L2 . Hence, Assumption 2.4 is satisfied. Moreover, it has been shown
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
403
in [19] that there exists a piecewise linear harmonic extension u of the piecewise linear function g = s−1 j =−s+1 gj θj into the interior of the domain Ω. Thus, we have verified the second Assumption 2.1, too. So, the optimal extension can be computed via Algorithm 2.5. For the computation of u, n systems of linear algebraic equations with a linear combination Mj , j = 1, . . . , n, of the one-dimensional mass matrix and one-dimensional stiffness matrix with respect to the basis Ξ have to be solved, cf. Remark 2.7. Since Ξ is the basis of the usual hat-functions, the system matrices Mj are tridiagonal matrices. Thus, one system solve requires O(m) operations and the total cost is O(mn). The approach can be extended to the 3D-case, where ω is a surface in R3 , [21]. 4. Numerical examples In this section, we compute the smallest possible constant cE in (2.7) for the proposed extension operator in the case of the p-version of the finite element method in two dimensions. In all experiments, the wavelet basis described in [7] is used with two types of scalings, an explicit one given by Theorem 3.3 and a implicit one given by the numbers dj(i) = βj,i , i = 1, 2, j = 1, . . . , n. In the first experiment, we use the numbers dj(i) , i = 1, 2, j = 1, . . . , n, defined via (3.7) and (3.9) instead of the numbers βj,i (2.16). This is a slight modification of Algorithm 2.5, cf. Remark 2.8. In the second example, we take dj(i) = βj,i , i = 1, 2, j = 1, . . . , n. Using (3.2) and (3.8), these entries can explicitly be computed in O(p2 ) operations. Table 1 displays the quality of the extension of Algorithm 2.5, i.e., Remark 2.8, in the first example. Moreover, the best constants c1,1 , c1,2 , c2,1 and c2,2 in (2.12) and (2.13) are computed. For all p, the constant cE2 lies between 1.4 and 1.9 and is bounded by a constant independent of p. Table 2 displays the constants cE , ci,j , i, j = 1, 2, for the second example. Now, the constants cE and c1,1 , c2,1 and c2,2 are closer to one. The only difference between the two proposed versions of extension is the type of diagonal scaling of the basis Ψ in L2 (ω) and H 1 (ω), i.e., the choice of the parameters dj(i) . So, the scaling factors dj(i) > 0 in Table 1 Smallest possible constant of (2.7) for the p-version extension operator p 3 7 15 31 2 in (2.7) cE 1.45 1.89 1.85 1.79 c1,1 0.22 0.20 0.18 0.17 0.94 0.37 0.25 0.21 c1,2 2.23 2.35 2.39 2.41 c2,1 3.94 3.94 3.42 3.23 c2,2
63 1.76 0.17 0.19 2.42 3.15
Table 2 Smallest possible constant of (2.7) for the p-version extension operator with implicitly given βj,i p 3 7 15 31 63 2 in (2.7) cE 1.00 1.13 1.26 1.34 1.37 c1,1 0.36 0.36 0.36 0.36 0.36 1.00 0.33 0.25 0.21 0.19 c1,2 1.63 1.68 1.77 1.83 1.88 c2,1 1.00 1.96 2.16 2.27 2.35 c2,2
127 1.75 0.16 0.18 2.43 3.11
127 1.39 0.36 0.18 1.92 2.40
404
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
order to solve dj(1) L2,I + dj(i) H1,I , cf. the proof of Lemma 3.5, are different between the first and second experiment. Since these systems are solved by a direct method and dj(i) > 0, the arithmetical cost for the computation of the extension itself is the same. For the second experiment, an additional initialization has to be done, i.e., the entries dj(i) = βj,i have to computed which is not required for the first experiment. Usually, the extension is embedded in a domain decomposition preconditioner C for A (1.2). Since the quality of the extension is not so good for the first experiment, the number of the pcg iterations may be higher than using the extension of the second experiment. The additional cost for the second experiment, i.e., the computation of the entries dj(i) = βj,i , can be done before starting the pcg method. Moreover, the arithmetical cost in order to compute βj,i is comparable with the cost for the computation of the extension one times. Summarizing, the method which uses dj(i) = βj,i , i.e., the second experiment should be preferred.
References [1] M. Ainsworth, A preconditioner based on domain decomposition for h–p finite element approximation on quasi-uniform meshes, SIAM J. Numer. Anal. 33 (4) (1996) 1358–1376. [2] M. Ainsworth, B. Gao, An additive Schwarz preconditioner for p-version boundary element approximation of the hypersingular operator in three dimensions, Numer. Math. 85 (3) (2000) 343–366. [3] I. Babuška, A. Craig, J. Mandel, J. Pitkäranta, Efficient preconditioning for the p-version finite element method in two dimensions, SIAM J. Numer. Anal. 28 (3) (1991) 624–661. [4] S. Beuchler, Preconditioning for the p-version of the FEM by bilinear elements, Technical Report SFB393 01-17, Technische Universität Chemnitz, May 2001. [5] S. Beuchler, Multi-grid solver for the inner problem in domain decomposition methods for p-FEM, SIAM J. Numer. Anal. 40 (3) (2002) 928–944. [6] S. Beuchler, AMLI preconditioner for the p-version of the FEM, Numer. Linear Algebra Appl. 10 (8) (2003) 721–732. [7] S. Beuchler, Optimal preconditioners for the p-version of the FEM, Technical Report SFB393 03-03, Technische Universität Chemnitz, March 2003. [8] S. Beuchler, R. Schneider, C. Schwab, Multiresolution weighted norm equivalences and applications, Numer. Math. 98 (1) (2004) 67–97. [9] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, I, Math. Comp. 47 (175) (1986) 103–134. [10] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, II, Math. Comp. 49 (179) (1987) 1–16. [11] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, III, Math. Comp. 51 (184) (1988) 415–430. [12] J. Bramble, J. Pasciak, A. Schatz, The construction of preconditioners for elliptic problems by substructuring, IV, Math. Comp. 53 (187) (1989) 1–24. [13] J. Bramble, J. Pasciak, J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (191) (1991) 1–22. [14] A. Cohen, I. Daubechies, P. Vial, Biorthogonal spline wavelets on the interval—stability and moment conditions, Appl. Comput. Harm. Anal. 6 (2) (1998) 132–196. [15] W. Dahmen, Wavelet and multiscale methods for operator equations, Acta Numerica 6 (1997) 55–228. [16] B. Guo, W. Cao, A preconditioner for the hp-version of the finite element method in two dimensions, Numer. Math. 75 (1996) 59–77. [17] B. Guo, W. Cao, An iterative and parallel solver based on domain decomposition for the hp-version of the finite element method, J. Comput. Appl. Math. 83 (1997) 71–85. [18] G. Haase, U. Langer, A. Meyer, The approximate Dirichlet domain decomposition method. Part I: An algebraic approach, Computing 47 (1991) 137–151.
S. Beuchler, J. Schöberl / Applied Numerical Mathematics 54 (2005) 391–405
405
[19] G. Haase, S.V. Nepomnyaschikh, Explicit extension operators on hierarchical grids, East–West J. Numer. Math. 5 (4) (1998) 231–284. [20] W. Hackbusch, Multigrid Methods and Applications, Springer, Heidelberg, 1985. [21] H. Harbrecht, R. Schneider, Biorthogonal wavelet bases for the boundary element method, Math. Nachr. 269–270 (1) (2004) 167–188. [22] S.A. Ivanov, V.G. Korneev, On the preconditioning in the domain decomposition technique for the p-version finite element method, Part I, Technical Report SPC 95-35, Technische Universität Chemnitz-Zwickau, December 1995. [23] S.A. Ivanov, V.G. Korneev, On the preconditioning in the domain decomposition technique for the p-version finite element method, Part II, Technical Report SPC 95-36, Technische Universität Chemnitz-Zwickau, December 1995. [24] S. Jensen, V.G. Korneev, On domain decomposition preconditioning in the hierarchical p-version of the finite element method, Comput. Methods Appl. Mech. Engrg. 150 (1–4) (1997) 215–238. [25] G.M. Karniadakis, S.J. Sherwin, Spectral/HP Element Methods for CFD, Oxford University Press, Oxford, 1999. [26] V.G. Korneev, Poqti optimalny metod rexeni zadaq Dirihle na podoblasth dekompozicii ierarhiqesko hp-versii, Differencialnye Uravneni 37 (7) (2001) 1–15. [27] V. Korneev, U. Langer, L. Xanthis, On fast domain decomposition methods solving procedures for hp-discretizations of 3d elliptic problems, Technical Report 03-07, Universität Linz, March 2003. [28] P.L. Lions, On the Schwarz alternating method I, in: R. Glowinski, G.H. Golub, G.A. Meurant, J. Périaux (Eds.), Proc. 1st Int. Symp. on Domain Decomposition Methods, SIAM, Philadelphia, PA, 1988, pp. 1–42. [29] P.L. Lions, On the Schwarz alternating method, II, in: Proc. 1st Internat. Symp. on Domain Decomposition Methods, Los Angeles, 1989, pp. 47–70. [30] A.M. Matsokin, S.V. Nepomnyaschikh, The Schwarz alternation method in a subspace, Izv. Vyssh. Uchebn. Zaved. Mat. 29 (10) (1985) 61–66. [31] A.M. Matsokin, S.V. Nepomnyaschikh, Norms in the space of traces of mesh functions, Sov. J. Numer. Anal. Math. Modelling 3 (3) (1988) 199–216. [32] S.V. Nepomnyaschikh, Method of splitting into subspaces for solving elliptic boundary value problems in Complex-form domains, Sov. J. Numer. Anal. Math. Modelling 6 (1991) 151–168. [33] R. Schneider, Multiskalen- und Wavelet Matrixkompression, Teubner, Leipzig, 1998. [34] C. Schwab, p- and hp-Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics, Clarendon, Oxford, 1998. [35] X. Zhang, Multilevel Schwarz methods, Numer. Math. 63 (1992) 521–539.