Numerische

Numer. Math. (****) ***:* **

Mathematik

DOI 10.1007/s00211-011-0401-4

Optimal multilevel methods for graded bisection grids

Long Chen Ricardo H. Nochetto Jinchao Xu

Received: 13 October 2010 / Revised: 15 June 2011 / Published online: 21 July 2011

Springer-Verlag 2011

Abstract We design and analyze optimal additive and multiplicative multilevel

methods for solving H 1 problems on graded grids obtained by bisection. We deal

with economical local smoothers: after a global smoothing in the nest mesh, local

smoothing for each added node during the re nement needs to be performed only for

three vertices - the new vertex and its two parent vertices. We show that our methods

lead to optimal complexity for any dimensions and polynomial degree. The theory

hinges on a new decomposition of bisection grids in any dimension, which is of inde-

pendent interest and yields a corresponding decomposition of spaces. We use the latter

to bridge the gap between graded and quasi-uniform grids, for which the multilevel

theory is well-established.

65M55 65N55 65N22 65F10

Mathematics Subject Classi cation (2000)

L. Chen was supported in part by NSF Grant DMS-0505454, DMS-0811272, and in part by 2010 2011

UC Irvine Academic Senate Council on Research, Computing and Libraries (CORCL). R. H. Nochetto

was supported in part by NSF Grant DMS-0505454 and DMS-0807811. J. Xu was supported in part by

NSF DMS-0609727, DMS 0915153, NSFC-10528102 and Alexander von Humboldt Research Award for

Senior US Scientists.

L. Chen (B )

Department of Mathematics, University of California at Irvine, Irvine, CA 92697, USA

e-mail: abpviw@r.postjobfree.com

R. H. Nochetto

Department of Mathematics, University of Maryland, College Park, MD 20742, USA

e-mail: abpviw@r.postjobfree.com

J. Xu

The School of Mathematical Science, Peking University, Beijing, China

J. Xu

Department of Mathematics, Pennsylvania State University, University Park, PA 16801, USA

e-mail: abpviw@r.postjobfree.com

123

2 L. Chen et al.

1 Introduction

Adaptive methods are now widely used in scienti c and engineering computation to

optimize the relation between accuracy and computational labor (degrees of freedom).

Standard adaptive nite element methods (AFEM) based on local mesh re nement

can be written as loops of the form

ESTIMATE MARK REFINE. (1.1)

SOLVE

The module ESTIMATE determines a posteriori error estimators; we refer to [58]. The

module MARK selects elements with largest error indicators and is critical for con-

vergence and optimality of AFEM. Neither of these two procedures plays a role in

the present discussion. The module REFINE re nes all marked elements and perhaps a

few more to keep mesh conformity. Of all possible re nement strategies, we are inter-

ested in bisection, a popular, elegant, and effective procedure for re nement in any

dimension [3,8,32,33,48,49,54,55]. Our goal is to design optimal multilevel solvers

that constitute the core of procedure SOLVE, and analyze them within the framework

of highly graded meshes created by bisection, from now on called bisection grids.

It is important to realize that having optimal solvers in SOLVE is crucial for the the-

ory and practice of AFEM. Convergence of adaptive loops (1.1) for dimension d > 1

started with the seminal work of D r er [26]. This was followed by Morin et al.

[40,41], who realized the role of data oscillation, and Mekchay and Nochetto [36],

who proved a contraction property for AFEM for general elliptic PDEs. More recently,

Binev et al. [10] proved quasi-optimal cardinality for a modi ed algorithm including

coarsening. Stevenson [53] was able to remove coarsening that still needs an arti cial

inner loop. The most standard AFEM has been later examined by Casc n et al. [18],

who have proved a contraction property and quasi-optimal cardinality. We refer to

Nochetto et al. [42] for an introduction to the theory of adaptive nite element meth-

ods. The gap to obtaining optimal complexity is precisely having optimal solvers and

storage for adaptive bisection grids the topic of this paper.

We consider a nested family of nite element spaces obtained by local mesh re ne-

ment:

V0 V1 V J = V .

A standard multilevel method contains a smoothing step on the spaces V j, j =

0, . . ., J . For graded grids obtained by AFEM, it is possible that V j results from

V j 1 by just adding few, say one, basis function. Thus smoothing on both V j and

V j 1 leads to a lot of redundancy. If we let N be the number of unknowns in the

nest space V, then the complexity of smoothing can be as bad as O( N 2 ). To achieve

optimal complexity O( N ), the smoothing in each space V j must be restricted to the

new unknowns and their neighbors. Such methods are referred to as local multilevel

methods [4,17]. Performing the smoothing only on the newly added nodes, the most

extreme choice, gives rise to the hierarchical basis (HB) method [6,66].

Since the literature on local multilevel methods is abundant, we restrict ourselves to

describing the papers most relevant to graded meshes. Brandt [17] proposed the multi-

level adaptive technique (MLAT) and further studied it in [4]. McCormick et al. [34,35]

123

Optimal multilevel methods for graded bisection grids 3

developed the fast adaptive composite grid (FAC) method, which requires exact solvers

on subdomains that are partitioned by uniform grids; hence mesh adaptivity is achieved

via superposition of tensor-product rectangular grids. Rivara [47] and Mitchell [37,38]

developed local multigrid methods on adaptive triangular grids obtained by longest

edge bisection and newest vertex bisection, respectively, for d = 2. Also for d = 2,

Bank et al. [7] proposed the red green re nement strategy, which was implemented

in the well-known piecewise linear triangular multigrid software package (PLTMG)

of Bank [5]. For the resulting locally re ned grids, Bank et al. [6] developed HB

multigrid, which is a variant of the HB preconditioner developed earlier by Yser-

entant [65]. They proved that the hierarchical basis methods are nearly optimal (up

to a logarithmic factor on the number of elements) for d = 2, and suboptimal for

d > 2. Bramble et al. [15] proposed the BPX preconditioner on both quasi-uniform

grids and locally re ned grids, and showed that it is nearly optimal for any dimension

d . Oswald [44] was able to remove the logarithmic factor, and thus proved optimal

complexity of BPX preconditioner and established a similar result on locally re ned

grids in [45]. Bramble and Pasciak [14] proved the optimality of multilevel algorithms

including BPX preconditioner and V-cycle multigrid methods on quasi-uniform and

locally re ned grids. Dahmen and Kunoth [24] proved optimal complexity of BPX,

for graded meshes created by red green re nement for d = 2; see also Bornemann

and Yserentant [13] for a simpler approach using a K-functor. Griebel [27,30] devel-

oped multilevel methods on adaptive sparse grids. Xu [61] introduced a uni ed abstract

framework for the various multilevel methods based on the method of subspace correc-

tions, which is the approach we pursue in this paper. More recently, Wu and Chen [60]

analyzed multigrid methods using Gauss-Seidel type smoothers on bisection grids gen-

erated by newest vertex bisection for d = 2. Finally, Aksoylu and Holst [2] extended

the optimality results of [24] to a practical, local red green re nement procedure for

d = 3 due to Bornemann et al. [11].

We now summarize our four main contributions and place them in context.

First we present a novel decomposition of bisection grids. Roughly speaking, for

any triangulation T N constructed from T0 by N bisections, we can write

T N = T0 + B, B = (b1, b2, . . ., b N ) (1.2)

where B denotes an ordered sequence of N elementary bisections bi . Each bi is

restricted to a local region, the star of the newly created vertex, and the corre-

sponding local grid is thus quasi-uniform. This decomposition induces a space

decomposition {Vi }iN 1 of the underlying subspace of continuous piecewise linear

=

functions over T N . Moreover, this decomposition serves as a general bridge to

transfer results from quasi-uniform grids to graded bisection grids and has some

intrinsic interest. For example, it is this geometric structure of bisection grids that

motivates the new ef cient implementation of multilevel methods developed by

Chen and Zhang [22] for d = 2 and by Chen [19] for d = 3, which avoids dealing

with the tree structure of the mesh and hinges on coarsening from the nest mesh

T N to nd (1.2). Such a grid decomposition may or may not coincide with the one

giving rise to T N via AFEM, but it does not matter.

123

4 L. Chen et al.

We stress that the popularity of multilevel methods among practitioners is somehow

hampered by their complicated data structures to store and access the hierarchical

grid structure. Ef cient multilevel algorithms developed in [22] and [19] require

only minimal bookkeeping (recording the two parent nodes for each new bisection

node) and a simple data structure (only the nest grid must be stored instead of

the whole re nement tree). They rely on a clever coarsening algorithm for bisec-

tion grids and exhibit both optimal storage and complexity. This paper provides a

theoretical basis for these methods and establishes their optimality.

Second, we introduce and analyze economical local smoothers which reduce the

complexity of the resulting multilevel methods. In fact, besides a smoothing in the

nest grid T N, the local smoothing (or local relaxation) associated with each bisec-

tion bi is just performed for the newly created vertex and its two parent vertices.

This implies that dim Vi = 3, whence the total complexity is proportional to the

size of the linear system with a relatively small constant.

From the algorithmic point of view, our algorithm is different from the traditional

local multigrid, which requires smoothing for the new vertex and all neighboring

vertices (and degrees of freedom for quadratic or higher order elements) but no

additional smoothing in the nest grid. In one iteration of V-cycle or BPX pre-

conditioner, our algorithm requires less operations than the traditional one. For

example, for linear elements, treating the cost of the smoothing on one node as

a unit, our algorithm requires 4 N operations while the traditional one may need

k N operations, where k is the number of neighboring nodes surrounding a node.

For bisection grids, the average of k is 5 for d = 2 and could be as high as 10

for d = 3. Furthermore our algorithm is easier to implement, especially for higher

order elements. The smoothing on the nest grid can be easily realized using Jacobi

or Gauss-Seidel methods since the matrix is given while the multilevel smooth-

ing only involves linear elements and three points. Note that the new vertex and

two parents vertices are minimal information needed, in any geometric multigrid

method, to construct the restriction and prolongation operators.

Our three-point smoother, inspired by an idea of Stevenson [52] for wavelets, can

be also thought of as an economical way to stabilize the hierarchical basis (HB)

methods for d > 2, which otherwise are known to be suboptimal. Other stabil-

ization attempts for the HB methods can be found in [56,57]. We stress that the

smoothing in the nest grid, which follows from the principle of auxiliary space

method [62], plays an important role in the stabilization.

Third, we provide an analysis without the so-called nested re nement assumption:

=, (1.3)

J 1

J 0

with each subdomain j being made of all the new elements created at the j -th

level which, therefore, were not present earlier. This also implies that all elements

contained in subdomain j possess the same generation j and a comparable size.

The grid corresponding to j is thus quasi-uniform, and local multigrid methods

can be analyzed using a truncated L 2 -projection [14,15].

A natural question then arises: how can we apply existing theories to bisection grids

which may not obey the nested re nement assumption (1.3)? For fully additive

123

Optimal multilevel methods for graded bisection grids 5

multilevel methods, e.g. the original BPX preconditioner, the ordering does not

matter in the implementation and analysis and thus we can always assume (1.3)

holds. For multiplicative methods (e.g. V-cycle multigrid or additive multilevel

method with multiplicative smoothers), further work needs to be done to apply the

existing theories.

One approach is to use the relation between the additive and the multiplica-

tive method [16,23,29]. Roughly speaking, if additive preconditioners lead to a

uniformly bounded condition number for the preconditioned system, then multipli-

cative methods converge with a rate depending at most on the number of levels J .

For quasi-uniform grids, J log h log N which is an acceptable factor in

practice. For bisection grids, however, the level J could be O( N ) in the worst

scenario and thus this estimate is not optimal.

When the tree structure of the local mesh re nement is available, one could recon-

struct a virtual local mesh re nement hierarchy by grouping all elements with the

same generation into one level such that the assumption (1.3) holds [11,12,31] and

implement multigrid algorithms on this virtual nested re nement. However, these

levels increase dynamically within AFEM and must be updated for every loop

(1.1). Consequently, reconstructing a virtual re nement hierarchy entails imple-

menting suitable bookkeeping data structures which might compromise optimal

storage and thus optimality.

The new algorithms [19,22], as well as [1], show that multilevel methods retain

optimality even when the nested re nement assumption (1.3) is violated or ele-

ments with disparate sizes are grouped together into one re nement patch. This

paper provides an alternative approach to analyze these more exible algorithms.

Four, we provide a uni ed framework for analysis of multilevel methods on graded

bisection grids, which is valid for any dimension d, any polynomial degree, and

minimal regularity. We should point out that Wu and Chen [60] have analyzed

multigrid methods for bisection grids and d = 2 without the nested re nement

assumption (1.3). Their proof of uniform convergence relies on the speci c geo-

metric structure of bisection grids for d = 2, and its extension to d > 2 seems

rather dif cult. Our approach below is conceptually simpler than [60], applies to

any dimension d as well as general smoothers, rather than just Gauss-Seidel, and

extends to BPX preconditioners. Our analysis, carried out in Sect. 2.4, hinges on

three basic properties: the contraction property (2.6) of the local (inexact) smoother,

the stability bound (1.6) of the subspace decomposition, and the strengthened

Cauchy Schwarz inequality (1.7). The proofs of (1.6) and (1.7) are the core of this

paper and are given in Sect. 4. They heavily rely on the decomposition of bisection

grids discussed earlier in Sect. 3.

In the rest of the introduction, we brie y present the model problem and outline

our approach. Let Rd, d 2 be a polyhedral domain, and consider the Dirichlet

form

a (u, v) := u v dx .

123

6 L. Chen et al.

Given a (graded) triangulation T of, we choose the nite element space V := {v

H0 : v Pm, for all T }, where Pm is the space of polynomial of degree

1

m, with an integer m 1, and let u V be the nite element solution of the second

order elliptic equation

( Au, v) := a (u, v) = f, v for all v V, (1.4)

where f H 1, A : V V is the discrete Laplacian, and, is the duality

pair. For ease of exposition, we let be partitioned exactly into a bisection grid T,

which is shape regular and conforming, and consider the Laplace operator.

Our analysis can be generalized to variable coef cients with small variations, in

which case the contraction factor will depend on the variations of the coef cients.

For variable coef cients, possibly with large jumps across interelement boundaries,

we refer to [21]: using local multigrid as a preconditioner of the conjugate gradient

method yields a robust method with respect to both the mesh size and the size of jump

discontinuities of the coef cients.

We now brie y discuss our approach to multilevel methods. Let { p } p be the

canonical basis functions of the space V and let V p = span{ p } with dim V p = 1

for p ; thus V = span{ p } p . Let Vi V with dim Vi = 3 be the space of

piecewise linear functions spanned by the newest vertex added by each elementary

bisection bi and its two parents vertices, for i = 1, . . ., N, and let V0 be the coarsest

space of piecewise linear elements over T0 . We then have the space decomposition

N

V= Vp + Vi, (1.5)

p i =0

for which we shall prove the following two key properties:

Stable Decomposition: For any v V, there exist v p V p, p, and vi Vi, i =

0, . . ., N such that v = p v p + iN 0 vi and

=

N

h 2 v p h i 2 vi

+ v 2 .

2 2

(1.6)

p 1

p i =0

Strengthened Cauchy Schwarz (SCS) Inequality: For any u i, vi Vi, i =

0, . . ., N, we have

1/2 1/2

N N N N

a (u i, v j ) u i 2 vi 2 . (1.7)

1 1

i =0 j =i +1 i =0 i =1

Hereafter h p or h i represent local meshsizes corresponding to V p or Vi, respectively.

With the help of (1.6) and (1.7), derived in Sect. 4, we are able to obtain optimal

multilevel methods including BPX preconditioner and V-cycle multigrid methods for

123

Optimal multilevel methods for graded bisection grids 7

solving the algebraic system (1.4) over graded bisection grids. We prove convergence

of these methods in Sect. 2.4.

We use standard Sobolev space notation: denotes the L 2 -norm and 1 the

1 -semi-norm, which is a norm on H 1 . We write x y to indicate x C y, with

H 0

constant C independent of problem size N and functions v V, as well as x y to

mean x y and y x .

The rest of this paper is organized as follows. In Sect. 2, we review the subspace

correction method and provide abstract convergence analysis based on three assump-

tions: (1.6), (1.7), and (2.6) below. In Sect. 3, we discuss bisection methods and

present the crucial decomposition of bisection grids. In Sect. 4, we rst obtain a space

decomposition based on the decomposition of bisection grids and next prove (1.6)

and (1.7). Finally, in Sect. 5 we summarize optimal complexity results for both local

BPX-preconditioner and V-cycle multigrid, which are valid for inexact local solvers

which induce the contraction (2.6).

2 The method of subspace corrections

Discretization of partial differential equations often leads to linear algebraic equations

of the form

Au = f, (2.1)

where A R N N is a sparse matrix and f R N . In this section, we give some gen-

eral and basic results that will be used in later sections to construct ef cient multilevel

iterative methods (such as multigrid methods) for (2.1) resulting from nite element

discretizations of elliptic partial differential equations. The presentation in this section

follows closely to Xu [61] with simpli ed analysis.

2.1 Iterative methods

A basic linear iterative method for Au = f can be written in the following form

u k +1 = u k + B ( f Au k ),

starting from an initial guess u 0 V ; B is called iterator. If A = (ai j ) R N N is

split into diagonal, lower and upper triangular parts, namely A = D + L + U, then

two classical examples are the Jacobi method B = D 1 and the Gauss-Seidel method

B = ( D + L ) 1 .

The art of constructing ef cient iterative methods lies on the design of B which

captures the essential information of A 1 and its action is easily computable. In this

context the notion of ef ciency entails two essential requirements:

One iteration requires a computational effort proportional to the number of

unknowns.

The rate of convergence is well below 1 and independent of the number of unknowns.

123

8 L. Chen et al.

The approximate inverse B, when it is SPD, can be used as a preconditioner for

Conjugate Gradient (CG) method. The resulting method, known as preconditioned

conjugate gradient method (PCG), admits the following error estimate in terms of the

condition number ( B A) = max ( B A)/ min ( B A)

k

u uk ( B A) 1

A

2 (k 1);

u u0 ( B A) + 1

A

B is called preconditioner. A good preconditioner should have the properties that the

action of B is easy to compute and that ( B A) is significantly smaller than ( A).

2.2 Space decomposition and method of subspace corrections

In the spirit of divide and conquer, we decompose the space V = iJ=0 Vi as the

summation of subspaces Vi V ; {Vi }iJ=0 is called a space decomposition of V .

Since iJ=0 Vi is not necessarily a direct sum, decompositions of u V of the form

u = iJ=0 u i are in general not unique. The original problem (2.1) can thus be split

into sub-problems in each Vi with smaller size which are relatively easier to solve.

Throughout this paper, we use the following operators, for i = 0, 1, . . ., J :

Q i : V Vi the projection in the inner product ;

Ii : Vi V the natural inclusion which is often called prolongation;

Pi : V Vi the projection in the inner product A = ( A, );

Ai : Vi Vi the restriction of A to the subspace Vi ;

Ri : Vi Vi an approximation of Ai 1 (often known as smoother);

Ti : V Vi Ti = Ri Q i A = Ri Ai Pi .

It is easy to verify the relation Q i A = Ai Pi and Q i = Iit with ( Iit u, vi ) := (u, Ii vi ).

The operator Iit is often called restriction. If Ri = Ai 1, then we have an exact local

solver and Ti = Pi . With slightly abused notation, we still use Ti to denote the restric-

tion Ti Vi : Vi Vi and Ti 1 = (Ti Vi ) 1 : Vi Vi .

For a given residual r V, we let ri = Q i r = Iit r denote the restriction of the

residual to the subspace Vi and solve the residual equation Ai ei = ri in Vi approxi-

mately

ei = Ri ri .

Subspace corrections ei are assembled to yield a correction in the space V, thereby

giving rise to the so-called method of subspace corrections. There are two basic ways

to assemble subspace corrections.

Parallel Subspace Correction (PSC) This method performs the correction on each

subspace in parallel. In operator form, it reads

u k +1 = u k + B ( f Au k ), (2.2)

123

Optimal multilevel methods for graded bisection grids 9

where

J

B= Ii Ri Iit . (2.3)

i =0

The subspace correction is ei = Ii Ri Iit ( f Au k ), and the correction in V is e =

J

ei . The error equation reads

i =0

J J

u u k +1 = I A (u u k ) = I Ti (u u k );

Ii Ri Iit

i =0 i =0

Successive Subspace Correction (SSC) This method performs the correction in a suc-

cessive way. In operator form, it reads

v 0 = u k, v i +1 = v i + Ii Ri Iit ( f Av i ), i = 0, . . ., J, u k +1 = v J +1, (2.4)

and the corresponding error equation is

J J

u u k +1 = ( I Ii Ri Iit A) (u u k ) = ( I Ti ) (u u k );

i =0 i =0

in the notation iJ=0 ai, we assume there is a built-in ordering from i = 0 to J,

i.e., iJ=0 ai = a0 a1 . . . a J . Therefore, PSC is an additive method whereas SSC is a

multiplicative method.

As a trivial example, we consider the space decomposition R J = iJ=1 span{ei }.

In this case, if we use exact (one dimensional) subspace solvers, the resulting SSC is

just the Gauss-Seidel method and the PSC is just the Jacobi method. More complicated

and effective examples, including multigrid methods and multilevel preconditioners,

will be discussed later on.

2.3 Sharp convergence identities

The analysis of parallel subspace correction methods relies on the following identity

which is well known in the literature [28,59,61,64].

Theorem 2.1 (Identity for PSC) If Ri is SPD on Vi for i = 0, . . ., J, then B de ned

by (2.3) is also SPD on V . Furthermore

J

( Ri 1 vi, vi ).

( B 1 v, v) = inf (2.5)

J

i =0 vi =v i =0

On the other hand, the analysis of Successive subspace correction methods hinges

on an identity of Xu and Zikatanov [64] to be described below. First we assume that

123

10 L. Chen et al.

each subspace smoother Ri induces a convergent iteration, i.e. the error operator I Ti

is a contraction.

(T) Contraction of Subspace Error Operator There exists i

When we choose exact local solvers, i.e., Ri = Ai 1 and consequently Ti = Pi for

i = 0, . . ., J, (T) holds with = 0. Therefore we have a more concise formulation

for such choice [64].

Corollary 2.3 (Identity of SSC with exact solver) One has the following identity

2

J

1

( I Pi ) =1,

1 + c0

i =0 A

where

2

J N

c0 = sup vj .

inf Pi (2.9)

J

v A =1 i =0 vi =v i =0 j =i +1 A

123

Optimal multilevel methods for graded bisection grids 11

2.4 Convergence analysis

We now present a convergence analysis based on three assumptions: (T) on Ti and

the two ones below on the space decomposition. The analysis here is adapted from

Xu [61] and simpli ed by using the XZ identity.

(A1) Stable Decomposition For any v V, there exists a decomposition v =

J

i =0 vi, vi Vi, i = 0, . . ., J such that

J

vi K1 v A.

2 2

(2.10)

A

i =0

(A2) Strengthened Cauchy Schwarz (SCS) Inequality For any u i, vi Vi, i = 0, . . ., J

1/2 1/2

J J J J

(u i, v j ) A K 2 vi 2 .

ui 2 (2.11)

A A

i =0 j =i +1 i =0 i =0

Theorem 2.4 (Multilevel preconditioning) Let V = iJ=0 Vi be a space decomposi-

tion satisfying assumptions (A1) and (A2), and let Ri be SPDs for i = 0, . . ., J such

that

( Ri 1 u i, u i ) K 3 u i

K4 1 ui A.

2 2

(2.12)

A

Then B de ned by (2.3) is SPD and

( B A) (1 + 2 K 2 ) K 1 K 3 K 4 . (2.13)

J

Proof Let v = i =0 vi be a decomposition satisfying (2.10). It follows from the

identity (2.5), and the definitions (2.10) of K 1 and (2.12) of K 3, that

J J

( Ri 1 vi, vi ) K 3

( B 1 v, v) vi K1 K3 v = K 1 K 3 ( Av, v),

2 2

A A

i =0 i =0

which implies

min ( B A) ( K 1 K 3 ) 1 . (2.14)

123

12 L. Chen et al.

J

For any decomposition v = i =0 vi, in view of (2.11) and (2.12), we have

J J J

( Av, v) (vi, vi ) A + 2 (vi, v j ) A

i =0 i =0 j =i +1

J J

( Ri 1 vi, vi ).

(1 + 2 K 2 ) vi (1 + 2 K 2 ) K 4

2

A

i =0 i =0

Taking the infimum and using again (2.5), we get

J

( Ri 1 vi, vi ) = (1 + 2 K 2 ) K 4 ( B 1 v, v),

( Av, v) (1 + 2 K 2 ) K 4 inf

J

k =0 vi =v i =0

which implies

max ( B A) (1 + 2 K 2 ) K 4 . (2.15)

The estimate (2.13) then follows from (2.14) and (2.15).

1

Lemma 2.5 (Estimate of eigenvalues of T i ) If (T) holds, then T i is non-singular

and

1

1 1

1 min (T i ) max (T i ) . (2.16)

1 2

Proof The estimates follow easily by (2.6) and the definition of T i, which satis es

I T i = ( I Ti )( I Ti ) = ( I Ti ) ( I Ti ). (2.17)

We omit the details.

Theorem 2.6 (Convergence of SSC) Let V = iJ=0 Vi be a space decomposition

satisfying assumptions (A1) and (A2), and let the subspace smoothers Ti satisfy (T).

We then have

2

J

1 2

( I Ti ) 1 .

2 K 1 (1 + (1 + )2 K 2 )

2

i =0 A

Proof We shall give an upper bound of the constant K in Theorem 2.2 by choosing a

stable decomposition v = i vi satisfying (2.10). By the inequality 2ab a 2 + b2,

we have

1 1 1

2(T i vi, Ti wi ) A (T i vi, vi ) A + (T i Ti wi, Ti wi ) A, (2.18)

123

Optimal multilevel methods for graded bisection grids 13

where w j = j >i v j . Therefore we only need to estimate the two terms on the right

hand side of (2.18).

1

Using the eigenvalue estimate of T i in Lemma 2.5, together with (2.10), we arrive

at

J J

K1

1 1

(T i vi, vi ) A max (T i ) vi v A.

2 2

(2.19)

A

1 2

i =0 i =0

To estimate the second term in (2.18), we use the SCS estimate (2.11) to get

J J J

1 1

(T i Ti wi, Ti wi ) A = (Ti T i Ti wi, v j ) A

i =0 i =0 j =i +1

1/2 1/2

J J

1

Ti T i Ti wi 2

K2 vi 2 .

A A

i =0 i =0

1 1

= Ti 1 + and T i = max (T i ) 1

Since Ti, we deduce

Ai Ai Ai 1 2

1+

1 1

Ti T i Ti wi Ti

Ti wi wi A.

Ti

A Ai Ai Ai A

1

Now using the SCS estimate (2.11) again, we get

1/2 1/2

J J J J J

(wi, wi ) A = (wi, v j ) A K 2 wi 2 vi 2,

A A

i =0 i =0 j =i +1 i =0 i =0

which leads to

J J

wi K2 vi A.

2 2 2

A

i =0 i =0

Consequently,

J J

1+ 1+

1

(T i Ti wi, Ti wi ) A K 2 vi K1 K2 v A.

2 2 2 2

(2.20)

A

1 1

i =0 i =0

Inserting (2.18), (2.19), and (2.18) into (2.8), we get the upper bound of K

2 K 1 1 + (1 + )2 K 2 /(1 2 ). Finally, the desired contraction estimate follows

from Theorem 2.2.

When we use exact local solvers Ri = Ai 1, we have a simpler proof and sharper

estimate.

123

14 L. Chen et al.

Corollary 2.7 Let the space decomposition satisfy (A1) and (A2), and let Ri = Ai 1

for all i . Then

2

J

1

( I Pi ) 1 .

1 + K1 K2

2

i =0 A

J

Proof We apply (2.11) with u i = Pi j =i +1 v j to obtain

J J J J J

u i, Pi vj =

= (u i, v j ) A

2

ui A

i =0 i =0 j =i +1 i =0 j =i +1

A

1/2 1/2

J J

K2 vi 2 .

ui 2

A A

i =0 i =0

J

Consequently, if v = k =0 vk is a stable decomposition satisfying (2.10), we get

2

J J J J

vj = K2 vi K1 K2 v A,

2 2 2 2 2

Pi ui A A

i =0 j =i +1 i =0 i =0

A

which implies c0 K 2 K 1 . The desired result then follows from Corollary 2.3.

2

3 Bisection methods

We now discuss bisection methods for simplicial grids for d 2, following [42,53],

and present a novel decomposition of conforming meshes obtained by bisection. We

do not discuss the alternative re nement method, called regular re nement, which

divides one simplex into 2d children; see [7,24] for d = 2 and [1,11] for d = 3.

3.1 Bisection rules

Given a simplex, we assign one of its edges as the re nement edge of . Starting

from an initial triangulation T0, a bisection method consists of the following rules:

R1. Assign re nement edges for each element T0 (initial labeling);

R2. Divide a simplex into two simplices by joining the midpoint of its re nement

edge with its vertices other than those in the re nement edge (bisection);

R3. Assign re nement edges to the two children of a bisected simplex (labeling).

There are several bisection methods proposed for d 3 [3,8,32,33,46,54], which

generalize the newest vertex bisection [37] and longest edge bisection [48] for d = 2.

We now give a mathematical description based on Kossaczky [32], Traxler [55], and

Stevenson [54]. For each simplex, rules R1-3 associate a unique re nement edge e.

123

Optimal multilevel methods for graded bisection grids 15

The pair (, e) is called labeled simplex, and (T, L) := {(, e) : T } is called

labeled triangulation. A d -simplex is a set of d + 1 ordered vertices {xi }i =0 and

d

type t :

= { x 0, x 1, . . ., x d }, t {0, 1, . . ., d }.

We let e = x0 xd be the re nement edge, and let x = 2 (x0 + xd ) be the midpoint of e.

1

The two children 1, 2 of are the simplices obtained by joining x with the vertices

of other than x0, xd . Ordering the vertices of the children, or equivalently labeling

them, is a crucial process that includes R2-3. We consider the following rule R3:

1 := {x0, x, x1, . . ., xt, xt +1, . . ., xd 1 }(t +1)mod d,

(3.1)

2 := {xd, x, x1, . . ., xt, xd 1, . . ., xt +1 }(t +1)mod d,

with the convention that arrows point in the direction of increasing indices and

{x1, . . ., x0 } =, {xd, . . ., xd 1 } = . For d = 2 rule R3 does not depend on the

element type and we get for = {x0, x1, x2 } the two children 1 = {x0, x, x1 }, 2 =

{x2, x, x1 }. Moreover, the re nement edge of the two children is opposite to the new

vertex x, whence this procedure coincides with the newest vertex bisection method.

We refer to the survey [42, Sect. 4] for a discussion for d 2, and stress that once

rule R1 is settled, then the subsequent labeled grids are uniquely de ned.

For a labeled triangulation (T, L), and T, a bisection b : {(, e)}

{( 1, e1 ), ( 2, e2 )} is a map that encodes the above procedure. We next de ne the

formal addition as follows:

T + b := (T, L, e)} {( 1, e1 ), ( 2, e2 )}.

For an ordered sequence of bisections B = (b 1, b 2, . . ., b N ), we de ne

T + B := ((T + b 1 ) + b 2 ) + + b N,

whenever the addition is well de ned (i.e. i should exists in the previous labeled

triangulation). These additions are a convenient mathematical description of bisection

on triangulations.

Given an initial grid T0 of and rules R1-3, we de ne the sets

G(T0 ) = {T : there exists a bisection sequence B such that T = T0 + B },

T(T0 ) = {T G(T0 ) : T is conforming}.

Therefore G(T0 ) contains all (possibly nonconforming) grids obtained from T0 using

the bisection method, which are uniquely de ned once the rules R1-3 have been set,

whereas T(T0 ) is the subset of conforming grids.

123

16 L. Chen et al.

It is essential for the discussion to de ne the sequence of uniformly re ned meshes

{T k } 0 by:

k=

T k := T k 1 + {b : T k 1 }, for k 1,

with T 0 := T0 . This means that T k is obtained by bisecting all elements in T k 1 only

once. Note that T k G(T0 ) but not necessarily in T(T0 ).

We thus consider bisection methods which satisfy the following two assumptions:

( B 1) Shape Regularity: G(T0 ) is shape regular.

( B 2) Conformity of Uniform Re nement: T k T(T0 ), i.e. T k is conforming, for all

k 0.

With the speci c rule (3.1), due to [32,54,55], we see that the type t increases by 1

and the vertex ordering changes with t, which in turn implies that after d recurrent

bisections of a simplex all its edges are bisected. This leads to a xed number of

similarity classes of elements, depending only on T0, and thus B1 holds for any d . We

refer to [42] for a thorough discussion.

We recall that for d = 2, rule (3.1) reduces to the newest vertex bisection, in which

case Sewell [51] showed that all the descendants of a triangle in T0 fall into four

similarity classes and hence (B1) holds. Note that (B2) may not hold for an arbitrary

rule R1, namely the re nement edge for elements in the initial triangulation cannot be

selected freely. Mitchell [37] came up with a rule R1 for which (B2) holds. He proved

the existence of such initial labeling scheme (so-called compatible initial labeling),

and Biedl et al. [9] gave an optimal O( N ) algorithm to nd a compatible initial labeling

for a triangulation with N elements. In summary, for d = 2, newest vertex bisection

with compatible initial labeling is a bisection method which satis es (B1) and (B2).

Enforcing (B2) for d > 2 requires also a labeling of the initial mesh T0, for which

there is no constructive procedure. The algorithms proposed by Kossaczk [32] for

d = 3 and Stevenson [54] for d 3 enforce such initial labeling upon further

re ning every element of the initial triangulation, which deteriorates the shape regu-

larity. Although (B2) imposes a severe restriction on the initial labeling, we emphasize

that it is also used to prove the optimal cardinality of adaptive nite element meth-

ods [18,42,53]. Finding conditions weaker than (B2) is a challenging open problem.

3.2 Compatible bisections

We denote by N (T ) the set of vertices of the triangulation T and by E (T ) the set of

all edges of T . By convention, all simplices are closed sets. For a vertex x N (T )

or an edge e E (T ), we de ne the rst ring (or the star) of x or e to be

Rx = { T x }, Re = { T e },

and the local patch of x or e as x = Rx, and e = Re . Note that x and

e are subsets of, while Rx and Re are subsets of T which can be thought of as

triangulations of x and e, respectively. We indicate with # S the cardinality of a set S .

Given a labeled triangulation (T, L), an edge e E (T ) is called a compatible edge

if e is the re nement edge of for all Re . For a compatible edge e, the ring Re

123

Optimal multilevel methods for graded bisection grids 17

Fig. 1 Two compatible bisections for d = 2. Left interior edge; right boundary edge. The edge with

boldface is the compatible re nement edge, and the dash-line represents the bisection

(b)

(a)

Fig. 2 A compatible bisection for d = 3: the edge e (in bold) is the re nement edge of all elements in the

patch e . Connecting the midpoint p of e to the other vertices bisects each element of the compatible ring

Re and keeps the mesh conforming without spreading re nement outside e . This is an atomic operation

is called a compatible ring, and the patch e is called a compatible patch. Let x be

the midpoint of a compatible edge e and Rx be the ring of x in T + {b : Re }.

A compatible bisection is a mapping be : Re Rx . We then de ne the addition

T + be := T + {b : Re } = T \Re Rx .

Note that if T is conforming, then T + be is conforming for a compatible bisection

be, whence compatible bisections preserve conformity of triangulations and are thus a

fundamental concept both in theory and practice. For a compatible bisection sequence

B = (bi )iJ=0, the addition T0 + B is de ned recursively Ti = Ti 1 + bi for 1 i J .

In two dimensions, a compatible bisection be has only two possible con gurations;

see Fig. 1. The rst one corresponds to bisecting an interior compatible edge, in which

case the patch e is a quadrilateral. The second case corresponds to bisecting a bound-

ary edge, which is always compatible, and e is a triangle. In three dimensions, the

con guration of compatible bisections depends on the initial labeling; see Fig. 2 for a

simple case.

The bisection of paired triangles was rst introduced by Mitchell for dimension d =

2 [37,38]. The idea was generalized by Kossaczk [32] to d = 3, and Maubach [33] and

Stevenson [54] to d 2. In the aforementioned references, ef cient recursive comple-

tion procedures of bisection methods are introduced based on compatible bisections.

We use them to characterize the conforming mesh obtained by bisection methods.

3.3 Decomposition of bisection grids

We now present a decomposition of meshes in T(T0 ) using compatible bisections,

which will be instrumental later.

123

18 L. Chen et al.

(a) (b) (c)

(d) (e) (f)

Fig. 3 Decomposition of a bisection grid for d = 2: Each frame displays a mesh Ti +k = Ti +

(bi +1, . . ., bi +k ) obtained from Ti by a sequence of compatible bisections (b j )ij+k +1 using the longest

Contact this candidate