TR2005-873 A BDDC ALGORITHM FOR PROBLEMS WITH MORTAR ...

Report 3 Downloads 25 Views
TR2005-873

A BDDC ALGORITHM FOR PROBLEMS WITH MORTAR DISCRETIZATION ∗ HYEA HYUN KIM †, MAKSYMILIAN DRYJA ‡, AND OLOF B. WIDLUND

§

September 4, 2005 Abstract. A BDDC (balancing domain decomposition by constraints) algorithm is developed for elliptic problems with mortar discretizations for geometrically non-conforming partitions in both two and three spatial dimensions. The coarse component of the preconditioner is defined in terms of one mortar constraint for each edge/face which is an intersection of the boundaries of a pair of subdomains. A condition number bound of the ˘ ¯ form C maxi (1 + log(Hi /hi ))3 is established. In geometrically conforming cases, the bound can be improved ˘ ¯ to C maxi (1 + log(Hi /hi ))2 . This estimate is also valid in the geometrically nonconforming case under an additional assumption on the ratio of mesh sizes and jumps of the coefficients. This BDDC preconditioner is also shown to be closely related to the Neumann-Dirichlet preconditioner for the FETI–DP algorithms of [9, 11] and it is shown that the eigenvalues of the BDDC and FETI–DP methods are the same except possibly for an eigenvalue equal to 1. Key words. BDDC, FETI–DP, mortar methods, preconditioner AMS subject classifications. 65N30, 65N55

1. Introduction. This study focuses on a scalable BDDC algorithm for solving linear systems arising from mortar finite element discretizations of elliptic problems. A BDDC method was first introduced by Dohrmann [4] as an improvement of the balancing NeumannNeumann method and using different coarse finite element spaces. The coarse space consists of a weighted sum of functions each of which minimizes the local discrete energy norm with certain constraints on the subdomain interfaces; continuity of the solutions at vertices, or average or momentum matching condition on solutions over edges/faces are considered in [4, 16, 17, 19, 20]. The resulting coarse problem then gives a more local coupling between the subdomains than for the older balancing methods and more freedom in choosing the constraints to improve the convergence. An additional advantage is that all linear systems will have positive definite, symmetric matrices at least for conforming finite element problems. The constraints on the coarse finite element space are basically the same as those of a ∗ The first author

was supported in part by the Applied Mathematical Sciences Program of the U.S. Department of

Energy under contract DE-FG02-00ER25053 and in part by the Post-doctoral Fellowship Program of Korea Science and Engineering Foundation (KOSEF). The second author was supported in part by Polish Sciences Foundation under grant 2P03A00524. The third author was supported in part by the U.S. Department of Energy under contract DE-FC02-01ER25482. † Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, N.Y. 10012. U.S.A. Electronic mail address: [email protected]. URL: http://cims.nyu.edu/˜hhk2. ‡ Department of Mathematics, Warsaw University, Banach 2, 02-097 Warsaw, Poland. Electronic mail address: [email protected]. § Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, N.Y. 10012. U.S.A. Electronic mail address: [email protected]. URL: http://cs.nyu.edu/cs/faculty/widlund/index.html. 1

2

KIM, DRYJA AND WIDLUND

FETI–DP algorithm. In a FETI–DP algorithm, a linear system formulated for a set of dual variables is solved after eliminating the primal unknowns related to the primal constraints, given by average matching condition over edges/faces or continuity of the solutions at vertices. The resulting linear system, in itself, contains a coarse problem while its preconditioner is built only from subdomain problems. In a BDDC method, a linear system of the primal variables is solved iteratively with a preconditioner that has both coarse and subdomain components. This provides BDDC methods with more flexibility, allowing for the use of inexact coarse problems. Thus, an inexact coarse problem can be introduced by applying the BDDC method recursively to the coarse problem; see Tu [22, 23]. The use of inexact local problems for the BDDC preconditioners has also been considered by Li and Widlund [18]. Recently the BDDC methods have been shown to be closely related to the FETI–DP methods. A condition number bound of the BDDC operator was first given by Mandel and Dohrmann [19]. They proved a C(1 + log(H/h))2 bound that is comparable to that for the FETI–DP methods. Further, Mandel, Dohrmann, and Tezaur [20] showed that the eigenvalues of the FETI–DP and BDDC operators are the same except possibly for eigenvalues equal to 0 and 1. Recently, a new formulation of the BDDC method was given by Li and Widlund [17]. They introduced a change of variables as well as an average operator for the BDDC method based on the jump operator used in [15] in the analysis of FETI–DP methods. The change of variables greatly simplifies the analysis; it has also led to a successful and robust implementation of FETI–DP algorithms [12, 13]. In this paper, we will first describe a BDDC algorithm with a mortar discretization and a change of variables. Primal constraints on edges/faces are introduced. We consider quite general geometrically non-conforming partitions and the second generation of the mortar method as well as the dual basis mortar methods. A preconditioner is then proposed which uses a cer tain weight matrix D, that leads to the condition number bound: C maxi (1 + log(Hi /hi ))3 . Section 4 is devoted to proving the condition number bound in terms of a bound of an average

operator ED in a certain norm. The algorithm can also be applied to a geometrically conform ing partition and then gives a better bound: C maxi (1 + log(Hi /hi ))2 . The same bound can be established for geometrically non-conforming partitions with an additional assumption on the mesh sizes. In Section 5, we show that the preconditioner proposed for our BDDC algorithm is closely connected to the Neumann-Dirichlet preconditioner of the FETI–DP algorithm given in [9, 11]. By establishing connections between the average and jump operators, the spectra of the BDDC and FETI–DP algorithms are then shown to be the same except possibly for an eigenvalue equal to 1. This approach was used by Li and Widlund [17] and provided a simpler proof for the condition number bound of the BDDC algorithm. Our BDDC algorithm is also applicable to elasticity problems employing a preconditioner that is closely connected to the Neumann-Dirichlet preconditioner of the FETI–DP formulation developed in [10].

A BDDC METHOD WITH MORTAR DISCRETIZATION

3

In the final section, numerical results show that the FETI–DP and BDDC algorithms perform very similarly when the same set of primal constraints are selected. Throughout this paper, C denotes a generic constant that does not depend on any mesh parameters and coefficients of the elliptic problems. We note that this paper originated from two projects developed separately by the first and the second authors; the contribution of the third author began with a suggestion that a theory could be developed for the geometrically non-conforming case. 2. Finite element spaces and mortar matching constraints. 2.1. A model problem and mortar methods. We consider a model elliptic problem in a polygonal/polyhedral domain Ω ⊂ R2 (R3 ): find u ∈ H01 (Ω) such that (2.1)

−∇ · (ρ(x)∇u) = f (x) u=0

∀x ∈ Ω,

on ∂Ω,

where ρ(x) ≥ ρ0 > 0 and f (x) ∈ L2 (Ω). Let Ω be partitioned into disjoint polygonal/polyhedral subdomains Ω=

N [

Ωi .

i=1

We assume that the partition can be geometrically non-conforming, see discussion below, and that ρ(x) = ρi , x ∈ Ωi for some positive constant ρi . We denote by Xi the P1 -conforming finite element space on a quasi-uniform triangulation Ti of each subdomain Ωi . The Ti might not align across subdomain interfaces. The space Wi is the trace space of Xi on ∂Ωi . We then introduce the product spaces X :=

N Y

Xi ,

W :=

N Y

Wi .

i=1

i=1

For functions in these spaces, we will impose the mortar matching condition across the interfaces using suitable Lagrange multiplier spaces. In a geometrically non-conforming partition, the intersection of the boundaries of neighboring subdomains can be only a part of a edge/face of a subdomain. Let us define the entire interface by 

Γ=

[ ij



∂Ωi ∩ ∂Ωj  \ ∂Ω.

Among the subdomain edges/faces, we select nonmortar edges/faces Fl for which [ l

F l = Γ,

Fl ∩ Fk = ∅, l 6= k.

4

KIM, DRYJA AND WIDLUND

Since the subdomain partition can be geometrically non-conforming, a single nonmortar edge/face Fl ⊂ ∂Ωi may intersect several subdomain boundaries ∂Ωj . This provides Fl with a partition Fl =

[

F ij ,

Fij = ∂Ωi ∩ ∂Ωj .

j

A dual or a standard Lagrange multiplier space Ml is given for each nonmortar edge/face Fl . ◦

We require that the space Ml has the same dimension as the space W (Fl ) := Wi |Fl ∩ H01 (Fl ) and that it contains the constant functions. Constructions of such Lagrange multiplier spaces were introduced in [1, 3] for standard Lagrange multiplier spaces and in [24, 25] for dual Lagrange multiplier spaces; see also [8]. For (w1 , · · · , wN ) ∈ W , we define φ ∈ L2 (Fl ) by φ = wj on Fij ⊂ Fl . The mortar matching condition for the geometrically non-conforming partition is given by Z (wi − φ)λ ds = 0, ∀λ ∈ Ml , ∀Fl . (2.2) Fl

We write its matrix representation as (2.3)

N X

B (i) w(i) = 0,

i=1

with w(i) a vector representation of wi using nodal basis functions. We further define the ◦

following product spaces by gathering the spaces Ml and W (Fl ) given on each nonmortar edges/faces: Y Y ◦ (2.4) M= Ml , Wn = W (Fl ). The mortar finite element method for problem (2.1) is to approximate the solution by

Galerkin’s method in the mortar finite element space b := {v ∈ X : v satisfies the mortar matching condition (2.2)} . X

2.2. Finite element spaces with a change of variables. In this subsection, we introduce a change of variables for the unknowns in the space W . This change of variables is based on the primal constraints that will be imposed in our BDDC algorithm. In mortar discretizations, we may consider the following sets of primal constraints; vertex constraints, vertex and edge average constraints, or edge average constraints only for two spatial dimensions, and vertex constraints and face average constraints, or face average constraints only for three spatial dimensions. We note that vertex constraints are appropriate only for the first generation of the mortar method. In order to reduce the number of primal constraints, we can select only some edges/faces or some vertices as primal where the primal constraints will be imposed. Such choices have been considered for the FETI–DP algorithms and conforming finite elements in [14] and for mortar finite elements in [10].

A BDDC METHOD WITH MORTAR DISCRETIZATION

5

In our BDDC formulation, we will introduce certain primal constraints over edges/faces that are selected from the mortar matching constraints (2.2). We consider {ψij,k }k , the basis functions in Ml that are supported in F ij , and define ψij =

X

ψij,k .

k

We assume that at least one such basis function ψij,k exists for each Fij ⊂ Fl . We now introduce the following primal constraints for (w1 , · · · , wN ) ∈ W over each edge/face Fij Z

(2.5)

(wi − wj )ψij ds = 0,

Fij

and define (2.6)

f = {w ∈ W : w satisfies the primal constraints (2.5)} . W

c ⊂ W f ⊂ W , where W c is the restriction of X b to Γ. For the case of a geoNote that W

metrically conforming partition, i.e., when Fij is a full face of two subdomains, the above constraints are edge/face average matching condition because ψij = 1. In addition to the

above constraints, vertex constraints can be considered for the first generation mortars if the partition is geometrically conforming. We now introduce a change of variables, following Li and Widlund [17], based on the primal constraints and in the two dimensional case. This approach can also be extended to the three dimensional case without any difficulty. We recall that F ⊂ ∂Ωi is a nonmortar edge/face and that {Fij }j is a partition of F given by Fij = F ∩ ∂Ωj , a mortar edge/face of Ωj . We denote by {zk }lk=1 the unknowns of wi ∈ Wi at the nodes in F related to the Lagrange multipliers {ψij,k } and by {vk }pk=1 the unknowns at the remaining nodes in F . We will now define a transform that retains the unknowns {vk }pk=1 and changes {zk }lk=1 into {b zk }lk=1 so that for a fixed m, chosen arbitrarily, zbm satisfies

R

F zbm = Rij

Let R

Fij

e hk = R

φek ψij ds

Fij

ψij ds

wi ψij ds

Fij

,

ψij ds

R

.

F hk = Rij

φk ψij ds

Fij

ψij ds

,

where φek and φk are the nodal basis functions of the unknowns vk and zk , respectively. For a simpler presentation, we assume that p = 2 but the following can be generalized to any p.

6

KIM, DRYJA AND WIDLUND

We then consider the following transform TFij : 

v1









0

0

···

0

0

0

···

1

0

···

0

0

0

···

0 .. .

1 .. .

··· .. .

0 .. .

A .. .

0 .. .

··· .. .

0

0

···

1

A

0

···

c2

r1

···

rm−1

0 .. .

0 .. .

··· .. .

0 0 zbl     0 v1     0  v2      1      zb1  .  .  .  .  .  .       ..   zb .  m−1      .   . = A  .  zbm +  zb0  ,   . zbm+1  .   .  ..     .   ..    .  .     ..  .    ..   .     ..  .  ..      zbl 1

0

···

v1

1

      v2   0  v2        zb   0  z   1    1   .  .  .   .  .  .   .  .  .       zm−1  = TFij zbm−1  =  0            zbm  c1  zm          zb z  m+1   0  m+1   .  .  .   .  .  .   .  .  .  zl

where

A rm+1

···

0 .. .

A .. .

1 .. .

··· .. .

0

A

0

···

0



v1



  0   v2      0   zb1    ..   ..   .  .      0  zbm−1    rl   zbm      0  zbm+1    ..   ..   .  .  1

zbl

zb0 = c1 v1 + c2 v2 + r1 zb1 + · · · + rm−1 zbm−1 + rm+1 zbm+1 + · · · + rl zbl , R

F

A = Pijl

ψij ds

k=1

hk

, c1 = −

e e h1 h2 hk , c2 = − , rk = − , k 6= m. hm hm hm

We see that this transform satisfies the requirements stated above. The transform TFij can be applied to each face Fij ⊂ F independently. For the case when an edge F ⊂ ∂Ωi is a mortar edge, there exists a Ωj across the interface with a nonmortar side. We then consider Fij = F ∩ ∂Ωj . In this case, the unknowns {zk }lk=1 are related to the nodes in F with its basis functions supported in F ij and the remaining unknowns in F are denoted by {vk }pk=1 . The transform TFij is then defined for these unknowns as before.

7

A BDDC METHOD WITH MORTAR DISCRETIZATION

ci → Wi By gathering the transforms TFij of all F ⊂ ∂Ωi , we get a transform T (i) : W

of the form

(i)

T

(i)

=

(i)

Trr

Trc

0

I

!

,

where c and r stand for the unknowns retained by the transform and the remaining unknowns, ci denotes the space of new unknowns. With this set of new unknowns, the respectively, and W

local stiffness matrix, the mortar matching matrix, and the local force vector are written as t Sb(i) = T (i) S (i) T (i) ,

b (i) = B (i) T (i) , B

t

gb(i) = T (i) g (i) .

The unknowns zbm , the unknowns for the averages over the edges, are the primal varif in (2.6) can be represented as ables. With this set of new variables, the space W f = W∆ ⊕ WΠ , W

(2.7)

where W∆ consists of functions with a zero value at the primal variables and WΠ consists (i)

of functions with a zero value at the other variables. We denote by RΠ the restriction of the primal unknowns uΠ ∈ WΠ to the subdomain Ωi . By using the set of new unknowns, the mortar matching condition (2.3) can be written as B∆ w∆ + BΠ wΠ = 0.

(2.8) Here

  (1) (N ) B∆ = B∆ , · · · B∆ (i)

BΠ =

N X

(i)

(i)

BΠ RΠ ,

i=1

(i)

b (i) with columns corresponding to the primal variwhere BΠ and B∆ are submatrices of B ables and the remaining unknowns, respectively.

f will be imposed by using Furthermore the mortar matching condition on functions in W

non-redundant Lagrange multipliers. We select the non-redundant Lagrange multipliers as

follows. From the bases of Ml , we eliminate one basis element among {ψij,k }k for each

Fij ⊂ Fl and denote the reduced Lagrange multiplier space by M l . The non-redundant Lagrange multiplier space is then defined as M=

Y

M l.

l

f by using the nonThe mortar matching condition (2.2) is imposed on the space W

redundant Lagrange multipliers λ ∈ M . To simplify the notation, we use the same notation as in (2.8) for this case, i.e.,

B∆ w∆ + BΠ wΠ = 0.

8

KIM, DRYJA AND WIDLUND

The space W∆ can be split into W∆ = W∆,n ⊕ W∆,m , where n and m denote unknowns at nonmortar edges/faces (interior) and the remaining unknowns, respectively. The above equation can be written as Bn wn + Bm wm + BΠ wΠ = 0.

(2.9)

After the change of variables, we order the local Schur complement matrix and the local Schur complement vector into Sb(i) =

Π∆

and define a matrix and vectors by

(2.10)

Se =

S∆∆

S∆Π

SΠ∆

SΠΠ

(i) Sb∆∆ (i) Sb

!

,

where

(2.11)

(i) Sb∆Π (i) Sb ΠΠ

N X i=1

(i)

,

gb(i) =

 (1) g∆ b  .  .  g∆ =   . , (N ) g∆ b 

  b(i) S∆∆ = diagN i=1 S∆∆ ,  (1) t b(1) SΠ∆ = RΠ SΠ∆ · · · SΠΠ =

!

t

(i) (i) (i) RΠ SbΠΠ RΠ .

(N ) t (N ) RΠ SbΠ∆

gb∆

(i)

gbΠ

gΠ =

!

N X i=1



,

(i) t (i)

RΠ gbΠ ,

t S∆Π = SΠ∆ ,

3. A BDDC algorithm for the mortar discretizations. In this section, we formulate a BDDC operator for the elliptic problem described in Section 2.1. We consider the same finite element space and subdomain partition as in Section 2.1 and use the unknowns after the change of variables introduced in Section 2.2. We will omit the hats for the transformed matrices to simplify the notation. We recall the mortar matching condition (2.9). Since the matrix Bn is invertible, we solve (2.9) for wn wn = −Bn−1 (Bm wm + BΠ wΠ ). We then define the matrix

(3.1)



 RΓ =  

−Bn−1 Bm

−Bn−1 BΠ

I

0

0

I



 , 

A BDDC METHOD WITH MORTAR DISCRETIZATION

9

t t t t t t which maps (wm , wΠ ) into a vector (wnt , wm , wΠ ) that satisfies the mortar matching con-

dition (2.9). Let us define the mortar finite element space by n o c= w∈W f : (wn , wm , wΠ ) satisfies (2.9) . W

In the BDDC method, we approximate the solution of the elliptic problem in the mortar finite c and obtain the following discrete problem: element space W ! ! gm wm t t e , = RΓ (3.2) RΓ SRΓ gΠ wΠ where gm is the component of the vector g∆ in (2.10) other than the nonmortar part.

We now introduce a coarse finite element space based on the primal constraints so as to solve (3.2) efficiently. In each subdomain, we solve the following problem ! ! ! (i) (i) (i) S∆∆ S∆Π Ψ∆ 0 (3.3) = (i) (i) (i) (i) (i) , SΠ∆ SΠΠ IΠ FΠΠ IΠ (i)

where the matrix IΠ is the identity matrix of a dimension equal to the number of primal variables of Ωi . We then obtain (i)

Ψ

(i)

=

Ψ∆

(i)



!

(i)

(i)

(i)

−(S∆∆ )−1 S∆Π IΠ

=

(i)



!

and also (i) −1

(i)

(i)

(i)

FΠΠ = SΠΠ − SΠ∆ S∆∆ (i)

(i)

S∆Π .

(i)

Let RΠ : WΠ → WΠ restrict the global primal variables to the subdomain Ωi . From Ψ(i) , we construct the coarse finite element space Ψ as follows:   (1) Ψ(1) RΠ   .. . Ψ= .   (N ) Ψ(N ) RΠ

Each column ψ of the matrix Ψ is related to a primal variable. Since, the vector ψ ∈ W has t t t the same values at the primal variables, we take ψ = (ψ∆ , ψΠ ) from the vector ψ and define

a matrix Ψ with the columns ψ. We then obtain (3.4)

t Ψ = RΠ −

N X

(i)

(i)

(i)

(i)

(R∆ )t (S∆∆ )−1 S∆Π RΠ ,

i=1

(i)

(i)

t where RΠ : WΠ → W∆ × WΠ and (R∆ )t : W∆ → W∆ × WΠ are zero extensions.

Let us now define

(3.5)



 RD,Γ =  



Dnn Dmm DΠΠ

  RΓ , 

10

KIM, DRYJA AND WIDLUND

where the matrices Dnn , Dmm and DΠΠ will be specified later. We then propose the following preconditioner M −1 for the problem (3.2) ( ! ) −1 S∆∆ 0 −1 t t −1 t (3.6) M = RD,Γ + Ψ(Ψ SΨ) Ψ RD,Γ , 0 0 where   S = diagi S (i)

and S∆∆ is given in (2.11). We will show that

Ψt SΨ = FΠΠ , where −1 FΠΠ = SΠΠ − SΠ∆ S∆∆ S∆Π =

N X i=1

  (i) (i) (i) (i) (i) (i) (RΠ )t SΠΠ − SΠ∆ (S∆∆ )−1 S∆Π RΠ .

From the definition of Ψ, we have t

Ψ SΨ =

N X

(i)

(i)

(RΠ )t (Ψ(i) )t S (i) Ψ(i) RΠ ,

i=1

and from (3.3), we obtain Ψt SΨ =

(3.7)

N X

(i)

(i)

(i)

(RΠ )t FΠΠ RΠ = FΠΠ .

i=1

Using the block Cholesky decomposition of Se as in Li and Widlund [17] and above, see

also (2.10), we have

−1 S∆∆

Se−1 =

0 +

t RΠ

! 0 0



N X

(i) (i) (i) (i) (R∆ )t (S∆∆ )−1 S∆Π RΠ

i=1

t RΠ −

N X

!

(i) (i) (i) (i) (R∆ )t (S∆∆ )−1 S∆Π RΠ

i=1

By combining the above equation with (3.4) and (3.7), we obtain ! −1 S∆∆ 0 t −1 e S = + Ψ(Ψt SΨ)−1 Ψ . 0 0

−1 FΠΠ

!t

.

Therefore, the BDDC operator, see (3.2), with the preconditioner M −1 in (3.6) can be written

as (3.8)

t e Γ. BDDC = RD,Γ Se−1 RD,Γ RΓt SR

A BDDC METHOD WITH MORTAR DISCRETIZATION

11

4. Condition number analysis using a bound on ED . In this section, we will estimate the condition number of the BDDC operator by using an approach introduced in [16]. A bound for the average operator ED in a certain norm is central in the analysis. We recall the definitions of RΓ and RD,Γ in (3.1) and (3.5), respectively. The operator ED is defined as t ED = RΓ RD,Γ ,

(4.1)

where the weight matrix D will be chosen so that t (P1) RΓt RD,Γ = RD,Γ RΓ = I (

(P2) |ED w|2Se ≤ C max i

Hi 1 + log hi

3 )

|w|2Se.

e wi. We then have Here |w|2Se = hSw, ! ! t wm −Bm (Bnt )−1 Dnn zn + Dmm wm t = , RΓ RD,Γ t wΠ −BΠ (Bnt )−1 Dnn zn + DΠΠ wΠ where zn = −Bn−1 (Bm wm + BΠ wΠ ). In order to satisfy property (P1), the weight matrix D is chosen so that (4.2)

Dnn = 0,

Dmm = I,

DΠΠ = I.

R EMARK 4.1. The weights above lead to an operator ED of the form     −Bn−1 (Bm wm + BΠ wΠ ) wn        ED  wm  wm  =  wΠ wΠ

that does not involve any averages across the interfaces in contrast to the average operator considered for conforming finite elements. We will still call ED the average operator just borrowing the name from the conforming case. We will now show that the average operator ED satisfies property (P2) with the weight matrix D just given. As a preparation, we need to establish an estimate for the mortar prof in the H 1/2 (F )-norm. For an edge/face F ⊂ ∂Ωi , the space jection of a function w in W 00

1/2

H00 (F ) consists of functions for which the zero extension to the whole boundary ∂Ωi belongs to the Sobolev space H 1/2 (∂Ωi ). It is equipped with the norm Z w(x)2 2 2 kwkH 1/2 (F ) = |w|H 1/2 (F ) + ds(x). 00 F dist(x, ∂F ) This norm has the well-known property (4.3)

c|w| e H 1/2 (∂Ωi ) ≤ kwkH 1/2 (F ) ≤ C|w| e H 1/2 (∂Ωi ) , 00

12

KIM, DRYJA AND WIDLUND

where w e is the zero extension of w to ∂Ωi \ F ; see [7, Lemma1.3.2.6].

We recall that the subdomain Ωj intersect the subdomain Ωi along Fij ⊂ F where F

is a nonmortar edge/face in ∂Ωi and that φ = wj on Fij . We then have φ ∈ H 1/2−ǫ (F ),

0 < ǫ ≤ 1/2 and the following estimate; see Proposition 3.2 in [2]. L EMMA 4.2. Assume that Ωi and Ωj are scaled by the diameter Hi of the Ωi . For φ ∈ H 1/2−ǫ (F ) and 0 < ǫ ≤ 1/2, we have ǫkφk2H 1/2−ǫ (F ) ≤ C

X

kwj k21/2,∂Ωj .

j

We need the following assumption on the coefficients of the elliptic problem. A SSUMPTION 4.3. The coefficients satisfy ρi ≤ Cρj where Ωi and Ωj are the nonmortar side and the mortar side of the common set Fij = ∂Ωi ∩ ∂Ωj . For any set A ⊂ ∂Ωi and wi ∈ Wi , we define a nodal value interpolant IA (wi ) ∈ Wi as ( wi (x) x ∈ A ∩ Vi , (4.4) IA (wi )(x) = 0 at the other nodes. Here Vi denotes the set of nodes in the finite element space Wi . Let F ⊂ ∂Ωi be a nonmortar edge/face. We denote by I(F ) the set containing the indices of the subdomains that intersect F , and by πF the mortar projection given on the edge/face F . We now provide the following bound for functions v ∈ L2 (F ) and with πF v = 0 on ∂F . f satisfies L EMMA 4.4. With Assumption 4.3 on the ρi , w = (w1 , · · · , wN ) ∈ W 3 X  Hi hS (k) wk , wk i, ρi kπF (φ − wi )k2H 1/2 (F ) ≤ C 1 + log hi 00 k∈I(F )

where F ⊂ ∂Ωi is a nonmortar edge/face. Proof. For any function v(x) ∈ L2 (F ) or L2 (Ωl ), let us define vb(x) = v(Hi x),

b l, x ∈ Fb or Ω

b l denote the dilated sets. From the definition where Hi is the diameter of the Ωi , and Fb and Ω

of the mortar projection, we see that

π\ b, F (v) = πF bv

(4.5)

where πFb denotes the mortar projection based on the finite element space dilated by Hi . We now consider kπF (φ − wi )k2H 1/2 (F ) ≤ 2kπF (φ)k2H 1/2 (F ) + 2kπF (wi )k2H 1/2 (F ) . 00

00

00

13

A BDDC METHOD WITH MORTAR DISCRETIZATION

Let φe = wj − cij ,

where

R

w fi = wi − cij Fij

cij = R

wi ψij ds

Fij

ψij ds

on Fij = ∂Ωj ∩ ∂Ωi ,

R

F = Rij

wj ψij ds

Fij

ψij ds

.

e w We then have φ, ei ∈ H 1/2−ǫ (F ) for 0 < ǫ ≤ 1/2 and φ − wi = φe − w ei , and can thus e replace φ and wi above by φ and w ei , respectively. By using a scaling argument and the identity (4.5), we have 2 kπF (φ)k2H 1/2 (F ) = Hid−2 kπ\ F (φ)kH 1/2 (F b) 00

00

= ≤

(4.6)

b 2 1/2 Hid−2 kπFb (φ)k b) H (F

2Hid−2



00

 2 b 2 1/2 b kπFb (φb − Qφ)k + kπ (Q φ)k 1/2 b b) b) , F H H (F (F 00

00

where Qφb is the L2 -projection of φb on the finite element space Wi (Fb), i.e., the dilated finite element space provided for the nonmortar edge/face Fb .

From an inverse inequality, the continuity of πFb in L2 (Fb ), the approximation property of Q for a function φb ∈ H 1/2−ǫ (Fb ), Lemma 4.2, and a scaling argument, we obtain b 2 1/2 b −1 b b 2 kπFb (φb − Qφ)k b) b) ≤ C hi kφ − QφkL2 (F H (F 00

b 1−2ǫ kφk b 2 1/2−ǫ ≤ Cb h−1 b) i hi H (F X kw bj k21,∂ Ω ≤ Cb h−2ǫ ǫ−1 b . i j

j

b Replacing φb with φe in the above estimate and using a Poincar´e inequality and a scaling argu-

ment, we find

b 2 1/2 b −2ǫ ǫ−1 kπFb (φb − Qφ)k b ) ≤ C hi H (F 00

(4.7)

We now estimate

(4.8)

X j

|w bnj |21,∂ Ω b

≤ CHi2−db h−2ǫ ǫ−1 i

j

X

|wj |21/2,∂Ωj .

j

  2

b + Qφb − I b (Qφ) b b 2 1/2 I (Q φ) kπFb (Qφ)k = π

1/2 b b b b) F F F H00 (F H00 (F )   −1 b 2 1/2 b b − I b (Qφ)k b 22 ≤ C kIFb (Qφ)k + h kQ φ b) , i b) F L (F H (F 00

where IFb (wi ) is the nodal value interpolant described in (4.4). Here, we have used that πFb is 1/2 a bounded map in H (Fb ) as well as in L2 (Fb), and also used an inverse inequality. By using 00

14

KIM, DRYJA AND WIDLUND

Lemma 4.24 in [21], an inverse inequality, the stability of Q in H 1/2−ǫ (Fb ), Lemma 4.2, a

Poincar´e inequality, and a scaling argument, we obtain  2 Hi b 2 1/2 b 2 1/2 kIFb (Qφ)k ≤ C 1 + log kQφk b) b) H (F H00 (F hi 2  Hi b −2ǫ b 2 hi kQφkH 1/2−ǫ (Fb) ≤ C 1 + log hi 2  Hi b −2ǫ b 2 hi kφkH 1/2−ǫ (Fb) ≤ C 1 + log hi 2  Hi b −2ǫ −1 X kw bj k21/2,∂ Ω hi ǫ ≤ C 1 + log bj hi j 2  X Hi b (4.9) |wj |21/2,∂Ωj . h−2ǫ ǫ−1 ≤ CHi2−d 1 + log i hi j

b has nonzero value only at the nodes on the boundary of Fb. We note that Qφb − IFb (Qφ)

In two dimensions, by using Lemma 4.15 in [21], we obtain

b 22 b b 2 kQφb − IFb (Qφ)k b) ≤ C hi kQφk∞,F b L (F   Hi b b 2 1/2 ≤ C hi 1 + log kQφk b) H (F hi

and in three dimensions, by using Lemma 4.17 in [21], we also obtain

b b 2 b 22 kQφb − IFb (Qφ)k b) ≤ C hi kQφkL2 (∂ F b) L (F   Hi b 2 1/2 . b kQφk ≤ C hi 1 + log b) H (F hi

b 2 The same estimate, as before, for the term kQφk b) gives H 1/2 (F

(4.10)

  Hi b −2ǫ −1 X 2−d b 2 b b |wj |21/2,∂Ωj . hi ǫ kQφ − IFb (Qφ)kL2 (Fb) ≤ CHi hi 1 + log hi j

Combining (4.6) with (4.7)-(4.10) results in 2  X ρi Hi 2 b hS (j) wj , wj i. h−2ǫ ǫ−1 ρi kπF (φ)kH 1/2 (F ) ≤ C 1 + log i hi ρ 00 j j

The desired bound follows by letting ǫ = 1/(2|logb hi |) and using the assumption that ρi /ρj ≤ −2ǫ C. We note that b hi = hi /Hi and log(b h ) = 1. The same analysis applied to kπF (wi )k2 1/2 i

gives

3  Hi hS (i) wi , wi i. ρi kπF (wi )k2H 1/2 (F ) ≤ C 1 + log hi 00

H00 (F )

A BDDC METHOD WITH MORTAR DISCRETIZATION

15

R EMARK 4.5. For the geometrically conforming case, Lemma 4.4 is valid with a factor (1 + log(Hi /hi ))2 using the same analysis as above. Therefore, in this case, we obtain a better condition number estimate; see also Remark 4.8 below. The estimate improves the one in [9, 11] by using the projection Q; we do not need the assumption on the mesh sizes hi and hj hj ≤C hi



ρj ρi



for some 0 ≤ γ ≤ 1,

considered in [9, 11] where Ωi is the nonmortar side and Ωj is the mortar side. With the help of Lemma 4.4, we can establish property (P2) for the operator ED . L EMMA 4.6. With Assumption 4.3, the operator ED satisfies

|ED w|2Se

≤ C max i

(

Hi 1 + log hi

3 )

|w|2Se.

Proof. Using the weight matrix D in (4.2), the average operator ED in (4.1) is given by 





 wn − Bn−1 (Bn wn + Bm wm + BΠ wΠ )       , ED  wm wm  =   wΠ wΠ wn

see Remark 4.1. Let

zn = wn − Bn−1 (Bn wn + Bm wm + BΠ wΠ ), and construct zi by restricting the unknowns (zn , wm , wΠ ) to the subdomain Ωi . Similarly, we construct wi from (wn , wm , wΠ ). We note that (w1 , · · · , wN ) satisfies the primal conc , i.e., z satisfies the mortar straints on the edges/faces. By definition, z = (z1 , · · · , zN ) ∈ W matching condition, and each zi is of the form

zi = wi +

X

(i)

EF πF (φ − wi ),

F ⊂∂Ωi

(i)

where F is a nonmortar edge/face in ∂Ωi , EF is the zero extension of functions defined on

16

KIM, DRYJA AND WIDLUND

F to all of ∂Ωi \ F , and φ = wj on Fij (:= ∂Ωj ∩ ∂Ωi ) ⊂ F . We then obtain |ED w|2Se =

N X

hS (i) zi , zi i

i=1

≤C

N X

hS

(i)

wi , wi i +

i=1

≤C

N X i=1

X

hS

(i)

(i) EF πF (φ



(i) wi ), EF πF (φ

− wi )i

F ⊂∂Ωi

hS (i) wi , wi i + (

N X X

i=1 F ⊂∂Ωi

!

ρi kπF (φ − wi )k2H 1/2 (F ) 00

3 ) X N Hi 1 + log ≤ C max hS (i) wi , wi i i hi i=1 ( 3 ) Hi e wi. 1 + log hSw, ≤ C max i hi

Here we have used that hS (i) wi , wi i ≃ ρi |wi |2H 1/2 (∂Ωi ) , the relation in (4.3), and Lemma 4.4. By using the properties (P1) and (P2), we can show the following condition number bound of the BDDC operator (3.8). A similar proof is given in Li and Widlund [16] in an analysis of a BDDC algorithm for the Stokes problem with conforming meshes. T HEOREM 4.7. With Assumption 4.3, we have the condition number bound ( 3 ) Hi . 1 + log κ(BDDC ) ≤ C max i hi Proof. We let

and we then have

t M −1 = RD,Γ Se−1 RD,Γ ,

e Γ, Sb = RΓt SR

b BDDC = M −1 S.

We will now provide a lower bound by proving

b b. hu, uiSb ≤ hu, M −1 Sui S

b From property (P1), Rt RD,Γ = Rt RΓ = I, we obtain u = Let w = Se−1 RD,Γ Su. Γ D,Γ −1 t e b S R Sw. We then consider Γ

b hu, uiSb = ut Su

e = ut RΓt Sw

= hw, RΓ uiSe 1/2 S 1/2 1/2 hw, wiSe hu, uiSb .

1/2 S

≤ hw, wi e hRΓ u, RΓ ui e ≤

A BDDC METHOD WITH MORTAR DISCRETIZATION

17

Here we have used the Cauchy-Schwarz inequality. Squaring and cancelling a common factor, we obtain hu, uiSb ≤ hw, wiSe . By combining the above estimate with b b t Se−1 SeSe−1 RD,Γ Su hw, wiSe = ut SR D,Γ t b b = hu, RD,Γ Se−1 RD,Γ Sui S

(4.11) we obtain the desired lower bound.

b b, = hu, M −1 Sui S

We will now find an upper bound by proving ( 3 ) Hi 1/2 −1 b −1 b 1/2 1 + log hM Su, M Sui b ≤ C max hu, ui b . S S i hi

We consider

b M −1 Sui b b = hRt w, Rt wi b hM −1 Su, D,Γ D,Γ S S

t t = hRΓ RD,Γ w, RΓ RD,Γ wiSe

= |ED w|2Se (

Hi 1 + log hi

≤ C max i

3 )

|w|2Se.

The last inequality follows from Lemma 4.6. Combining the above estimate with (4.11), we obtain b M −1 Sui b b ≤ C max hM −1 Su, S i

(

Hi 1 + log hi

3 )

b b. hu, M −1 Sui S

b b, the desired upper By applying the Cauchy-Schwarz inequality to the term hu, M −1 Sui S bound follows.

R EMARK 4.8. The analysis above can be modified for the geometrically conforming

case and leads to the condition number bound κ(BDDC ) ≤ C max i

(

Hi 1 + log hi

2 )

,

when Assumption 4.3 holds; see Remark 4.5. R EMARK 4.9. For a geometrically non-conforming partition, the number of primal constraints tends to be bigger than for a conforming partition in case only edge/face constraints are used. We note that there have been several previous studies which explore the possibility of selecting primal constraints for only some of the edges/faces; see [10, 14, 15].

18

KIM, DRYJA AND WIDLUND

Ωi F

1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1

00 11 11 00 00 11

Ωl

11 00 00 11 00 11 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1

00 11 11 00 0 00 1 11 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 00 11 0 1 00 11 0 1 00 11 0 1 0 1 00 11 0 1 00 1 11 0 00 11 0 1 0 1 0 1 00 11 0 1 00 11 0 1 00 11 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 00 11 111 000 0 1 00 11 00 11 00 11 111111111111111 000000000000000 00 11 00 11 11 00 00 11 00 11

i Fil

l

Fil

k

Fik

Fiki

Fi l

Fik Ωk

F IG . 1. Geometrically non-conforming partition: white circles (nodes in Vjij (⊂ ∂Ωj ), j = l, k or Viij (⊂ ∂Ωi ), j = l, k), black circles (nodes in Njij (⊂ ∂Ωj ), j = l, k or Niij (⊂ ∂Ωi ), j = l, k), each faces F , Fij , i and F j for j = l, k are described. Fij ij

We will now provide a better estimate for geometrically non-conforming partitions under an assumption on the meshes that is considered in [9, 11]. We will prove our result only for the two-dimensional case; in three dimensions there are some additional technical difficulties. We conjecture that the result also holds in that case. A SSUMPTION 4.10. The mesh sizes hi and hj satisfy  γ hj ρj ≤C for some 0 ≤ γ ≤ 1, hi ρi where Ωi is the nonmortar side and Ωj is the mortar side. f satisfies L EMMA 4.11. With Assumptions 4.3 and 4.10, (w1 , · · · , wN ) ∈ W ) ( 2 X Hk 2 hS (k) wk , wk i, 1 + log ρi kπF (φ − wi )kH 1/2 (F ) ≤ C max hk k∈I(F ) 00 k∈I(F )

where F ⊂ ∂Ωi is a nonmortar edge/face and I(F ) is the set of the indices of the subdomains that intersect F . Proof. We consider the case in Figure 1. The nonmortar edge F ⊂ ∂Ωi is partitioned f into Fil and Fik and φ is given by wj on Fij , j = l, k. Since the function (w1 , · · · , wN ) ∈ W satisfies the primal constraints, we have Z Z (φ − wi )ψij ds = and we then define

R

Fij

Let

(wj − wi )ψij = 0,

j = l, k,

Fij

Fij

cij = R

wi ψij ds

Fij

ψij ds

R

F = Rij

wj ψij ds

Fij

ψij ds

,

w ei = wi − cij and φe = wj − cij on Fij ,

j = l, k.

j = l, k.

19

A BDDC METHOD WITH MORTAR DISCRETIZATION

We note that φ − wi = φe − w ei .

Let Viij be the set of nodes in Vi with nodal basis functions supported in Fij . We denote

by Fiji the union of these supports. We note that Fiji ⊂ Fij . The set Vjij and Fijj are defined

similarly; see Figure 1. Let Niij be the set of nodes in Vi \ Viij with nodal basis functions with support that intersects Fij . The set Njij is defined similarly. In general, we may assume that the number of nodes in each set Njij and Niij is bounded uniformly with respect to the mesh parameters. We consider (4.12) kπF (φ − wi )k2H 1/2 (F ) 00



X X

 = IF j (wj − cij ) − IFiji (wi − cij )

πF ij

j=l,k j=l,k

 2

e  +φ − ei + IF j (wj − cij ) − w IFiji (wi − cij )

ij

j=l,k j=l,k X

X

1/2 H00 (F )

1/2

Since the first two terms in the above equation are in H00 (F ), the continuity of the 1/2

mortar projection in H00 (F ) and Lemma 4.24 in [21] give kπF (IF j (wj − cij ))k2H 1/2 (F ) ≤ CkIF j (wj − cij )k2H 1/2 (F ) ij

ij

00

00

cij )k2H 1/2 (F j ) 00 ij

= CkIF j (wj − ij 2  Hj kwj − cij k2H 1/2 (∂Ωj ) , ≤ C 1 + log hj

(4.13) and (4.14)

 2 Hi kπF (IFiji (wi − cij ))k2H 1/2 (F ) ≤ C 1 + log kwi − cij k2H 1/2 (∂Ωi ) . hi 00

We now bound the third term in (4.12) as follows: X kπF (φe − IF j (wj − cij ))k2H 1/2 (F ) ij

00

j=l,k

e ≤ Ch−1 i kπF (φ −

(4.15)

e ≤ Ch−1 i kφ −

≤ Ch−1 i

X

X

j=l,k

X

IF j (wj − cij ))k2L2 (F ) ij

IF j (wj − cij )k2L2 (F ) ij

j=l,k

kwj − cij − IF j (wj − cij )k2L2 (Fij ) ij

j=l,k

≤ Ch−1 i

X

hj kwj − cij k2L∞ (∂Ωj )

j=l,k

(4.16)



Ch−1 i

X

j=l,k

hj



Hj 1 + log hj



kwj − cij k2H 1/2 (∂Ωj ) .

.

20

KIM, DRYJA AND WIDLUND

We have used an inverse inequality, the continuity of πF in L2 (F ) and Lemma 4.15 in [21]. The expression in (4.15) has nonzero values at the nodes in Njij . By using the fact that the number of nodes in Njij is bounded independently of any mesh parameters (at most three in Figure 1), we have kwj − cij − IF j (wj − cij )k2L2 (Fij ) ≤ Chj kwj − cij k2L∞ (∂Ωj ) . ij

Similarly, we have a bound for the last term in (4.12): X kπF (w ei − IFiji (wi − cij ))k2H 1/2 (F ) 00

j=l,k

(4.17)

  Hi X ≤ C 1 + log kwi − cij k2H 1/2 (∂Ωi ) . hi j=l,k

As a result, (4.12), (4.13), (4.14), (4.16), and (4.17) give ( 2 ) H k 1 + log kπF (φ − wi )k2H 1/2 (F ) ≤ C max hk k∈I(F ) 00     X hj kwj − cij k2H 1/2 (∂Ωj ) . kwi − cij k2H 1/2 (∂Ωi ) + 1 + hi j=l,k

A Poincar´e inequality can be applied to the functions wi − cij and wj − cij and this replaces the norms by semi-norms. By using the relation ρi |wi |2H 1/2 (∂Ωi ) ≃ hS (i) wi , wi i, we obtain the following bound ρi kπF (φ − wi )k2H 1/2 (F ) 00

(

2 ) Hk ≤ C max 1 + log hk k∈I(F )     X hj ρi (j) hS (i) wi , wi i + 1 + hS wj , wj i . h i ρj j=l,k

By using Assumptions 4.3 and 4.10, we then have   1−γ !  ρi h j ρi ≤ C. ≤C 1+ 1+ h i ρj ρj Therefore the required bound holds with a constant C which does not depend further on the mesh parameters and the jumps of the coefficients. By using Lemma 4.11 and the same analysis as in Theorem 4.4, we obtain a better condition number bound for the geometrically non-conforming case. T HEOREM 4.12. For a geometrically non-conforming subdomain partition and in two dimensions, the BDDC operator satisfies κ(BDDC ) ≤ C max i

when Assumptions 4.3 and 4.10 hold.

(

Hi 1 + log hi

2 )

,

A BDDC METHOD WITH MORTAR DISCRETIZATION

21

5. A connection between FETI–DP and BDDC methods. In this section, we will show that the BDDC algorithm developed in the previous sections is closely connected to the FETI–DP algorithm developed by the first author in [9, 10] and by her jointly with Lee in [11]. These two algorithms will be shown to share the same spectra except possibly for an eigenvalue equal to 1. A study comparing the spectra of the BDDC algorithm to that of the FETI–DP algorithm was carried out by Mandel, Dohrmann and Tezaur [20] for conforming finite elements. They showed that these two algorithms have the same set of eigenvalues except possibly for eigenvalues equal to 0 and 1. Recently, a quite simple proof of this fact was given by Li and Widlund [17]. They formulated the BDDC operators as well as the FETI–DP operators using a change of variables and introducing certain projections and average operators. These projections and average operators provide an important connection between the FETI–DP and the BDDC operators. We first formulate a FETI–DP operator with the change of variables introduced in Section 2.2. We then show that the FETI–DP operator has essentially the same spectrum as the BDDC operator by establishing several properties of the projections and average operators that are used in the analysis by Li and Widlund [17]. After the change of variables, the linear system considered in the FETI–DP formulation is given by

(5.1)



S∆∆

  SΠ∆  B∆

S∆Π SΠΠ BΠ

t B∆



u∆





g∆



    t     BΠ   uΠ  =  gΠ  , λ

0

0

where the matrices S∆∆ , S∆Π , SΠ∆ , and SΠΠ are defined in (2.11) and the matrices B∆ and BΠ are obtained from the mortar matching condition (2.8). We note that the subscripts Π and ∆ stand for the unknowns or submatrices related to the primal variables and the remaining part, respectively, and that λ ∈ M , the non-redundant Lagrange multiplier space. After eliminating the unknowns u∆ and uΠ , we obtain an equation for λ ∈ M : BΓ Se−1 BΓt λ = d,

(5.2) where (5.3)

 BΓ = B∆



BΠ ,

Se =

S∆∆

S∆Π

SΠ∆

SΠΠ

!

,

and d is also the result of Gaussian elimination.

We will now express the Neumann-Dirichlet preconditioner considered in [9, 10, 11] using the change of variables. We recall the space Wn , defined in (2.4) and then define ) ( Z fn := wn ∈ Wn : wn ψij ds = 0, ∀Fij , ∀F ⊂ ∂Ωi , ∀i , W Fij

22

KIM, DRYJA AND WIDLUND

where F ⊂ ∂Ωi are nonmortar edges/faces with the partition {Fij }j . For the geometrically fn consists of functions with zero average on each nonmortar conforming case, the space W

edge/face F because ψij = 1.

−1 The Neumann-Dirichlet preconditioner MDP is defined by

hMDP λ, λi = max

(5.4)

fn wn ∈W

hBE(wn ), λi2 , hSE(wn ), E(wn )i

 B (1) · · · B (N ) . Here we consider λ ∈ M , a non-redundant Lagrange multiplier space, and hence the morwhere E(wn ) is the zero extension of wn to all the interfaces and B =



tar matching matrix B has one less row for each nonmortar edge/face than in the original formulation in [9, 10, 11]. We recall the space W∆ given in (2.7) and note that it can be split into W∆ = W∆,n ⊕ W∆,m , where n and m denote the unknowns of the nonmortar edges/faces and mortar edges/faces, respectively. The vectors in these spaces are represented by the unknowns after the change of fn except that the bases are different. variables. The space W∆,n is then identical to W By using the change of variables, (5.4) can be written as hMDP λ, λi =

(5.5)

max

w∆,n ∈W∆,n

where  b= B b (1) B

···

b E(w b ∆,n ), λi2 hB , b ∆,n ), E(w b ∆,n )i hSb E(w

 b (N ) , B

Sb = diag(Sb(i) ).

b act on the new unknowns w(i) and w(i) that result from the change of The matrices Sb and B Π ∆ b ∆,n ) = (w1 , · · · , wN ) is given by variables. The extension E(w (i)

wi =

w∆

(i)



!

,

(i)

(i)

where w∆ is zero on the mortar edges/faces, w∆ is equal to w∆,n on the nonmortar edges/faces, (i)

and wΠ is zero. The formula (5.5) can be written as (5.6)

hMDP λ, λi =

max

w∆,n ∈W∆,n

hBn w∆,n , λi2 . hSnn w∆,n , w∆,n i

Here the matrices Bn and Snn are submatrices of B∆ and S∆∆ in (5.1) corresponding to the nonmortar part. We see that Snn : W∆,n → W∆,n and Bn : W∆,n → M are invertible. The maximum in (5.6) occurs when Snn w∆,n = Bnt λ and hence it follows that −1 MDP = (Bnt )−1 Snn Bn−1 .

A BDDC METHOD WITH MORTAR DISCRETIZATION

23

Further this matrix can be written as −1 t e Σ,Γ , MDP = BΣ,Γ SB

(5.7) where



 t BΣ,Γ = 

Σnn Σmm

with the weights given by

Σnn = (Bnt Bn )−1 ,



Bnt



   B t    m t BΠ ΣΠΠ

Σmm = 0,

ΣΠΠ = 0.

Here the matrix Bm is a submatrix of B∆ corresponding to the unknowns of the mortar part. −1 Therefore the FETI–DP operator with the Neumann-Dirichlet preconditioner MDP is

given by −1 e t BΓ Se−1 B t , MDP FDP = BΣ,Γ SB Γ Σ,Γ

while the preconditioned BDDC operator is given by

t e Γ. BDDC = RD,Γ Se−1 RD,Γ RΓt SR

Let us now define the following jump and average operators t PΣ = BΣ,Γ BΓ ,

t ED = RΓ RD,Γ .

The following results are provided in [17, Section 5]. T HEOREM 5.1. Assume that PΣ and ED satisfy 1. ED + PΣ = I, 2 2. ED = ED , PΣ2 = PΣ ,

3. ED PΣ = PΣ ED = 0. −1 Then the operators MDP FDP and BDDC have the same eigenvalues except possibly for the

eigenvalue equal to 1. We will now show that the assumptions in Theorem 5.1 hold for the operators PΣ and f by using the unknowns wn , wm , and wΠ : ED . We express the space W  f = (wt , wt , wt )t : ∀wn , wm , wΠ , W n m Π

and we recall the mortar finite element space

c = {w ∈ W f : Bm wm + BΠ wΠ + Bn wn = 0}. W

f. We note that PΣ and ED are operators defined on the space W L EMMA 5.2. The operators PΣ and ED satisfy

24

KIM, DRYJA AND WIDLUND

1. ED + PΣ = I, 2 2. ED = ED , PΣ2 = PΣ ,

3. ED PΣ = PΣ ED = 0. Proof. From Σmm = 0, ΣΠΠ = 0, Σnn = (Bnt Bn )−1 , Dmm = I, DΠΠ = I, Dnn = 0, we have 

 Bn−1 (Bm wm + BΠ wΠ + Bn wn )   , PΣ w =  0   0   −Bn−1 (Bm wm + BΠ wΠ )   . ED w =  wm   wΠ

Hence,

ED + PΣ = I. 2 c and ED w = w for all We will now show that ED = ED . Since Range(ED ) ∈ W c , we obtain w∈W

This implies that

f. ED (ED w) = ED w for all w ∈ W 2 ED = ED .

(5.8)

2 From ED + PΣ = I and ED = ED , we have

ED (ED + PΣ ) = ED and therefore ED PΣ = 0. c and Range(ED ) ∈ W c , we can show that Moreover, from PΣ w = 0 for all w ∈ W PΣ ED = 0.

To show that PΣ2 = PΣ , we consider PΣ (ED + PΣ ) = PΣ ,

A BDDC METHOD WITH MORTAR DISCRETIZATION

25

and from PΣ ED = 0, we obtain PΣ2 = PΣ .

R EMARK 5.3. Other FETI–DP preconditioners in two dimensions with different weights   Σnn    Σ= Σmm   ΣΠΠ

have been developed and shown to give a condition number bound  C max (1 + log(Hi /hi ))2 i

for some geometrically conforming cases with nonzero weights Σmm and ΣΠΠ ; see [6, 5]. We have not found a weight matrix D that results in ED + PΣ = I for such a choice of Σ. 6. Numerical tests. In this section, we discuss numerical tests which compare the efficiency of the BDDC method and that of the FETI–DP method. For Ω = [0, 1]2 , we solve the elliptic problem with the exact solution u(x, y) = sin(πx)(1 − y)y; −∆u = f in Ω, u = 0 on ∂Ω. We have carried out experiments for both matching and non-matching grids employing the mortar matching conditions across the interfaces. The CG (Conjugate Gradient) iteration continues until the residual norm has been reduced by a factor 10−6 . The domain Ω is divided into square subdomains as in Figure 2. For matching grids, we introduce uniform meshes with n nodes on each horizontal and vertical edge. To make the meshes non-matching across subdomain interfaces, we generate triangulations in each subdomain in the following way: for each subdomain, we choose n random quasi-uniform nodes on each horizontal and vertical edges. From these nodes, we generate nonuniform structured grids in each subdomain. Since we choose the same number of quasi-uniform nodes n for all subdomains, the mesh sizes of neighboring subdomains are comparable. First, we compare the two algorithms with the matching grids employing the mortar matching condition and primal constraints at the vertices. In Table 1, we divide Ω into N = 4 × 4 subdomains (see Figure 2) and increase the number of nodes n. We compute L2 - and H 1 -errors between the exact solution and the solution of the iterative method, the number of CG iterations, and the minimum and the maximum eigenvalues of the BDDC and the FETI– DP operators. For the H 1 -error, we compute the broken H 1 -norm based on the subdomain partition. Table 2 shows the numerical results when we fix n − 1 = 4 and increase N , the

26

KIM, DRYJA AND WIDLUND

Ω33 Ωij Ω01 Ω00

Ω10

F IG . 2. Partition of subdomains when N = 4 × 4

number of subdomains. For N = 8 × 8, 16 × 16 and 32 × 32, we divide Ω into square subdomains in the same manner as for N = 4 × 4. We see that both methods gives the same accuracy. The minimum eigenvalues of the BDDC operator are always equal to 1 while those of the FETI–DP operator are greater than 1. The maximum eigenvalues of both operators are almost the same. In Table 3 and 4, we perform the same computations for non-matching grids. The results shows similar patterns for the minimum and maximum eigenvalues as for matching grids except that the minimum eigenvalues of FETI–DP operator converge to 1 when the number of nodes increases; see Table 3. From the numerical results, we see that the BDDC operator always has the minimum eigenvalue 1 while the FETI–DP operator has all its eigenvalues greater than 1 and that these operators have almost the same maximum eigenvalues. Generally, we can conclude that the two algorithms perform quite similarly. TABLE 1 (Matching grids) Comparison of FETI–DP and BDDC methods when n increases with a fixed number of subdomain N = 4 × 4 −1 MDP FDP h

h

BDDC

n−1

ku − u k0

ku − u k1

Iter

λmin

λmax

Iter

λmin

λmax

4

4.1293e-4

5.7497e-2

10

1.43

4.01

11

1.00

4.01

8

1.0399e-4

2.8798e-2

12

1.35

5.64

13

1.00

5.64

16

2.6057e-5

1.4405e-2

14

1.31

7.64

15

1.00

7.64

32

6.5183e-6

7.2036e-3

15

1.31

1.00e+1

16

1.00

1.00e+1

64

1.6315e-6

3.6019e-3

16

1.35

1.27e+1

18

1.00

1.27e+1

REFERENCES [1] F. B. B ELGACEM AND Y. M ADAY, The mortar element method for three dimensional finite elements, M2 AN Math. Model. Numer. Anal., 31 (1997), pp. 289–302.

27

A BDDC METHOD WITH MORTAR DISCRETIZATION TABLE 2

(Matching grids) Comparison of FETI–DP and BDDC methods when N increases with a fixed number of nodes (n − 1 = 4) in each subdomain −1 MDP FDP h

h

BDDC

N

ku − u k0

ku − u k1

Iter

λmin

λmax

Iter

λmin

λmax

4×4

4.1293e-4

5.7497e-2

10

1.43

4.01

11

1.00

4.01

8×8

1.0399e-4

2.8798e-2

11

1.39

4.21

12

1.00

4.21

16 × 16

2.6057e-5

1.4405e-2

11

1.39

4.20

12

1.00

4.26

32 × 32

6.5183e-6

7.2036e-3

11

1.41

4.20

12

1.00

4.27

TABLE 3 (Non-matching grids) Comparison of FETI–DP and BDDC methods when n increases with a fixed number of subdomain N = 4 × 4 −1 MDP FDP h

h

BDDC

n−1

ku − u k0

ku − u k1

Iter

λmin

λmax

Iter

λmin

λmax

4

5.0850e-4

6.0126e-2

10

1.40

4.09

12

1.00

4.09

8

1.2865e-4

3.0128e-2

13

1.01

5.72

15

1.00

5.72

16

3.2231e-5

1.5072e-2

15

1.00

7.72

16

1.00

7.72

32

8.0621e-6

7.5374e-3

16

1.01

1.00e+1

17

1.00

1.00e+1

64

2.0134e-6

3.7688e-3

17

1.01

1.28e+1

19

1.00

1.28e+1

[2] C. B ERNARDI

AND

Y. M ADAY, Spectral, spectral element and mortar element methods, tech. rep., Berlin,

2001. [3] C. B ERNARDI , Y. M ADAY,

AND

A. T. PATERA, A new nonconforming approach to domain decomposition:

the mortar element method, in Nonlinear partial differential equations and their applications. Coll`ege de France Seminar, Vol. XI (Paris, 1989–1991), vol. 299 of Pitman Res. Notes Math. Ser., Longman Sci. Tech., Harlow, 1994, pp. 13–51. [4] C. R. D OHRMANN, A preconditioner for substructuring based on constrained energy minimization, SIAM J. Sci. Comput., 25 (2003), pp. 246–258. [5] N. D OKEVA , M. D RYJA ,

AND

W. P ROSKUROWSKI, A FETI-DP preconditioner with special scaling for

mortar discretization of elliptic problems with discontinuous coefficients, to appear in SIAM J. Numer. Anal., (2005). [6] M. D RYJA AND O. B. W IDLUND, A generalized FETI-DP method for a mortar discretization of elliptic problems, in Domain decomposition methods in science and engineering, Natl. Auton. Univ. Mex., M´exico, 2003, pp. 27–38. Proceedings of the 14th international conference in domain decomposition methods. [7] P. G RISVARD, Elliptic problems in nonsmooth domains, vol. 24 of Monographs and Studies in Mathematics, Pitman (Advanced Publishing Program), Boston, MA, 1985. [8] C. K IM , R. D. L AZAROV, J. E. PASCIAK , AND P. S. VASSILEVSKI, Multiplier spaces for the mortar finite element method in three dimensions, SIAM J. Numer. Anal., 39 (2001), pp. 519–538. [9] H. H. K IM, A preconditioner for the FETI-DP formulation with mortar methods in three dimensions, tech. rep., KAIST, 2004. [10]

, A FETI-DP formulation of three dimensional elasticity problems with mortar discretization, tech. rep., New York University, 2005.

28

KIM, DRYJA AND WIDLUND TABLE 4 (Non-matching grids) Comparison of FETI–DP and BDDC methods when N increases with a fixed number of

nodes (n − 1 = 4) in each subdomain −1 MDP FDP h

h

BDDC

N

ku − u k0

ku − u k1

Iter

λmin

λmax

Iter

λmin

λmax

4×4

5.0850e-4

6.0126e-2

10

1.40

4.09

12

1.00

4.09

8×8

1.1744e-4

2.9900e-2

11

1.37

4.41

12

1.00

4.41

16 × 16

2.9743e-5

1.4980e-2

12

1.32

4.49

13

1.00

4.49

32 × 32

7.4317e-6

7.4917e-3

12

1.30

4.57

13

1.00

4.62

[11] H. H. K IM

AND

C.-O. L EE, A preconditioner for the FETI-DP formulation with mortar methods in two

dimensions, SIAM J. Numer. Anal., 42 (2005), pp. 2159–2175. [12] A. K LAWONN

AND

O. R HEINBACH, A parallel implementation of Dual-Primal FETI methods for three

dimensional linear elasticity using a transformation of basis, tech. rep., Technical Report SM-E-601, University of Essen, 2005. [13]

, Robust FETI-DP methods for heterogeneous three dimensional elasticity problems, tech. rep., Preprint, 2005.

[14] A. K LAWONN AND O. B. W IDLUND, Dual-Primal FETI methods for linear elasticity, tech. rep., New York University, 2004. [15] A. K LAWONN , O. B. W IDLUND , AND M. D RYJA, Dual-primal FETI methods for three-dimensional elliptic problems with heterogeneous coefficients, SIAM J. Numer. Anal., 40 (2002), pp. 159–179. [16] J. L I

AND

O. W IDLUND, BDDC algorithms for incompressible Stokes equations, tech. rep., New York Uni-

versity, 2005. [17] [18]

, FETI-DP, BDDC, and block Cholesky methods, (2005). , On the use of inexact subdomain solvers for BDDC algorithms, tech. rep., New York University, 2005.

[19] J. M ANDEL AND C. R. D OHRMANN, Convergence of a balancing domain decomposition by constraints and energy minimization, Numer. Linear Algebra Appl., 10 (2003), pp. 639–659. [20] J. M ANDEL , C. R. D OHRMANN , AND R. T EZAUR, An algebraic theory for Primal and Dual substrucruting methods by constraints, Appl. Numer. Math., 54 (2005), pp. 167–193. [21] A. T OSELLI AND O. W IDLUND, Domain decomposition methods—algorithms and theory, vol. 34 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 2005. [22] X. T U, Three-level BDDC in two dimensions, tech. rep., New York University, 2004. [23]

, Three-level BDDC in three dimensions, tech. rep., New York University, 2005.

[24] B. I. W OHLMUTH, A mortar finite element method using dual spaces for the Lagrange multiplier, SIAM J. Numer. Anal., 38 (2000), pp. 989–1012. [25]

, Discretization methods and iterative solvers based on domain decomposition, vol. 17 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2001.