Citation
Analysis of lagrange multiplier based domain decomposition

Material Information

Title:
Analysis of lagrange multiplier based domain decomposition
Creator:
Tezaur, Radek
Publication Date:
Language:
English
Physical Description:
113 leaves : ; 28 cm

Subjects

Subjects / Keywords:
Finite element method ( lcsh )
Boundary value problems ( lcsh )
Neumann problem ( lcsh )
Lagrangian functions ( lcsh )
Boundary value problems ( fast )
Finite element method ( fast )
Lagrangian functions ( fast )
Neumann problem ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 106-113).
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Radek Tezaur.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
40337054 ( OCLC )
ocm40337054
Classification:
LD1190.L622 1998d .T49 ( lcc )

Downloads

This item has the following downloads:


Full Text
ANALYSIS OF LAGRANGE MULTIPLIER BASED
DOMAIN DECOMPOSITION
by
Radek Tezaur
M.S., Univerzita Karlova, Praha, Czech Republic, 1993
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1998


This thesis for the Doctor of Philosophy
degree by
Eadek Tezaur
has been approved
by
Jan Mandel
Charbel Farhat
Leo Franca
Andrew Knyazev
Tom Russell
Date


Tezaur, Eadek (Ph. D., Applied Mathematics)
Analysis of Lagrange multiplier based domain decomposition
Thesis directed by Professor Jan Mandel
ABSTRACT
The convergence of a substructuring iterative method with Lagrange
multipliers known as Finite Element Tearing and Interconnecting (FETI) meth-
od is analyzed in this thesis. This method, originally proposed by Farhat and
Roux, decomposes finite element discretization of an elliptic boundary value
problem into Neumann problems on the subdomains, plus a coarse problem
for the subdomain null space components. For linear conforming elements
and preconditioning by Dirichlet problems on the subdomains, the asymptotic
bound on the condition number (7(1 +log(H/h))'7, where 7 = 2 or 3, is proved
for a second order problem, h denoting the characteristic element size and H
the size of subdomains. A similar method proposed by Park is shown to be
equivalent to FETI with a special choice of some components and the bound
(7(1 + log (H/h))2 on the condition number is established. Next, the origi-
nal FETI method is generalized to fourth order plate bending problems. The
main idea there is to enforce continuity of the transversal displacement held at
the subdomain crosspoints throughout the preconditioned conjugate gradient
iterations. The resulting method is shown to have a condition number that
does not increase with the number of subdomains, and again grows at most
111


poly-logarithmically with the number of elements per subdomain; the condi-
tion number is bounded by (7(1 + log(H/h))3. These optimal properties hold
for numerous plate bending elements that are used in practice including the
HCT, DKT, and a class of non-locking elements for the Reissner-Mindlin plate
models. The theoretical results are confirmed by numerical experiments.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed __________________________________________
Jan Mandel


DEDICATION
Dedicated to my parents.


ACKNOWLEDGEMENT
Foremost I would like to thank my advisor, Jan Mandel, for his encourage-
ment in preparation of this thesis, his support, guidance, and time spent in
discussions throughout my studies at the University of Colorado at Denver. I
would also like to extend my thanks to Patrick Le Tallec for useful suggestions
in early stages of this work and Charbel Farhat for computational results.


CONTENTS
Chapter
1. Introduction................................................. 1
2. Overview of Related Domain
Decomposition Methods........................................ 7
2.1 Finite Element Approximation, Model Problem and PCG 7
2.2 Abstract Schwarz Methods............................... 12
2.3 Overlapping Domain Decomposition....................... 15
2.4 Substructuring Methods................................. 18
2.4.1 Poincare-Steklov Operators and Schur Complement . 18
2.4.2 Various Preconditioners for the Reduced Problem ... 25
2.4.3 The Neumann-Neumann and Balancing Domain
Decomposition Methods................................ 27
2.4.4 The Neumann-Neumann Domain Decomposition
Algorithm for Plates and Shells ..................... 29
2.4.5 Lagrange Multipliers and Poincare-Steklov
Operators ........................................... 31
2.5 A Two Level Method Based on Smoothed
Prolongation........................................... 33
3. Finite Element Tearing and
Interconnecting (Derivation and Extensions)................. 37
3.1 Original FETI.......................................... 37
3.1.1 Problem Setting and Assumptions...................... 38
3.1.2 Derivation of the Dual Problem ...................... 42
vii


3.1.3 Saddle Point Formulation
45
3.1.4 Algebraic Properties of the Original FETI Method ... 46
3.2 Generalized FETI........................................ 48
3.2.1 Duality Derivation of the Generalized FETI.............. 50
3.2.2 Derivation of Generalized FETI Algorithm
without Duality........................................ 52
3.2.3 Method Selection for Plate Problems and
Other Generalizations.................................. 57
3.3 The Dirichlet Preconditioner............................ 58
3.4 Interface Formulation of Park et al..................... 59
4. Convergence Analysis....................................... 63
4.1 Preliminaries........................................... 64
4.1.1 Estimates for Pl/Ql Elements............................ 67
4.1.2 Estimates for HCT Elements.............................. 71
4.2 Abstract Analysis Framework............................. 77
4.3 Analysis of the Original FETI .......................... 81
4.3.1 Assumptions............................................. 81
4.3.2 Discrete Norm Bounds ................................... 83
4.3.3 Condition Number Estimate............................... 89
4.4 Analysis of the Method by Park.......................... 90
4.5 Convergence Estimates for Plate Bending................. 92
4.5.1 Assumptions............................................. 93
4.5.2 Discrete Norm Bounds ................................... 95
4.5.3 Condition Number Estimate............................... 99
5. Computational Results ..................................... 101
viii


References
106
IX


1. Introduction
In this thesis, the Finite Element Tearing and Interconnecting (FETI)
method is analyzed. It is one of domain decomposition methods for solving
large systems of linear equations arising from finite element discretizations of
elliptic differential equations.
Elliptic equations often arise in modeling physical phenomena. The
Laplace equation,
A0 = div V0 = / in U. (1.1)
models electrostatic interactions and many other potential problems. Phys-
ical quantities are often governed by systems of equations, as in the case of
linear elasticity which is modeled by a system of equations for the unknown
displacement vector u
div In the isotropic case, the stress tensor a can be written as a(u) = 2jie(u) +
A tr (e(u))S, where A and /i are called Lame constants, 8 is the Kronecker ten-
sor, and e(u) is the strain, e(u) = l/2(Vtt + (Vit)T). An equation describing
a plate bending problem can be obtained from the linear elasticity equation by
a limit process that considers thickness of the plate to be infinitesimally small.
In Kirchhoff-Love model, this yields a fourth order elliptic equation.
In order to set up a well-posed problem, the equations need to be
complemented by boundary conditions. A boundary value problem may be
1


reformulated in weak (variational) form [56]. Two basic kinds of boundary
conditions are recognized. Essential boundary conditions, such as a Dirichlet
boundary condition (f> = 0 on dQ for (1.1), need to be explicitly imposed in the
weak form of the problem. Natural boundary conditions, such as a Neumann
boundary condition dip/dn = 0 dQ for (1.1), are naturally incorporated into
the weak form.
The process of discretization may be based on a Galerkin approxima-
tion. In that case, an approximate solution is sought in a finite-dimensional
subspace of the space in which the weak form is posed. Finding the Galerkin
approximation then requires solving a system of linear equations for the coeffi-
cients of the solution relative to a basis of the subspace. Finite element spaces
are widely used in this context. They are generated by basis functions that are
usually polynomial on each element of a triangulation of the domain and the
supports of basis functions have only small overlaps.
The matrix of the discrete system has several special properties. It is
symmetric and positive definite. Also, due to limited overlaps of finite element
basis functions, it is sparse; the number of nonzero entries in a row of the ma-
trix is a lot smaller than total number of entries in the row. The size and the
condition number of this matrix are affected by fineness of the triangulation.
For quasi-uniform meshes for which the finite element discretization is charac-
terized by the mesh size h of the triangulation, the condition number of the
matrix is of the order p for second order elliptic problems and p for fourth
order elliptic problems.
Two classes of methods for solving the linear system exist, direct and
iterative. Direct methods, usually based on a variant of Gaussian elimination,
2


are common in everyday engineering practice. They are fairly robust; their dis-
advantage lies in their memory requirements and their speed. Gaussian elimi-
nation leads to fill-in and the destruction of the sparse structure of the matrix.
The fill-in can be reduced by re-numbering the variables, but for problems in
3D this helps only to a certain extent. In general, the number of operations
required to solve the system arising from the finite element discretization is
proportional to the square of the number of unknowns or worse.
Since systems arising from discretizations of problems in engineering
practice constantly push available computer resources to the limits, iterative
methods are often preferable for large-scale problems. Their main disadvan-
tage is probably their lack of robustness when too wide a class of problems is
considered. On the other hand, the memory requirements of iterative methods
are typically much smaller than those of the direct methods. Usually only a
small multiple of the number of nonzero entries of the original matrix needs to
be stored. The number of operations can be as low as a multiple of the number
of unknowns, but often the computational cost is not known in advance. The
speed of convergence depends on the condition number of the matrix of the
system. Since the condition number deteriorates with decreasing mesh size,
so does the speed of convergence of simple iterative methods. It is therefore
desirable to construct methods that overcome bad conditioning of the matrix
and that are, if possible, independent of other singular perturbations such as
inhomogeneities or bad Poisson ratios in the linear elasticity problem or, in the
case of plate bending model, the problems arising from the thickness of a plate
approaching zero and the model becoming a fourth order problem. Due to the
elliptic nature of the underlying problem, this sort of stability is often possible
3


to achieve and various algorithms have been proposed.
In this thesis, we will study one of domain decomposition methods.
This class of methods has gained enormous popularity in the last decade. Fol-
lowing a divide-and-conquer idea, domain decomposition methods divide the
original problem into a number of smaller subproblems. Sometimes such a divi-
sion arises from breaking up a complicated geometry. In many cases, though,
it is entirely artificial. The subproblems are easier to solve because of their
smaller size and often parallelism can be exploited. This is especially impor-
tant with the onset of parallel computing.
Domain decomposition methods can be seen from two different points
of view. They may arise from separation of a physical domain into regions,
where a problem is modeled by separate partial differential equations, with the
interfaces between the subdomains being handled by various conditions, such
as continuity. The opposite approach is to see domain decomposition methods
purely as methods for solving large algebraic linear systems arising from the
discretization of PDEs. In this context, the large system is subdivided into
smaller problems, whose solutions can be used to produce a preconditioner for
the large system.
Every domain decomposition algorithm involves two principal issues.
It breaks up the original problem into subproblems, that are solved by some
known method, and it resolves interactions of local solutions, usually by the
means of an iterative method. Many domain decomposition methods have
been developed, originally, for the case of two subdomains. As the complexity
of problems of interest demands splitting the problem into small enough sub-
problems, a multi-domain case proves to be of more importance. Early works
4


demonstrate independence of proposed algorithms on the characteristic mesh
size. However, it has been shown that to achieve independence of the number
of subdomains, a coarse space needs to be introduced. Solving a small coarse
space problem distributes the information about the solution globally, thus
resolving global characteristics of the solution not visible by the subdomains.
In the next chapter, we give an overview of several well-know domain
decomposition algorithms. In particular, we will describe abstract Schwarz
methods, overlapping and substructuring methods and their relationship to
the finite element tearing and interconnecting method (FETI) [27, 32], the
main subject of this thesis. We also briefly discuss a domain decomposition
method based on smoothed prolongation [68, 64],
The remaining chapters, except for Sections 3.4 and 4.4, are based
on papers [50, 51, 53] by Mandel and Tezaur, and by Mandel, Tezaur and
Farhat. In Chapter 3, the derivation of the original FETI method is shown.
FETI tears the computational domain into non-overlapping subdomains and
enforces intersubdomain continuity via Lagrange multipliers applied at the
subdomain interfaces. The Lagrange multipliers are used as the unknowns
and FETI formulation is obtained by solving the saddle point problem for the
Lagrangian. Consistent treatment of subdomain singularities leads to a small
coarse problem which is solved in each iteration of the preconditioned conjugate
gradient method iteration. In Chapter 3, we also present a generalization of
the original FETI algorithm. We describe its particular application to a plate
bending problem.
Chapter 4 is concerned with analyzing the original and generalized
FETI methods. We show that the condition number of the preconditioned
5


FETI method is bounded independently of the number of subdomains and
poly-logarithmically in terms of subdomain size. That is, the condition num-
ber is bounded by (7(1 + log^)3, where H is the characteristic subdomain size
and h the mesh size. We present a complete analysis for decompositions with
crosspoints in 2D and edges and crosspoints in 3D for second order elliptic
problems. Furthermore, we demonstrate that the choice of the Lagrange mul-
tipliers given by Park [59, 58] assures that the condition number estimate can
be improved to (7(1 + log y)2. Finally, we also show that the generalized FETI
converges poly-logarithmically for a biharmonic problem.
In order to illustrate the potential of the generalized FETI method,
Chapter 5 summarizes some computational results. We consider the plate
bending problem on a unit square and demonstrate that the condition number
is almost independent of the number of subdomains and the size of the problem,
as predicted by the theoretical analysis.
6


2. Overview of Related Domain
Decomposition Methods
2.1 Finite Element Approximation, Model Problem and PCG
In this section, we are going to summarize some concepts and algo-
rithms used throughout this chapter and the rest of the thesis. We will consider
the following model problem. Let fi be a bounded Lipschitz domain in IE2 or
IE3. We will study the second-order elliptic equation with Dirichlet boundary
condition
with the matrix [a^-] symmetric, uniformly positive definite and bounded on
fi. The corresponding bilinear form is then given by
and is defined for all v and w in the Sobolev space if1(Q) (the space of gen-
eralized functions with square integrable first derivatives). Let us denote the
L2(0) inner product
Au = f in fi
(2,1)
u = 0 on dVt
where
7


By the divergence theorem (integration by parts), the problem (2.1) can be
written in weak (variational) form: Find u E Hq(Q) such that
a(u,v) = (f,v) (2.2)
for all v in V = Hq(Q), where Hq is the completion of smooth functions with
support in Q with respect to the norm in if1(Q). Since the bilinear form a(u, v)
is symmetric, continuous and elliptic (coercive) [56], the problem has a unique
solution by the Lax-Milgram theorem.
Let us now use the standard Galerkin approximation. We look for
an approximate solution of (2.2) in a finite-dimensional subspace 14 (Q) of the
space V. The Galerkin approximation is the solution of the following problem:
Find uh E Vh (fi) such that
a(uh,vh) = ( f, vh), (2.3)
for all vh E 14 (Q). This problem leads to a system of linear equation when a
basis for Vh is chosen. Let i be a basis for 14(0). Assuming
N
uh =
i=1
(2.3) becomes
N
YaiVuVj)ui =
i=1
j = 1,...,N.
Matrix K = [a((pi,(pj)\ is called the stiffness or Gramm matrix and
the right hand side vector [(/, (pj)\ is called the load vector. The vector of
unknowns, [iq], is also called the vector of degrees of freedom.
For a wide class of approximation spaces 14, uh is a good approxima-
tion of u. In the case of finite element methods, the construction of suitable
8


spaces Vh relies on triangulation, which splits the domain fi into small disjoint
regions of simple geometric shape, such as triangles or quadrangles in IE2 or
tetrahedrons or hexahedrons in IR3; h refers to the characteristic size of these
regions. Under certain assumptions which prevent the triangulation from de-
generating (e.g, the size of the angles in triangles is bounded from below so as
they do not become too sharp), the finer the triangulation, the closer a finite
element Galerkin solution to the exact solution, [14],
Functions in a finite element space Vh usually arise from a polynomial
interpolation on the elements of the triangulation. Every polynomial defined on
a given region is uniquely determined by its values and perhaps also the values
of its derivatives at some nodal points, usually the vertices of the region. Each
function in Vh is thus determined by a set of values at nodal points, so called
degrees of freedom. Simple examples of finite element spaces include spaces
formed by continuous functions linear on triangle regions in IR2 or tetrahedrons
in IR3 or bilinear functions on quadrangles. These spaces are being referred
to as Pi and Qi, respectively, and their entries are uniquely determined by
their function values at the vertices of the triangulation. Another type of finite
elements we will encounter later in this thesis is HCT elements, which are C1
continuous and determined by the function values and the values of the first
derivatives at the vertices of triangular regions in IR2.
The standard basis for a finite element space is the one in which each
of the basis functions has exactly one degree of freedom equal to 1 and the rest
0. The unknowns in the linear system arising from the discretization are di-
rectly the degrees of freedom of the Galerkin approximation (hence the name).
The overlaps of the supports of the basis functions are small, which causes
9


the stiffness matrix to be sparse. Due to the properties of the bilinear form,
the stiffness matrix is also symmetric and positive definite. Preconditioned
conjugate gradient (PCG) method is therefore often considered for solving the
system. A detailed description and analysis of PCG (and other methods) can
be found, e.g., in Golub and Van Loan, [38]. Here we give only a brief summary.
The conjugate gradient method is based on the observation that the
solution of Ax = b is the only minimum of
Q(x)
xTb.
A sequence of vectors which converges to the solution of the system can be
obtained by setting
xk+l = xk + akPk>
where Pk is a suitably chosen direction and vector + a^p^ minimizes Q on
the line Xk + apk- In PCG, vectors pu are obtained by A-orthogonalization of
residua, which leads to the following algorithm:
Algorithm 1 (Conjugate Gradients) Given x0, set
r0 = b Axq.
For k = 1,... do:
Pk
Pk
k
Xk
rl-in-i
n-1 + -p-Pk-1 (Pi = r0)
Pk-1
0k
plApk
Xk-1 + O-kPk
Xk Xk-i otkApk
10


Among iteration methods the conjugate gradient method is rather exceptional.
Assuming roundoff errors are not present, it gives the exact solution of the
system after a finite number of iteration, which equals at most the number of
unknowns. For large systems, however, performing so many iterations would
be impractical. Fortunately, much smaller number of iterations is often needed
to get a good enough approximate solution.
This is not true if the condition number of A is large. The performance
then can be improved by introducing a preconditioner M = ETE, where E is
some nonsingular matrix. We point out that M defined in this way is symmetric
and positive definite. Preconditioner M must have two essential properties:
matrix E^TAE^1 is better conditioned than A,
system Md = e can be easily solved.
E^T denotes the transpose of the inverse of E.
The preconditioned conjugate gradient (PCG) method is the conju-
gate gradient method applied to the system
E-TAE-1{Ex) = E~Tb.
This leads to the following algorithm.
Algorithm 2 (Preconditioned Conjugate Gradients) Given
xq, set
r0 = (b- Axo).
For i = 1,... do:
zk-i = M lrk-1
A = aT-iAc-i
11


Pk = zk^i + -p-Pk-1 (pi = Zq)
Pk-l
_ (h
ak jr-
PkAPk
xk = xk-1 + akPk
xk xk 1 Otk-Apk
We point out that matrix M 1 is not formed explicitly and, instead, a
linear system with matrix M is solved in each iteration of PCG method. Stan-
dard preconditioners used in practice include diagonal scaling and incomplete
factorizations (e.g., Cholesky) [38].
2.2 Abstract Schwarz Methods
Abstract Schwarz methods take their name after Herman Schwarz, a
German mathematician of the past century, who used an alternating method
to construct harmonic functions on regions with non-smooth boundaries. They
provide an abstract framework which allows analysis of many different domain
decomposition methods. Here we briefly describe the general abstract algo-
rithm and refer to literature for details on analysis and practical/non-practical
choices of particular methods. This presentation is based on [6] and [23].
Let V be a Hilbert space with an inner product (,). We also con-
sider the inner product a(-, ) implied by the bilinear form of the variational
formulation of our problem: Find u V so that
a(u,v) = (f,v) Vt> e V.
Let Vi, % = 0,..., n be closed subspaces of V that form a decomposi-
tion of V, i.e.
V
W + W + W + ... + fy
12


In addition to the inner product implied by the the bilinear form a, we consider
the inner products a,i(-,-) implied by some bilinear forms a* defined on V* x
Vi, i = 0,... ,n. The Schwarz additive method can be then written as follows.
Algorithm 3 (Additive Schwarz) Given an initial approximation
u, start with k = 0 and do:
1. Find iVi E V such that
adwi, Vi) = (/, Vi) a(uk, i>*), e V.
2. Deline the next iterate as
n
uk+1 = uk + Yswi-
i=0
3. Set k = k + 1 and go to 1.
The Schwarz multiplicative method is then defined by the following
algorithm.
Algorithm 4 (Multiplicative Schwarz) Given an initial approxi-
mation u, starting with k = 0 do:
1. For / = 0.....n.
(i) find iVi E Vt so that
Oi(m, Vi) = (/, vt) a{uk+l/{n+l),Vi), Mvi E Vi}
(ii) set
yk+(i+l)/(n+l) yk+i/(n+1) +
2. Set k = k + 1 and go to step 1.
The algorithms above can be combined to form hybrid methods. They
are usually used as preconditioners in the conjugate gradient method. Various
13


choices of the subspaces \\ yield a wide variety of different algorithms. Since
the multiplicative method above is not symmetric, it is usually symmetrized
by considering the sequence of subspaces in forward and reverse order; that is
F0, F,.....F, Vn, V. V ,. .... F- F0.
Each sub-step of an additive and a multiplicative Schwarz method
can be interpreted as an approximate projection of the error. Denoting P, the
operator corresponding to the i-th sub-step of an additive method, that is the
operator that satisfies
a,i(PiW,v) = a(w,v) Vu e F,
one step of the additive method satisfies
n
i=0
where ek = u uk is the error of the k-th approximation. Similarly, the
reduction of the error of one step of the multiplicative method is governed by
the equation
ek+1 = (/ Pn)...(/ P\){I Po)ek.
We refer to [25] and references included there for abstract analysis.
Here we only summarize the main results. They are due to a number of authors:
Dryja and Widlund [25, 24]; Nepomnyaschikh [55]; Bramble, Pasciak, Wang,
Xu [9, 70]; Lions [44]; Bjprstad and Mandel [6], and others.
Let Cq be the constant such that for all u E V, there exists a repre-
sentation a = Y'i u 11 > ui ^ Vi, and
n
^2ai(ui,Ui) < C0 a(u, u).
i=0
14


Let B = [by] be the matrix of strengthened Cauchy-Schwarz coefficients,
Ia(vi,Vj)\ < bija(vi,vi)1/2a(vj,vj)1/2 e V^vj G Vj,i,j = l,...,n.
Furthermore, let u be the constant such that
a(ui, Ui) ^ ui a,i{u,h Ui) 'tfu-i G Vj, i 0,..., n.
Then the smallest eigenvalue of the operator Y,Pi of the abstract additive
method used a preconditioner in the conjugate gradients method is bounded
by 1/Co and the largest by u(p(B) +1). The bound on the condition number of
the symmetric multiplicative method used as a preconditioner in the conjugate
gradients method is C0(l + 2ujp(B))/(2 cD), where Cj = max(l,o;).
The space V0 is usually regarded as a coarse space. It has been shown
that without the presence of the coarse space, the condition number of an
abstract Schwarz method grows with the number of subspaces. The coarse
space ensures a mechanism of global exchange of information otherwise lacking
in decompositions involving a large number of subdomains.
2.3 Overlapping Domain Decomposition
The first known domain decomposition method is due to Schwarz [63]
and is known as the Schwarz alternating method. It divides the domain fi,
on which the problem is defined, into two overlapping regions fii and fi2 and
performs a multiplicative Schwarz type algorithm.
Let us consider the model problem (2.1)
Au = f in fi,
u = 0 on
15


Let fii and fi2 be overlapping subdomains of fi, fii U fi2 = fh Furthermore,
let I'i = (){11 fl ilj and T2 = dfl2 H (-h-
Algorithm 5 (Alternating Schwarz) Given an initial guess m2 for
the values on fi2, for k = 1, 2, ... do:
1. Solve Auk = / in f^!,
uk = 0 on \ r
uk = |n on Tj
2. Solve Auk = / in fl2,
uk = 0 on cKl2 \ T:
uk _ r,.k\ i|r2 on r2.
To be precise, f|ri and if2|r2 stand in the algorithm above for the
trace of u\ on Id and uk on r2, respectively, rather than the simple restrictions.
Approximate solutions of the model system can be assembled as
uk = nt in tt2,
uk = uk in fii \ fl2.
The variational formulation of the algorithm above is due to Lions
[44], It is the multiplicative abstract Schwarz algorithm using the subspaces
Vi and V2 of V = Hi (fi) formed by functions that vanish outside of fii and
fl2, respectively. The bilinear forms ai and a2 are simple restrictions of a.
16


Choosing some other bilinear forms ai and a2 can model the situation when
the problems in step 1. and 2. of the algorithm above were to be solved only
approximately. An additive version of the algorithm has been described by
Dryja [20] and Matsokin and Nepomnyaschikh [54],
In the discrete, multidomain case, overlapping methods often start
with a non-overlapping partition of Q into subdomains (substructures) fij,
i = 1,,NS, the diameter of the subdomains being of order H. Then, each
subdomain fij is extended to a larger region fi'. It is assumed that neither
Vti nor fi' cuts across any of the elements. The overlaps of subdomains fi',
i = 1,..., Ns are said to be generous, if the distance between boundaries of fij
and fi' is greater than some fixed fraction of H. The subspaces Vt, i = 1,..., n
of the abstract Schwarz method are then chosen as restrictions of the finite
element space 14(fi) to the subregions fi'. That is,
\; = C,n //>( Approximate or exact solvers then can be chosen on the subspaces and addi-
tive, multiplicative or hybrid two-level Schwarz methods can be used. As far
as parallelization is concerned, additive versions clearly gain an edge, because
they can be easily parallelized. Multiplicative versions always involve sequen-
tial steps. Some parallelization is possible when mutually non-overlapping
subdomains fi' are grouped together. This process is often called coloring of
subdomains.
Without a coarse space, the condition number of such methods is of
the order 1 /H2 [24], where H is the diameter of the subdomains. Dryja and
17


Widlund in [24] propose to use the space of continuous piecewise linear func-
tions on the coarse mesh defined by the subdomain Q* as the coarse space V0.
They prove that, in the case of generous overlap, the condition number relevant
for the conjugate gradient iteration is uniformly bounded. Mandel shows in
[47] that a similar method by Cowsar [16], which uses discrete harmonic ex-
tensions determined by piecewise constant values on interfaces of subdomains,
leads to a poly-logarithmic bound on the condition number.
In their later paper [21], Dryja and Widlund mention, that their nu-
merical experiments indicate that the convergence rate is often satisfactory
even for small overlaps. Running time is often smallest when the overlaps are
at minimum. The number of conjugate gradient iterations is higher in such
a case, because the condition number deteriorates, but this is being compen-
sated for by the fact that the local problems are smaller and therefore cheaper
to solve. Also, they are better conditioned and, if they are solved by iteration
solvers, the rate of convergence is faster.
2.4 Substructuring Methods
2.4.1 Poincare-Steklov Operators and Schur Complement
Substructuring methods borrow their name form structural engineer-
ing. It was in the context of structural engineering that several substructuring
algorithms have been pioneered. We shall concentrate on substructuring meth-
ods on interfaces of subdomains.
Let us consider the model problem (2.1) and let fi be divided into two
non-overlapping subdomains fii and fi2 (Figure 2.4.1). Our goal is to solve the
model problem only on the subdomains. Let T = dQ\ fi <9fi2 be the subdomain
18


interface. It is well known that the solution u E Hq(Q) exists for all f E L2(0).
The solution satisfies continuity conditions on the subdomain interface: conti-
nuity of fluxes of the solution and continuity of the solution. Let us introduce
a few notations. For k = 1,2, we define the space Hjj(Qk) as the subspace of
such that the functions in H}j(Vtk) vanish on dVtkndVt. Let Hpl(Q,k) be
the dual of H}j(Vtk). Let jk be the trace operator mapping functions in H}j(Vtk)
to their traces on T. Let hH2(T) be the fractional order Sobolev space on T
consisting of traces of functions in HlD(Vtk) and let (hIq2{T))' denote its dual.
Using the continuity of fluxes, we will split the problem into two
subproblems, k = 1,2:
Auk = fk in ttk
n ([aijjVuk) = (-1 )k+l(f on T
uk = 0 on dkl fl dQk.
The vector n is the normal vector to T oriented, for example, from fii to fi2-
In variational terms, these are the problems of finding uk E = 1,2
19


such that
a.k(uk,Vk)= g*Vk+ fkVk, Vufe 6 HlD(Qk), (2.4)
JV J
where
ak{v,w) = Y, /
in
i,j=1
dv dw
hik dxi dxj
a,i
At this point it would be possible to design an iterative method working on the
subproblems using the continuity conditions. Usually, however, we first reduce
the problem to a problem on the subdomain interface using Poincare-Steklov
operators.
First, we will look for the unknown Neumann data g* on T. We define
the Poincare-Steklov operators Qk : (hIq2{T))' ifdt)2(r),£; = 1,2 by
Qk9* = IkUki (2.5)
where, for g* e is the solution of (2.4) with fk = 0. Such uk is
the harmonic function satisfying the Neumann condition given by g*. In other
words, the Poincare-Steklov operator maps the Neumann boundary condition
into the corresponding Dirichlet boundary condition
duk
Qk ' ----^ TWfc-
on
Furthermore, we define Rk : H]J1(Qk) Hq12(T), k = 1,2 by the equation
Rkfk
where for fk e L2(f2fe), uk is the solution of (2.4) with g* = 0. In terms of the
Poincare-Steklov operators, the problem is to find the solution g* such that
(Qi + Q2)g* = R2f2-Rifi
(2.6)
20


That is to find the Neumann data g* on T such that the traces of the solutions
Uk, k = 1, 2 of (2.4) coincide on T.
The second possibility is to enforce continuity of the solution on the
boundary a priori. This yields a dual formulation to (2.6). Consider the
Dirichlet problems, k = 1,2,
Auk = fk in klk
uk = 9 on r
uk = 0 on dn n dttk
We will look for the unknown Dirichlet data g on T such that the fluxes are
continuous for the solutions Uk,k = 1,2 of the problems above. Using the
Poincare-Steklov operators, this problem can be written as
(Or1 + orOs = or'-Ri/i + Q^Rih- (2.7)
The equations above and iterative methods for solving them have
been studied, for example, by Agoshkov [2] and in a mixed-method setting
by Glowinski and Wheller [36]. The paper by Bakhvalov and Knyazev [4] is
concerned with highly discontinuous coefficients between the subdomains.
As the case of two subdomains is hardly of practical importance, we
conclude this section by introducing a multi-domain discrete analog of (2.7).
An analog of (2.6) is described in the next chapter.
Let us recall now the concept of Schur complement.
Definition 6 (Schur complement) Let A be an n x n matrix and a, fi C
{l,...,n} be index sets. We denote A(a, (3) the submatrix that lies in the
rows of A indexed by a and the columns indexed by /3, and A{a!, (Â¥) the sub-
matrix given by deleting the rows indexed by a and the columns given by (3.
21


Let A(a, a) be nonsingular. Then the matrix
A(ex, ex ) A(ex, ex^jA^ex^ o;) A^cXj ex )
is called the Schur complement of A(a, a) in A.
We note that the Schur complement is the matrix that arises when
eliminating the unknowns x{a) from the equation Ax = b. Then,
b(ex') A(ex', ex)A(ex, a:)_1&(a:) =
(A(ex ex ) A(ex cx^)A(cX] ex'} A^ex^ ex }}x(ex ).
We also observe that the Schur complement S of A(ex, a) in A has the property
(Sy,y)
inf (Ax,x).
xmn,x(r)=y
(2.8)
Let Q be a domain in IR2 or IR3 decomposed into Ns non-overlapping
subdomains fii, , &ns- Let iq be the vector of degrees of freedom for
subdomain Qj corresponding to a conforming hnite element discretization of
the second order elliptic problem (2.1) defined on fi, such that each subdomain
is a union of some of the elements. Let u,. !\ ,. and /*, be the vector of degrees of
freedom, the local stiffness matrix, and the load vectors, respectively, associated
with the subdomain Q*.
Let Li denote the zero-one assembly matrix mapping the subdomain
degrees of freedom iq into global degrees of freedom u, that is iq = Lju. The
stiffness matrix then is
Ns
K = jTLlKlLj
i=1
and the load vector
Ns
/ = ££./.
i=1
22


The problem to be solved,
Ku = /,
can be reduced to an interface problem by splitting the degrees of freedom
into interface and interior degrees of freedom. Let us assume that the interior
degrees of freedom are listed first. The subdomain stiffness matrices Kt and
the restriction matrices Lt can be split accordingly:
Ui
Ui Ki kf
Ki =
Ui ki k
Li \Li, Li
By eliminating the interior degrees of freedom, we obtain the Schur comple-
ments of Ki in Ki
~ ' (2.9)
Si = Ki^ KjKS K'
The problem then reduces to the interface problem
Su = f
for the interface unknowns u with the global Schur complement
Ns
s = y,l1s1zj
i=1
and the interface right hand side
Ns
f = Y^Li{fi- KiKf'fi).
i=1
In most cases, explicit computation of Si would be too expensive, but
for gradient-descent methods, only evaluation of the action of Si is necessary.
This evaluation can be performed efficiently by factorizing and it corre-
sponds to solving a Dirichlet problem on every subdomain. The subdomain
problems can be solved in parallel.
23


Let us demonstrate that the action of S) can be evaluated by solving
a Dirichlet problem on the subdomain Q*. Consider the discretized problem
on the subdomain
Ki kT Ui 0
Ki ki Ui fi
where the Dirichlet data tq is known. From the first equation, we find -tq =
l\; ll\jd,. This the solution of the Dirichlet problem given by the data tq.
Then, we substitute into the second equation and obtain
{Ki ~ kikr'kDiii = SiUi =
Many preconditioners for the reduced problem are based on compu-
tations of the action of Sf1. This can be evaluated by solving a Neumann
problem on the subdomain, for unknown tq and -tq and given rq,
k kT i Ui 0
k k Ui 9i
where rq corresponds to the Neumann data. Then, tq = Sil9i- Indeed, from
the first equation, we find -tq = lI\ jH,. Substituting into the second
equation, we find
{Ki ~ kik^kT)ui = SiUi = g,.
We note that the Schur complement Si is the discrete analog of the
inverse of the Poincare-Steklov operator Q^1 : u|r |^, and the inverse of the
Schur complement is the analog of Qit : u|r (2.5).
One of the advantages of the reduced problem is that it improves the
condition number of the original problem. For our model problem, assuming
24


the Poincare-Friedrichs inequality \\u/||^ 2 n. < cp\v/l\\2Cl. holds, the subdomain
stiffness matrix Kt satisfies
r\U,h\\2 < < r\Uih\\2
HI ui 111,2,Hi 22: ai 22: Ly\\ai 111,2,0,)
where u/ is the finite element function corresponding to the vector of degrees
of freedom iq. Then, the subdomain Schur complement S) satisfies (cf. Lemma
26 and Lemma 33)
c\\u
h ||2
i 11-5,2,(90,
< U
TSm, < C\\uh
t 11 ^ ,2,(902;
This implies that, for a second order problem, the condition number of the
reduced problem is of order 1/h for triangulations of characteristic mesh size
h, while the condition number of the original problem is of the order 1/h2.
2.4.2 Various Preconditioners for the Reduced Problem
Instead of solving the reduced problem directly, iterative methods
such as the preconditioned conjugate gradient method can be applied in which
only matrix vector products are required. Since each iteration is quite ex-
pensive (it requires solving a Dirichlet problem on each subdomain), efficient
preconditioning is important to keep the number of iterations small.
Diagonal or block diagonal preconditioning of the reduced problem
requires knowing the diagonal of S. Instead of computing it directly, Chan
[13] proposes a boundary probing technique to construct an approximate
diagonal by evaluating actions of S on carefully selected vectors.
For geometrically simple subdomain in 2D, Bramble, Pasciak and
Schatz in [10] propose a method for preconditioning the original problem by
25


first splitting the functions in the finite element space into discrete harmon-
ics on each subdomain and functions which vanish on the boundary of each
subdomain. This is equivalent to eliminating the interior degrees of freedom
as in the reduced system. Then the space of discrete harmonic functions is
divided into functions linear on interfaces and functions whose values are zero
at subdomain crosspoints. The preconditioner is then based on a bilinear form
defined on the space splitting. The condition number of the method is shown
to be bounded by (7(1 + log H/h)2. A similar idea is exploited in [3] for the
p-version of finite elements and a bound (7(1 + logp)2 is obtained.
A different algorithm [8] by Bramble, Pasciak and Schatz is described
for the case of two subdomains fii and fi2 separated by the interface T. They
propose to precondition the original problem in the following way: First split
the local solution

homogeneous equation with the zero boundary condition on <9fi2 and a discrete
harmonic function A Dirichlet problem is solved on fi2 to obtain ipP. This is followed by a so-
lution of a mixed Neumann-Dirichlet problem on fii using the Neumann data
obtained by solving the problem on fi2. Finally another Dirichlet problem
is solved to obtain ipH, a discrete harmonic function, that satisfies the corre-
sponding homogeneous equation. This is sometimes called Neumann-Dirichlet
decomposition. This algorithm has also been described in [24] and it is in fact
preconditioning the reduced system by either Af1 or A71. The inverses of the
Schur complements need not be explicitly computed; the action corresponds to
solving a Neumann (or mixed Neumann-Dirichlet) problem on one subdomain
and then extending the solution to the other subdomain by solving a Dirichlet
26


problem with the Dirichlet data found from the Neumann problem. The idea of
the Neumann-Dirichlet preconditioner can be extended to multiple subdomain
case when there is a red-black ordering of the subdomains.
Several methods solving an equation with the global Schur comple-
ment are proposed by Mandel [47]. The methods are set in the space of all dis-
crete harmonic functions on the union of subdomain interfaces T = U*50* 5Q.
The methods are hybrid Schwarz methods using the coarse space Vo in a mul-
tiplicative fashion and other spaces are treated additively. A method of this
kind can also be interpreted as a two level variational multigrid. The methods
use V* that are associated with subdomains or globs. A glob is a vertex, an
edge that does not contain its endpoints, or a face of a subdomain interface.
Piecewise linear functions defined on the subdomain triangulation or a piece-
wise constant space based on glob-wise averaging of subdomain values is used
as the coarse space. Poly-logarithmic bounds are obtained.
Many other choices of the coarse space are discussed by Dryja, Smith
and Widlund [23]. Their exhaustive investigation comprises vertex based coarse
spaces, based on piecewise linear functions on substructures used as elements,
wire basket algorithms, that use averages on substructures, and face based
algorithms.
2.4.3 The Neumann-Neumann and Balancing Domain
Decomposition Methods
The balancing domain decomposition (BDD) is based on the so called
Neumann-Neumann preconditioner that preconditions the reduced problem
by a weighted sum of inverses of the local Schur complements. It is called
27


Neumann-Neumann because it corresponds to solving Neumann problems only
on subdomains (as opposed to the Dirichlet-Neumann preconditioner, that uses
both Dirichlet and Neumann problems). Since this is a preconditioner in a way
dual to the method that is the subject of this thesis, we describe it in a little
more detail. To describe the preconditioner, weight matrices Dh i = 1,..., Ns
satisfying the decomposition of unity
Ns
i=1
are used. A simple choice for D-i is a diagonal matrix with the diagonal elements
being the reciprocals of the number of subdomains the degree of freedom is
associated with.
The Neumann-Neumann method used as a preconditioner of the prob-
lem is as follows:
Given the residual r, distribute it to subdomains r* = DfLJr, solve the local
problems S;ti; = and average the results \l lr = l-jlh'U [7, 43].
The drawback of the Neumann-Neumann preconditioner is that it
lacks a mechanism of distributing the error globally. This has been resolved
by adding a coarse space to the problem. The resulting method is called the
balancing domain decomposition. Let Z;t be the matrix with linearly indepen-
dent columns that generate the kernel of l\ ,. Im Z% = Ker K,i. If Kt is regular,
Zi is a void matrix. The balancing preconditioner is as follows:
1. Balance the original residual r by solving the auxiliary problems
Ns
ZjDjZf (r S £ L3D3Z3v3) = 0, i = I.......Y
j=i 2
2. Distribute the balanced residual to the subdomains and find a solution of
28


the local problems
Ns
SiUi = Dj LJ(r S ^2 LjDjZjVj), i = l,...,Ns
3=1
3. Balance by solving the auxiliary problems
Na
4. Average the result on the interface
Na
M 1r = Y, I-,!),(, + Zm)
i=1
The BDD is proven to be independent of number of subdomains,
and its condition number is independent of jumps of coefficients between sub-
domains with appropriate weight matrices D*. It can be written also as an
abstract additive Schwarz algorithm [24, 41].
2.4.4 The Neumann-Neumann Domain Decomposition
Algorithm for Plates and Shells
Let us introduce the Kirchhoff-Love model of plate bending following
[41]. We consider a plate occupying a domain in IR2, which is clamped on the
part of boundary <9fi0o and simply supported on <9fi0 \ <9^oo- The Kirchhoff-
Love model plate model characterizes the vertical displacement u E V of the
plate as the solution of the variational problem
where the bilinear form on the right hand side is symmetric, continuous and
coercive
a(u, v) = F(v), Vv E V,
(2.10)
29


the functional on the right hand side is given by
F{v) = f fv + [ ma^T + f 90'
Jn Jdn-dQoo on Jan-an,,
The space of kinematically admissible fields is
f)v
V = {ti6 H2(Q),v = 0 on 8Qq,v = = 0 on <9fi0o}-
on
In the definitions above / is the density of vertical forces, mg the density
of flexion moments applied on the part of the boundary where the plate is
free to rotate, and g is the density of vertical boundary loading. The symbol
denotes a tensor product, e(9) = |(V0 + (V0)T) is the curvature tensor,
9(u) = Vu represents the in-plane notation of the plate and K is the plate
flexural stiffness. K is symmetric, elliptic and continuous in the sense that
e(9(u)) : K : e(9(u)) > ct3\e(9(u))\2
e{9{u)) : K : e(9(v)) < Ct3\£(9(u))\\£(9(v))\,
where t is the plate thickness. For a simple case of an isotropic plate made of
an homogeneous material with Young modulus E and Poisson coefficient u, it
is given by
Ef3
£(9(u)) : K : e(9(u)) = -----r-((l v)V2u : V2v + vAuAv).
12(1 v2)
Since V is a subspace of H2(Q), the appropriate finite element spaces are C1
continuous. Examples of such elements include the discrete Kirchhoff triangle
(DKT) and HCT elements.
The Neumann-Neumann and BDD preconditioner do not perform well
for the plate bending problem [41]. The condition number estimate of BDD is
30


based on the estimate [45]
k < sup
Y^iVs II rTy^iVs r.n, 112
2^=1 IWj 2si=l J-'iJ-Siai\\Sj
E
^ lb,.II2
i=l ll'OllS,
: u.
Ker Si, SjUi Im Z,:
(2,11)
At crosspoints of subdomains, the vectors in the estimate are constructed from
contributions several subdomains different from those that share an edge lead-
ing to the crosspoint. This leads to a discontinuity that is inappropriate for
a fourth order problem. The BDD method for plates [42] avoids this prob-
lem arising at subdomain crosspoints by enhancing the coarse space of the
balancing domain decomposition algorithm. The coarse space is again Im Z%,
where
Zi \%il........%irii i Vil j Uirrii ]
{xn,...,Xirii} is a basis for Ker Si and for each crosspoint j = 1,,nij of
the subdomain yij is the solution of the problem Sit/ij = e%j, with the
vector corresponding to the unit normal load at the crosspoint j. With this
choice, the normal displacement component of the vectors Ui coming out of
the coarse space problem is zero since SiUi yij implies that Ui Sit/ij = %.
Then, the supremum (2.11) is taken over functions with zero at endpoints of
the edges, which makes it possible to prove a poly-logarithmic bound. For
another approach, imposing zeros at crosspoint directly, see [41].
2.4.5 Lagrange Multipliers and Poincare-Steklov
Operators
Following [19], we show how Lagrange multiplier approach to en-
forcing solution continuity is related to interface formulations using Poincare-
Steklov operators defined in Section 2.4.1. The discrete multi-domain case will
be treated in the next chapter.
31


We consider the problem from Section 2.4.1 and the notation intro-
duced there. We reformulate our problem as a constrained minimization prob-
lem: Find the solution (u\,u2) E Hjj(Qi) x H}j(Vt2) that minimizes
subject to the condition 7iU\ = 72m2. Then, for each (Ui,u2,g*) E 1) x
The Lagrange multiplier A* is used to to enforce the continuity of the solution
on the boundary. The rest of the Lagrangian is the usual quadratic functional
implied by the weak form of the problem. The critical points (-17, u2, A*) of the
saddle point problem of A now must satisfy the variational equality
This is the equation (2.6). It shows that solving the Lagrange multiplier for-
mulation is equivalent to finding Neumann interface data on T (cf. [18]).
The paper [19] uses the Lagrange formulation to introduce finite ele-
ment spaces of Lagrangians of small dimension per interface for regular meshes.
This can reduce the size of the problem substantially, but it is restricted to reg-
ular meshes. The space of Lagrangians can be chosen as the restriction of the
finite elements space on the subdomains to the interface. This is similar to the
approach taken by FETI as explained in the next chapter.
HlD(Vt2) x (#oo2(r))', we define the Lagrangian
for all (t>i, v2, n*) E if|)(fii) x HlD(Vt2) x (H^fT))'. This problem has a unique
solution [19] which is the solution of the problem
(Q1 + Q2)A* = Rah R1.f1
32


2.5 A Two Level Method Based on Smoothed
Prolongation
The two-level method described in [68, 66] develops a simple abstract
framework based on the concept of smoothed tentative prolongator introduced
in [65]. The tentative prolongator is derived from a system of nonoverlapping
subdomains. As opposed to the previously described methods, the union of all
subdomains Q*, i = 1,..., Ns does not cover whole fi. Instead, there is a layer
one element wide between each two subdomains.
Our algorithm can be written as a variational two-level multigrid
with a special choice of components. We will first describe components of the
method and abstract assumptions that ensure coarse space size independent
convergence. Let us consider the system of linear algebraic equations
Ku = /,
where K is an n x n symmetric positive definite stiffness matrix arising from
a discretization of a second order elliptic PDE, for example (2.1).
Interpolation from the coarse space to the fine space is represented by
an operator MP, composed from a tentative prolongator P and prolongator
smoother M. Let us denote
Km = M2K, M' = 1- -Km,
P
where u e (0, 2) is a given constant and p(KM) is an estimated upper bound
for p(KM), the spectral radius of KM The following assumption specifies
requirements on M, P and p(Km)- are needed for proving coarse space size
independent convergence.
33


Assumption 7 There exist positive constants Ci, C2 independent of n and
Ns and a positive constant Cd(Ns,ti) such that
(1) For every u E IE", there exists v E IE^" such that
||m Pv\\ < ClCD{Ns,v)p-ll2{K)\\u\\K. (2.12)
(2) The prolongator smoother M. is symmetric, commutes with K, p(M) < 1,
and
p(KM) < p(KM) < ClC^(Ns,n)p(K). (2.13)
We also need to define a pre-smoother SM and post-smoother SM';
SM(u,f) = Mu + Nf,
and analogously for M'. The matrix N is chosen so the smoother is consistent
with the system Ku = f. The smoothers SM and SM' are relaxation operators
used to smooth the approximate solution corrected by an interpolated coarse
level error. The algorithm can now be written as follows.
Algorithm 8 Given an initial guess u, repeat until convergence:
(1) u t Sm(u, /),
(2) solve (PtMKMP)w = PTM(Ku /),
(3) u { u MPw,
(4) u { (u, f).
Post process u SM(u, /).
The tentative prolongator P may be defined as P = Dl/2P.j where
D = diag(A) and
Pij = 1, if the node corresponding to tp belongs to subdomain Qj,
= 0, otherwise.
34


The most practical choice of M is a polynomial in K. The choice described
below is a polynomial in K that attempts to minimize p(Km) given an upper
bound on p(K) and a chosen degree of the polynomial. For certain degrees,
this polynomial can be found explicitly by the following recursive algorithm:
Algorithm 9 Let p be the estimate of p(K) satisfying p(K) < p < Cpp(K),
with a given constant Cp. Set pi = A0 = K and for i = 1, 2,... do:
(1) Ai = (I
(2) d/, = 1 It y(/ Ipj 1.1()
Notice that deg(Mj) < For a 2D problem, we choose the prolongator
smoother M = Mk, where
deg( Mk i) > qNlj2 > deg(Mk), (2.14)
where q e (0,1] is a given parameter and Nes the average average number of
degrees of freedom per subdomain, and Mk are the polynomials constructed
by the algorithm above.
With the choices of the tentative prolongator and prolongator smooth-
er described above, the convergence of our method can be shown to be indepen-
dent of the meshsize, Nes, inhomogeneities between subdomains, and boundary
conditions [68]. It thus overcomes one of the main disadvantages of standard
two-level multigrid methods, suffering from the dependence of the rate of con-
vergence on the size of the coarse space. On the other hand, its computational
complexity is lower in comparison with domain decomposition methods using
direct local solvers. This is true even for simple direct solvers as internal solvers
in our method (which are preferred, as iterative solvers tend to have negative
impact on the robustness of the method).
35


Application of the method to solid problems is treated in [64] and
similar results are obtained. The method is provably robust with respect to
jumps in coefficients [67].
36


3. Finite Element Tearing and
Interconnecting (Derivation and Extensions)
3.1 Original FETI
The FETI (Finite Element Tearing and Interconnecting) method is
a domain decomposition algorithm derived from a hybrid variational principle
and designed for the iterative solution of systems of equations arising from the
finite element discretization of self-adjoint elliptic partial differential equations.
It was developed in [28, 27, 33], and also discussed in detail in monograph [34],
In this method, a given spatial domain is torn into non-overlapping
subdomains where an incomplete solution of the primary held is first evalu-
ated using a direct solver. Next, intersubdomain held continuity is enforced
via Lagrange multipliers applied at the subdomain interfaces. This gluing
phase generates a smaller size symmetric dual problem where the unknowns
are the Lagrange multipliers, and which is best solved by a preconditioned con-
jugate gradient (PCG) algorithm. This idea is related to the hctitious domain
method where the Lagrange multipliers enforce boundary conditions as in Dinh
et al. [17].
In contrast with other related domain decomposition methods using
Lagrange multipliers as unknowns [36, 62], the FETI method distinguishes itself
with the treatment of the null spaces of the subdomain stiffness matrices (rigid
37


body modes) associated with the so-called floating subdomains, i.e., subdo-
mains without a sufficient number of essential boundary conditions to prevent
the local stiffness matrix from being singular. Resolving the rigid body modes
leads to a small coarse problem that is solved in each PCG iteration. This
is an added complication, but also a blessing. Farhat, Mandel, and Roux [32]
have shown numerically, and proved for the FETI method without precondi-
tioning, that the auxiliary problem plays the role of a coarse problem, namely,
it causes the condition number to be bounded independently of the number of
subdomains. In [26], Farhat, Chen, and Mandel extended to time-dependent
problems, which lack the naturally occurring coarse problem.
The FETI method is in a sense dual to the Neumann-Neumann meth-
od with a coarse problem, developed by Mandel under the name Balancing
Domain Decomposition [45] (BDD), which we described in Chapter 2 and which
is based on an earlier method of de Roeck and Le Tallec [61]. A modified
method was analyzed by Dryja and Widlund [25]. It should be noted that
while the underlying ideas of FETI and BDD are in a way dual, FETI is not
the BDD method applied to the dual problem.
3.1.1 Problem Setting and Assumptions
Let us first introduce some notation used throughout this chapter and
the rest of this thesis. For u, v E IR", the inner product (u,v) = uTv serves also
as duality pairing. The £2 norm is denoted ||u|| = (u/u)1/2. For a symmetric
positive semidefinite matrix A, ||u||u = (/kpu)1/2 the induced seminorm. This
is a norm if A is positive definite. The superscript + denotes pseudoinverse,
defined as follows.
38


Definition 10 (Pseudoinverse) Let A be a linear operator. The
pseudoinverse A+ is any linear operator such that if a e Im A then AA+a = a.
In general, a pseudoinverse is not unique. Our algorithms will be
invariant to a specific choice of the pseudoinverse. If A is a symmetric oper-
ator on a finite dimensional space, A+ can be chosen to be also symmetric.
Considering the spectral decomposition of A,
A = y a-ly-cj, Ava = ava, v^va = 1, (3.1)
a
a pseudoinverse A+ can be chosen as
V = E Wit
tx^O a
For A positive semidehnite, we denote
Aa = y aavay.
a> 0
In particular, with this notation we have, A = A1!2A1!2 and A+ = A^1I2A^1I2.
We point out that Ker Aa = Ker A for any real a.
We refer to Section 2.2 and the model problem (2.1) for explanation
of some of the terminology used next.
Let Q be a domain in IE2 or IE3 decomposed into Ns non-overlapping
subdomains fii, fi2, , &ns- Let Ui be the vector of degrees of freedom for
subdomain Qj corresponding to a conforming finite element discretization of an
elliptic problem (e.g. linear elasticity, plate bending) defined on fi, such that
each subdomain is a union of some of the elements. Let iq, /\(. and /*, be the
vector of degrees of freedom, the local stiffness matrix, and the load vectors,
39


respectively, associated with the subdomain fV We will use the block notation
U\ h Ad 0 . . 0
u = u2 , / = h , K = 0 : ^ to . 0
uNs In, 0 0 .
Depending on the imposed boundary conditions and the location of the subdo-
main, the local stiffness matrix K,i is positive definite or positive semideftnite.
A subdomain without sufficient essential boundary conditions to prevent the
subdomain stiffness matrix Ad from being singular is called a floating subdo-
main. Let Zi be the matrix with linearly independent columns that generate
the kernel of Im Z% = Ker AT*. If Kt is regular, Zt is a void matrix. Denote
0 . .. 0
0 ^2 .. 0
o o ZN,
thus,
Im Z = Ker K, Ker Z = {0}.
A single mesh point x E Q has several degrees of freedom associated
with it if it lies on the interface (intersection of the boundaries) of two or
more subdomains, see Figure 3.1. Let B be a given matrix such that Bu =
0 expresses the condition that for each mesh node shared by two or more
subdomains the values of the degrees of freedom associated with that node
coincide. We denote by W the space of all vectors of degrees of freedom, and
40


Figure 3.1. Model domain: the mesh node x shared by subdo-
mains Qj and Qj and the degrees of freedom uirn and Ujn associ-
ated with the mesh node x.
41


by A the space of the vectors of values of the continuity constraint; thus,
K : W -> W, B : W -> A.
The problem to be solved is the minimization of the energy of the
system subject to intersubdomain continuity conditions
£{u) = -utKu fTu > min subject to Bu = 0, u E W. (3.2)
We assume that the global structure is not floating, that is, the solu-
tion of (3.2) is unique. From (3.2), this is equivalent to the assumption that
Ker K fi Ker B = {0}.
(3.3)
The following notation will be needed to write the FETI algorithm
concisely:
G = BZ,
F = BK+B1',
d = BK+f, (3.4)
e = ZTf,
P = / G(GtG)-1Gt.
Note that P is an £2 orthogonal projection on (Im G)^ = Ker GT.
Lemma 12 justifies the expression for P by showing that GTG is invertible.
3.1.2 Derivation of the Dual Problem
The method, as originally derived by Farhat and Roux [27], introduces
Lagrange multipliers to enforce the continuity of the solution. Solving the
42


problem (3.2) leads to the system of equations
Ku + BT A = /
(3.5)
Bu
0
A solution u of the first equation in (3.5) exists if and only if
f Bt A e Im K.
(3.6)
It must then have the form
u = K+(f BtX) + Za.
(3.7)
where a is to be determined. Substituting u from (3.7) into the second equation
of (3.5) yields
Now (3.8) multiplied by P together with (3.6) show that A satisfies the system
of equations
where we have used notation (3.4).
The first of the equations (3.9) is solved by a projected preconditioned
conjugate gradient method using an initial approximation Ao that satisfies the
second equation. The conjugate gradient algorithm requires evaluating only
the actions of PF. Most of the computational work in /' = BK+BT is concen-
trated in the action of K+. Since K+ is subdomain block diagonal, its action
can be computed in parallel and involves solving subdomain problems only.
Application of P requires solving a small coarse problem. For a scalar prob-
lem, the size of this problem is less than the number of subdomains Ns. For a
BK+(f BTX) + BZa = 0.
(3.8)
P(FA -d) = 0
GT A = e,
(3.9)
43


linear elasticity problem, the size of the problem is the number of rigid body
modes, which again does not exceed a small multiple of A'.,. A preconditioned
projected CG for the equation PFA = d using a symmetric preconditioner D
can be written as follows.
Algorithm 11 (FETI) Given an initial A0, compute the initial es-
timate
A0 = G(GTG)-1e + PA0
and the initial residual
fo = P(F\q d).
Repeat for k = 1, 2,... until convergence:
Zk 1 F^'kl
yk-1 = Pzk-1
Cfe = rl_xyk-i
Cfe
Pk = Vk-i + tPk-1 (pi = yo)
1
Uk PkPFPk
A k = Afe_i + VkPk,
i'k = rk ukPFpk
The choice of the preconditioner is discussed in Section 3.3.
Lagrange multipliers A in FETI can be seen as interface forces and
moments in the physical system. From (3.7) and the definition of F in (3.4),
the residual P(FA d) = Bu has the interpretation of jumps of the values
of degrees of freedom, between subdomains. The condition / BTA Ker K
44


means that the action of the loads and intersubdomain forces and moments
does not excite rigid body motions.
3.1.3 Saddle Point Formulation
Let us explain in more detail how enforcing the continuity using La-
grange multipliers yields the system (3.9). This approach is presented by
Farhat, Mandel and Roux [32], and it will be later used to derive a gener-
alized FETI method.
For the Lagrangian of the minimization problem (3.2),
£( n. A) = -vTKu f1 n + A7 /if/, n e W, A e A,
we solve the dual problem: find A* such that
C(A*) = maxC(A)
7 AeA v 7
max inf C(u, A).
AeA uew
By a direct computation,
(3.10)
^oo if / BTA / Ker K,
C( A)= r (3.11)
^|(iL+(/ BTA), / BTA) otherwise.
The dual problem (3.10) is thus equivalent to maximizing C(A) on the
admissible set
A = {A e A | c(A) > "x,}.
The space of admissible increments is
{Ai A2 | Ai E Al, A2 E Al} = {/i A |Brn Ker K} = Ker GT. (3.12)
At the maximum of C(A), A e A. the derivative of C, DC (A; /i), is zero in all
directions in fj, e Ker GT:
DC (A; /i) = 0 V/i G Ker GT.
45


By a straightforward computation, this becomes
(~BK+Bt\ + BK+f, n) = 0, V/i E Ker GT. (3.13)
In order to express (3.13) as a linear equation in the space Ker GT, we use the
l2 orthogonal projection defined by (3.4). Then for /i E Ker GT,
(~BK+BTX + BK+f,n) = (~BK+BTX + BK+f,Pn)
= {P(^BK+BT\ + BK+f),n)
Therefore, the dual problem (3.10) is equivalent to the linear equation in
Ker GT for the unknown /i,
H 6 Ker Gt, P(^BK+Bt(/i + A0) + BK+f) = 0, (3.14)
where A0 is an arbitrary starting feasible solution, that is, A0 E A. Denoting
A = A0 + /i, we obtain the system (3.9).
3.1.4 Algebraic Properties of the Original FETI Method
Lemma 12 (GrG) 1 exists.
Proof. Let Ga = BZa = 0. Then Za E Ker B. and Za G Ker K by the
definition of Z. It follows from (3.3) that Za = 0, hence a = 0, since Z was
assumed to be of full rank. See also [34, Theorem 5.4],
Theorem 13 The solution A of (3.9) is unique up to addition of a vector from
Ker Bt. Any solution A of (3.9) yields the same solution u of the minimization
problem (3.2), using (3.7) with a = (GTG)^1GT(d FA).
Proof. The relation between A and a follows by a direct computation. To prove
uniqueness of A, it is sufficient to show that
Ker PF n Ker G1 = Ker F n Ker G1 = Ker B1. (3.15)
46


First, BTA = 0 implies FA = 0 and GTA = 0, so
Ker Bt C Ker F fi Ker GT C Ker PF fi Ker GT.
Conversely, assume GTA = 0 and BIX = 0. Then, since G1 = ZTBT, the
first equation implies PTA Im Z = Ker F. Thus BTA e Im F. From the
dehnition of P, GTA = 0 is equivalent to PA = A. From PFA = 0, we obtain
0 = A7 PFA = A7 / 'A = (B'XfFB'X.
Since F+ is positive semidehnite, this implies that BTA e Ker F+ = Ker F.
Together, PTA e Ker F fi Im F = {0}.
Non-uniqueness of the multipliers A corresponds to redundant inter-
subdomain continuity constraints, which occur naturally at crosspoints of more
than two subdomains.
Define the space of the Lagrange multipliers as the factorspace
A = A/ Ker Br.
Since for any A e Ker BT, PA, FA, and GTA are the zero vectors, the operators
P, F, and GT induce operators on A, which will be denoted by the same
symbols. To avoid confusion, all null spaces will refer to the space A, not the
factorspace A. For example, the null space of the induced operator GT on A
will be denoted by Ker GT/ Ker BT.
It is easy to see that F is symmetric and positive semi-definite and
hence so is PFP. The next lemma shows that the operator PF restricted
to Ker GT/ Ker BT, which coincides with the restriction of PFP, is positive
definite.
47


Lemma 14 The operator PF is symmetric and positive definite on the fac-
torspace Ker GT/ Ker BT.
Proof. For u, v E Ker GT, we have from the definition of P,
(PFPu, v) = (PFu, v) = (Fa. r).
which proves that PF is symmetric positive semidefinite on Ker GT. To prove
that PF is nonsingular on Ker GT/ Ker BT, let PFu = 0 and GTu = 0. Since
F = BKBt and GT = ZTBT by definition, it follows that BTu E Ker K and
BTu Ker K, hence BTu = 0.
The original FETI method is therefore the method of preconditioned
conjugate gradients in the factorspace A for the operator equation, equivalent
to (3.9),
PFX = Pd, X E XQ + Ker GT, (3.16)
where A0 is an initial approximation to the conjugate gradient method chosen
so that GTA0 = e, and all search directions are in the space Ker GT C A.
3.2 Generalized FETI
There have been several extensions to the original method. Extension
to time-dependent problems was done in [26]. Generalization to plate-bending
problems is discussed in [51, 53] and, in detail, in this thesis.
For the original method for plate bending problems, the condition
number was observed to grow fast with the number of elements per subdo-
main [32], This is caused by the fact that plate bending is a fourth order
problem, while the FETI domain decomposition method tears the approxi-
mate solution at subdomain crosspoints, which is suitable only for second order
problems.
48


The limitation of method can be cured to extend the FETI method-
ology to obtain a non-overlapping domain decomposition method for plate
bending problems. This new method has the properties one usually looks for in
iterative substructuring methods: the condition number can be bounded inde-
pendently of the number of subdomains, and it grows only poly-logarithmically
with the number of elements per subdomain. The computational cost per it-
eration is proportional to the solution of a boundary value problem in each
subdomain, plus the solution of a sparse coarse problem with only few vari-
ables per subdomain. Such methods are commonly referred to as scalable and
quasi-optimal, though, of course, for very large number of subdomains, the
solution of the coarse problem would dominate.
The key idea of our method is to enforce the continuity of the ap-
proximate solution at the subdomain crosspoints throughout the iterations by
adding the corresponding Lagrange multipliers to the coarse problem. A simi-
lar idea was employed in the Balancing Domain Decomposition (BDD) method
for plates [42], where approximate continuity at crosspoints is enforced by
adding new basis functions to the original coarse space [46, 49] in order to keep
the energy of the approximate solution minimal with respect to displacements
that are solutions for point loads at the subdomain crosspoints (cf. Section
2.4.4). The distinguishing features of both the present method and the method
from [42] is that they are non-overlapping and work for standard finite elements
used in everyday engineering practice.
For other domain decomposition methods for the biharmonic equation
and plate bending see, for example, [12, 71]. Extensions to shells, implemen-
tation issues, and further computational can be found in [29, 31].
49


3.2.1 Duality Derivation of the Generalized FETI
In this section, we present a derivation and formulation [52] of the
generalized FETI, based on the concept of coarse optimality of a dual objective
function.
Preserving the condition GTA& = e throughout the iterations of the
Algorithm 11 can be interpreted as enforcing that every A& be optimal with
respect to all possible increments of the form Got:
GTXu = e C(Afc) > C(Afc Ga), Va, (3.17)
where C is defined by (3.11).
The key to the generalization of the FETI method to plate bending
problems is to make all Ak optimal in more directions. Let C be some given
matrix with the same number of rows as G. Each column of C will give rise to
an additional variable in the coarse problem. We shall satisfy in each iteration
the coarse optimality property
C{X) > C(X Ga C(3), Va,/3. (3.18)
with A = A*,. To satisfy this property, consider an auxiliary problem: For a
given A find a and /3 so that
C(A) max, A = A Ga (3.19)
Since we only need solutions satisfying C(A) > ^oo, we consider the maximiza-
tion problem (3.19) along with the constraint
GtA = e (3.20)
50


Introducing new Lagrange multipliers fx for (3.20), we get that a, /3,
and n solve the saddle point problem
inf sup £(A Ga Cfi, n) = sup inf £(A Ga Cfi, n)
a,f) ft, ft, a,f)
where
C(A, fi) = -G'FX + ATd + ,,r(Gr\- e)
From the optimality conditions
d£(\,n)_ <9£(A ,h)
da d i
d£(\,n)
dfi
A = A Ga C/3,
we obtain that (3.19) is equivalent to the block linear system
where
a GTF GTd
P = CTF A- CTd
n GT e
(3.21)
rfF T?ri Kjr F Lx rfF T?ri \jr r Ls GTG GT 0
T jjs s-i Ls r \jr CTFC CTG = CT 0
GTG GTC 0 0 GT
F I G c 0
I 0 0 0 G
(3.22)
The solution A of (3.19) is unique up to the addition of a vector in Ker F fi
Ker Gr. and we write it as
( GTF GTd \
A = A ^ G C 0 M+ CTF A - CTd
\ GT e )
(3.23)
51


Note that since (3.19) has a solution for any A, so does (3.21), hence
GtF
Im CTF C Im M. (3.24)
GT
Further, if A is coarse optimal, then A + 8 is also coarse optimal if and only if
GTF
G C 0 M+ cTF
GT
The generalized FETI algorithm is thus obtained from Algorithm 11 by pre-
scribing the initial approximation A0 by
S e Ker F Pi Ker GT.
(3.25)
Aq Aq
G C 0
M+
/ GTF GTd \
CTF Xq CTd
V Gt e /
and by replacing the step yu-i = Pzk-i by
Uk-i = zk-i [G, C, 0}M+
cf. Algorithm 20 in the next section.
GTF
CtF
GT
Zk-U
3.2.2 Derivation of Generalized FETI Algorithm
without Duality
The generalized FETI method can be also obtained by forcing the
iterates to satisfy also a weighted residual condition [53]. That is, we require
throughout the iterations that
CtP(FX -d) = 0, GtA = e, (3.26)
52


where C is some given matrix. Search directions that preserve (3.26) form the
space
V = {A 6 A | GTX = 0, CtPFX = 0}. (3.27)
The corresponding space of residuals is
V = {u E Im B | Gtu = 0, CTu = 0}. (3.28)
Denote the associated subspaces of the factorspace A as
V = {u + Ker BT \ u E V} = V/ Ker Br, V = {u + Ker BT \ uE V}.
Note that V and V are isomorphic: since V C Im B and Im BflKer BT = {0},
each class of V contains exactly one element of V.
We will need several properties of the spaces V and V.
Lemma 15
V = PFV'.
Proof. From the dehnition, clearly PFV C V, hence PFV C V. To show that
V C PFV, let u + Ker BT E V. Since PF is a bijection on Ker GT/ Ker BT
by Lemma 14, there is A e Ker GT / Ker BT such that PF A = u + Ker BT. It
follows that A e V.
Lemma 16 The space V is the dual of V with the duality pairing (, ).
Proof. Any A e V dehnes a linear functional A' on V by \'(v) = (X,v). Let A'
be an arbitrary linear functional on V From Lemma 15, it follows that A 'PF
is a linear functional on V, and, from Riesz representation theorem, there is a
unique A e V such that X'(PFv) = (A, v) for all v E V. From Lemma 14, PF
is a bijection on V, and it follows that the mapping between a linear functional
A' and its representation A e V is an isomorphism.
53


Lemma 17
A = V'@V~.
Proof. From Lemma 14, PF is symmetric positive definite on Ker GT/ Ker BT.
Hence, (PF)-1 is also symmetric, positive definite on Ker GT/ Ker BT, and,
using Lemma 15, it follows that
A = V' (V')-(FF)-1 = V' ((PF)V)-(^)-1 = V' V~, (3.29)
which was to be proved.
We now define a projection operator Q by
Q: A^A, Q2 = Q, Im Q = V'. Ker Q = F-, (3.30)
and compute the matrix representation of Q. This representation defines also
an operator on A, which will be denoted by the same symbol Q.
Lemma 18 The projector Q is given by the formula
Q = I
G C 0
M+
GTF
CTF
GT
where M is defined by (3.22).
Proof. Let A,A 6 A, Q(A + Ker Br) = A + Ker BT. Since A e Vr, from the
definition of P, there exists /i such that PF A = FA + G/i. From the definition
of V, GTA = 0 and CTPFX = 0, hence
+ GTGn = 0,
+ CTGn = 0,
GtF\
CTFX
GT X
54
0.
(3.31)


From the definition of Q, we have A A V, so for any liClmB, GTu = 0
and CiTu = 0 implies (A A, u) = 0. Consequently, there exist a and /3 such
that for all aelmB,
(A A, u) = (Ga,u) + (C/3,u),
which implies that
A = A + Ga + C/3 + j, j e (Im B)~ = Ker BT (3.32)
Since A e V, substituting (3.32) into (3.31) gives that a, /3, /i satisfy the linear
system
GTF(X + Ga + C/3) + G'G/i = 0
CiTF(X + Ga + C0) + CTGn = 0 (3.33)
G>(\ + Ga + C/3) = 0.
On the other hand, it is easy to see from (3.31) and the definitions of Q and
V that if a, (3, // satisfy (3.33), then A = A + Ga + C/3 £ QA + Ker If
Similarly as in the case of the original FETI, we need to find an initial
approximation Ao satisfying certain conditions, in this case (3.26).
Lemma 19 For any A0 e A, the system of equations
GtF(Xq + Ga + Cfi) + G'G/i = G'd
CtF(Xq + Ga + C/3) + CTGn = G'd (3.34)
GT(A0 + Ga + Cfi) = e
has a solution n. i. //. and
A0 = A0 + Ga + Cfi (3.35)
satisfies
CtP(FXq -d) = 0, GtA0 = e. (3.36)
55


Proof. As in the preceding proof, A0 satisfies (3.36) if and only if there is a /i
such that
GtFXq + GTGn = GTd
CtFXq + (,rG,> = ("d (3.37)
GTA0 = e
Substituting Ao from (3.35) into (3.37) yields the system (3.34), which can be
written as
d G C 0
X =
G(GTG)-le 1 o o Cl i
U J
From the factorization (3.22), it follows that Ker M = Ker X. Using symmetry
of M, we have Im M = (Ker M)- = (Ker .V) = Im A7: hence, (3.34) has a
solution.
The generalized FETI method is the conjugate gradients method for
the operator PF : V V, preconditioned by QDQT : V > V, where the
D : A > A is a given operator symmetric on V. Since QT is also a projection
and Im QT = (Ker Q)^ = V, the application of QT on V can be omitted,
and one obtains the following algorithm similar to the algorithm of the original
FETI (Algorithm 11).
Algorithm 20 (Generalized FETI) Given an initial A0, compute
the initial Ao from (3.34) and (3.35), and compute the initial residual
A) = P(FXa ^ d).
Repeat for k = 1, 2,... until convergence:
Zk-i = Drk-1
56


Vk1 QZk-1
Cfe = rl_
Cfe
Pk = Vk-1 + 7------Pfc-i (Pi = Po)
4fe-i
^ PkPFPk
Afe = Afe_i + VkPk
>'k = i'k i VkPFpk
The choice of the preconditioner is discussed in Section 3.3.
3.2.3 Method Selection for Plate Problems and
Other Generalizations
We choose the columns of matrix C, which appears in the description
of the generalized FETI method above, as vectors with a one at the position of a
Lagrange multiplier that enforces the continuity of the transversal displacement
at a crosspoint, and zeros everywhere else. By a crosspoint we understand an
interface node adjacent to at least three subdomains or to two subdomains and
the complement of fi. For a more precise formulation, see Section 4.5.1.
A similar idea can be exploited for shell problems. The continuity of
the normal displacement needs to be enforced. To avoid finding normals and
the added complexity of enforcing continuity of the normal displacement, one
may enforce continuity of all displacement degrees of freedom. This, however,
increases the size of the coarse space. Other possibilities and computational
issues are discussed in [30].
Farhat et al. in [26] studies the case of time-dependent problems
where the subdomain stiffness matrices K% are perturbed by the addition of
57


a multiple of the subdomain mass matrix, thus making the new local matrix
positive definite. Consequently, all matrices Z% are void and the natural coarse
problem is lost in time-dependent applications. The methodology developed
in [26] for reintroducing a coarsening operator in the FETI algorithm for dy-
namics problems is a special case of the present generalization where C is
taken to be the matrix G before the perturbation, that is, C = [BiZf\ where
the columns of Z% are the basis of the kernel of the local stiffness matrix of
the subdomain fV The reason why the preconditioner works for the dynamics
problems is quite different from that for the plate bending.
3.3 The Dirichlet Preconditioner
Let us decompose the space of all degrees of freedom W into the
space of internal degrees of freedom and the degrees of freedom on subdomain
interfaces,
W = W x W.
In the corresponding block notation,
B = [0 B], B : W - A,
since B has nonzero entries for interface degrees of freedom only. Also,
Z
_ 5
Z
Ker Bt = Ker BT.
and we have
G = BZ = BZ,
58


Let S be the Schur complement of K obtained by elimination of degrees of
freedom internal to all subdomains. Then
F = BK+Bt = BS Br (3.38)
and Ker S = Im Z. Evaluation of the matrix-vector product S+u reduces to
the solution of independent Neumann problems on all subdomains, cf. Section
2.4.1. Inspired by (3.38), we choose D = BSBT, giving the preconditioner
QD = QBSBt. (3.39)
This preconditioner is called the Dirichlet preconditioner, since evaluating the
matrix-vector product Sr is equivalent to solving independent Dirichlet prob-
lems on all subdomains, cf. Section 2.4.1.
3.4 Interface Formulation of Park et al.
In [59, 58], K.C. Park has developed a similar substructuring method.
We will show that this method can essentially be written as FETI with a special
choice of the interface continuity operator B.
The method augments the Lagrangian by introducing another vari-
able ug, the global vector of degrees of freedom, on the interfaces between
subdomains. This variable is redundant and is later eliminated, yielding an
analog of FETI. We will show that the method can be written as FETI and
discuss the resulting choice of B and its implications.
Let L be the subdomain assembly matrix, that is a zero-one matrix
mapping the local subdomain degrees of freedom to global ones. Using this
matrix, the continuity of the solution is enforced through the constraint
U = Lug.
59


The Lagrangian can be written as
C(u, ug, A) = -u Ku f1 u + X1 (u Lug).
Decomposing u = Za + /3 into a rigid body mode component and its comple-
ment (pure deformation component), we obtain
£(n. /3, a g. A) = 3tK/3 f'i fTZa + A T(Za + (3 Lug).
The stationarity condition leads to the system of equations
1 h^ O O 1 P f
I 0 Z -L A 0
0 ZT 0 0 a -zf
0 0 0 ug 0
Eliminating the deformation component, (3 = K+(f A), the system becomes
K+ Z A K+f
ZT 0 0 a = -Zf
-LT 0 0 ug 0
The first equation of the system is
K+A + Za Lug = K+f. (3.40)
We first eliminate the term Za by multiplying (3.40) by the orthogonal pro-
jection Pz onto (Im Z) = Ker Z1. Pz = I Z(Z' Z) lZ'. and obtain the
equation
PZK+A PzLug = Pzl\ f. (3.41)
60


equivalent to (3.40) holding for some a. Next, using the orthogonal projection
Pi onto (Im PZL)- = Ker (PZL)T,
Pi = i PzL{LTPzL)~lLTPz,
in a similar way,
P,,PZK+A = P,PzK+f.
This is the equation (15) in [58].
We may reverse the order of elimination of the terms, and we obtain
from (3.40)
PzPlK+X = PzPLK+f, (3.42)
where the orthogonal projections Pz and Pz are given by
PL = I L(LTLylLT
Pz = I^PlZ(ZtPlZ)-1ZtPl.
Since Im Pz C Im PL = Ker LT, we have v = PLv for v e Im (Pz). Thus,
I>i,l\ i' = PlK+Plv.
Now choosing
we may write
B = Pl,
(3.43)
PZ = I PlZ(ZtPlZ)-1ZtPl = I G(GtG)-1Gt = P,
where we have used the notation (3.4). This shows that the method with the
equation (3.42) can be written as FETI with the special choice of B (3.43).
61


Let us illustrate this choice. If a mesh node is on the interface of
two subdomains, then, omitting all other degrees of freedom and Lagrange
multipliers, the corresponding block of the assembly matrix L is
1
1
and the corresponding block of the projection PL is
1
2
1 -1
1 1
Similarly for a three node and four node interfaces, we get
2 -1 -1
i -1 2 -1 -
-1 -1 2
and
3 -1 -1 -1
1 3 -1 -1
1 -1 3 -1
1 -1 -1 3
respectively. The choice of B = PL has some interesting implications on the
analysis, which will be discussed in Section 4.4.
62


4. Convergence Analysis
In this chapter, we prove that for a second order problem the con-
dition number of the preconditioned FETI method is bounded independently
of the number of subdomains and poly-logarithmically in terms of subdomain
size, similarly as it is in the case for other optimal non-overlapping domain
decomposition methods [11, 23, 25, 49, 48]. We prove that the choice of the
interface continuity operator as given by Park et al. assures that the condition
number is bounded by (7(1 -l-log ^)2. Finally we also show that the generalized
FETI converges poly-logarithmically for the plate bending problem problem.
Analysis of domain decomposition methods typically demonstrates
spectral equivalence of the quadratic form that defines the problem in a varia-
tional setting and the quadratic form that defines the preconditioner, often by
way of P.L. Lions lemma [6, 23, 24, 44], Since the preconditioner in the FETI
method is quite complicated and is not defined in terms of a quadratic form,
we proceed differently. We find a bound on the norm of the product of the
system operator and the preconditioner, so as to bound the maximal eigen-
value, and a bound on the inverse, to bound the minimal eigenvalue. Related
analyses were previously done for methods without crosspoints between the
subdomains, or done formally in functional spaces, cf., for example, Glowin-
ski and Wheeler [37] and Bakhvalov and Knyazev [4], We present a complete
analysis in terms of upper and lower bound on the preconditioned operator
for decompositions with crosspoints in 2D and edges and crosspoints in 3D for
63


second order elliptic problems. We show that the condition number is bounded
by (7(1 + log ^)7, where 7 = 2, 3, h is the characteristic mesh size, and H the
diameter of the subdomains.
4.1 Preliminaries
In this section, we present some results that will be used in analysis
of FETI methods. The results are mostly concerned with estimates for finite
element functions used to discretize our model problems.
We will assume that a domain fi is divided into a set of nonoverlapping
subdomains fij, i = 1,..., Ns, Cl = Cli U... U ClNs The subdomains are assumed
to be shape regular of diameter H according to the definition below. We will
formulate all results for one of the subdomains fij and assume that constants
in the estimates do not depend on the index of the subdomain.
Definition 21 A subdomain C IRd is said to be shape-regular of diameter
0(H) if it can be generated from a reference domain (square or cube) Cl of unit
diameter by a mapping 7q such that fij = 77(f)). The mapping is assumed to
satisfy
m\ where <977 is the Jacobian of the mapping, ||.|| is the Euclidean IRd matrix
norm, and cs is a positive constant.
We will be using the Sobolev spaces IF/,.,,. For p = 2, the spaces
are Hilbert spaces and we denote them Hk = Definitions of Sobolev
spaces can be found, for example, in [1], For k ^ 1/2,1, we use || ||fejPjn and
I \k,p,n to denote the standard Sobolev norms and seminorms of functions in
WkjP. Following [22], we define the scaled Sobolev norm for a scalar function
64


u E H^Qi)
II II2 | |2 , ^ || 112
Pill,2,0, P11,2,0* + ff2 11 U\10,2,0,!
where the Sobolev seminorm is defined by
l?,2,o, = / W^u{x)fdx.
We define the scaled Sobolev norm for a scalar function u E Hl!2(dVtj)
II II2 _ | 12 . 1 || 112
11 111/2,2,ao, P11/2,2,50* + jj\\ 110,2,50,!
where the Sobolev seminorm is defined
u
2
1/2,2,50,
an, Jan
by
| u(x) u(y) \2
Ik y\\d
dxdy.
Here || || is the Euclidean norm in IRd. We note that the space Hl/2(dVtj) is
the space of traces of functions in the space Hl(Vtj).
If u = (tii, U2) is a vector function, then we define, for example,
111,2,0, lll||l,2,0, + ll2||lj2,0,-
Other norms and seminorms of vector functions are defined analogously.
We present two variants of the Poincare-Friedrichs inequality. We
prove the first for the sake of completeness and to demonstrate the technique
that is used to prove its variants for finite element functions. We note that
//,. /. = (). I_denotes the space of polynomials of degree at most k.
Lemma 22 Let Q be a continuous linear functional on Hll2(dVt,i) such that for
all u E Pq, G{u) = § implies u = 0. Then there exists a constant c independent
of H such that for all u E Hll2{dVti) u E Ker Q
u
4,2,an, < CIU,2,an,-
65


Proof. We prove the theorem for a reference domain of diameter 1 by contra-
diction. The result in the scaled norms then follows.
If the inequality is not true, there exists a sequence {n} C Ker Q
such that
ll'Un||i,2,an = 1 and \Un\l.2,dh 0-
Due to the compact imbedding of Hll2(dVt) in L2(d£i), there exists a subse-
quence of {} that converges in L2(dQ). Since |n|i 2an the subsequence
{unk} is Cauchy in Hll2(dVt). Therefore it converges in Hll2(dVt) to some
u E Hll2(dVt). The continuity of the norm and the seminorm on Hll2(dVt)
implies
\\uk2,2,dh = 1 and Mf,2,ah= - (4-2)
Therefore, by the second equation in (4.2), u = k almost everywhere, where
k is some real number. Let us assume without loss of generality that u = k
everywhere; that is u E P0. Since Q is continuous Q(unk) Q(u) = 0. This
yields, by the assumption, u = 0 which contradicts the first equation in (4.2).

The proof of this Poincare-Friedrichs inequality in scaled norms is
similar to the proof of the previous lemma.
Lemma 23 Let Q be a continuous linear functional on Hl(Vt,j) such that for
all u E Pq, Q{u) = 0 implies u = 0. Then there exists a constant c independent
of H such that for all u E if1(Qi), E Ker Q
IM|i,2A < c\u\ij2,ni-
66


4.1.1 Estimates for Pl/Ql Elements
Let 1 (Qj) be a conforming finite element space of PI or Q1 elements
[14] satisfying the usual regularity and inverse properties and possibly some
essential boundary conditions. That is, for example in IE2, we assume that
where each element K of the triangulation Th,i is a triangle or a rectangle.
Furthermore, for all K E %i
c\d(lC) < h < C2p(JC), (4.3)
where d(K) is the diameter of /C, and p(K) the diameter of the circle inscribed
in K. Then h is called the characteristic mesh size. A vertex of an element
K E %,i will be referred to as a mesh point, a nodal point, or just a vertex.
If essential boundary conditions are prescribed on E C we assume that
/i(£) > ccp(dQi), where p(-) denotes the measure. We note that functions in
V^l(Vtj) are continuous.
As in the previous chapter, the corresponding space of vectors of
degrees of freedom is denoted Wy. Let IPi : Wy V[l(Vtri) denote the linear
one-to-one transformation that maps a vector of degrees of freedom to the
corresponding finite element function. This transformation is often called finite
element interpolation.
We decompose a vector of degrees of freedom U E Wy into internal
and boundary degrees of freedom assuming the boundary degrees of freedom
are listed last. That is, in block notation, U = [UT, UT]T, where U is the vector
of boundary degrees of freedom. We denote Wy the space of boundary degrees
67


of freedom and define the matrix T] so that U = Till. Then, IpiU dehnes a
function on dift which depends on U only. By abuse of notation, we will write
IpiU = Ipi[0T, UT}T on Oil,.
We summarize some well known results and inequalities in a form
suitable for our purposes. The next lemma summarizes the fact that the H1!2
norm of a zero extension of a piece of a function can be bounded by the
norm of the function.
Lemma 24 Let E C be a vertex, edge, or face (if d = 3) of subdomain
fij. A face is understood not to contain adjacent edges, and an edge does not
contains its endpoints. For for all z E V^l(dVtft.j define w E V^l(dVtj) by
w(x) = z{x) on all nodes of triangulation x E E, w(x) = 0 on all other nodes
of <9fV Then
\w\
5,2 ,an,

I Alio,2,an,;
where
(3 = 1 if d = 2 and E is a vertex, or d = 3 and E is an edge or a vertex;
ft = 2 if d = 2 and E is an edge, or d = 3 and E is a face.
Proof. The inequality for d = 2 was proved in [49, 48]. The case when d = 3
follows from Lemmas 4.1 and 4.2 in [11] if E is an edge or a vertex, and Lemma
4.3 in [11] if E is a face. Cf. also [23].
The following lemma can be proved by using Lemmas 4.1 and 4.2
in [11] and estimates the H1!2 seminorm of a spike function.
Lemma 25 There exists a constant c independent of h and H such that for
all u E u(x) = 0 for all mesh points x E dift, x ^ xq,
cIMIo, 00,50,
68


The following lemma shows the equivalence of the discrete seminorm
defined by the Schur complement and the if1/2 seminorm. In the next section
we will apply this result to the stiffness matrix K,h which is why we the same
symbol here. This result is standard [10, 69] and it is proved here for the sake
of completeness.
Lemma 26 Let /\, be a symmetric matrix that satisfies
c\IpiU
2
1,2,0,
< {KjUjU) < C\IpiU\f 2n
VU 6 Wp
Let U = TjU and the matrix be decomposed so that
KiU
1 E-h 1 U
1 1 U
Let Si be the Schur complement of K in Kj. Then, there exist constants c0
and CQ independent of h and H such that
cQ\IpiU\
1,2,50,
< (SiU,u) < cQ\iP1u\i
2,50,
MU e Wi
Proof. From the dehnition of the Schur complement (Definition 6) and the
property (2.8), it follows that
{SjU, U) = (KU, U) (f\ lI\'r. kTU) = inf {KiU, U)
U=[UT,UT}T
Hence, from the assumption by invariance of the | 11,2,0, seminorm to adding
a constant and by the discrete extension theorem [22] for the scaled norms, we
obtain
(SiU,U) U=[UT>UT]T 2 Ll
69


where k E JR. Choosing k
yields the equivalence of the scaled
norm to the seminorm by the Poincare-Friedrichs inequality (Lemma 22). On
the other hand, by the trace theorem, invariance of the seminorms to adding
a constant, and the Poincare-Friedrichs inequality (Lemma 23), we have
where ct and cp are constants arising from the trace and Poincare-Friedrichs
inequality respectively.
The following theorem shows a variant of the Poincare-Friedrichs in-
equality. We note that the dimension of W] implicitly depends on h and for a
uniform mesh, the lemma follows from Lemma 22.
Lemma 27 Let S) be the Schur complement from Lemma 26. There exists a
constant C independent of H and h such that for all U E U). U Ker S)
Proof. If there is an essential boundary condition imposed on then
Ker Si = {0} and the statement follows from the Poincare-Friedrichs inequal-
ity Lemma 22. Otherwise, Ker Si is spanned by the vector of ones. Then the
condition U Ker S) implies XqLi C, = 0 for (.' = [Ui,..., Un]t- Following
along the lines of the proof of Lemma 22, we obtain a sequence of vectors
Un E Wi = IR Vi"i such that
< c inf
U=[UT,UT}T
IpiU\i,2,nt < (StU,U)
IpiU\\ 1,2,50, ^ 011 pi U | i '
{Ip\Un) o J~i } u in L2{dtt)
and
(4.4)
70


Therefore, u = k almost everywhere, where k E IR. Thus, (IPiUn) o Jri > k in
L,2(d£l). Then, from the quasi-uniformity of the triangulation, we obtain
/N(n)
fe
\(IP1Un) O Ti ^ k\
cJ^\\U'-k\\*>J^['£(U?-k)
'N(n)
Since Un Ker S), it holds that U? = 0. Thus
,J=i
ctf <\\{IPlUn)oTt^kt2dh ->0.
This is possible only if k = 0 which is a contradiction with (4.4).
4.1.2 Estimates for HCT Elements
Let V^CT(Vtj) be the finite element space of HCT elements satisfying
the usual regularity and inverse properties and possibly some essential bound-
ary conditions. That is, we assume that
(-h C/,,
where each element /C of the triangulation Thji is a triangle. Furthermore,
condition (4.3) is satisfied and h denotes the characteristic mesh size. We note
that functions in V^CT(Vtj) are C\ continuous and the space V^CT(Vtj) can be
written as
rHCT
'h
(^
{v E Cll(Qi) : V/C E %i, v\£ E P?,{K.j) for all subtriangles K3 of K,
dv
|aka. E Pi{ap(ij) for all sides {a^aj) of the triangle K.}
dn
On each triangle, a function v in Vf^CT(Qi) is determined by the values v(a,i)
and the values of its derivatives J^(a*) at the vertices of the triangle. The
71


internal vertex of subtriangles is not arbitrary, but it is determined to guarantee
the unisolvence of the element. This is as much detail as we will need in our
considerations. We refer to [14] for more details about the HCT element.
The corresponding space of vectors of degrees of freedom is denoted
Wi. Let I pi : Wi VlfICT(£li) denote the linear one-to-one transformation
that maps a vector of degrees of freedom to the corresponding finite element
function. As in the previous section, we decompose a vector of degrees of free-
dom U into internal and boundary degrees of freedom assuming the boundary
degrees of freedom are listed last. That is, in block notation, U = [UT, UT]T,
where U is the vector of boundary degrees of freedom. We denote Wi the space
of boundary degrees of freedom and define the matrix T] so that U = TjU.
Then, IhctU defines a function on which depends on U only. By abuse of
notation, we will write IhctU = Ihct[Ot,Ut}t on cKh.
We summarize here some well known results and inequalities in a
form suitable for our purposes. The following lemma estimates the norm of a
spike function and is proved in [42],
Lemma 28 Let x be a vertex of a subdomain Q*. For u E V^CT(Vtj) such
that u{x) = 0, deline z E \ '/J*1 (il,) by z{x) = n(.r) = 0, Vz(x) = Vn(.r). and
z(y) = 0, Vz(y) = 0 at all other nodes y of <9fThen
11^111,2,50, < C (X + 1Sx) +
The following estimate of the trace norm of the extension by zero is
proved as in [10, Lemma 3.5].
Lemma 29 Let u E VlfCT(Qi). Then there exists a constant C such that if
72


supp u fl dQi is contained in a segment a of dQi of length r, then
+ C ^1 + log ^
The following estimate is a modification of the previous lemma.
Lemma 30 Let u E VlfCT(Q,i) and supp-u PidQi be contained in a segment of
dQi and X\ and 2/1 be its endpoints. Let x2 and y2 be the nodal points next to
x\ and |/i respectively in the segment. Let (x2,y2) denote the segment between
x2 and y2. Let the length of this segment be r > h. Assume that the function
Proof. By the previous lemma, the inequality holds with (x2,y2) replaced by
the segment (xi,yi). By the definition of the if1/2 seminorm, we have
the last two terms are bounded by c2||'u||o,oo,02,?/2)- ^
We will also need a straightforward extension of the discrete Sobolev
inequality of Dryja [22] to piecewise polynomial functions of order p > 1 [57].
Lemma 31 For every p > 1 there exists a constant C = C(p) such that for
every u continuous on a C dQi of length r such that u E Pp on the side of
every element,
73


The following lemma establishes a useful inequality between discrete
and Sobolev norms.
Lemma 32 Let u E Vlf1CT(Q.i) and u(x0) = 0 for some x0 E Let U =
Ihctu and U = Till. Then, for all such U,
h\\Uf < c(l + H2)\\Vu\\l2jmi,
where c is independent of h and H.
Proof. Let E be an element edge on and x\ and x2 be its endpoints. Each
component of Vu is a polynomial of degree at most two on E. Since all norms
in a finite dimensional space are equivalent, it holds that
/i(||V-u(:ri)||2 + II2) < Ci min Ilf lino p
feP2XP2,f(xi)=Vu(xi),f(x2)=Vu(x2)
< Ci||V||gj2jS
Summing over all edges of the boundary, we obtain
where U0 denotes the displacement degrees of freedom of U. We show that
h\\UQ\\2 < c2ff211Vm||o2,5nj Since u(xo) = 0, we can write for any mesh point
x on the boundary
u{x) = / Vu(y) r(y)dy,
J(xo,x)
where (xq,x) C dQi is the part of the boundary between xq and x, and r is
the vector tangential to the boundary. Squaring, using the Cauchy inequality,
and considering that the length of (xq,x) is bounded by c2H, we obtain
u{x)2 < C2H||Vw 11^2,00,
74


The statement of the lemma then follows by summing over all mesh points on
the boundary.
The following lemma shows the equivalence of the discrete seminorm
defined by the Schur complement and the if1/2 seminorm of the gradient of
the corresponding finite element function and is an analog of Lemma 26 of the
previous section.
Lemma 33 Let /\( be a matrix that satisfies
c\IHCtU\12^ < (Kill, U) < C\IECtU\\%n. VU E Wi.
Let U = TjU and let the matrix Kt be decomposed so that
Kill
1 E-h 1 U
1 1 u
Let Si be the Schur complement of K in Ki. Then, there exist constants c0
and Cq independent of h and H such that
co\VIhctU\{ 2dQi < (SiU, U) < c0\viHCTu\i2mi VU E Wi
Proof. From the definition of the Schur complement it follows that
(SiU, U) = (KU, U) u=\uT ,iiT]T
Hence, by the discrete extension theorem [41], we obtain
(SiU, U) 4,2
U=[UT,UT)T
On the other hand, the trace theorem and invariance of the
to adding a linear function, and the Poincare-Friedrichs theorem
io0. seminorm
1 2 Li
|VIp\U\\2m_ < ctCp inf \IpiU\2,2,nr
2 r u=\UT UT]T
75



The following theorem shows another variant of the Poincare-Fried-
richs inequality. We again emphasize that the dimension of Wi implicitly de-
pends on h and is an analog of Lemma 27 of the previous section.
Lemma 34 Let Si be the Schur complement from Lemma 33. There exists a
constant C independent of H and h such that for all U E U',-. U Ker Si
Proof. Let us prove the case when Ker S is nontrivial. In the other case, the
statement follows from the standard Poincare-Friedrichs inequality. Following
along the lines of the proof of Lemma 22, we obtain a sequence of vectors
Un E Wi = IR Vi"i such that
Therefore, u = (k\, k2) E IE x IK almost everywhere. Thus, (VIhctU71) o Ti
corresponding to the partial derivative with respect to x, f/f, and the partial
derivative with respect to y, U2 Then, by the proof of Lemma 32, we obtain
VIiici C |^hiri f
(VIrctUn) o y u in L2(dtt) x L2(dtt)
and
(4.5)
(k\,k2) in L2(dVt) x L2(dVt). We decompose the vector Un into two vectors
(VIHCTUn) oT^{kuk2)
Since Ker Si, it holds that Xilo^ f L = 0 and Xil^ U2^ = 0. Thus
c{k\ + kl) < \\(VIHCTUn) o X, (h, k2)\\l2M 0.
76


This is possible only if (ki, k2) = (0, 0) which is a contradiction with (4.5).
4.2 Abstract Analysis Framework
This section develops a framework for estimating the condition num-
ber of our method. We estimate the minimum and maximum eigenvalues of
the matrix of the FETI formulation and the preconditioner.
The following simple lemma will be the basis of our estimates. It will
allow to reduce estimates of norms to estimates of boundedness and coercivity.
The proof follows a standard argument and it is presented for completeness
only.
Lemma 35 Let X be a Banach space, X' the dual of X. and A : X X' a
linear operator that is onto and satisfies the conditions
{//. Ax) = (x, Ay), Mx, y e X
ca\\x\\2x < (x,Ax) < CU||a;||!', Va; E X
(4.6)
(4.7)
with constants C \. ca > 0. Then
^ A. } 11
1
Proof. From (4.6)
sup .. ..
x£X 11 X\\ x
Ax Hx'
From (4.7) and the fact that A is onto, we obtain
concluding the proof.
77


Assume that the preconditioner of a FETI method can be written in
the form RDRT, where R is a projection, and RDRT is an isomorphism from
% onto Let % and its dual %' are subspaces of Ker GT. Furthermore,
let H = Im RT. The original FETI satisfies this with R = P and % =
Ker GT/ Ker BT. The generalized FETI satisfies the assumption with R = Q
and H = \ .
For the purpose of analysis, we equip the space % with the norm
\\v\\n = \\BTv\\s = (SBTv,BTv)1/2. (4.8)
Since BTv Ker S for v E Ti, (4.8) indeed defines a norm rather than only
seminorm. The dual space %' is equipped with the dual norm
|[A|[^' = sup . (4.9)
ven \\v\\n
Since Hi A. it immediately follows
(X,v)
A
w
sup
sup
(A, Bw)
(4.10)
vG'H 11B o||g w&Y, BwE'H \\^
The norm on % was chosen so that the preconditioner RD is trivially
coercive and bounded.
Lemma 36 For all v E %,
(v, RDv) =
Proof. Let v E % = Im RT. Since RT is a projection, we have by definition of
the preconditioner D,
(v, RDv) = (R r. BSB r) = (r. BSB r) = (B r.SB r) = \\v\\2n,
which was to be shown.
78


Coercivity and boundedness of the system operator PF on %' will be
estimated using the following lemma.
Lemma 37 For all A E ,
(A, FA)
sup
weW,w-Ker S
(A, Bw}2
IMI!
Proof. Let A 6 "W'. Then
sup
w&V
(A, Bw}2
IMI!
(A, FA) = (S+BtX, Bt\) = (S-1/2BtX, S-1/2BtX)
(S-1/2BTX,x)2
\x\
since
||5-1/2FtA||2 = sup
xew
{BTX,S-^2x)2 {BtX,S-V2x2)2
sup -----------------7-r- = sup ----------------1|--7-r----
xew. x=xt+x2 \\Xi+X2\\2 x2ew, X2-Ker S 11^21|
i£Ker S £2 Ker S
S^l!2xi = 0 and lldl2 = Ibill2 + lla^ll2. Now write any w E W as
w = wi + w2, Wi E Ker S, w2 = S 1^2x2 Ker S.
X E'H! implies that {BTX,Wi) = 0, hence (BTX,w2) = {BTX,w) = (X,Bw). It
follows that
(A, FA)
sup
W2&V, w2Ker S
(BtX,w2)2
(w2, Sw2)
sup
w&V
(A, Bw)2
il iT? 5
\m\s
which was to be proved.
It is well known [39] that after k iterations of the preconditioned
conjugate gradient method, the energy norm or the error |||e||| = (PFe, e)1!2
is reduced by a factor of at least 2((y/re 1)/(s/k+l))fe, where k is the condition
number. The condition number is given in our case by
k = k(RDPF)
Xmax(RDPF)
Amin (RDPF)
(4.11)
where Amax and Am;n are the maximum and minimum eigenvalue, respectively.
We are now ready to prove an abstract bound on k.
79


Theorem 38 Assume there exist constants (\. C2 such that
(i) for any A E %' and w E W such that Bw E H, there is w E W such that
(A, Bw) = (A, Bw), and \\w\\2s < C\\\BTBw\\2s]
(ii) for any X E H' and w E W, w Ker S, there is w E W such that Bw E %,
(A, Bw) = (A, Bw), and \\BTBw\\2s < C2\\w\\2s.
Then k(RDPF) < C\C2.
Proof. Lemma 35 applied to the operator RD together with Lemma 36 give
WQ-^Wn^n1 II (-^0 1|l n'-tu 1- (4-12)
From assumption (i), we have
(A, Bw}2
1
> sup
sup 2 _
wew INIs Li wew,Bwen
while assumption (ii) gives the converse bound
(A, Bw)2
\\BTBw\\2
(4.13)
(X,Bw)2 (X,Bw)2
sup rjp < C2 sup ,.=Tp ||2
weW.weKerS INiS w&Â¥,BweU \\B Bw\\s
(4.14)
Using (4.10) and Lemma 37, we see that the inequalities (4.13) and (4.14)
imply the inequality
i|A||^ < (A, PFX) < C2||A|&,, VA 6 W.
Ui
Applying Lemma 35 to the operator PF, we obtain
I^IIW-,H < C'2
|(PF)-1||^W (4.15)
From (4.12) and (4.15), we have
i.)I)PF\}(' < \\QD\\n^w\\PF\\w^n < C2
80


and
(QDPF) 1\\w^w < \\{QD) 1\\u'^h\\(PF) 1\\'h^w The result follows.
4.3 Analysis of the Original FETI
In this section, we use the abstract analysis framework from the pre-
vious section to prove a bound on the condition number of FETI for a second
elliptic problem, utilizing results from Section 4.1.
4.3.1 Assumptions
We need more specific assumptions in order to be able to prove a
bound on the condition number k. Let us recall the model problem (2.1).
Assume we are solving the boundary value problem
with a(x) a measurable function such that 0 < aQ < a(x) < cti a.e. in Q.
The domain Q is assumed to be divided into non-overlapping subdo-
mains fij, i = 1, ...,NS, that are shape regular of diameter H (Definition 21).
Assume that Vh(fi) is a conforming PI or Q1 finite element space on a triangu-
lation of fi, which satisfies the standard regularity and inverse assumptions (cf.
Section 4.1). In particular, we recall that the degrees of freedom are values at
nodes of the triangulation h denotes the characteristic element size. Each sub-
domain fij is assumed to be a union of some of the elements, and all functions
in 14(fi) are zero on <9fi.
Au = g in fi, u = 0 on dQ
where A is a second order elliptic operator
81


Let us define the restriction operator Rt : W Wq by the equation
RjW = iVi,
where tv = [tvj'tvl' ... Wi E U). i = 1,..., Ns.
We will assume that the interface continuity operator B is defined as
follows. Let wr(x) and ws{x) denote the pair of degrees of freedom correspond-
ing to the mesh node x E dVtr fi dVts and let (Bw)rs(x) be the entry of vector
Bw that corresponds to the mesh point x, and subdomains Qr and Qs. We
require each such (Bw)rs(x) to have form
(Bw)rs(x) = ars( wr(x) ws(x)), (4.16)
where ars is either 1 or -1. Note that ars does not depend on x, that is,
coefficients ars are uniform along edges (and faces, in 3D) between two sub-
domains. For each node x that belongs to the interface of two and more
subdomains, x E i = Si, s2, , sn, vector Bw contains n 1 entries,
{Bw)Sksk+1{x), k = 1,..., s-i. For an example of the definition of the values
of Bw for (si, 52,53) = (1,3,2) in 2D around a crosspoint, see Fig. 4.1. We
point out that B chosen in this way has full row rank, that is, this definition
implies that there are no redundant constraints in enforcing the continuity of
the solution at the nodes where more than two subdomains meet.
Only the improved estimate in statement 3 of Lemma 39 will require
the specific definition of this section. If redundant constraints are allowed, the
estimates, with the exceptions of the improved estimates, still hold which can
be shown by minor modifications of the proofs. See Section 4.4.
An additional connectivity assumption on the decomposition is need-
ed. Let Nrs, r,s = 1 ,...,NS, r ^ s be the number of interface conditions
82


between subdomains fir and Qs, Nrr the number of degrees of freedom of
subdomain ilr and N = max(v N^. We will assume that there exists constants
c and no independent of h and H such that for all r, s = 1,..., Ns,r ^ s, for
which Nrs > 0, there exists a sequence of indices {rj}, rQ = r, Tk = s, k < n0
such that
Nri_iri > cN \fi = 1,..., k.
Throughout the next section, c, C, C\, c2, c3, C4, c5 and c6 denote posi-
tive constants independent of H and h.
4.3.2 Discrete Norm Bounds
The stiffness matrices /\(. i = 1,...,NS satisfy the assumption of
Lemma 26. Therefore, using Lemma 26 for each subdomain and summing
over all subdomains, we obtain
Ns Ns
ci ^2 \IpiR,
i=1
W\\,2,dSli ^ llWlls
< C2 \IpiRiw\l 2 an- e 11
(4.17)
i=1
where the positive constants ci,c2 are independent of the characteristic mesh
size h and the subdomain diameter H. The constants ci,c2 may depend on
the regularity of the shape of subdomains.
The following lemmas verify the assumptions of Theorem 38.
Lemma 39 For all A e Ker GT and all w E W such that BTBw Ker S,
there exists w E W such that
(A, Bw) = (A, Bw) and \\w\\2s < C(1 + logH/h)a\\BTBw\\2s.
where a = 1, and a = 0 in the following special cases:
(1) |BBt = I, which means that there are no nodes shared by more than
two subdomains.
83


Figure 4.1. Definition of B
do,
dO-2
dili
fi3
84


(2) d = 2, and the matrix B has the following property: If w E Im B1', x is
a crosspoint (node shared by more than two subdomains), IpiRiw(x) =
IpiRiw(y) for all i such that x E and all nodes y that are adjacent
to x on then IPiRjw(x) = 0 for all i such that x E <9f
(3) d = 2, B is dehned as in Section 4.3.1 and all nodes in the triangulation
belong to either one, two, or an odd number of subdomains.
Proof. Let us first prove that, in the general case, we obtain a < 1. Let tv E W
and define w = BT(BBT)^1Bw. The triangle inequality shows that
\w\\s <
-BtBw +
2 s
\BT V
-BBt
2
Bw
(4.18)
Denote 2: = |BT(I (^BBT)^r)Bw. From the definition of B in
(4.16), 2 is zero at all nodes that belong to at most two subdomains. The
remaining nodes lie on crosspoints or edges (in the 3D case) of subdomains.
From the definition of B, at every such node x, z is a linear combination of the
entries of BTBw that correspond to the same node x and the coefficients of the
linear combinations are bounded only in terms of the number of subdomains
the node belongs to. Using Lemma 24 for the crosspoints of subdomains and
the equivalence (4.17), we obtain for the 2D case that
Ns
Ns
Mis < CZ - CX1 +!g(^A))E \\LPiRiBTBw\\%2,d^-
(4.19)
i=1
i=1
The Poincare inequality, Lemma 27, yields
Ne
Y. llpiR.BBw |i < \\BTBwfs
i=1 2
In the 3D case, the argument for subdomain crosspoints is the same. In ad-
dition, we note that the coefficients of the linear combination do not change
along a subdomain edge, so it remains to apply Lemma 24 on every edge.
85


Let us now turn to the special cases that give a = 0. If \BB7 = I,
we choose w = \BtBw = w, which proves the special case 1.
Now we prove special case 2. From the definition of the H1!2 norm,
the equivalence of norms (4.17), and the fact that IpiRiB7Bw is a piecewise
linear function, it follows that
Ns
\\BTBwfs > cYJ\IpiRiBTBw\2i!2mi (4.20)
i= 1
> c ^2 (lp\RiBTBw(x) Ip\RiBTBw(y)j
x crosspoint, x£dQ^
y adjacent to x,y£dQ^
For any crosspoint x, it follows from the assumption that for every tw elm BT,
T; (IpiRiw(x) IpiRiw(y))2 = 0 =>- ^2 (Ip\Riw{x))2 = 0.
i.d nt3x
y adjacent to x, yEdQ^
Consequently, by compactness, and since there are only finitely many different
numbers of subdomains sharing a crosspoint, for all w E Im BT
*52 (IpiRiw(x))2 < c (IpiRjw(x) IpiRjw(y))2 .
i,ant3x
y adjacent to x, yEdQ^
By summation over all crosspoints x, using Lemma 25 and (4.20), we get
MU < C\\BTu\\2s,
which concludes the proof of this case.
In order to prove case 3, we verify the assumptions of case 2. We for-
mulate the proof only for a crosspoint shared by three subdomains (Fig. 4.1).
The proof is similar for a different odd number of subdomains. Let w Im B7'.
Since IpiRiw(xy) IpiR\w{xa) = 0, and IpiRiw(xy) IpiRiw(x$) = 0,
we have IpiRiw(xa) = IpiRiw(xs). Similarly, we obtain IpiR2w(xa) =
86


IpiR2w(x7), and IpiR3w(xs) = IpiRzwix^). Moreover w E Im BT implies
IpiR\w{xa) = IpiR2w(xa), IP1R2w(xy) = //>i R:iir(.r. ). and IpiR$w{xs) =
IpiRiw{xs), which can be satisfied only if IpiRiw(xa) = IpiRiw(x$) = ... =
0.
Remark 40 In general, the exponent a = 1 in Lemma 39 cannot be improved.
To see that, let us consider the configuration with values of u and Bu in the
neighborhood of a crosspoint as in Fig. 4.2, which violate the assumptions of
special case 2. Extending the values of u in Fig. 4.2 to decay as log7(t/if),
7 < 1/2, where t is the distance from the crosspoint, we obtain a vector u E A
such that
\\BTu\\s ~ C, \\Ipiu\\ij2janinan2 ~ I logh/if|7.
If u = Bir. then on dQ\ fi dQ2, u = w2 wi, which gives
11,2,50x0502 < I^Plwl|i,2,50xn502 + 1^1^211,2,50x0502
< |^P1Wi|i,2,50i + |^PlW2|i)2,502
< ||w||s,
and therefore \\w\\s > C(y)| logh/H\J for all 7 < 1/2.
Lemma 41 Let A e Ker GT. Then for all w E W, w Ker S, there is a
w E W such that B1 Bir Ker S,
(A, Bw) = (A, Bw), and \\BTBw\\2s < C(1 + logH/h)2\\w\\2s
Proof. Let tv E W be arbitrary, and put Bw = PBw. Then
(A, Bw) = (A, Bw). (4-21)
Denote Bz = Bw Bir. z E Ker S. Then, since P is an orthogonal projection,
\\Bz\\ < \\Bw\\.
87


Figure 4.2. Counter-example
fil 1/2' 1/2 -1/2. -1 1/2 1/2 1/2
f o, > s
1 1 1 1
1/2 1/2 ^4 l/2i 1 -1/2 1/2 > 1/2 ^3

88


Now, from the definition of B and the Poincare inequality, Lemma 27,
we obtain
h\\Bw\\2 < 2/i||tc||2 < C!H\\w\\2s.
Also, since z E Ker S, it is constant on each using the connectivity as-
sumption, we have the following by Lemma 24:
\\BTBzfs < C!^-\\Bz\\2(1+logH/h)2.
H
Together this yields
\\BTBz\\2s < C( 1 + logH/h)2\\w\\2s.
By the definition of B, RiBTBw on U dVtj is a linear combination (with
bounded coefficients) of (a bounded number of) Wk from all dVtk adjacent to
dVtjUdVtj. From Lemma 24 and the Poincare-Friedrichs inequality Lemma 27,
we obtain
\\BTBw\\s < C'( 1 + log(H/hj)\\w\\s, Vie E W.
Finally, summarizing,
\\BtBw\\s < \\BTBw\\s + \\BTBz\\s < C(1 + logH/h)\\w\\s-
From this and (4.21), the statement of the lemma follows.
4.3.3 Condition Number Estimate
We are now ready to prove the final result. It follows from the ab-
stract estimate in Lemma 35 with the assumptions verified by Lemma 39 and
Lemma 41.
89


Theorem 42 Under the assumptions of section 4.3.1, the condition number
of the FETI method with the Dirichlet preconditioner satisfies
^max(PDPF)
A min(PDPF)
with 7 = 3, and 7 = 2 in the special cases listed in Lemma 39.
4.4 Analysis of the Method by Park
We consider the same assumptions as for the original FETI for second
order problems except for the choice of B. We consider B given by projection
matrix PL (3.43). The essential property of this matrix B used here is that
B = BtB = Bt, since B is an orthogonal projection.
We present two lemmas that verify assumptions of Theorem 38.
Lemma 43 For all A E Ker GT/ Ker W and all w E W such that Bw E
Ker GT/ Ker B1. there exists w E W such that
Proof. Since BT = B and IP = B. we may choose w = w.
Lemma 44 Let A e Ker GT/ Ker B1. Then for all w E W, w Ker S, there
is a w E W such that BTBw Ker S and
(A,Bw) = (A,Bw), \\BtBw\\2s < C( 1 + logH/h)2\\w\\2s
Proof. Let tv E W and Bw = PBw. Then
Let z E Ker S, Bz = Bw Bw. Then, since P is an orthogonal projection,
11 BUI < \\Bw\\.
(A, Bw) = (A, Bw), \\w\\s = \\BT Bw\\s
(A, Bw) = (A, Bw)
(4.22)
90


Now, from the definition of B and the Poincare inequality, we obtain
h\\Bw\\2 < 2/i||tc||2 < C!H\\w\\2s.
Also, since 2: e Ker S, it is constant on each dVt,h and using the connectivity
assumption, we have the following by Lemma 24
\\Bzfs < cb\Bzf{l + log if/h)2.
ri
Together this yields
Bz\|| < (7(1 + logH/h)2\\w\\2s
By the definition of B, (Bw)i on dQ,i U dQj is a linear combination (with
bounded coefficients) of (a bounded number of) wk from all dVtk adjacent to
dVti U dVtj. From Lemma 24 and the Poincare-Friedrichs inequality (Lemma
27),
\\Bw\\s < (7(1 + log(if/h))11w115, \/w e W.
Finally, summarizing,
1111^ < ll-Bm A + \\BtBz\\s < (7(1 + logif/h)||ic115.
From this, (4.22), and BTB = B, the result follows.
Theorem 45 Under the assumption of Section 4.3.1, the condition number
of the Parks variant of FETI method with the Dirichlet preconditioner (3.39)
satisfies
KtaxjQDPF)
A min(QDPF)
2
Proof. Lemmas 43 and 44 verify the assumptions (i) and (ii) of Theorem 38
with C\ = 1 and (72 = (7(1 + logH/h)2 respectively.
91


Full Text

PAGE 1

ANAL YSIS OF LA GRANGE MUL TIPLIER BASED DOMAIN DECOMPOSITION b y Radek T ezaur M.S., Univ erzita Karlo v a, Praha, Czec h Republic, 1993 A thesis submitted to the Univ ersit y of Colorado at Den v er in partial fulllmen t of the requiremen ts for the degree of Doctor of Philosoph y Applied Mathematics 1998

PAGE 2

This thesis for the Doctor of Philosoph y degree b y Radek T ezaur has been appro v ed b y Jan Mandel Charbel F arhat Leo F ranca Andrew Kn y azev T om Russell Date

PAGE 3

T ezaur, Radek (Ph. D., Applied Mathematics) Analysis of Lagrange m ultiplier based domain decomposition Thesis directed b y Professor Jan Mandel ABSTRA CT The con v ergence of a substructuring iterativ e method with Lagrange m ultipliers kno wn as Finite Elemen t T earing and In terconnecting (FETI) method is analyzed in this thesis. This method, originally proposed b y F arhat and Roux, decomposes nite elemen t discretization of an elliptic boundary v alue problem in to Neumann problems on the subdomains, plus a coarse problem for the subdomain n ull space componen ts. F or linear conforming elemen ts and preconditioning b y Diric hlet problems on the subdomains, the asymptotic bound on the condition n um ber C (1+log( H=h )) r where r = 2 or 3, is pro v ed for a second order problem, h denoting the c haracteristic elemen t size and H the size of subdomains. A similar method proposed b y P ark is sho wn to be equiv alen t to FETI with a special c hoice of some componen ts and the bound C (1 + log( H=h )) 2 on the condition n um ber is established. Next, the original FETI method is generalized to fourth order plate bending problems. The main idea there is to enforce con tin uit y of the transv ersal displacemen t eld at the subdomain crosspoin ts throughout the preconditioned conjugate gradien t iterations. The resulting method is sho wn to ha v e a condition n um ber that does not increase with the n um ber of subdomains, and again gro ws at most iii

PAGE 4

poly-logarithmicall y with the n um ber of elemen ts per subdomain; the condition n um ber is bounded b y C (1 + log( H=h )) 3 These optimal properties hold for n umerous plate bending elemen ts that are used in practice including the HCT, DKT, and a class of non-loc king elemen ts for the Reissner-Mindlin plate models. The theoretical results are conrmed b y n umerical experimen ts. This abstract accurately represen ts the con ten t of the candidate's thesis. I recommend its publication. Signed Jan Mandel iv

PAGE 5

DEDICA TION Dedicated to m y paren ts.

PAGE 6

A CKNO WLEDGEMENT F oremost I w ould lik e to thank m y advisor, Jan Mandel, for his encouragemen t in preparation of this thesis, his support, guidance, and time spen t in discussions throughout m y studies at the Univ ersit y of Colorado at Den v er. I w ould also lik e to extend m y thanks to P atric k Le T allec for useful suggestions in early stages of this w ork and Charbel F arhat for computational results.

PAGE 7

CONTENTS Chapter 1. In troduction . . . . . . . . . . . . . . 1 2. Ov erview of Related Domain Decomposition Methods . . . . . . . . . . . 7 2.1 Finite Elemen t Appro ximation, Model Problem and PCG 7 2.2 Abstract Sc h w arz Methods . . . . . . . . 12 2.3 Ov erlapping Domain Decomposition . . . . . . 15 2.4 Substructuring Methods . . . . . . . . . 18 2.4.1 P oincar e-Steklo v Operators and Sc h ur Complemen t . 18 2.4.2 V arious Preconditioners for the Reduced Problem . 25 2.4.3 The Neumann-Neumann and Balancing Domain Decomposition Methods . . . . . . . . . 27 2.4.4 The Neumann-Neumann Domain Decomposition Algorithm for Plates and Shells . . . . . . 29 2.4.5 Lagrange Multipliers and P oincar e-Steklo v Operators . . . . . . . . . . . . 31 2.5 A Tw o Lev el Method Based on Smoothed Prolongation . . . . . . . . . . . . 33 3. Finite Elemen t T earing and In terconnecting (Deriv ation and Extensions) . . . . . 37 3.1 Original FETI . . . . . . . . . . . . 37 3.1.1 Problem Setting and Assumptions . . . . . . 38 3.1.2 Deriv ation of the Dual Problem . . . . . . 42 vii

PAGE 8

3.1.3 Saddle P oin t F orm ulation . . . . . . . . 45 3.1.4 Algebraic Properties of the Original FETI Method . 46 3.2 Generalized FETI . . . . . . . . . . . 48 3.2.1 Dualit y Deriv ation of the Generalized FETI . . . 50 3.2.2 Deriv ation of Generalized FETI Algorithm without Dualit y . . . . . . . . . . . 52 3.2.3 Method Selection for Plate Problems and Other Generalizations . . . . . . . . . 57 3.3 The Diric hlet Preconditioner . . . . . . . . 58 3.4 In terface F orm ulation of P ark et al. . . . . . . 59 4. Con v ergence Analysis . . . . . . . . . . . 63 4.1 Preliminaries . . . . . . . . . . . . 64 4.1.1 Estimates for P1/Q1 Elemen ts . . . . . . . 67 4.1.2 Estimates for HCT Elemen ts . . . . . . . 71 4.2 Abstract Analysis F ramew ork . . . . . . . . 77 4.3 Analysis of the Original FETI . . . . . . . 81 4.3.1 Assumptions . . . . . . . . . . . . 81 4.3.2 Discrete Norm Bounds . . . . . . . . . 83 4.3.3 Condition Num ber Estimate . . . . . . . 89 4.4 Analysis of the Method b y P ark . . . . . . . 90 4.5 Con v ergence Estimates for Plate Bending . . . . 92 4.5.1 Assumptions . . . . . . . . . . . . 93 4.5.2 Discrete Norm Bounds . . . . . . . . . 95 4.5.3 Condition Num ber Estimate . . . . . . . 99 5. Computational Results . . . . . . . . . . . 101 viii

PAGE 9

References . . . . . . . . . . . . . . . . . 106 ix

PAGE 10

1. In troduction In this thesis, the Finite Elemen t T earing and In terconnecting (FETI) method is analyzed. It is one of domain decomposition methods for solving large systems of linear equations arising from nite elemen t discretizations of elliptic dieren tial equations. Elliptic equations often arise in modeling ph ysical phenomena. The Laplace equation, = div r = f in n ; (1.1) models electrostatic in teractions and man y other poten tial problems. Ph ysical quan tities are often go v erned b y systems of equations, as in the case of linear elasticit y whic h is modeled b y a system of equations for the unkno wn displacemen t v ector u div ( u ) = f : In the isotropic case, the stress tensor can be written as ( u ) = 2 ( u ) + tr ( ( u )) where and are called Lam e constan ts, is the Kronec k er tensor, and ( u ) is the strain, ( u ) = 1 = 2( r u + ( r u ) T ). An equation describing a plate bending problem can be obtained from the linear elasticit y equation b y a limit process that considers thic kness of the plate to be innitesimally small. In Kirc hho-Lo v e model, this yields a fourth order elliptic equation. In order to set up a w ell-posed problem, the equations need to be complemen ted b y boundary conditions. A boundary v alue problem ma y be 1

PAGE 11

reform ulated in w eak (v ariational) form [56]. Tw o basic kinds of boundary conditions are recognized. Essen tial boundary conditions, suc h as a Diric hlet boundary condition = 0 on @ n for (1.1), need to be explicitly imposed in the w eak form of the problem. Natural boundary conditions, suc h as a Neumann boundary condition @=@ n = 0 @ n for (1.1), are \naturally" incorporated in to the w eak form. The process of discretization ma y be based on a Galerkin appro ximation. In that case, an appro ximate solution is sough t in a nite-dimensional subspace of the space in whic h the w eak form is posed. Finding the Galerkin appro ximation then requires solving a system of linear equations for the coecien ts of the solution relativ e to a basis of the subspace. Finite elemen t spaces are widely used in this con text. They are generated b y basis functions that are usually polynomial on eac h elemen t of a triangulation of the domain and the supports of basis functions ha v e only small o v erlaps. The matrix of the discrete system has sev eral special properties. It is symmetric and positiv e denite. Also, due to limited o v erlaps of nite elemen t basis functions, it is sparse; the n um ber of nonzero en tries in a ro w of the matrix is a lot smaller than total n um ber of en tries in the ro w. The size and the condition n um ber of this matrix are aected b y neness of the triangulation. F or quasi-uniform meshes for whic h the nite elemen t discretization is c haracterized b y the mesh size h of the triangulation, the condition n um ber of the matrix is of the order 1 h 2 for second order elliptic problems and 1 h 4 for fourth order elliptic problems. Tw o classes of methods for solving the linear system exist, direct and iterativ e. Direct methods, usually based on a v arian t of Gaussian elimination, 2

PAGE 12

are common in ev eryda y engineering practice. They are fairly robust; their disadv an tage lies in their memory requiremen ts and their speed. Gaussian elimination leads to ll-in and the destruction of the sparse structure of the matrix. The ll-in can be reduced b y re-n um bering the v ariables, but for problems in 3D this helps only to a certain exten t. In general, the n um ber of operations required to solv e the system arising from the nite elemen t discretization is proportional to the square of the n um ber of unkno wns or w orse. Since systems arising from discretizations of problems in engineering practice constan tly push a v ailable computer resources to the limits, iterativ e methods are often preferable for large-scale problems. Their main disadv antage is probably their lac k of robustness when too wide a class of problems is considered. On the other hand, the memory requiremen ts of iterativ e methods are t ypically m uc h smaller than those of the direct methods. Usually only a small m ultiple of the n um ber of nonzero en tries of the original matrix needs to be stored. The n um ber of operations can be as lo w as a m ultiple of the n um ber of unkno wns, but often the computational cost is not kno wn in adv ance. The speed of con v ergence depends on the condition n um ber of the matrix of the system. Since the condition n um ber deteriorates with decreasing mesh size, so does the speed of con v ergence of simple iterativ e methods. It is therefore desirable to construct methods that o v ercome bad conditioning of the matrix and that are, if possible, independen t of other singular perturbations suc h as inhomogeneities or bad P oisson ratios in the linear elasticit y problem or, in the case of plate bending model, the problems arising from the thic kness of a plate approac hing zero and the model becoming a fourth order problem. Due to the elliptic nature of the underlying problem, this sort of stabilit y is often possible 3

PAGE 13

to ac hiev e and v arious algorithms ha v e been proposed. In this thesis, w e will study one of domain decomposition methods. This class of methods has gained enormous popularit y in the last decade. F ollo wing a divide-and-conquer idea, domain decomposition methods divide the original problem in to a n um ber of smaller subproblems. Sometimes suc h a division arises from breaking up a complicated geometry In man y cases, though, it is en tirely articial. The subproblems are easier to solv e because of their smaller size and often parallelism can be exploited. This is especially importan t with the onset of parallel computing. Domain decomposition methods can be seen from t w o dieren t poin ts of view. They ma y arise from separation of a ph ysical domain in to regions, where a problem is modeled b y separate partial dieren tial equations, with the in terfaces bet w een the subdomains being handled b y v arious conditions, suc h as con tin uit y The opposite approac h is to see domain decomposition methods purely as methods for solving large algebraic linear systems arising from the discretization of PDE's. In this con text, the large system is subdivided in to smaller problems, whose solutions can be used to produce a preconditioner for the large system. Ev ery domain decomposition algorithm in v olv es t w o principal issues. It breaks up the original problem in to subproblems, that are solv ed b y some kno wn method, and it resolv es in teractions of local solutions, usually b y the means of an iterativ e method. Man y domain decomposition methods ha v e been dev eloped, originally for the case of t w o subdomains. As the complexit y of problems of in terest demands splitting the problem in to small enough subproblems, a m ulti-domain case pro v es to be of more importance. Early w orks 4

PAGE 14

demonstrate independence of proposed algorithms on the c haracteristic mesh size. Ho w ev er, it has been sho wn that to ac hiev e independence of the n um ber of subdomains, a coarse space needs to be in troduced. Solving a small coarse space problem distributes the information about the solution globally th us resolving global c haracteristics of the solution not visible b y the subdomains. In the next c hapter, w e giv e an o v erview of sev eral w ell-kno w domain decomposition algorithms. In particular, w e will describe abstract Sc h w arz methods, o v erlapping and substructuring methods and their relationship to the nite elemen t tearing and in terconnecting method (FETI) [27, 32 ], the main subject of this thesis. W e also briery discuss a domain decomposition method based on smoothed prolongation [68, 64]. The remaining c hapters, except for Sections 3.4 and 4.4, are based on papers [50 51, 53] b y Mandel and T ezaur, and b y Mandel, T ezaur and F arhat. In Chapter 3, the deriv ation of the original FETI method is sho wn. FETI tears the computational domain in to non-o v erlapping subdomains and enforces in tersubdomain con tin uit y via Lagrange m ultipliers applied at the subdomain in terfaces. The Lagrange m ultipliers are used as the unkno wns and FETI form ulation is obtained b y solving the saddle poin t problem for the Lagrangian. Consisten t treatmen t of subdomain singularities leads to a small coarse problem whic h is solv ed in eac h iteration of the preconditioned conjugate gradien t method iteration. In Chapter 3, w e also presen t a generalization of the original FETI algorithm. W e describe its particular application to a plate bending problem. Chapter 4 is concerned with analyzing the original and generalized FETI methods. W e sho w that the condition n um ber of the preconditioned 5

PAGE 15

FETI method is bounded independen tly of the n um ber of subdomains and poly-logarithmicall y in terms of subdomain size. That is, the condition n umber is bounded b y C (1+log H h ) 3 where H is the c haracteristic subdomain size and h the mesh size. W e presen t a complete analysis for decompositions with crosspoin ts in 2D and edges and crosspoin ts in 3D for second order elliptic problems. F urthermore, w e demonstrate that the c hoice of the Lagrange m ultipliers giv en b y P ark [59, 58] assures that the condition n um ber estimate can be impro v ed to C (1+log H h ) 2 Finally w e also sho w that the generalized FETI con v erges poly-logarithmicall y for a biharmonic problem. In order to illustrate the poten tial of the generalized FETI method, Chapter 5 summarizes some computational results. W e consider the plate bending problem on a unit square and demonstrate that the condition n um ber is almost independen t of the n um ber of subdomains and the size of the problem, as predicted b y the theoretical analysis. 6

PAGE 16

2. Ov erview of Related Domain Decomposition Methods 2.1 Finite Elemen t Appro ximation, Model Problem and PCG In this section, w e are going to summarize some concepts and algorithms used throughout this c hapter and the rest of the thesis. W e will consider the follo wing model problem. Let n be a bounded Lipsc hitz domain in I R 2 or I R 3 W e will study the second-order elliptic equation with Diric hlet boundary condition A u = f in n u = 0 on @ n ; (2.1) where A u = d X i;j =1 @ @x i a ij @u @x j : with the matrix [ a ij ] symmetric, uniformly positiv e denite and bounded on n. The corresponding bilinear form is then giv en b y a ( v; w ) = d X i;j =1 Z n a ij @v @x i @w @x j and is dened for all v and w in the Sobolev space H 1 (n) (the space of generalized functions with square in tegrable rst deriv ativ es). Let us denote the L 2 (n) inner product ( v; w ) = Z n vw: 7

PAGE 17

By the div ergence theorem (in tegration b y parts), the problem (2.1) can be written in w eak (v ariational) form: Find u 2 H 1 0 (n) suc h that a ( u; v ) = ( f; v ) (2.2) for all v in V = H 1 0 (n), where H 1 0 is the completion of smooth functions with support in n with respect to the norm in H 1 (n). Since the bilinear form a ( u; v ) is symmetric, con tin uous and elliptic (coerciv e) [56], the problem has a unique solution b y the Lax-Milgram theorem. Let us no w use the standard Galerkin appro ximation. W e look for an appro ximate solution of (2.2) in a nite-dimensional subspace V h (n) of the space V The Galerkin appro ximation is the solution of the follo wing problem: Find u h 2 V h (n) suc h that a ( u h ; v h ) = ( f; v h ) ; (2.3) for all v h 2 V h (n). This problem leads to a system of linear equation when a basis for V h is c hosen. Let f g N i =1 be a basis for V h (n). Assuming u h = N X i =1 u i i ; (2.3) becomes N X i =1 a ( i ; j ) u i = ( f; j ) ; j = 1 ; : : :; N Matrix K = [ a ( i ; j )] is called the stiness or Gramm matrix and the righ t hand side v ector [( f; j )] is called the load v ector. The v ector of unkno wns, [ u i ], is also called the v ector of degrees of freedom. F or a wide class of appro ximation spaces V h u h is a good appro ximation of u In the case of nite elemen t methods, the construction of suitable 8

PAGE 18

spaces V h relies on triangulation, whic h splits the domain n in to small disjoin t regions of simple geometric shape, suc h as triangles or quadrangles in I R 2 or tetrahedrons or hexahedrons in I R 3 ; h refers to the c haracteristic size of these regions. Under certain assumptions whic h prev en t the triangulation from degenerating (e.g, the size of the angles in triangles is bounded from belo w so as they do not become too sharp), the ner the triangulation, the closer a nite elemen t Galerkin solution to the exact solution, [14]. F unctions in a nite elemen t space V h usually arise from a polynomial in terpolation on the elemen ts of the triangulation. Ev ery polynomial dened on a giv en region is uniquely determined b y its v alues and perhaps also the v alues of its deriv ativ es at some nodal poin ts, usually the v ertices of the region. Eac h function in V h is th us determined b y a set of v alues at nodal poin ts, so called degrees of freedom. Simple examples of nite elemen t spaces include spaces formed b y con tin uous functions linear on triangle regions in I R 2 or tetrahedrons in I R 3 or bilinear functions on quadrangles. These spaces are being referred to as P 1 and Q 1 respectiv ely and their en tries are uniquely determined b y their function v alues at the v ertices of the triangulation. Another t ype of nite elemen ts w e will encoun ter later in this thesis is HCT elemen ts, whic h are C 1 con tin uous and determined b y the function v alues and the v alues of the rst deriv ativ es at the v ertices of triangular regions in I R 2 The standard basis for a nite elemen t space is the one in whic h eac h of the basis functions has exactly one degree of freedom equal to 1 and the rest 0. The unkno wns in the linear system arising from the discretization are directly the degrees of freedom of the Galerkin appro ximation (hence the name). The o v erlaps of the supports of the basis functions are small, whic h causes 9

PAGE 19

the stiness matrix to be sparse. Due to the properties of the bilinear form, the stiness matrix is also symmetric and positiv e denite. Preconditioned conjugate gradien t (PCG) method is therefore often considered for solving the system. A detailed description and analysis of PCG (and other methods) can be found, e.g., in Golub and V an Loan, [38 ]. Here w e giv e only a brief summary The conjugate gradien t method is based on the observ ation that the solution of Ax = b is the only minim um of Q ( x ) = 1 2 x T Ax x T b: A sequence of v ectors whic h con v erges to the solution of the system can be obtained b y setting x k +1 = x k + k p k ; where p k is a suitably c hosen direction and v ector x k + k p k minimizes Q on the line x k + p k In PCG, v ectors p k are obtained b y A -orthogonalization of residua, whic h leads to the follo wing algorithm: Algorithm 1 (Conjugate Gradien ts) Giv en x 0 set r 0 = b Ax 0 : F or k = 1 ; : : : do: k = r T k 1 r k 1 p k = r k 1 + k k 1 p k 1 ( p 1 = r 0 ) k = k p T k Ap k x k = x k 1 + k p k r k = r k 1 k Ap k 10

PAGE 20

Among iteration methods the conjugate gradien t method is rather exceptional. Assuming roundo errors are not presen t, it giv es the exact solution of the system after a nite n um ber of iteration, whic h equals at most the n um ber of unkno wns. F or large systems, ho w ev er, performing so man y iterations w ould be impractical. F ortunately m uc h smaller n um ber of iterations is often needed to get a good enough appro ximate solution. This is not true if the condition n um ber of A is large. The performance then can be impro v ed b y in troducing a preconditioner M = E T E where E is some nonsingular matrix. W e poin t out that M dened in this w a y is symmetric and positiv e denite. Preconditioner M m ust ha v e t w o essen tial properties: matrix E T AE 1 is better conditioned than A system Md = e can be easily solv ed. E T denotes the transpose of the in v erse of E The preconditioned conjugate gradien t (PCG) method is the conjugate gradien t method applied to the system E T AE 1 ( Ex ) = E T b: This leads to the follo wing algorithm. Algorithm 2 (Preconditioned Conjugate Gradien ts) Giv en x 0 set r 0 = ( b Ax 0 ) : F or i = 1 ; : : : do: z k 1 = M 1 r k 1 k = r T k 1 z k 1 11

PAGE 21

p k = z k 1 + k k 1 p k 1 ( p 1 = z 0 ) k = k p T k Ap k x k = x k 1 + k p k r k = r k 1 k Ap k W e poin t out that matrix M 1 is not formed explicitly and, instead, a linear system with matrix M is solv ed in eac h iteration of PCG method. Standard preconditioners used in practice include diagonal scaling and incomplete factorizations (e.g., Cholesky) [38]. 2.2 Abstract Sc h w arz Methods Abstract Sc h w arz methods tak e their name after Herman Sc h w arz, a German mathematician of the past cen tury who used an alternating method to construct harmonic functions on regions with non-smooth boundaries. They pro vide an abstract framew ork whic h allo ws analysis of man y dieren t domain decomposition methods. Here w e briery describe the general abstract algorithm and refer to literature for details on analysis and practical/non-practical c hoices of particular methods. This presen tation is based on [6 ] and [23 ]. Let V be a Hilbert space with an inner product h ; i W e also consider the inner product a ( ; ) implied b y the bilinear form of the v ariational form ulation of our problem: Find u 2 V so that a ( u; v ) = h f; v i 8 v 2 V: Let V i ; i = 0 ; : : :; n be closed subspaces of V that form a decomposition of V i.e. V = V 0 + V 1 + V 2 + : : : + V n : 12

PAGE 22

In addition to the inner product implied b y the the bilinear form a w e consider the inner products a i ( ; ) implied b y some bilinear forms a i dened on V i V i ; i = 0 ; : : :; n The Sc h w arz additiv e method can be then written as follo ws. Algorithm 3 (Additiv e Sc h w arz) Giv en an initial appro ximation u 0 start with k = 0 and do: 1. Find w i 2 V i suc h that a i ( w i ; v i ) = h f; v i i a ( u k ; v i ) ; 8 v i 2 V i : 2. Dene the next iterate as u k +1 = u k + n X i =0 w i : 3. Set k = k + 1 and go to 1. The Sc h w arz m ultiplicativ e method is then dened b y the follo wing algorithm. Algorithm 4 (Multiplicativ e Sc h w arz) Giv en an initial appro ximation u 0 starting with k = 0 do: 1. F or i = 0 ; : : :; n (i) nd w i 2 V i so that a i ( w i ; v i ) = h f; v i i a ( u k + i= ( n +1) ; v i ) ; 8 v i 2 V i ; (ii) set u k +( i +1) = ( n +1) = u k + i= ( n +1) + w i : 2. Set k = k + 1 and go to step 1. The algorithms abo v e can be com bined to form h ybrid methods. They are usually used as preconditioners in the conjugate gradien t method. V arious 13

PAGE 23

c hoices of the subspaces V i yield a wide v ariet y of dieren t algorithms. Since the m ultiplicativ e method abo v e is not symmetric, it is usually symmetrized b y considering the sequence of subspaces in forw ard and rev erse order; that is V 0 ; V 1 ; : : :; V n 1 ; V n ; V n ; V n 1 ; : : :; V 1 ; V 0 Eac h sub-step of an additiv e and a m ultiplicativ e Sc h w arz method can be in terpreted as an appro ximate projection of the error. Denoting P i the operator corresponding to the i -th sub-step of an additiv e method, that is the operator that satises a i ( P i w; v ) = a ( w; v ) 8 v 2 V i ; one step of the additiv e method satises e k +1 = ( I n X i =0 P i ) e k ; where e k = u u k is the error of the k -th appro ximation. Similarly the reduction of the error of one step of the m ultiplicativ e method is go v erned b y the equation e k +1 = ( I P n ) : : : ( I P 1 )( I P 0 ) e k : W e refer to [25] and references included there for abstract analysis. Here w e only summarize the main results. They are due to a n um ber of authors: Dryja and Widlund [25 24 ]; Nepomn y asc hikh [55]; Bram ble, P asciak, W ang, Xu [9, 70]; Lions [44 ]; Bjrstad and Mandel [6 ], and others. Let C 0 be the constan t suc h that for all u 2 V there exists a represen tation u = P n i =0 u i ; u i 2 V i and n X i =0 a i ( u i ; u i ) C 0 a ( u; u ) : 14

PAGE 24

Let B = [ b ij ] be the matrix of strengthened Cauc h y-Sc h w arz coecien ts, j a ( v i ; v j ) j b ij a ( v i ; v i ) 1 = 2 a ( v j ; v j ) 1 = 2 8 v i 2 V i ; 8 v j 2 V j ; i; j = 1 ; : : :; n: F urthermore, let be the constan t suc h that a ( u i ; u i ) a i ( u i ; u i ) 8 u i 2 V i ; i = 0 ; : : :; n: Then the smallest eigen v alue of the operator P P i of the abstract additiv e method used a preconditioner in the conjugate gradien ts method is bounded b y 1 =C 0 and the largest b y ( ( B )+1). The bound on the condition n um ber of the symmetric m ultiplicativ e method used as a preconditioner in the conjugate gradien ts method is C 0 (1 + 2^ ( B )) = (2 ^ ), where ^ = max(1 ; ). The space V 0 is usually regarded as a coarse space. It has been sho wn that without the presence of the coarse space, the condition n um ber of an abstract Sc h w arz method gro ws with the n um ber of subspaces. The coarse space ensures a mec hanism of global exc hange of information otherwise lac king in decompositions in v olving a large n um ber of subdomains. 2.3 Ov erlapping Domain Decomposition The rst kno wn domain decomposition method is due to Sc h w arz [63] and is kno wn as the Sc h w arz alternating method. It divides the domain n, on whic h the problem is dened, in to t w o o v erlapping regions n 1 and n 2 and performs a m ultiplicativ e Sc h w arz t ype algorithm. Let us consider the model problem (2.1) A u = f in n ; u = 0 on @ n : 15

PAGE 25

Let n 1 and n 2 be o v erlapping subdomains of n, n 1 [ n 2 = n. F urthermore, let 1 = @ n 1 \ n 2 and 2 = @ n 2 \ n 1 Algorithm 5 (Alternating Sc h w arz) Giv en an initial guess u 0 2 for the v alues on n 2 for k = 1 ; 2 ; : : : do: 1. Solv e A u k 1 = f in n 1 ; u k 1 = 0 on @ n 1 n 1 ; u k 1 = u k 1 2 j 1 on 1 ; 2. Solv e A u k 2 = f in n 2 ; u k 2 = 0 on @ n 2 n 2 ; u k 2 = u k 1 j 2 on 2 : T o be precise, u k 1 j 1 and u k 2 j 2 stand in the algorithm abo v e for the trace of u k 1 on 1 and u k 2 on 2 respectiv ely rather than the simple restrictions. Appro ximate solutions of the model system can be assem bled as u k = u k 2 in n 2 ; u k = u k 1 in n 1 n n 2 : The v ariational form ulation of the algorithm abo v e is due to Lions [44]. It is the m ultiplicativ e abstract Sc h w arz algorithm using the subspaces V 1 and V 2 of V = H 1 0 (n) formed b y functions that v anish outside of n 1 and n 2 respectiv ely The bilinear forms a 1 and a 2 are simple restrictions of a 16

PAGE 26

Choosing some other bilinear forms a 1 and a 2 can model the situation when the problems in step 1. and 2. of the algorithm abo v e w ere to be solv ed only appro ximately An additiv e v ersion of the algorithm has been described b y Dryja [20] and Matsokin and Nepomn y asc hikh [54]. In the discrete, m ultidomain case, o v erlapping methods often start with a non-o v erlapping partition of n in to subdomains (substructures) n i i = 1 ; : : :; N s the diameter of the subdomains being of order H Then, eac h subdomain n i is extended to a larger region n 0 i It is assumed that neither n i nor n 0 i cuts across an y of the elemen ts. The o v erlaps of subdomains n 0 i i = 1 ; : : :; N s are said to be generous, if the distance bet w een boundaries of n i and n 0 i is greater than some xed fraction of H The subspaces V i ; i = 1 ; : : :; n of the abstract Sc h w arz method are then c hosen as restrictions of the nite elemen t space V h (n) to the subregions n 0 i That is, V i = V h \ H 1 0 (n 0 i ) ; i = 1 ; : : :; N s : Appro ximate or exact solv ers then can be c hosen on the subspaces and additiv e, m ultiplicativ e or h ybrid t w o-lev el Sc h w arz methods can be used. As far as parallelization is concerned, additiv e v ersions clearly gain an edge, because they can be easily parallelized. Multiplicativ e v ersions alw a ys in v olv e sequential steps. Some parallelization is possible when m utually non-o v erlapping subdomains n 0 i are grouped together. This process is often called coloring of subdomains. Without a coarse space, the condition n um ber of suc h methods is of the order 1 =H 2 [24], where H is the diameter of the subdomains. Dryja and 17

PAGE 27

Widlund in [24] propose to use the space of con tin uous piecewise linear functions on the coarse mesh dened b y the subdomain n i as the coarse space V 0 They pro v e that, in the case of generous o v erlap, the condition n um ber relev an t for the conjugate gradien t iteration is uniformly bounded. Mandel sho ws in [47] that a similar method b y Co wsar [16], whic h uses discrete harmonic extensions determined b y piecewise constan t v alues on in terfaces of subdomains, leads to a poly-logarithmic bound on the condition n um ber. In their later paper [21], Dryja and Widlund men tion, that their n umerical experimen ts indicate that the con v ergence rate is often satisfactory ev en for small o v erlaps. Running time is often smallest when the o v erlaps are at minim um. The n um ber of conjugate gradien t iterations is higher in suc h a case, because the condition n um ber deteriorates, but this is being compensated for b y the fact that the local problems are smaller and therefore c heaper to solv e. Also, they are better conditioned and, if they are solv ed b y iteration solv ers, the rate of con v ergence is faster. 2.4 Substructuring Methods 2.4.1 P oincar e-Steklo v Operators and Sc h ur Complemen t Substructuring methods borro w their name form structural engineering. It w as in the con text of structural engineering that sev eral substructuring algorithms ha v e been pioneered. W e shall concen trate on substructuring methods on in terfaces of subdomains. Let us consider the model problem (2.1) and let n be divided in to t w o non-o v erlapping subdomains n 1 and n 2 (Figure 2.4.1). Our goal is to solv e the model problem only on the subdomains. Let = @ n 1 \ @ n 2 be the subdomain 18

PAGE 28

W 2 W G 1 nFigure 2.1. Model problem in terface. It is w ell kno wn that the solution u 2 H 1 0 (n) exists for all f 2 L 2 (n). The solution satises con tin uit y conditions on the subdomain in terface: con tin uit y of ruxes of the solution and con tin uit y of the solution. Let us in troduce a few notations. F or k = 1 ; 2, w e dene the space H 1 D (n k ) as the subspace of H 1 (n k ) suc h that the functions in H 1 D (n k ) v anish on @ n k \ @ n. Let H 1 D (n k ) be the dual of H 1 D (n k ). Let r k be the trace operator mapping functions in H 1 D (n k ) to their traces on Let H 1 = 2 00 () be the fractional order Sobolev space on consisting of traces of functions in H 1 D (n k ) and let ( H 1 = 2 00 ()) 0 denote its dual. Using the con tin uit y of ruxes, w e will split the problem in to t w o subproblems, k = 1 ; 2: A u k = f k in n k n ([ a ij ] r u k ) = ( 1) k +1 g ? on u k = 0 on @ n \ @ n k : The v ector n is the normal v ector to orien ted, for example, from n 1 to n 2 In v ariational terms, these are the problems of nding u k 2 H 1 D (n k ) ; k = 1 ; 2 19

PAGE 29

suc h that a k ( u k ; v k ) = Z g ? v k + Z n k f k v k ; 8 v k 2 H 1 D (n k ) ; (2.4) where a k ( v; w ) = d X i;j =1 Z n k a ij @v @x i @w @x j : A t this poin t it w ould be possible to design an iterativ e method w orking on the subproblems using the con tin uit y conditions. Usually ho w ev er, w e rst reduce the problem to a problem on the subdomain in terface using P oincar e-Steklo v operators. First, w e will look for the unkno wn Neumann data g ? on W e dene the P oincar e-Steklo v operators Q k : ( H 1 = 2 00 ()) 0 H 1 = 2 00 () ; k = 1 ; 2 b y Q k g ? = r k u k ; (2.5) where, for g ? 2 ( H 1 = 2 00 ()) 0 u k is the solution of (2.4) with f k = 0. Suc h u k is the harmonic function satisfying the Neumann condition giv en b y g ? In other w ords, the P oincar e-Steklo v operator maps the Neumann boundary condition in to the corresponding Diric hlet boundary condition Q k : @u k @ n r k u k : F urthermore, w e dene R k : H 1 D (n k ) H 1 = 2 00 () ; k = 1 ; 2 b y the equation R k f k = r k u k ; where for f k 2 L 2 (n k ), u k is the solution of (2.4) with g ? = 0. In terms of the P oincar e-Steklo v operators, the problem is to nd the solution g ? suc h that ( Q 1 + Q 2 ) g ? = R 2 f 2 R 1 f 1 : (2.6) 20

PAGE 30

That is to nd the Neumann data g ? on suc h that the traces of the solutions u k ; k = 1 ; 2 of (2.4) coincide on The second possibilit y is to enforce con tin uit y of the solution on the boundary a priori. This yields a dual form ulation to (2.6). Consider the Diric hlet problems, k = 1 ; 2, A u k = f k in n k u k = g on u k = 0 on @ n \ @ n k : W e will look for the unkno wn Diric hlet data g on suc h that the ruxes are con tin uous for the solutions u k ; k = 1 ; 2 of the problems abo v e. Using the P oincar e-Steklo v operators, this problem can be written as ( Q 1 1 + Q 1 2 ) g = Q 1 1 R 1 f 1 + Q 1 2 R 2 f 2 : (2.7) The equations abo v e and iterativ e methods for solving them ha v e been studied, for example, b y Agoshk o v [2] and in a mixed-method setting b y Glo winski and Wheller [36]. The paper b y Bakh v alo v and Kn y azev [4] is concerned with highly discon tin uous coecien ts bet w een the subdomains. As the case of t w o subdomains is hardly of practical importance, w e conclude this section b y in troducing a m ulti-domain discrete analog of (2.7). An analog of (2.6) is described in the next c hapter. Let us recall no w the concept of Sc h ur complemen t. Denition 6 (Sc h ur complemen t) Let A be an n n matrix and ; f 1 ; : : :; n g be index sets. W e denote A ( ; ) the submatrix that lies in the ro ws of A indexed b y and the columns indexed b y and A ( 0 ; 0 ) the submatrix giv en b y deleting the ro ws indexed b y and the columns giv en b y 21

PAGE 31

Let A ( ; ) be nonsingular. Then the matrix A ( 0 ; 0 ) A ( 0 ; ) A ( ; ) 1 A ( ; 0 ) is called the Sc h ur complemen t of A ( ; ) in A W e note that the Sc h ur complemen t is the matrix that arises when eliminating the unkno wns x ( ) from the equation Ax = b Then, b ( 0 ) A ( 0 ; ) A ( ; ) 1 b ( ) = ( A ( 0 ; 0 ) A ( 0 ; ) A ( ; ) 1 A ( ; 0 )) x ( 0 ) : W e also observ e that the Sc h ur complemen t S of A ( ; ) in A has the propert y h Sy; y i = inf x 2 I R n ;x ( 0 )= y h Ax; x i : (2.8) Let n be a domain in I R 2 or I R 3 decomposed in to N s non-o v erlapping subdomains n 1 ; n 2 ; :::; n N s Let u i be the v ector of degrees of freedom for subdomain n i corresponding to a conforming nite elemen t discretization of the second order elliptic problem (2.1) dened on n, suc h that eac h subdomain is a union of some of the elemen ts. Let u i K i and f i be the v ector of degrees of freedom, the local stiness matrix, and the load v ectors, respectiv ely associated with the subdomain n i Let L i denote the zero-one assem bly matrix mapping the subdomain degrees of freedom u i in to global degrees of freedom u that is u i = L T i u The stiness matrix then is K = N s X i =1 L i K i L T i and the load v ector f = N s X i =1 L i f i : 22

PAGE 32

The problem to be solv ed, Ku = f; can be reduced to an in terface problem b y splitting the degrees of freedom in to in terface and in terior degrees of freedom. Let us assume that the in terior degrees of freedom are listed rst. The subdomain stiness matrices K i and the restriction matrices L i can be split accordingly: u i = 2 6 6 4 u i u i 3 7 7 5 ; K i = 2 6 6 4 K i ~ K T i ~ K i K i 3 7 7 5 ; L i = [ L i ; L i ] : By eliminating the in terior degrees of freedom, w e obtain the Sc h ur complemen ts of K i in K i S i = K i ~ K i K 1 i ~ K T i : (2.9) The problem then reduces to the in terface problem S u = f for the in terface unkno wns u with the global Sc h ur complemen t S = N s X i =1 L i S i L T i and the in terface righ t hand side f = N s X i =1 L i ( f i ~ K i K 1 i f i ) : In most cases, explicit computation of S i w ould be too expensiv e, but for gradien t-descen t methods, only ev aluation of the action of S i is necessary This ev aluation can be performed ecien tly b y factorizing K i and it corresponds to solving a Diric hlet problem on ev ery subdomain. The subdomain problems can be solv ed in parallel. 23

PAGE 33

Let us demonstrate that the action of S i can be ev aluated b y solving a Diric hlet problem on the subdomain n i Consider the discretized problem on the subdomain 2 6 6 4 K i ~ K T i ~ K i K i 3 7 7 5 2 6 6 4 u i u i 3 7 7 5 = 2 6 6 4 0 f i 3 7 7 5 ; where the Diric hlet data u i is kno wn. F rom the rst equation, w e nd u i = K 1 i ~ K T i u i This the solution of the Diric hlet problem giv en b y the data u i Then, w e substitute in to the second equation and obtain ( K i ~ K i K 1 i ~ K T i ) u i = S i u i = f i : Man y preconditioners for the reduced problem are based on computations of the action of S 1 i This can be ev aluated b y solving a Neumann problem on the subdomain, for unkno wn u i and u i and giv en g i 2 6 6 4 K i ~ K T i ~ K i K i 3 7 7 5 2 6 6 4 u i u i 3 7 7 5 = 2 6 6 4 0 g i 3 7 7 5 ; where g i corresponds to the Neumann data. Then, u i = S 1 i g i Indeed, from the rst equation, w e nd u i = K 1 i ~ K T i u i Substituting in to the second equation, w e nd ( K i ~ K i K 1 i ~ K T i ) u i = S i u i = g i : W e note that the Sc h ur complemen t S i is the discrete analog of the in v erse of the P oincar e-Steklo v operator Q 1 i : u j @u @ and the in v erse of the Sc h ur complemen t is the analog of Q i : @u @ u j (2.5). One of the adv an tages of the reduced problem is that it impro v es the condition n um ber of the original problem. F or our model problem, assuming 24

PAGE 34

the P oincar e-F riedric hs inequalit y k u h i k 2 1 ; 2 ; n i c p j u h i j 2 1 ; 2 ; n i holds, the subdomain stiness matrix K i satises c k u h i k 2 1 ; 2 ; n i u T i K i u i C k u h i k 2 1 ; 2 ; n i ; where u h i is the nite elemen t function corresponding to the v ector of degrees of freedom u i Then, the subdomain Sc h ur complemen t S i satises (cf. Lemma 26 and Lemma 33) c k u h i k 2 1 2 ; 2 ;@ n i u T i S i u i C k u h i k 2 1 2 ; 2 ;@ n i : This implies that, for a second order problem, the condition n um ber of the reduced problem is of order 1 =h for triangulations of c haracteristic mesh size h while the condition n um ber of the original problem is of the order 1 =h 2 2.4.2 V arious Preconditioners for the Reduced Problem Instead of solving the reduced problem directly iterativ e methods suc h as the preconditioned conjugate gradien t method can be applied in whic h only matrix v ector products are required. Since eac h iteration is quite expensiv e (it requires solving a Diric hlet problem on eac h subdomain), ecien t preconditioning is importan t to k eep the n um ber of iterations small. Diagonal or bloc k diagonal preconditioning of the reduced problem requires kno wing the diagonal of S Instead of computing it directly Chan [13] proposes a \boundary probing" tec hnique to construct an appro ximate diagonal b y ev aluating actions of S on carefully selected v ectors. F or geometrically simple subdomain in 2D, Bram ble, P asciak and Sc hatz in [10 ] propose a method for preconditioning the original problem b y 25

PAGE 35

rst splitting the functions in the nite elemen t space in to discrete harmonics on eac h subdomain and functions whic h v anish on the boundary of eac h subdomain. This is equiv alen t to eliminating the in terior degrees of freedom as in the reduced system. Then the space of discrete harmonic functions is divided in to functions linear on in terfaces and functions whose v alues are zero at subdomain crosspoin ts. The preconditioner is then based on a bilinear form dened on the space splitting. The condition n um ber of the method is sho wn to be bounded b y C (1 + log H=h ) 2 A similar idea is exploited in [3] for the p-v ersion of nite elemen ts and a bound C (1 + log p ) 2 is obtained. A dieren t algorithm [8] b y Bram ble, P asciak and Sc hatz is described for the case of t w o subdomains n 1 and n 2 separated b y the in terface They propose to precondition the original problem in the follo wing w a y: First split the local solution on n 2 in to the componen t P that satises the nonhomogeneous equation with the zero boundary condition on @ n 2 and a discrete harmonic function H that satises the corresponding homogeneous equation. A Diric hlet problem is solv ed on n 2 to obtain P This is follo w ed b y a solution of a mixed Neumann-Diric hlet problem on n 1 using the Neumann data obtained b y solving the problem on n 2 Finally another Diric hlet problem is solv ed to obtain H a discrete harmonic function, that satises the corresponding homogeneous equation. This is sometimes called Neumann-Diric hlet decomposition. This algorithm has also been described in [24 ] and it is in fact preconditioning the reduced system b y either S 1 1 or S 1 2 The in v erses of the Sc h ur complemen ts need not be explicitly computed; the action corresponds to solving a Neumann (or mixed Neumann-Diric hlet) problem on one subdomain and then extending the solution to the other subdomain b y solving a Diric hlet 26

PAGE 36

problem with the Diric hlet data found from the Neumann problem. The idea of the Neumann-Diric hlet preconditioner can be extended to m ultiple subdomain case when there is a red-blac k ordering of the subdomains. Sev eral methods solving an equation with the global Sc h ur complemen t are proposed b y Mandel [47]. The methods are set in the space of all discrete harmonic functions on the union of subdomain in terfaces = [ i @ n i @ n. The methods are h ybrid Sc h w arz methods using the coarse space V 0 in a m ultiplicativ e fashion and other spaces are treated additiv ely A method of this kind can also be in terpreted as a t w o lev el v ariational m ultigrid. The methods use V i that are associated with subdomains or globs. A glob is a v ertex, an edge that does not con tain its endpoin ts, or a face of a subdomain in terface. Piecewise linear functions dened on the subdomain triangulation or a piecewise constan t space based on glob-wise a v eraging of subdomain v alues is used as the coarse space. P oly-logarithmic bounds are obtained. Man y other c hoices of the coarse space are discussed b y Dryja, Smith and Widlund [23 ]. Their exhaustiv e in v estigation comprises v ertex based coarse spaces, based on piecewise linear functions on substructures used as elemen ts, wire bask et algorithms, that use a v erages on substructures, and face based algorithms. 2.4.3 The Neumann-Neumann and Balancing Domain Decomposition Methods The balancing domain decomposition (BDD) is based on the so called Neumann-Neumann preconditioner that preconditions the reduced problem b y a w eigh ted sum of in v erses of the local Sc h ur complemen ts. It is called 27

PAGE 37

Neumann-Neumann because it corresponds to solving Neumann problems only on subdomains (as opposed to the Diric hlet-Neumann preconditioner, that uses both Diric hlet and Neumann problems). Since this is a preconditioner in a w a y dual to the method that is the subject of this thesis, w e describe it in a little more detail. T o describe the preconditioner, w eigh t matrices D i ; i = 1 ; : : :; N s satisfying the decomposition of unit y I = N s X i =1 L i D i L T i are used. A simple c hoice for D i is a diagonal matrix with the diagonal elemen ts being the reciprocals of the n um ber of subdomains the degree of freedom is associated with. The Neumann-Neumann method used as a preconditioner of the problem is as follo ws: Giv en the residual r distribute it to subdomains r i = D T i L T i r solv e the local problems S i u i = r i and a v erage the results M 1 r = P N s i =1 L i D i u i [7, 43 ]. The dra wbac k of the Neumann-Neumann preconditioner is that it lac ks a mec hanism of distributing the error globally This has been resolv ed b y adding a coarse space to the problem. The resulting method is called the balancing domain decomposition. Let Z i be the matrix with linearly independen t columns that generate the k ernel of K i Im Z i = Ker K i If K i is regular, Z i is a v oid matrix. The balancing preconditioner is as follo ws: 1. Balance the original residual r b y solving the auxiliary problems Z T i D T i L T i ( r S N s X j =1 L j D j Z j v j ) = 0 ; i = 1 ; : : :; N s 2. Distribute the balanced residual to the subdomains and nd a solution of 28

PAGE 38

the local problems S i u i = D T i L T i ( r S N s X j =1 L j D j Z j v j ) ; i = 1 ; : : :; N s 3. Balance b y solving the auxiliary problems Z T i D T i L T i ( r S N s X j =1 L j D j ( u j + Z j w j ) = 0 ; i = 1 ; : : :; N s 4. Av erage the result on the in terface M 1 r = N s X i =1 L i D i ( u i + Z i w i ) The BDD is pro v en to be independen t of n um ber of subdomains, and its condition n um ber is independen t of jumps of coecien ts bet w een subdomains with appropriate w eigh t matrices D i It can be written also as an abstract additiv e Sc h w arz algorithm [24 41]. 2.4.4 The Neumann-Neumann Domain Decomposition Algorithm for Plates and Shells Let us in troduce the Kirc hho-Lo v e model of plate bending follo wing [41]. W e consider a plate occup ying a domain in I R 2 whic h is clamped on the part of boundary @ n 00 and simply supported on @ n 0 n @ n 00 The Kirc hhoLo v e model plate model c haracterizes the v ertical displacemen t u 2 V of the plate as the solution of the v ariational problem a ( u; v ) = F ( v ) ; 8 v 2 V; (2.10) where the bilinear form on the righ t hand side is symmetric, con tin uous and coerciv e a ( u; v ) = Z n ( ( u )) : K : ( ( v )) ; 29

PAGE 39

the functional on the righ t hand side is giv en b y F ( v ) = Z n fv + Z @ n @ n 00 m g @v @ n + Z @ n @ n 0 gv: The space of kinematically admissible elds is V = f v 2 H 2 (n) ; v = 0 on @ n 0 ; v = @v @ n = 0 on @ n 00 g : In the denitions abo v e f is the densit y of v ertical forces, m g the densit y of rexion momen ts applied on the part of the boundary where the plate is free to rotate, and g is the densit y of v ertical boundary loading. The sym bol ":" denotes a tensor product, ( ) = 1 2 ( r + ( r ) T ) is the curv ature tensor, ( u ) = r u represen ts the in-plane notation of the plate and K is the plate rexural stiness. K is symmetric, elliptic and con tin uous in the sense that ( ( u )) : K : ( ( u )) ct 3 j ( ( u )) j 2 ( ( u )) : K : ( ( v )) Ct 3 j ( ( u )) jj ( ( v )) j ; where t is the plate thic kness. F or a simple case of an isotropic plate made of an homogeneous material with Y oung modulus E and P oisson coecien t it is giv en b y ( ( u )) : K : ( ( u )) = Et 3 12(1 2 ) ((1 ) r 2 u : r 2 v + u v ) : Since V is a subspace of H 2 (n), the appropriate nite elemen t spaces are C 1 con tin uous. Examples of suc h elemen ts include the discrete Kirc hho triangle (DKT) and HCT elemen ts. The Neumann-Neumann and BDD preconditioner do not perform w ell for the plate bending problem [41]. The condition n um ber estimate of BDD is 30

PAGE 40

based on the estimate [45 ] sup 8 < : P N s j =1 k L T j P N s i =1 L i D i u i k 2 S j P N s i =1 k u i k 2 S i : u i ? Ker S i ; S i u i ? Im Z i 9 = ; (2.11) A t crosspoin ts of subdomains, the v ectors in the estimate are constructed from con tributions sev eral subdomains dieren t from those that share an edge leading to the crosspoin t. This leads to a discon tin uit y that is inappropriate for a fourth order problem. The BDD method for plates [42] a v oids this problem arising at subdomain crosspoin ts b y enhancing the coarse space of the balancing domain decomposition algorithm. The coarse space is again Im Z i where Z i = [ x i 1 ; : : :; x in i ; y i 1 ; : : :y im i ] : f x i 1 ; : : :; x in i g is a basis for Ker S i and for eac h crosspoin t j = 1 ; : : :; m i of the subdomain n i y ij is the solution of the problem S i y ij = e ij with e ij the v ector corresponding to the unit normal load at the crosspoin t j With this c hoice, the normal displacemen t componen t of the v ectors u i coming out of the coarse space problem is zero since S i u i ? y ij implies that u i ? S i y ij = e ij Then, the suprem um (2.11) is tak en o v er functions with zero at endpoin ts of the edges, whic h mak es it possible to pro v e a poly-logarithmic bound. F or another approac h, imposing zeros at crosspoin t directly see [41]. 2.4.5 Lagrange Multipliers and P oincar e-Steklo v Operators F ollo wing [19], w e sho w ho w Lagrange m ultiplier approac h to enforcing solution con tin uit y is related to in terface form ulations using P oincar eSteklo v operators dened in Section 2.4.1. The discrete m ulti-domain case will be treated in the next c hapter. 31

PAGE 41

W e consider the problem from Section 2.4.1 and the notation in troduced there. W e reform ulate our problem as a constrained minimization problem: Find the solution ( u 1 ; u 2 ) 2 H 1 D (n 1 ) H 1 D (n 2 ) that minimizes 1 2 2 X k =1 a k ( u k ; u k ) Z n k fv k subject to the condition r 1 u 1 = r 2 u 2 Then, for eac h ( u 1 ; u 2 ; g ? ) 2 H 1 D (n 1 ) H 1 D (n 2 ) ( H 1 = 2 00 ()) 0 w e dene the Lagrangian ( u 1 ; u 2 ; ? ) = 1 2 2 X k =1 a k ( u k ; u k ) Z n k fu k Z ? ( u 1 u 2 ) : The Lagrange m ultiplier ? is used to to enforce the con tin uit y of the solution on the boundary The rest of the Lagrangian is the usual quadratic functional implied b y the w eak form of the problem. The critical poin ts ( u 1 ; u 2 ; ? ) of the saddle poin t problem of no w m ust satisfy the v ariational equalit y 2 X i =1 a k ( u k ; v k ) Z ( ? ( v 1 v 2 ) + ? ( u 1 u 2 ) d = 2 X i =1 Z n i fv k ; for all ( v 1 ; v 2 ; ? ) 2 H 1 D (n 1 ) H 1 D (n 2 ) ( H 1 = 2 00 ()) 0 This problem has a unique solution [19] whic h is the solution of the problem ( Q 1 + Q 2 ) ? = R 2 f 2 R 1 f 1 This is the equation (2.6). It sho ws that solving the Lagrange m ultiplier form ulation is equiv alen t to nding Neumann in terface data on (cf. [18]). The paper [19] uses the Lagrange form ulation to in troduce nite elemen t spaces of Lagrangians of small dimension per in terface for regular meshes. This can reduce the size of the problem substan tially but it is restricted to regular meshes. The space of Lagrangians can be c hosen as the restriction of the nite elemen ts space on the subdomains to the in terface. This is similar to the approac h tak en b y FETI as explained in the next c hapter. 32

PAGE 42

2.5 A Tw o Lev el Method Based on Smoothed Prolongation The t w o-lev el method described in [68 66 ] dev elops a simple abstract framew ork based on the concept of smoothed ten tativ e prolongator in troduced in [65]. The ten tativ e prolongator is deriv ed from a system of nono v erlapping subdomains. As opposed to the previously described methods, the union of all subdomains n i i = 1 ; : : :; N s does not co v er whole n. Instead, there is a la y er one elemen t wide bet w een eac h t w o subdomains. Our algorithm can be written as a v ariational t w o-lev el m ultigrid with a special c hoice of componen ts. W e will rst describe componen ts of the method and abstract assumptions that ensure coarse space size independen t con v ergence. Let us consider the system of linear algebraic equations Ku = f; where K is an n n symmetric positiv e denite stiness matrix arising from a discretization of a second order elliptic PDE, for example (2.1). In terpolation from the coarse space to the ne space is represen ted b y an operator MP composed from a ten tativ e prolongator P and prolongator smoother M Let us denote K M = M 2 K; M 0 = I K M ; where 2 (0 ; 2) is a giv en constan t and ( K M ) is an estimated upper bound for ( K M ), the spectral radius of K M The follo wing assumption species requiremen ts on M P and ( K M ). are needed for pro ving coarse space size independen t con v ergence. 33

PAGE 43

Assumption 7 There exist positiv e constan ts C 1 C 2 independen t of n and N s and a positiv e constan t C D ( N s ; n ) suc h that (1) F or ev ery u 2 I R n there exists v 2 I R N s suc h that k u Pv k C 1 C D ( N s ; n ) 1 = 2 ( K ) k u k K : (2.12) (2) The prolongator smoother M is symmetric, comm utes with K ( M ) < 1, and ( K M ) ( K M ) C 2 2 C 2 D ( N s ; n ) ( K ) : (2.13) W e also need to dene a pre-smoother S M and post-smoother S M 0 ; S M ( u; f ) = Mu + Nf; and analogously for M 0 The matrix N is c hosen so the smoother is consisten t with the system Ku = f The smoothers S M and S M 0 are relaxation operators used to smooth the appro ximate solution corrected b y an in terpolated coarse lev el error. The algorithm can no w be written as follo ws. Algorithm 8 Giv en an initial guess u repeat un til con v ergence: (1) u S M ( u; f ), (2) solv e ( P T MKMP ) w = P T M ( Ku f ) ; (3) u u MPw (4) u S M 0 ( u; f ). P ost process u S M ( u; f ). The ten tativ e prolongator P ma y be dened as P = D 1 = 2 ~ P; where D = diag( K ) and ~ P ij = 1 ; if the node corresponding to u i belongs to subdomain n j ; = 0 ; otherwise. 34

PAGE 44

The most practical c hoice of M is a polynomial in K The c hoice described belo w is a polynomial in K that attempts to minimize ( K M ) giv en an upper bound on ( K ) and a c hosen degree of the polynomial. F or certain degrees, this polynomial can be found explicitly b y the follo wing recursiv e algorithm: Algorithm 9 Let ^ be the estimate of ( K ) satisfying ( K ) ^ C ( K ), with a giv en constan t C Set ^ i = ^ 9 i A 0 = K and for i = 1 ; 2 ; : : : do: (1) A i = ( I 4 3 ^ 1 i 1 A i 1 ) 2 A i 1 (2) M i = Q i 1 j =0 ( I 4 3 ^ 1 j A j ) Notice that deg( M i ) 3 i 1 2 F or a 2 D problem, w e c hoose the prolongator smoother M = M k where deg( M k +1 ) qN 1 = 2 es deg( M k ) ; (2.14) where q 2 (0 ; 1] is a giv en parameter and N es the a v erage a v erage n um ber of degrees of freedom per subdomain, and M k are the polynomials constructed b y the algorithm abo v e. With the c hoices of the ten tativ e prolongator and prolongator smoother described abo v e, the con v ergence of our method can be sho wn to be independen t of the meshsize, N es inhomogeneities bet w een subdomains, and boundary conditions [68 ]. It th us o v ercomes one of the main disadv an tages of standard t w o-lev el m ultigrid methods, suering from the dependence of the rate of conv ergence on the size of the coarse space. On the other hand, its computational complexit y is lo w er in comparison with domain decomposition methods using direct local solv ers. This is true ev en for simple direct solv ers as in ternal solv ers in our method (whic h are preferred, as iterativ e solv ers tend to ha v e negativ e impact on the robustness of the method). 35

PAGE 45

Application of the method to solid problems is treated in [64] and similar results are obtained. The method is pro v ably robust with respect to jumps in coecien ts [67]. 36

PAGE 46

3. Finite Elemen t T earing and In terconnecting (Deriv ation and Extensions) 3.1 Original FETI The FETI (Finite Elemen t T earing and In terconnecting) method is a domain decomposition algorithm deriv ed from a h ybrid v ariational principle and designed for the iterativ e solution of systems of equations arising from the nite elemen t discretization of self-adjoin t elliptic partial dieren tial equations. It w as dev eloped in [28 27 33], and also discussed in detail in monograph [34]. In this method, a giv en spatial domain is \torn" in to non-overlapping subdomains where an incomplete solution of the primary eld is rst ev aluated using a direct solv er. Next, in tersubdomain eld con tin uit y is enforced via Lagrange m ultipliers applied at the subdomain in terfaces. This \gluing" phase generates a smaller size symmetric dual problem where the unkno wns are the Lagrange m ultipliers, and whic h is best solv ed b y a preconditioned conjugate gradien t (PCG) algorithm. This idea is related to the ctitious domain method where the Lagrange m ultipliers enforce boundary conditions as in Dinh et al. [17]. In con trast with other related domain decomposition methods using Lagrange m ultipliers as unkno wns [36 62], the FETI method distinguishes itself with the treatmen t of the n ull spaces of the subdomain stiness matrices (rigid 37

PAGE 47

body modes) associated with the so-called roating subdomains, i.e., subdomains without a sucien t n um ber of essen tial boundary conditions to prev en t the local stiness matrix from being singular. Resolving the rigid body modes leads to a small \coarse" problem that is solv ed in eac h PCG iteration. This is an added complication, but also a blessing. F arhat, Mandel, and Roux [32] ha v e sho wn n umerically and pro v ed for the FETI method without preconditioning, that the auxiliary problem pla ys the role of a coarse problem, namely it causes the condition n um ber to be bounded independen tly of the n um ber of subdomains. In [26], F arhat, Chen, and Mandel extended to time-dependen t problems, whic h lac k the naturally occurring coarse problem. The FETI method is in a sense dual to the Neumann-Neumann method with a coarse problem, dev eloped b y Mandel under the name Balancing Domain Decomposition [45 ] (BDD), whic h w e described in Chapter 2 and whic h is based on an earlier method of de Roec k and Le T allec [61 ]. A modied method w as analyzed b y Dryja and Widlund [25]. It should be noted that while the underlying ideas of FETI and BDD are in a w a y dual, FETI is not the BDD method applied to the dual problem. 3.1.1 Problem Setting and Assumptions Let us rst in troduce some notation used throughout this c hapter and the rest of this thesis. F or u; v 2 I R n the inner product h u; v i = u T v serv es also as dualit y pairing. The ` 2 norm is denoted k u k = h u; u i 1 = 2 F or a symmetric positiv e semidenite matrix A k u k A = h Au; u i 1 = 2 the induced seminorm. This is a norm if A is positiv e denite. The superscript + denotes pseudoin v erse, dened as follo ws. 38

PAGE 48

Denition 10 (Pseudoin v erse) Let A be a linear operator. The pseudoin v erse A + is an y linear operator suc h that if a 2 Im A then AA + a = a In general, a pseudoin v erse is not unique. Our algorithms will be in v arian t to a specic c hoice of the pseudoin v erse. If A is a symmetric operator on a nite dimensional space, A + can be c hosen to be also symmetric. Considering the spectral decomposition of A A = X v v T ; Av = v ; v T v = 1 ; (3.1) a pseudoin v erse A + can be c hosen as A + = X 6 =0 1 v v T : F or A positiv e semidenite, w e denote A = X > 0 v v T : In particular, with this notation w e ha v e, A = A 1 = 2 A 1 = 2 and A + = A 1 = 2 A 1 = 2 W e poin t out that Ker A = Ker A for an y real W e refer to Section 2.2 and the model problem (2.1) for explanation of some of the terminology used next. Let n be a domain in I R 2 or I R 3 decomposed in to N s non-o v erlapping subdomains n 1 ; n 2 ; :::; n N s Let u i be the v ector of degrees of freedom for subdomain n i corresponding to a conforming nite elemen t discretization of an elliptic problem (e.g. linear elasticit y plate bending) dened on n, suc h that eac h subdomain is a union of some of the elemen ts. Let u i K i and f i be the v ector of degrees of freedom, the local stiness matrix, and the load v ectors, 39

PAGE 49

respectiv ely associated with the subdomain n i W e will use the bloc k notation u = 2 6 6 6 6 6 6 6 6 6 6 4 u 1 u 2 . u N s 3 7 7 7 7 7 7 7 7 7 7 5 ; f = 2 6 6 6 6 6 6 6 6 6 6 4 f 1 f 2 . f N s 3 7 7 7 7 7 7 7 7 7 7 5 ; K = 2 6 6 6 6 6 6 6 6 6 6 4 K 1 0 : : : 0 0 K 2 : : : 0 : : : : : : : : : : : : 0 0 : : : K N s 3 7 7 7 7 7 7 7 7 7 7 5 : Depending on the imposed boundary conditions and the location of the subdomain, the local stiness matrix K i is positiv e denite or positiv e semidenite. A subdomain without sucien t essen tial boundary conditions to prev en t the subdomain stiness matrix K i from being singular is called a ro ating sub domain Let Z i be the matrix with linearly independen t columns that generate the k ernel of K i Im Z i = Ker K i If K i is regular, Z i is a v oid matrix. Denote Z = 2 6 6 6 6 6 6 6 6 6 6 4 Z 1 0 : : : 0 0 Z 2 : : : 0 : : : : : : : : : : : : 0 0 : : : Z N s 3 7 7 7 7 7 7 7 7 7 7 5 ; th us, Im Z = Ker K; Ker Z = f 0 g : A single mesh poin t x 2 n has sev eral degrees of freedom associated with it if it lies on the in terface (in tersection of the boundaries) of t w o or more subdomains, see Figure 3.1. Let B be a giv en matrix suc h that Bu = 0 expresses the condition that for eac h mesh node shared b y t w o or more subdomains the v alues of the degrees of freedom associated with that node coincide. W e denote b y W the space of all v ectors of degrees of freedom, and 40

PAGE 50

W W j i l l u u i m j n xFigure 3.1. Model domain: the mesh node x shared b y subdomains n i and n j and the degrees of freedom u i m and u j n associated with the mesh node x 41

PAGE 51

b y the space of the v ectors of v alues of the con tin uit y constrain t; th us, K : W W ; B : W : The problem to be solv ed is the minimization of the energy of the system subject to in tersubdomain con tin uit y conditions E ( u ) = 1 2 u T Ku f T u min subject to Bu = 0 ; u 2 W : (3.2) W e assume that the global structure is not roating, that is, the solution of (3.2) is unique. F rom (3.2), this is equiv alen t to the assumption that Ker K \ Ker B = f 0 g : (3.3) The follo wing notation will be needed to write the FETI algorithm concisely: G = BZ; F = BK + B T ; d = BK + f; (3.4) e = Z T f; P = I G ( G T G ) 1 G T : Note that P is an ` 2 orthogonal projection on (Im G ) ? = Ker G T Lemma 12 justies the expression for P b y sho wing that G T G is in v ertible. 3.1.2 Deriv ation of the Dual Problem The method, as originally deriv ed b y F arhat and Roux [27 ], in troduces Lagrange m ultipliers to enforce the con tin uit y of the solution. Solving the 42

PAGE 52

problem (3.2) leads to the system of equations Ku + B T = f Bu = 0 (3.5) A solution u of the rst equation in (3.5) exists if and only if f B T 2 Im K: (3.6) It m ust then ha v e the form u = K + ( f B T ) + Z: (3.7) where is to be determined. Substituting u from (3.7) in to the second equation of (3.5) yields BK + ( f B T ) + BZ = 0 : (3.8) No w (3.8) m ultiplied b y P together with (3.6) sho w that satises the system of equations P ( F d ) = 0 G T = e; (3.9) where w e ha v e used notation (3.4). The rst of the equations (3.9) is solv ed b y a projected preconditioned conjugate gradien t method using an initial appro ximation 0 that satises the second equation. The conjugate gradien t algorithm requires ev aluating only the actions of PF Most of the computational w ork in F = BK + B T is concentrated in the action of K + Since K + is subdomain bloc k diagonal, its action can be computed in parallel and in v olv es solving subdomain problems only Application of P requires solving a small coarse problem. F or a scalar problem, the size of this problem is less than the n um ber of subdomains N s F or a 43

PAGE 53

linear elasticit y problem, the size of the problem is the n um ber of rigid body modes, whic h again does not exceed a small m ultiple of N s A preconditioned projected CG for the equation PF = d using a symmetric preconditioner D can be written as follo ws. Algorithm 11 (FETI) Giv en an initial 0 compute the initial estimate 0 = G ( G T G ) 1 e + P 0 and the initial residual r 0 = P ( F 0 d ) : Repeat for k = 1 ; 2 ; : : : un til con v ergence: z k 1 = Dr k 1 y k 1 = Pz k 1 k = r T k 1 y k 1 p k = y k 1 + k k 1 p k 1 ( p 1 = y 0 ) k = k p T k PFp k k = k 1 + k p k r k = r k 1 k PFp k The c hoice of the preconditioner is discussed in Section 3.3. Lagrange m ultipliers in FETI can be seen as interfac e for c es and moments in the ph ysical system. F rom (3.7) and the denition of F in (3.4), the residual P ( F d ) = Bu has the in terpretation of jumps of the values of de gr e es of fr e e dom bet w een subdomains. The condition f B T ? Ker K 44

PAGE 54

means that the action of the loads and in tersubdomain forces and momen ts does not excite rigid body motions. 3.1.3 Saddle P oin t F orm ulation Let us explain in more detail ho w enforcing the con tin uit y using Lagrange m ultipliers yields the system (3.9). This approac h is presen ted b y F arhat, Mandel and Roux [32], and it will be later used to deriv e a generalized FETI method. F or the Lagrangian of the minimization problem (3.2), L ( u; ) = 1 2 u T Ku f T u + T Bu; u 2 W ; 2 ; w e solv e the dual problem: nd suc h that C ( ) = max 2 C ( ) max 2 inf u 2 W L ( u; ) : (3.10) By a direct computation, C ( ) = 8 > > < > > : 1 if f B T 6? Ker K; 1 2 h K + ( f B T ) ; f B T i otherwise : (3.11) The dual problem (3.10) is th us equiv alen t to maximizing C ( ) on the admissible set A = f 2 j C ( ) > 1g : The space of admissible incremen ts is f 1 2 j 1 2 A ; 2 2 Ag = f 2 j B T ? Ker K g = Ker G T : (3.12) A t the maxim um of C ( ), 2 A the deriv ativ e of C D C ( ; ), is zero in all directions in 2 Ker G T : D C ( ; ) = 0 8 2 Ker G T : 45

PAGE 55

By a straigh tforw ard computation, this becomes h BK + B T + BK + f; i = 0 ; 8 2 Ker G T : (3.13) In order to express (3.13) as a linear equation in the space Ker G T w e use the l 2 orthogonal projection dened b y (3.4). Then for 2 Ker G T h BK + B T + BK + f; i = h BK + B T + BK + f; P i = h P ( BK + B T + BK + f ) ; i Therefore, the dual problem (3.10) is equiv alen t to the linear equation in Ker G T for the unkno wn 2 Ker G T ; P ( BK + B T ( + 0 ) + BK + f ) = 0 ; (3.14) where 0 is an arbitrary starting feasible solution, that is, 0 2 A Denoting = 0 + w e obtain the system (3.9). 3.1.4 Algebraic Properties of the Original FETI Method Lemma 12 ( G T G ) 1 exists. Pr o of. Let G = BZ = 0. Then Z 2 Ker B and Z 2 Ker K b y the denition of Z It follo ws from (3.3) that Z = 0, hence = 0, since Z w as assumed to be of full rank. See also [34, Theorem 5.4]. Theorem 13 The solution of (3.9) is unique up to addition of a v ector from Ker B T An y solution of (3.9) yields the same solution u of the minimization problem (3.2), using (3.7) with = ( G T G ) 1 G T ( d F ). Pr o of. The relation bet w een and follo ws b y a direct computation. T o pro v e uniqueness of it is sucien t to sho w that Ker PF \ Ker G T = Ker F \ Ker G T = Ker B T : (3.15) 46

PAGE 56

First, B T = 0 implies F = 0 and G T = 0, so Ker B T Ker F \ Ker G T Ker PF \ Ker G T : Con v ersely assume G T = 0 and PF = 0. Then, since G T = Z T B T the rst equation implies B T ? Im Z = Ker K: Th us B T 2 Im K F rom the denition of P G T = 0 is equiv alen t to P = F rom PF = 0, w e obtain 0 = T PF = T F = ( B T ) T K + B T : Since K + is positiv e semidenite, this implies that B T 2 Ker K + = Ker K T ogether, B T 2 Ker K \ Im K = f 0 g : Non-uniqueness of the m ultipliers corresponds to redundan t in tersubdomain con tin uit y constrain ts, whic h occur naturally at crosspoin ts of more than t w o subdomains. Dene the space of the Lagrange m ultipliers as the factorspace ~ = = Ker B T : Since for an y 2 Ker B T P F and G T are the zero v ectors, the operators P F and G T induce operators on ~ whic h will be denoted b y the same sym bols. T o a v oid confusion, all n ull spaces will refer to the space not the factorspace ~ F or example, the n ull space of the induced operator G T on ~ will be denoted b y Ker G T = Ker B T It is easy to see that F is symmetric and positiv e semi-denite and hence so is PFP The next lemma sho ws that the operator PF restricted to Ker G T = Ker B T whic h coincides with the restriction of PFP is positiv e denite. 47

PAGE 57

Lemma 14 The operator PF is symmetric and positiv e denite on the factorspace Ker G T = Ker B T Pr o of. F or u; v 2 Ker G T w e ha v e from the denition of P h PFPu; v i = h PFu; v i = h Fu; v i ; whic h pro v es that PF is symmetric positiv e semidenite on Ker G T T o pro v e that PF is nonsingular on Ker G T = Ker B T let PFu = 0 and G T u = 0. Since F = BKB T and G T = Z T B T b y denition, it follo ws that B T u 2 Ker K and B T u ? Ker K hence B T u = 0. The original FETI method is therefore the method of preconditioned conjugate gradien ts in the factorspace ~ for the operator equation, equiv alen t to (3.9), PF = Pd; 2 0 + Ker G T ; (3.16) where 0 is an initial appro ximation to the conjugate gradien t method c hosen so that G T 0 = e and all searc h directions are in the space Ker G T 3.2 Generalized FETI There ha v e been sev eral extensions to the original method. Extension to time-dependen t problems w as done in [26]. Generalization to plate-bending problems is discussed in [51, 53] and, in detail, in this thesis. F or the original method for plate bending problems, the condition n um ber w as observ ed to gro w fast with the n um ber of elemen ts per subdomain [32]. This is caused b y the fact that plate bending is a fourth order problem, while the FETI domain decomposition method \tears" the appro ximate solution at subdomain crosspoin ts, whic h is suitable only for second order problems. 48

PAGE 58

The limitation of method can be cured to extend the FETI methodology to obtain a non-o v erlapping domain decomposition method for plate bending problems. This new method has the properties one usually looks for in iterativ e substructuring methods: the condition n um ber can be bounded independen tly of the n um ber of subdomains, and it gro ws only poly-logarithmicall y with the n um ber of elemen ts per subdomain. The computational cost per iteration is proportional to the solution of a boundary v alue problem in eac h subdomain, plus the solution of a sparse c o arse pr oblem with only few v ariables per subdomain. Suc h methods are commonly referred to as sc alable and quasi-optimal though, of course, for v ery large n um ber of subdomains, the solution of the coarse problem w ould dominate. The k ey idea of our method is to enforce the con tin uit y of the appro ximate solution at the subdomain crosspoin ts throughout the iterations b y adding the corresponding Lagrange m ultipliers to the coarse problem. A similar idea w as emplo y ed in the Balancing Domain Decomposition (BDD) method for plates [42], where appro ximate con tin uit y at crosspoin ts is enforced b y adding new basis functions to the original coarse space [46 49 ] in order to k eep the energy of the appro ximate solution minimal with respect to displacemen ts that are solutions for poin t loads at the subdomain crosspoin ts (cf. Section 2.4.4). The distinguishing features of both the presen t method and the method from [42] is that they are non-overlapping and w ork for standard nite elemen ts used in ev eryda y engineering practice. F or other domain decomposition methods for the biharmonic equation and plate bending see, for example, [12, 71]. Extensions to shells, implementation issues, and further computational can be found in [29, 31 ]. 49

PAGE 59

3.2.1 Dualit y Deriv ation of the Generalized FETI In this section, w e presen t a deriv ation and form ulation [52 ] of the generalized FETI, based on the concept of coarse optimalit y of a dual objectiv e function. Preserving the condition G T k = e throughout the iterations of the Algorithm 11 can be in terpreted as enforcing that ev ery k be optimal with respect to all possible incremen ts of the form G : G T k = e ( ) C ( k ) C ( k G ) ; 8 ; (3.17) where C is dened b y (3.11). The k ey to the generalization of the FETI method to plate bending problems is to mak e all k optimal in more directions. Let C be some giv en matrix with the same n um ber of ro ws as G Eac h column of C will giv e rise to an additional v ariable in the coarse problem. W e shall satisfy in eac h iteration the c o arse optimality pr op erty C ( ) C ( G C ) ; 8 ; : (3.18) with = k T o satisfy this propert y consider an auxiliary problem: F or a giv en ~ nd and so that C ( ) max ; = ~ G C (3.19) Since w e only need solutions satisfying C ( ) > 1 w e consider the maximization problem (3.19) along with the constrain t G T = e (3.20) 50

PAGE 60

In troducing new Lagrange m ultipliers for (3.20), w e get that , and solv e the saddle poin t problem inf ; sup ~ L ( ~ G C; ) = sup inf ; ~ L ( ~ G C; ) where ~ L ( ; ) = 1 2 T F + T d + T ( G T e ) F rom the optimalit y conditions @ ~ L ( ; ) @ = 0 ; @ ~ L ( ; ) @ = 0 ; @ ~ L ( ; ) @ = 0 ; = ~ G C; w e obtain that (3.19) is equiv alen t to the bloc k linear system M 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 ~ 2 6 6 6 6 6 6 4 G T d C T d e 3 7 7 7 7 7 7 5 (3.21) where M = 2 6 6 6 6 6 6 4 G T FG G T FC G T G C T FG C T FC C T G G T G G T C 0 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 G T 0 C T 0 0 G T 3 7 7 7 7 7 7 5 2 6 6 4 F I I 0 3 7 7 5 2 6 6 4 G C 0 0 0 G 3 7 7 5 : (3.22) The solution of (3.19) is unique up to the addition of a v ector in Ker F \ Ker G T and w e write it as = ~ G C 0 M + 0 B B B B B B @ 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 ~ 2 6 6 6 6 6 6 4 G T d C T d e 3 7 7 7 7 7 7 5 1 C C C C C C A : (3.23) 51

PAGE 61

Note that since (3.19) has a solution for an y ~ so does (3.21), hence Im 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 Im M: (3.24) F urther, if is coarse optimal, then + is also coarse optimal if and only if G C 0 M + 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 2 Ker F \ Ker G T : (3.25) The generalized FETI algorithm is th us obtained from Algorithm 11 b y prescribing the initial appro ximation 0 b y 0 = 0 G C 0 M + 0 B B B B B B @ 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 0 2 6 6 6 6 6 6 4 G T d C T d e 3 7 7 7 7 7 7 5 1 C C C C C C A and b y replacing the step y k 1 = Pz k 1 b y y k 1 = z k 1 [ G; C; 0] M + 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 z k 1 ; cf. Algorithm 20 in the next section. 3.2.2 Deriv ation of Generalized FETI Algorithm without Dualit y The generalized FETI method can be also obtained b y forcing the iterates to satisfy also a w eigh ted residual condition [53]. That is, w e require throughout the iterations that C T P ( F d ) = 0 ; G T = e; (3.26) 52

PAGE 62

where C is some giv en matrix. Searc h directions that preserv e (3.26) form the space V 0 = f 2 j G T = 0 ; C T PF = 0 g : (3.27) The corresponding space of residuals is V = f u 2 Im B j G T u = 0 ; C T u = 0 g : (3.28) Denote the associated subspaces of the factorspace ~ as ~ V 0 = f u + Ker B T j u 2 V 0 g = V 0 = Ker B T ; ~ V = f u + Ker B T j u 2 V g : Note that V and ~ V are isomorphic: since V Im B and Im B \ Ker B T = f 0 g eac h class of ~ V con tains exactly one elemen t of V W e will need sev eral properties of the spaces ~ V and ~ V 0 Lemma 15 ~ V = PF ~ V 0 : Pr o of. F rom the denition, clearly PFV 0 V hence PF ~ V 0 ~ V T o sho w that ~ V PF ~ V 0 let u + Ker B T 2 ~ V Since PF is a bijection on Ker G T = Ker B T b y Lemma 14, there is ~ 2 Ker G T = Ker B T suc h that PF ~ = u + Ker B T It follo ws that ~ 2 V 0 Lemma 16 The space ~ V 0 is the dual of ~ V with the dualit y pairing h ; i Pr o of. An y 2 ~ V 0 denes a linear functional 0 on ~ V b y 0 ( v ) = h ; v i Let 0 be an arbitrary linear functional on ~ V F rom Lemma 15, it follo ws that 0 PF is a linear functional on ~ V 0 and, from Riesz represen tation theorem, there is a unique 2 ~ V 0 suc h that 0 ( PFv ) = h ; v i for all v 2 ~ V 0 F rom Lemma 14, PF is a bijection on ~ V and it follo ws that the mapping bet w een a linear functional 0 and its represen tation 2 ~ V 0 is an isomorphism. 53

PAGE 63

Lemma 17 ~ = ~ V 0 ~ V ? : Pr o of. F rom Lemma 14, PF is symmetric positiv e denite on Ker G T = Ker B T Hence, ( PF ) 1 is also symmetric, positiv e denite on Ker G T = Ker B T and, using Lemma 15, it follo ws that ~ = ~ V 0 ( ~ V 0 ) ? ( PF ) 1 = ~ V 0 (( PF ) ~ V ) ? ( PF ) 1 = ~ V 0 ~ V ? ; (3.29) whic h w as to be pro v ed. W e no w dene a projection operator Q b y Q : ~ ~ ; Q 2 = Q; Im Q = ~ V 0 ; Ker Q = ~ V ? ; (3.30) and compute the matrix represen tation of Q This represen tation denes also an operator on whic h will be denoted b y the same sym bol Q Lemma 18 The projector Q is giv en b y the form ula Q = I G C 0 M + 2 6 6 6 6 6 6 4 G T F C T F G T 3 7 7 7 7 7 7 5 ; where M is dened b y (3.22). Pr o of. Let ; 2 Q ( + Ker B T ) = + Ker B T Since 2 V 0 from the denition of P there exists suc h that PF = F + G F rom the denition of V 0 G T = 0 and C T PF = 0, hence G T F + G T G = 0 ; C T F + C T G = 0 ; G T = 0 : (3.31) 54

PAGE 64

F rom the denition of Q w e ha v e ? V so for an y u 2 Im B G T u = 0 and C T u = 0 implies h ; u i = 0. Consequen tly there exist and suc h that for all u 2 Im B h ; u i = h G; u i + h C; u i ; whic h implies that = + G + C + r; r 2 (Im B ) ? = Ker B T (3.32) Since 2 V 0 substituting (3.32) in to (3.31) giv es that , satisfy the linear system G T F ( + G + C ) + G T G = 0 C T F ( + G + C ) + C T G = 0 G T ( + G + C ) = 0 : (3.33) On the other hand, it is easy to see from (3.31) and the denitions of Q and V that if , satisfy (3.33), then = + G + C 2 Q + Ker B T Similarly as in the case of the original FETI, w e need to nd an initial appro ximation 0 satisfying certain conditions, in this case (3.26). Lemma 19 F or an y 0 2 the system of equations G T F ( 0 + G + C ) + G T G = G T d C T F ( 0 + G + C ) + C T G = C T d G T ( 0 + G + C ) = e (3.34) has a solution ; ; and 0 = 0 + G + C (3.35) satises C T P ( F 0 d ) = 0 ; G T 0 = e: (3.36) 55

PAGE 65

Pr o of. As in the preceding proof, 0 satises (3.36) if and only if there is a suc h that G T F 0 + G T G = G T d C T F 0 + C T G = C T d G T 0 = e (3.37) Substituting 0 from (3.35) in to (3.37) yields the system (3.34), whic h can be written as M 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 = X T 2 6 6 4 d G ( G T G ) 1 e 3 7 7 5 ; X = 2 6 6 4 G C 0 0 0 G 3 7 7 5 : F rom the factorization (3.22), it follo ws that Ker M = Ker X Using symmetry of M w e ha v e Im M = (Ker M ) ? = (Ker X ) ? = Im X T ; hence, (3.34) has a solution. The generalized FETI method is the conjugate gradien ts method for the operator PF : V 0 V preconditioned b y QDQ T : V V 0 where the D : is a giv en operator symmetric on V Since Q T is also a projection and Im Q T = (Ker Q ) ? = V the application of Q T on V can be omitted, and one obtains the follo wing algorithm similar to the algorithm of the original FETI (Algorithm 11). Algorithm 20 (Generalized FETI) Giv en an initial 0 compute the initial 0 from (3.34) and (3.35), and compute the initial residual r 0 = P ( F 0 d ) : Repeat for k = 1 ; 2 ; : : : un til con v ergence: z k 1 = Dr k 1 56

PAGE 66

y k 1 = Qz k 1 k = r T k 1 y k 1 p k = y k 1 + k k 1 p k 1 ( p 1 = y 0 ) k = k p T k PFp k k = k 1 + k p k r k = r k 1 k PFp k The c hoice of the preconditioner is discussed in Section 3.3. 3.2.3 Method Selection for Plate Problems and Other Generalizations W e c hoose the columns of matrix C whic h appears in the description of the generalized FETI method abo v e, as v ectors with a one at the position of a Lagrange m ultiplier that enforces the con tin uit y of the transv ersal displacemen t at a crosspoin t, and zeros ev erywhere else. By a crosspoin t w e understand an in terface node adjacen t to at least three subdomains or to t w o subdomains and the complemen t of n. F or a more precise form ulation, see Section 4.5.1. A similar idea can be exploited for shell problems. The con tin uit y of the normal displacemen t needs to be enforced. T o a v oid nding normals and the added complexit y of enforcing con tin uit y of the normal displacemen t, one ma y enforce con tin uit y of all displacemen t degrees of freedom. This, ho w ev er, increases the size of the coarse space. Other possibilities and computational issues are discussed in [30]. F arhat et al. in [26 ] studies the case of time-dependen t problems where the subdomain stiness matrices K i are perturbed b y the addition of 57

PAGE 67

a m ultiple of the subdomain mass matrix, th us making the new local matrix positiv e denite. Consequen tly all matrices Z i are v oid and the natural coarse problem is lost in time-dependen t applications. The methodology dev eloped in [26] for rein troducing a coarsening operator in the FETI algorithm for dynamics problems is a special case of the presen t generalization where C is tak en to be the matrix G b efor e the p erturb ation that is, C = [ B i ~ Z i ] where the columns of ~ Z i are the basis of the k ernel of the local stiness matrix of the subdomain n i The reason wh y the preconditioner w orks for the dynamics problems is quite dieren t from that for the plate bending. 3.3 The Diric hlet Preconditioner Let us decompose the space of all degrees of freedom W in to the space of in ternal degrees of freedom and the degrees of freedom on subdomain in terfaces, W = W W: In the corresponding bloc k notation, B = [ 0 B ] ; B : W ; since B has nonzero en tries for in terface degrees of freedom only Also, Z = 2 6 6 4 Z Z 3 7 7 5 ; and w e ha v e G = BZ = B Z; Ker B T = Ker B T : 58

PAGE 68

Let S be the Sc h ur complemen t of K obtained b y elimination of degrees of freedom in ternal to all subdomains. Then F = BK + B T = BS + B T (3.38) and Ker S = Im Z Ev aluation of the matrix-v ector product S + u reduces to the solution of independen t Neumann pr oblems on all subdomains, cf. Section 2.4.1. Inspired b y (3.38), w e c hoose D = BS B T giving the preconditioner QD = Q BS B T : (3.39) This preconditioner is called the Dirichlet pr e c onditioner since ev aluating the matrix-v ector product Sr is equiv alen t to solving independen t Dirichlet pr oblems on all subdomains, cf. Section 2.4.1. 3.4 In terface F orm ulation of P ark et al. In [59, 58 ], K.C. P ark has dev eloped a similar substructuring method. W e will sho w that this method can essen tially be written as FETI with a special c hoice of the in terface con tin uit y operator B The method augmen ts the Lagrangian b y in troducing another v ariable u g the global v ector of degrees of freedom, on the in terfaces bet w een subdomains. This v ariable is redundan t and is later eliminated, yielding an analog of FETI. W e will sho w that the method can be written as FETI and discuss the resulting c hoice of B and its implications. Let L be the subdomain assem bly matrix, that is a zero-one matrix mapping the local subdomain degrees of freedom to global ones. Using this matrix, the con tin uit y of the solution is enforced through the constrain t u = Lu g : 59

PAGE 69

The Lagrangian can be written as L ( u; u g ; ) = 1 2 u T Ku f T u + T ( u Lu g ) : Decomposing u = Z + in to a rigid body mode componen t and its complemen t (pure deformation componen t), w e obtain L ( ; ; u g ; ) = 1 2 T K f T f T Z + T ( Z + Lu g ) : The stationarit y condition leads to the system of equations 2 6 6 6 6 6 6 6 6 6 6 4 K I 0 0 I 0 Z L 0 Z T 0 0 0 L 0 0 3 7 7 7 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 6 6 6 4 u g 3 7 7 7 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 6 6 6 4 f 0 Zf 0 3 7 7 7 7 7 7 7 7 7 7 5 Eliminating the deformation componen t, = K + ( f ), the system becomes 2 6 6 6 6 6 6 4 K + Z L Z T 0 0 L T 0 0 2 6 6 6 6 6 6 4 2 6 6 6 6 6 6 4 u g 3 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 4 K + f Zf 0 3 7 7 7 7 7 7 5 The rst equation of the system is K + + Z Lu g = K + f: (3.40) W e rst eliminate the term Z b y m ultiplying (3.40) b y the orthogonal projection P Z on to (Im Z ) ? = Ker Z T P Z = I Z ( Z T Z ) 1 Z T and obtain the equation P Z K + P Z Lu g = P Z K + f; (3.41) 60

PAGE 70

equiv alen t to (3.40) holding for some Next, using the orthogonal projection P l on to (Im P Z L ) ? = Ker( P Z L ) T P l = I P Z L ( L T P Z L ) 1 L T P Z ; in a similar w a y P l P Z K + = P l P Z K + f: This is the equation (15) in [58 ]. W e ma y rev erse the order of elimination of the terms, and w e obtain from (3.40) P z P L K + = P z P L K + f; (3.42) where the orthogonal projections P L and P z are giv en b y P L = I L ( L T L ) 1 L T P z = I P L Z ( Z T P L Z ) 1 Z T P L : Since Im P z Im P L = Ker L T w e ha v e = P L for 2 Im( P z ). Th us, P L K + = P L K + P L : No w c hoosing B = P L ; (3.43) w e ma y write P z = I P L Z ( Z T P L Z ) 1 Z T P L = I G ( G T G ) 1 G T = P; where w e ha v e used the notation (3.4). This sho ws that the method with the equation (3.42) can be written as FETI with the special c hoice of B (3.43). 61

PAGE 71

Let us illustrate this c hoice. If a mesh node is on the in terface of t w o subdomains, then, omitting all other degrees of freedom and Lagrange m ultipliers, the corresponding bloc k of the assem bly matrix L is 2 6 6 4 1 1 3 7 7 5 and the corresponding bloc k of the projection P L is 1 2 2 6 6 4 1 1 1 1 3 7 7 5 : Similarly for a three node and four node in terfaces, w e get 1 3 2 6 6 6 6 6 6 4 2 1 1 1 2 1 1 1 2 3 7 7 7 7 7 7 5 ; and 1 4 2 6 6 6 6 6 6 6 6 6 6 4 3 1 1 1 1 3 1 1 1 1 3 1 1 1 1 3 3 7 7 7 7 7 7 7 7 7 7 5 ; respectiv ely The c hoice of B = P L has some in teresting implications on the analysis, whic h will be discussed in Section 4.4. 62

PAGE 72

4. Con v ergence Analysis In this c hapter, w e pro v e that for a second order problem the condition n um ber of the preconditioned FETI method is bounded independen tly of the n um ber of subdomains and poly-logarithmical ly in terms of subdomain size, similarly as it is in the case for other optimal non-o v erlapping domain decomposition methods [11 23 25 49 48 ]. W e pro v e that the c hoice of the in terface con tin uit y operator as giv en b y P ark et al. assures that the condition n um ber is bounded b y C (1+log H h ) 2 Finally w e also sho w that the generalized FETI con v erges poly-logarithmicall y for the plate bending problem problem. Analysis of domain decomposition methods t ypically demonstrates spectral equiv alence of the quadratic form that denes the problem in a v ariational setting and the quadratic form that denes the preconditioner, often b y w a y of P .L. Lions lemma [6, 23, 24 44 ]. Since the preconditioner in the FETI method is quite complicated and is not dened in terms of a quadratic form, w e proceed dieren tly W e nd a bound on the norm of the product of the system operator and the preconditioner, so as to bound the maximal eigenv alue, and a bound on the in v erse, to bound the minimal eigen v alue. Related analyses w ere previously done for methods without crosspoin ts bet w een the subdomains, or done formally in functional spaces, cf., for example, Glo winski and Wheeler [37 ] and Bakh v alo v and Kn y azev [4]. W e presen t a complete analysis in terms of upper and lo w er bound on the preconditioned operator for decompositions with crosspoin ts in 2D and edges and crosspoin ts in 3D for 63

PAGE 73

second order elliptic problems. W e sho w that the condition n um ber is bounded b y C (1 + log H h ) r where r = 2 ; 3, h is the c haracteristic mesh size, and H the diameter of the subdomains. 4.1 Preliminaries In this section, w e presen t some results that will be used in analysis of FETI methods. The results are mostly concerned with estimates for nite elemen t functions used to discretize our model problems. W e will assume that a domain n is divided in to a set of nono v erlapping subdomains n i ; i = 1 ; : : :; N s n = n 1 [ : : : [ n N s The subdomains are assumed to be shape regular of diameter H according to the denition belo w. W e will form ulate all results for one of the subdomains n i and assume that constan ts in the estimates do not depend on the index of the subdomain. Denition 21 A subdomain n i I R d is said to be shape-regular of diameter O ( H ) if it can be generated from a reference domain (square or cube) ^ n of unit diameter b y a mapping F i suc h that n i = F i ( ^ n). The mapping is assumed to satisfy k @ F i k c s H; k @ F 1 i k c s H 1 (4.1) where @ F i is the Jacobian of the mapping, k : k is the Euclidean I R d matrix norm, and c s is a positiv e constan t. W e will be using the Sobolev spaces W k;p F or p = 2, the spaces are Hilbert spaces and w e denote them H k W k; 2 Denitions of Sobolev spaces can be found, for example, in [1]. F or k 6 = 1 = 2 ; 1, w e use k k k;p; n and j j k;p; n to denote the standard Sobolev norms and seminorms of functions in W k;p F ollo wing [22 ], w e dene the scaled Sobolev norm for a scalar function 64

PAGE 74

u 2 H 1 (n i ) k u k 2 1 ; 2 ; n i = j u j 2 1 ; 2 ; n i + 1 H 2 k u k 2 0 ; 2 ; n i ; where the Sobolev seminorm is dened b y j u j 2 1 ; 2 ; n i = Z n i kr u ( x ) k 2 dx: W e dene the scaled Sobolev norm for a scalar function u 2 H 1 = 2 ( @ n i ) k u k 2 1 = 2 ; 2 ;@ n i = j u j 2 1 = 2 ; 2 ;@ n i + 1 H k u k 2 0 ; 2 ;@ n i ; where the Sobolev seminorm is dened b y j u j 2 1 = 2 ; 2 ;@ n i = Z @ n i Z @ n i j u ( x ) u ( y ) j 2 k x y k d dxdy: Here k k is the Euclidean norm in I R d W e note that the space H 1 = 2 ( @ n i ) is the space of traces of functions in the space H 1 (n i ). If u = ( u 1 ; u 2 ) is a v ector function, then w e dene, for example, k u k 2 1 ; 2 ; n i = k u 1 k 2 1 ; 2 ; n i + k u 2 k 2 1 ; 2 ; n i : Other norms and seminorms of v ector functions are dened analogously W e presen t t w o v arian ts of the P oincar e-F riedric hs inequalit y W e pro v e the rst for the sak e of completeness and to demonstrate the tec hnique that is used to pro v e its v arian ts for nite elemen t functions. W e note that P k ; k = 0 ; 1 ; : : : denotes the space of polynomials of degree at most k Lemma 22 Let G be a con tin uous linear functional on H 1 = 2 ( @ n i ) suc h that for all u 2 P 0 G ( u ) = 0 implies u = 0. Then there exists a constan t c independen t of H suc h that for all u 2 H 1 = 2 ( @ n i ) ; u 2 Ker G k u k 1 2 ; 2 ;@ n i c j u j 1 2 ; 2 ;@ n i : 65

PAGE 75

Pr o of. W e pro v e the theorem for a reference domain of diameter 1 b y con tradiction. The result in the scaled norms then follo ws. If the inequalit y is not true, there exists a sequence f u n g Ker G suc h that k u n k 1 2 ; 2 ;@ ^ n = 1 and j u n j 1 2 ; 2 ;@ ^ n 0 : Due to the compact im bedding of H 1 = 2 ( @ ^ n) in L 2 ( @ ^ n), there exists a subsequence of f u n g that con v erges in L 2 ( @ ^ n). Since j u n j 1 2 ; 2 ;@ ^ n 0, the subsequence f u n k g is Cauc h y in H 1 = 2 ( @ ^ n). Therefore it con v erges in H 1 = 2 ( @ ^ n) to some u 2 H 1 = 2 ( @ ^ n). The con tin uit y of the norm and the seminorm on H 1 = 2 ( @ ^ n) implies k u k 1 2 ; 2 ;@ ^ n = 1 and j u j 1 2 ; 2 ;@ ^ n = 0 : (4.2) Therefore, b y the second equation in (4.2), u = k almost ev erywhere, where k is some real n um ber. Let us assume without loss of generalit y that u = k ev erywhere; that is u 2 P 0 Since G is con tin uous G ( u n k ) G ( u ) = 0. This yields, b y the assumption, u = 0 whic h con tradicts the rst equation in (4.2). The proof of this P oincar e-F riedric hs inequalit y in scaled norms is similar to the proof of the previous lemma. Lemma 23 Let G be a con tin uous linear functional on H 1 (n i ) suc h that for all u 2 P 0 G ( u ) = 0 implies u = 0. Then there exists a constan t c independen t of H suc h that for all u 2 H 1 (n i ) ; u 2 Ker G k u k 1 ; 2 ; n i c j u j 1 ; 2 ; n i : 66

PAGE 76

4.1.1 Estimates for P1/Q1 Elemen ts Let V P 1 h (n i ) be a conforming nite elemen t space of P1 or Q1 elemen ts [14] satisfying the usual regularit y and in v erse properties and possibly some essen tial boundary conditions. That is, for example in I R 2 w e assume that n i = [ K2T h;i K ; where eac h elemen t K of the triangulation T h;i is a triangle or a rectangle. F urthermore, for all K 2 T h;i c 1 d ( K ) h c 2 ( K ) ; (4.3) where d ( K ) is the diameter of K and ( K ) the diameter of the circle inscribed in K Then h is called the c haracteristic mesh size. A v ertex of an elemen t K 2 T h;i will be referred to as a mesh poin t, a nodal poin t, or just a v ertex. If essen tial boundary conditions are prescribed on @ n i w e assume that () c c ( @ n i ), where ( ) denotes the measure. W e note that functions in V P 1 h (n i ) are con tin uous. As in the previous c hapter, the corresponding space of v ectors of degrees of freedom is denoted W i Let I P 1 : W i V P 1 h (n i ) denote the linear one-to-one transformation that maps a v ector of degrees of freedom to the corresponding nite elemen t function. This transformation is often called nite elemen t in terpolation. W e decompose a v ector of degrees of freedom U 2 W i in to in ternal and boundary degrees of freedom assuming the boundary degrees of freedom are listed last. That is, in bloc k notation, U = [ U T ; U T ] T where U is the v ector of boundary degrees of freedom. W e denote W i the space of boundary degrees 67

PAGE 77

of freedom and dene the matrix T i so that U = T i U Then, I P 1 U denes a function on @ n i whic h depends on U only By abuse of notation, w e will write I P 1 U = I P 1 [0 T ; U T ] T on @ n i W e summarize some w ell kno wn results and inequalities in a form suitable for our purposes. The next lemma summarizes the fact that the H 1 = 2 norm of a zero extension of a "piece" of a function can be bounded b y the norm of the function. Lemma 24 Let @ n i be a v ertex, edge, or face (if d = 3) of subdomain n i A face is understood not to con tain adjacen t edges, and an edge does not con tains its endpoin ts. F or for all z 2 V P 1 h ( @ n i ), dene w 2 V P 1 h ( @ n i ) b y w ( x ) = z ( x ) on all nodes of triangulation x 2 w ( x ) = 0 on all other nodes of @ n i Then k w k 2 1 2 ; 2 ;@ n i C (1 + log H h ) ( j z j 2 1 2 ; 2 ;@ n i + 1 H k z k 2 0 ; 2 ;@ n i ) ; where = 1 if d = 2 and is a v ertex, or d = 3 and is an edge or a v ertex; = 2 if d = 2 and is an edge, or d = 3 and is a face. Pr o of. The inequalit y for d = 2 w as pro v ed in [49, 48]. The case when d = 3 follo ws from Lemmas 4.1 and 4.2 in [11] if is an edge or a v ertex, and Lemma 4.3 in [11] if is a face. Cf. also [23]. The follo wing lemma can be pro v ed b y using Lemmas 4.1 and 4.2 in [11 ] and estimates the H 1 = 2 seminorm of a spik e function. Lemma 25 There exists a constan t c independen t of h and H suc h that for all u 2 V P 1 h (n i ), u ( x ) = 0 for all mesh poin ts x 2 @ n i x 6 = x 0 j u j 2 1 2 ; 2 ;@ n i c k u k 2 0 ; 1 ;@ n i 68

PAGE 78

The follo wing lemma sho ws the equiv alence of the discrete seminorm dened b y the Sc h ur complemen t and the H 1 = 2 seminorm. In the next section w e will apply this result to the stiness matrix K i whic h is wh y w e the same sym bol here. This result is standard [10 69 ] and it is pro v ed here for the sak e of completeness. Lemma 26 Let K i be a symmetric matrix that satises c j I P 1 U j 2 1 ; 2 ; n i h K i U; U i C j I P 1 U j 2 1 ; 2 ; n i 8 U 2 W i : Let U = T i U and the matrix K i be decomposed so that K i U = 2 6 6 4 K ~ K T ~ K K 3 7 7 5 2 6 6 4 U U 3 7 7 5 : Let S i be the Sc h ur complemen t of K in K i Then, there exist constan ts c 0 and C 0 independen t of h and H suc h that c 0 j I P 1 U j 2 1 2 ; 2 ;@ n i h S i U; U i C 0 j I P 1 U j 2 1 2 ; 2 ;@ n i 8 U 2 W i : Pr o of. F rom the denition of the Sc h ur complemen t (Denition 6) and the propert y (2.8), it follo ws that h S i U; U i = h K U; U i h K 1 ~ K T U; ~ K T U i = inf U =[ U T ; U T ] T h K i U; U i Hence, from the assumption b y in v ariance of the j j 1 ; 2 ; n i seminorm to adding a constan t and b y the discrete extension theorem [22] for the scaled norms, w e obtain h S i U; U i C inf U =[ U T ; U T ] T j I P 1 U + k j 1 ; 2 ; n i Cc e k I P 1 U + k k 1 2 ; 2 ;@ n i ; 69

PAGE 79

where k 2 I R. Choosing k = R @ n i I P 1 Udx R @ n i dx yields the equiv alence of the scaled norm to the seminorm b y the P oincar e-F riedric hs inequalit y (Lemma 22). On the other hand, b y the trace theorem, in v ariance of the seminorms to adding a constan t, and the P oincar e-F riedric hs inequalit y (Lemma 23), w e ha v e 1 c t c p c j I P 1 U j 2 1 2 ; 2 ;@ n i c inf U =[ U T ; U T ] T j I P 1 U j 1 ; 2 ; n i h S i U; U i ; where c t and c p are constan ts arising from the trace and P oincar e-F riedric hs inequalit y respectiv ely The follo wing theorem sho ws a v arian t of the P oincar e-F riedric hs inequalit y W e note that the dimension of W i implicitly depends on h and for a uniform mesh, the lemma follo ws from Lemma 22. Lemma 27 Let S i be the Sc h ur complemen t from Lemma 26. There exists a constan t C independen t of H and h suc h that for all U 2 W i U ? Ker S i k I P 1 U k 1 2 ; 2 ;@ n i C j I P 1 U j 1 2 ; 2 ;@ n i : Pr o of. If there is an essen tial boundary condition imposed on @ n i then Ker S i = f 0 g and the statemen t follo ws from the P oincar e-F riedric hs inequalit y Lemma 22. Otherwise, Ker S i is spanned b y the v ector of ones. Then the condition U ? Ker S i implies P N j =1 U j = 0 for U = [ U 1 ; : : :; U N ] T F ollo wing along the lines of the proof of Lemma 22, w e obtain a sequence of v ectors U n 2 W i I R N ( n ) suc h that ( I P 1 U n ) F i u in L 2 ( @ ^ n) ; and k u k 1 2 ; 2 ;@ ^ n = 1 and j u j 1 2 ; 2 ;@ ^ n = 0 (4.4) 70

PAGE 80

Therefore, u = k almost ev erywhere, where k 2 I R. Th us, ( I P 1 U n ) F i k in L 2 ( @ ^ n). Then, from the quasi-uniformit y of the triangulation, w e obtain k ( I P 1 U n ) F i k k 2 0 ; 2 ;@ ^ n = c 1 N ( n ) k U n k k 2 c N ( n ) 2 0 @ N ( n ) X j =1 ( U n j k ) 1 A 2 Since U n ? Ker S i it holds that P N ( n ) j =1 U n j = 0. Th us, ck 2 k ( I P 1 U n ) F i k k 2 0 ; 2 ;@ ^ n 0 : This is possible only if k = 0 whic h is a con tradiction with (4.4). 4.1.2 Estimates for HCT Elemen ts Let V HCT h (n i ) be the nite elemen t space of HCT elemen ts satisfying the usual regularit y and in v erse properties and possibly some essen tial boundary conditions. That is, w e assume that n i = [ K 2T h;i K ; where eac h elemen t K of the triangulation T h;i is a triangle. F urthermore, condition (4.3) is satised and h denotes the c haracteristic mesh size. W e note that functions in V HCT h (n i ) are C 1 con tin uous and the space V HCT h (n i ) can be written as V HCT h (n i ) = f v 2 C 1 (n i ) : 8K 2 T h;i ; v j K j 2 P 3 ( K j ) for all subtriangles K j of K ; @v @n j a k a j 2 P 1 ( a k a j ) for all sides ( a k ; a j ) of the triangle Kg On eac h triangle, a function v in V HCT h (n i ) is determined b y the v alues v ( a i ) and the v alues of its deriv ativ es @v @x j ( a i ) at the v ertices of the triangle. The 71

PAGE 81

in ternal v ertex of subtriangles is not arbitrary but it is determined to guaran tee the unisolv ence of the elemen t. This is as m uc h detail as w e will need in our considerations. W e refer to [14] for more details about the HCT elemen t. The corresponding space of v ectors of degrees of freedom is denoted W i Let I P 1 : W i V HCT h (n i ) denote the linear one-to-one transformation that maps a v ector of degrees of freedom to the corresponding nite elemen t function. As in the previous section, w e decompose a v ector of degrees of freedom U in to in ternal and boundary degrees of freedom assuming the boundary degrees of freedom are listed last. That is, in bloc k notation, U = [ U T ; U T ] T where U is the v ector of boundary degrees of freedom. W e denote W i the space of boundary degrees of freedom and dene the matrix T i so that U = T i U Then, I HCT U denes a function on @ n i whic h depends on U only By abuse of notation, w e will write I HCT U = I HCT [0 T ; U T ] T on @ n i W e summarize here some w ell kno wn results and inequalities in a form suitable for our purposes. The follo wing lemma estimates the norm of a \spik e" function and is pro v ed in [42 ]. Lemma 28 Let x be a v ertex of a subdomain n i F or u 2 V HCT h (n i ) suc h that u ( x ) = 0, dene z 2 V HCT h (n i ) b y z ( x ) = u ( x ) = 0 ; r z ( x ) = r u ( x ), and z ( y ) = 0 ; r z ( y ) = 0 at all other nodes y of @ n i Then kr z k 2 1 2 ; 2 ;@ n i C 1 + log H h kr u k 2 1 2 ; 2 ;@ n i + 1 H kr u k 2 0 ; 2 ;@ n i : The follo wing estimate of the trace norm of the extension b y zero is pro v ed as in [10 Lemma 3.5]. Lemma 29 Let u 2 V HCT h (n i ). Then there exists a constan t C suc h that if 72

PAGE 82

supp u \ @ n i is con tained in a segmen t of @ n i of length then j u j 2 1 2 ; 2 ;@ n i j u j 2 1 2 ; 2 ; + C 1 + log h k u k 2 0 ; 1 ; : The follo wing estimate is a modication of the previous lemma. Lemma 30 Let u 2 V HCT h (n i ) and supp u \ @ n i be con tained in a segmen t of @ n i and x 1 and y 1 be its endpoin ts. Let x 2 and y 2 be the nodal poin ts next to x 1 and y 1 respectiv ely in the segmen t. Let ( x 2 ; y 2 ) denote the segmen t bet w een x 2 and y 2 Let the length of this segmen t be h Assume that the function u satises the condition k u k 0 ; 1 ; ( x 1 ;y 1 ) c k u k 0 ; 1 ; ( x 2 ;y 2 ) Then j u j 2 1 2 ; 2 ;@ n i j u j 2 1 2 ; 2 ; ( x 2 ;y 2 ) + C 1 + log h k u k 2 0 ; 1 ; ( x 2 ;y 2 ) : Pr o of. By the previous lemma, the inequalit y holds with ( x 2 ; y 2 ) replaced b y the segmen t ( x 1 ; y 1 ). By the denition of the H 1 = 2 seminorm, w e ha v e j u j 2 1 2 ; 2 ; ( x 1 ;y 1 ) = j u j 2 1 2 ; 2 ; ( x 2 ;y 2 ) + Z ( x 1 ;x 2) Z ( x 1 ;y 2) j u ( x ) u ( y ) j 2 k x y k 2 dxdy + Z ( y 1 ;y 2) Z ( x 1 ;y 2) j u ( x ) u ( y ) j 2 k x y k 2 dxdy: Using the fact that j u ( x ) u ( y ) j min f c 1 c h k u k 0 ; 1 ; ( x 2 ;y 2 ) k x y k ; 2 c k u k 0 ; 1 ; ( x 2 ;y 2 ) g ; the last t w o terms are bounded b y c 2 k u k 0 ; 1 ; ( x 2 ;y 2 ) W e will also need a straigh tforw ard extension of the discrete Sobolev inequalit y of Dryja [22 ] to piecewise polynomial functions of order p > 1 [57]. Lemma 31 F or ev ery p 1 there exists a constan t C = C ( p ) suc h that for ev ery u con tin uous on @ n i of length suc h that u 2 P p on the side of ev ery elemen t, k u k 2 0 ; 1 ; C 1 + log h j u j 2 1 2 ; 2 ; + 1 k u k 2 0 ; 2 ; : 73

PAGE 83

The follo wing lemma establishes a useful inequalit y bet w een discrete and Sobolev norms. Lemma 32 Let u 2 V HCT h (n i ) and u ( x 0 ) = 0 for some x 0 2 @ n i Let U = I 1 HCT u and U = T i U Then, for all suc h U h k U k 2 c (1 + H 2 ) kr u k 2 0 ; 2 ;@ n i ; where c is independen t of h and H Pr o of. Let E be an elemen t edge on @ n i and x 1 and x 2 be its endpoin ts. Eac h componen t of r u is a polynomial of degree at most t w o on E Since all norms in a nite dimensional space are equiv alen t, it holds that h ( kr u ( x 1 ) k 2 + kr u ( x 2 ) k 2 ) c 1 min f 2 P 2 P 2 ;f ( x 1 )= r u ( x 1 ) ;f ( x 2 )= r u ( x 2 ) k f k 2 0 ; 2 ;E c 1 kr u k 2 0 ; 2 ;E Summing o v er all edges of the boundary w e obtain h k U k 2 c 1 2 kr u k 2 0 ; 2 ;@ n i + h k U 0 k 2 ; where U 0 denotes the displacemen t degrees of freedom of U W e sho w that h k U 0 k 2 < c 2 H 2 kr u k 2 0 ; 2 ;@ n i Since u ( x 0 ) = 0, w e can write for an y mesh poin t x on the boundary u ( x ) = Z ( x 0 ;x ) r u ( y ) ( y ) dy; where ( x 0 ; x ) @ n i is the part of the boundary bet w een x 0 and x and is the v ector tangen tial to the boundary Squaring, using the Cauc h y inequalit y and considering that the length of ( x 0 ; x ) is bounded b y c 2 H w e obtain u ( x ) 2 c 2 H kr u k 2 0 ; 2 ;@ n i : 74

PAGE 84

The statemen t of the lemma then follo ws b y summing o v er all mesh poin ts on the boundary The follo wing lemma sho ws the equiv alence of the discrete seminorm dened b y the Sc h ur complemen t and the H 1 = 2 seminorm of the gradien t of the corresponding nite elemen t function and is an analog of Lemma 26 of the previous section. Lemma 33 Let K i be a matrix that satises c j I HCT U j 2 2 ; 2 ; n i h K i U; U i C j I HCT U j 2 2 ; 2 ; n i 8 U 2 W i : Let U = T i U and let the matrix K i be decomposed so that K i U = 2 6 6 4 K ~ K T ~ K K 3 7 7 5 2 6 6 4 U U 3 7 7 5 : Let S i be the Sc h ur complemen t of K in K i Then, there exist constan ts c 0 and C 0 independen t of h and H suc h that c 0 jr I HCT U j 2 1 2 ; 2 ;@ n i h S i U; U i C 0 jr I HCT U j 2 1 2 ; 2 ;@ n i 8 U 2 W i : Pr o of. F rom the denition of the Sc h ur complemen t it follo ws that h S i U; U i = h K U; U i h K 1 ~ K T U; ~ K T U i = inf U =[ U T ; U T ] T h K i U; U i Hence, b y the discrete extension theorem [41], w e obtain h S i U; U i C inf U =[ U T ; U T ] T j I P 1 U j 2 ; 2 ; n i Cc e jr I P 1 U j 1 2 ; 2 ;@ n i : On the other hand, the trace theorem and in v ariance of the j j 1 2 ; 2 ; n i seminorm to adding a linear function, and the P oincar e-F riedric hs theorem, jr I P 1 U j 2 1 2 ; 2 ;@ n i c t c p inf U =[ U T ; U T ] T j I P 1 U j 2 ; 2 ; n i : 75

PAGE 85

The follo wing theorem sho ws another v arian t of the P oincar e-F riedric hs inequalit y W e again emphasize that the dimension of W i implicitly depends on h and is an analog of Lemma 27 of the previous section. Lemma 34 Let S i be the Sc h ur complemen t from Lemma 33. There exists a constan t C independen t of H and h suc h that for all U 2 W i U ? Ker S i kr I HCT U k 1 2 ; 2 ;@ n i C jr I HCT U j 1 2 ; 2 ;@ n i : Pr o of. Let us pro v e the case when Ker S is non trivial. In the other case, the statemen t follo ws from the standard P oincar e-F riedric hs inequalit y F ollo wing along the lines of the proof of Lemma 22, w e obtain a sequence of v ectors U n 2 W i I R N ( n ) suc h that ( r I HCT U n ) F i u in L 2 ( @ ^ n) L 2 ( @ ^ n) ; and k u k 1 2 ; 2 ;@ ^ n = 1 and j u j 1 2 ; 2 ;@ ^ n = 0 (4.5) Therefore, u = ( k 1 ; k 2 ) 2 I R I R almost ev erywhere. Th us, ( r I HCT U n ) F i ( k 1 ; k 2 ) in L 2 ( @ ^ n) L 2 ( @ ^ n). W e decompose the v ector U n in to t w o v ectors corresponding to the partial deriv ativ e with respect to x U n 1 and the partial deriv ativ e with respect to y U n 2 Then, b y the proof of Lemma 32, w e obtain k ( r I HCT U n ) F i ( k 1 ; k 2 ) k 2 0 ; 2 ;@ ^ n c 1 N ( n ) ( k U n 1 k 1 k 2 + k U n 2 k 2 k 2 ) c N ( n ) 2 0 @ N ( n ) X i =1 ( U n 1 ;i k 1 ) 1 A 2 + c N ( n ) 2 0 @ N ( n ) X i =1 ( U n 2 ;i k 2 ) 1 A 2 : Since U n ? Ker S i it holds that P N ( n ) i =1 U n 1 ;i = 0 and P N ( n ) i =1 U n 2 ;i = 0. Th us, c ( k 2 1 + k 2 2 ) k ( r I HCT U n ) F i ( k 1 ; k 2 ) k 2 0 ; 2 ;@ ^ n 0 : 76

PAGE 86

This is possible only if ( k 1 ; k 2 ) = (0 ; 0) whic h is a con tradiction with (4.5). 4.2 Abstract Analysis F ramew ork This section dev elops a framew ork for estimating the condition n umber of our method. W e estimate the minim um and maxim um eigen v alues of the matrix of the FETI form ulation and the preconditioner. The follo wing simple lemma will be the basis of our estimates. It will allo w to reduce estimates of norms to estimates of boundedness and coercivit y The proof follo ws a standard argumen t and it is presen ted for completeness only Lemma 35 Let X be a Banac h space, X 0 the dual of X and A : X X 0 a linear operator that is on to and satises the conditions h y; Ax i = h x; Ay i ; 8 x; y 2 X (4.6) c A k x k 2 X h x; Ax i C A k x k 2 X ; 8 x 2 X (4.7) with constan ts C A ; c A > 0. Then k A k X X 0 C A ; k A 1 k X 0 X 1 c A : Pr o of. F rom (4.6), k A k X X 0 = sup x 2 X k Ax k X 0 k x k X = sup x; ~ x 2 X h x; A ~ x i k x k X k ~ x k X = sup x 2 X h x; Ax i k x k 2 X C A : F rom (4.7) and the fact that A is on to, w e obtain 1 k A 1 k X 0 X = inf x 2 X k Ax k X 0 k x k X = inf x 2 X sup ~ x 2 X h x; A ~ x i k x k X k ~ x k X inf x 2 X h x; Ax i k x k 2 X c A ; concluding the proof. 77

PAGE 87

Assume that the preconditioner of a FETI method can be written in the form RDR T where R is a projection, and RDR T is an isomorphism from H on to H 0 Let H and its dual H 0 are subspaces of Ker G T F urthermore, let H = Im R T The original FETI satises this with R = P and H = Ker G T = Ker B T The generalized FETI satises the assumption with R = Q and H = ~ V F or the purpose of analysis, w e equip the space H with the norm k v k H = k B T v k S = h S B T v; B T v i 1 = 2 : (4.8) Since B T v ? Ker S for v 2 H (4.8) indeed denes a norm rather than only seminorm. The dual space H 0 is equipped with the dual norm k k H 0 = sup v 2H h ; v i k v k H : (4.9) Since H ~ it immediately follo ws k k H 0 = sup v 2H h ; v i k B T v k S = sup w 2 W; Bw 2H h ; Bw i k B T Bw k S : (4.10) The norm on H w as c hosen so that the preconditioner RD is trivially coerciv e and bounded. Lemma 36 F or all v 2 H h v; RDv i = k v k 2 H : Pr o of. Let v 2 H = Im R T Since R T is a projection, w e ha v e b y denition of the preconditioner D h v; RDv i = h R T v; BS B T v i = h v; BS B T v i = h B T v; S B T v i = k v k 2 H ; whic h w as to be sho wn. 78

PAGE 88

Coercivit y and boundedness of the system operator PF on H 0 will be estimated using the follo wing lemma. Lemma 37 F or all 2 H 0 h ; F i = sup w 2 W;w ? Ker S h ; Bw i 2 k w k 2 S = sup w 2 W h ; Bw i 2 k w k 2 S : Pr o of. Let 2 H 0 Then h ; F i = h S + B T ; B T i = h S 1 = 2 B T ; S 1 = 2 B T i = k S 1 = 2 B T k 2 = sup x 2 W h S 1 = 2 B T ; x i 2 k x k 2 = sup x 2 W; x = x 1 + x 2 x 1 2 Ker S x 2 ? Ker S h B T ; S 1 = 2 x i 2 k x 1 + x 2 k 2 = sup x 2 2 W; x 2 ? Ker S h B T ; S 1 = 2 x 2 i 2 k x 2 k 2 since S 1 = 2 x 1 = 0 and k x k 2 = k x 1 k 2 + k x 2 k 2 No w write an y w 2 W as w = w 1 + w 2 ; w 1 2 Ker S; w 2 = S 1 = 2 x 2 ? Ker S: 2 H 0 implies that h B T ; w 1 i = 0, hence h B T ; w 2 i = h B T ; w i = h ; Bw i It follo ws that h ; F i = sup w 2 2 W; w 2 ? Ker S h B T ; w 2 i 2 h w 2 ; Sw 2 i = sup w 2 W h ; Bw i 2 k w k 2 S ; whic h w as to be pro v ed. It is w ell kno wn [39] that after k iterations of the preconditioned conjugate gradien t method, the energy norm or the error jjj e jjj = h PFe; e i 1 = 2 is reduced b y a factor of at least 2(( p 1) = ( p +1)) k ; where is the condition n um ber. The condition n um ber is giv en in our case b y = ( RDPF ) = max ( RDPF ) min ( RDPF ) ; (4.11) where max and min are the maxim um and minim um eigen v alue, respectiv ely W e are no w ready to pro v e an abstract bound on 79

PAGE 89

Theorem 38 Assume there exist constan ts C 1 C 2 suc h that (i) for an y 2 H 0 and w 2 W suc h that Bw 2 H there is ~ w 2 W suc h that h ; B ~ w i = h ; Bw i ; and k ~ w k 2 S C 1 k B T Bw k 2 S ; (ii) for an y 2 H 0 and w 2 W w ? Ker S there is ~ w 2 W suc h that B ~ w 2 H h ; B ~ w i = h ; Bw i ; and k B T B ~ w k 2 S C 2 k w k 2 S : Then ( RDPF ) C 1 C 2 Pr o of. Lemma 35 applied to the operator RD together with Lemma 36 giv e k QD k 2 H!H 0 1 ; k ( RD ) 1 k 2 H 0 !H 1 : (4.12) F rom assumption (i), w e ha v e sup w 2 W h ; Bw i 2 k w k 2 S 1 C 1 sup w 2 W; Bw 2H h ; Bw i 2 k B T Bw k 2 S (4.13) while assumption (ii) giv es the con v erse bound sup w 2 W;w 2 Ker S h ; Bw i 2 k w k 2 S C 2 sup w 2 W; Bw 2H h ; Bw i 2 k B T Bw k 2 S : (4.14) Using (4.10) and Lemma 37, w e see that the inequalities (4.13) and (4.14) imply the inequalit y 1 C 1 k k 2 H 0 h ; PF i C 2 k k 2 H 0 ; 8 2 H 0 : Applying Lemma 35 to the operator PF w e obtain k PF k 2 H 0 !H C 2 ; k ( PF ) 1 k 2 H!H 0 C 1 : (4.15) F rom (4.12) and (4.15), w e ha v e k QDPF k H 0 !H 0 k QD k H!H 0 k PF k H 0 !H C 2 80

PAGE 90

and k ( QDPF ) 1 k H 0 !H 0 k ( QD ) 1 k H 0 !H k ( PF ) 1 k H!H 0 C 1 : The result follo ws. 4.3 Analysis of the Original FETI In this section, w e use the abstract analysis framew ork from the previous section to pro v e a bound on the condition n um ber of FETI for a second elliptic problem, utilizing results from Section 4.1. 4.3.1 Assumptions W e need more specic assumptions in order to be able to pro v e a bound on the condition n um ber Let us recall the model problem (2.1). Assume w e are solving the boundary v alue problem A u = g in n ; u = 0 on @ n ; where A is a second order elliptic operator A v = d X i;j =1 @ @x i ( x ) @v ( x ) @x j ; with ( x ) a measurable function suc h that 0 < 0 ( x ) 1 a.e. in n. The domain n is assumed to be divided in to non-o v erlapping subdomains n i i = 1 ; :::; N s that are shape regular of diameter H (Denition 21). Assume that V h (n) is a conforming P1 or Q1 nite elemen t space on a triangulation of n, whic h satises the standard regularit y and in v erse assumptions (cf. Section 4.1). In particular, w e recall that the degrees of freedom are v alues at nodes of the triangulation h denotes the c haracteristic elemen t size. Eac h subdomain n i is assumed to be a union of some of the elemen ts, and all functions in V h (n) are zero on @ n. 81

PAGE 91

Let us dene the restriction operator R i : W W i b y the equation R i w = w i ; where w = [ w T 1 w T 2 : : :w T N s ] T w i 2 W i ; i = 1 ; : : :; N s W e will assume that the in terface con tin uit y operator B is dened as follo ws. Let w r ( x ) and w s ( x ) denote the pair of degrees of freedom corresponding to the mesh node x 2 @ n r \ @ n s and let ( B w ) rs ( x ) be the en try of v ector B w that corresponds to the mesh poin t x and subdomains n r and n s W e require eac h suc h ( B w ) rs ( x ) to ha v e form ( Bw ) rs ( x ) = rs ( w r ( x ) w s ( x )) ; (4.16) where rs is either 1 or -1. Note that rs does not depend on x that is, coecien ts rs are uniform along edges (and faces, in 3D) bet w een t w o subdomains. F or eac h node x that belongs to the in terface of t w o and more subdomains, x 2 T @ n i ; i = s 1 ; s 2 ; : : :; s n v ector Bw con tains n 1 en tries, ( Bw ) s k s k +1 ( x ) ; k = 1 ; : : :; s n 1 F or an example of the denition of the v alues of Bw for ( s 1 ; s 2 ; s 3 ) = (1 ; 3 ; 2) in 2D around a crosspoin t, see Fig. 4.1. W e poin t out that B c hosen in this w a y has full ro w rank, that is, this denition implies that there are no redundan t constrain ts in enforcing the con tin uit y of the solution at the nodes where more than t w o subdomains meet. Only the impro v ed estimate in statemen t 3 of Lemma 39 will require the specic denition of this section. If redundan t constrain ts are allo w ed, the estimates, with the exceptions of the impro v ed estimates, still hold whic h can be sho wn b y minor modications of the proofs. See Section 4.4. An additional connectivit y assumption on the decomposition is needed. Let N rs ; r ; s = 1 ; : : :; N s ; r 6 = s be the n um ber of in terface conditions 82

PAGE 92

bet w een subdomains n r and n s N rr the n um ber of degrees of freedom of subdomain n r and N = max N s i =1 N ii W e will assume that there exists constan ts c and n 0 independen t of h and H suc h that for all r ; s = 1 ; : : :; N s ; r 6 = s for whic h N rs > 0, there exists a sequence of indices f r i g r 0 = r r k = s k < n 0 suc h that N r i 1 r i cN 8 i = 1 ; : : :; k: Throughout the next section, c; C; c 1 ; c 2 ; c 3 ; c 4 ; c 5 and c 6 denote positiv e constan ts independen t of H and h 4.3.2 Discrete Norm Bounds The stiness matrices K i ; i = 1 ; : : :; N s satisfy the assumption of Lemma 26. Therefore, using Lemma 26 for eac h subdomain and summing o v er all subdomains, w e obtain c 1 N s X i =1 j I P 1 R i w j 2 1 2 ; 2 ;@ n i k w k 2 S c 2 N s X i =1 j I P 1 R i w j 2 1 2 ; 2 ;@ n i 8 w 2 W ; (4.17) where the positiv e constan ts c 1 ; c 2 are independen t of the c haracteristic mesh size h and the subdomain diameter H The constan ts c 1 ; c 2 ma y depend on the regularit y of the shape of subdomains. The follo wing lemmas v erify the assumptions of Theorem 38. Lemma 39 F or all 2 Ker G T and all w 2 W suc h that B T Bw ? Ker S there exists ~ w 2 W suc h that h ; B ~ w i = h ; Bw i and k ~ w k 2 S C (1 + log H=h ) k B T Bw k 2 S : where = 1, and = 0 in the follo wing special cases: (1) 1 2 B B T = I whic h means that there are no nodes shared b y more than t w o subdomains. 83

PAGE 93

Figure 4.1. Denition of B A A A A A A A v v v v v v v v v r r r r r n 1 n 2 n 3 @ n 3 @ n 2 @ n 1 w 3 ( x ) w 3 ( x ) w 3 ( x r ) w 2 ( x r ) w 1 ( x ) w 1 ( x ) w 2 ( x ) w 2 ( x ) w 1 ( x ) w 1 ( x ) w 3 ( x ) w 1 ( x ) w 3 ( x ) w 2 ( x ) w 3 ( x ) w 2 ( x r ) w 3 ( x r ) w 1 ( x ) w 2 ( x ) 84

PAGE 94

(2) d = 2, and the matrix B has the follo wing propert y: If w 2 Im B T x is a crosspoin t (node shared b y more than t w o subdomains), I P 1 R i w ( x ) = I P 1 R i w ( y ) for all i suc h that x 2 @ n i and all nodes y that are adjacen t to x on @ n i then I P 1 R i w ( x ) = 0 for all i suc h that x 2 @ n i (3) d = 2, B is dened as in Section 4.3.1 and all nodes in the triangulation belong to either one, t w o, or an odd n um ber of subdomains. Pr o of. Let us rst pro v e that, in the general case, w e obtain 1. Let w 2 W and dene ~ w = B T ( B B T ) 1 Bw The triangle inequalit y sho ws that k ~ w k S r r r r 1 2 B T Bw r r r r S + r r r r r 1 2 B T I 1 2 B B T 1 Bw r r r r r S : (4.18) Denote z = 1 2 B T ( I ( 1 2 B B T ) 1 ) Bw F rom the denition of B in (4.16), z is zero at all nodes that belong to at most t w o subdomains. The remaining nodes lie on crosspoin ts or edges (in the 3D case) of subdomains. F rom the denition of B at ev ery suc h node x z is a linear com bination of the en tries of B T Bw that correspond to the same node x and the coecien ts of the linear com binations are bounded only in terms of the n um ber of subdomains the node belongs to. Using Lemma 24 for the crosspoin ts of subdomains and the equiv alence (4.17), w e obtain for the 2D case that k z k 2 S c N s X i =1 j I P 1 R i z j 2 1 2 ; 2 ;@ n i C (1 + log( H=h )) N s X i =1 k I P 1 R i B T Bw k 2 1 2 ; 2 ;@ n i : (4.19) The P oincar e inequalit y Lemma 27, yields N s X i =1 k I P 1 R i B T Bw k 2 1 2 ; 2 ;@ n i k B T Bw k 2 S In the 3D case, the argumen t for subdomain crosspoin ts is the same. In addition, w e note that the coecien ts of the linear com bination do not c hange along a subdomain edge, so it remains to apply Lemma 24 on ev ery edge. 85

PAGE 95

Let us no w turn to the special cases that giv e = 0. If 1 2 B B T = I w e c hoose ~ w = 1 2 B T Bw = w whic h pro v es the special case 1. No w w e pro v e special case 2. F rom the denition of the H 1 = 2 norm, the equiv alence of norms (4.17), and the fact that I P 1 R i B T Bw is a piecewise linear function, it follo ws that k B T Bw k 2 S c N s X i =1 j I P 1 R i B T Bw j 2 1 2 ; 2 ;@ n i (4.20) c X x crosspoin t ;x 2 @ n i y adjacen t to x;y 2 @ n i I P 1 R i B T Bw ( x ) I P 1 R i B T Bw ( y ) 2 : F or an y crosspoin t x it follo ws from the assumption that for ev ery w 2 Im B T X i;@ n i 3 x y adjacen t to x; y 2 @ n i ( I P 1 R i w ( x ) I P 1 R i w ( y )) 2 = 0 ) X i;@ n i 3 x ( I P 1 R i w ( x )) 2 = 0 : Consequen tly b y compactness, and since there are only nitely man y dieren t n um bers of subdomains sharing a crosspoin t, for all w 2 Im B T X i;@ n i 3 x ( I P 1 R i w ( x )) 2 C X i;@ n i 3 x y adjacen t to x; y 2 @ n i ( I P 1 R i w ( x ) I P 1 R i w ( y )) 2 : By summation o v er all crosspoin ts x using Lemma 25 and (4.20), w e get k z k 2 S C k B T u k 2 S ; whic h concludes the proof of this case. In order to pro v e case 3, w e v erify the assumptions of case 2. W e form ulate the proof only for a crosspoin t shared b y three subdomains (Fig. 4.1). The proof is similar for a dieren t odd n um ber of subdomains. Let w 2 Im B T Since I P 1 R 1 w ( x ) I P 1 R 1 w ( x ) = 0, and I P 1 R 1 w ( x ) I P 1 R 1 w ( x ) = 0, w e ha v e I P 1 R 1 w ( x ) = I P 1 R 1 w ( x ). Similarly w e obtain I P 1 R 2 w ( x ) = 86

PAGE 96

I P 1 R 2 w ( x r ), and I P 1 R 3 w ( x ) = I P 1 R 3 w ( x r ). Moreo v er w 2 Im B T implies I P 1 R 1 w ( x ) = I P 1 R 2 w ( x ), I P 1 R 2 w ( x r ) = I P 1 R 3 w ( x r ), and I P 1 R 3 w ( x ) = I P 1 R 1 w ( x ), whic h can be satised only if I P 1 R 1 w ( x ) = I P 1 R 1 w ( x ) = ::: = 0. Remark 40 In general, the exponen t = 1 in Lemma 39 cannot be impro v ed. T o see that, let us consider the conguration with v alues of u and Bu in the neigh borhood of a crosspoin t as in Fig. 4.2, whic h violate the assumptions of special case 2. Extending the v alues of u in Fig. 4.2 to deca y as log r ( t=H ), r < 1 = 2, where t is the distance from the crosspoin t, w e obtain a v ector u 2 suc h that k B T u k S C; k I P 1 u k 1 2 ; 2 ;@ n 1 \ @ n 2 j log h=H j r : If u = Bw then on @ n 1 \ @ n 2 u = w 2 w 1 whic h giv es j I P 1 u j 1 2 ; 2 ;@ n 1 \ @ n 2 j I P 1 w 1 j 1 2 ; 2 ;@ n 1 \ @ n 2 + j I P 1 w 2 j 1 2 ; 2 ;@ n 1 \ @ n 2 j I P 1 w 1 j 1 2 ; 2 ;@ n 1 + j I P 1 w 2 j 1 2 ; 2 ;@ n 2 k w k S ; and therefore k w k S C ( r ) j log h=H j r for all r < 1 = 2. Lemma 41 Let 2 Ker G T Then for all w 2 W w ? Ker S there is a ~ w 2 W suc h that B T B ~ w ? Ker S h ; Bw i = h ; B ~ w i ; and k B T B ~ w k 2 S C (1 + log H=h ) 2 k w k 2 S Pr o of. Let w 2 W be arbitrary and put B ~ w = PBw Then h ; Bw i = h ; B ~ w i : (4.21) Denote Bz = Bw B ~ w z 2 Ker S Then, since P is an orthogonal projection, k Bz k k Bw k : 87

PAGE 97

Figure 4.2. Coun ter-example t t t t t t t t t t t t q q q q q q q n 1 n 2 n 3 n 4 1 1 1 1 1 0 1 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 1 = 2 88

PAGE 98

No w, from the denition of B and the P oincar e inequalit y Lemma 27, w e obtain h k Bw k 2 2 h k w k 2 CH k w k 2 S : Also, since z 2 Ker S it is constan t on eac h @ n i using the connectivit y assumption, w e ha v e the follo wing b y Lemma 24: k B T Bz k 2 S C h H k Bz k 2 (1 + log H=h ) 2 : T ogether this yields k B T Bz k 2 S C (1 + log H=h ) 2 k w k 2 S : By the denition of B R i B T Bw on @ n i [ @ n j is a linear com bination (with bounded coecien ts) of (a bounded n um ber of) w k from all @ n k adjacen t to @ n i [ @ n j F rom Lemma 24 and the P oincar e-F riedric hs inequalit y Lemma 27, w e obtain k B T Bw k S C (1 + log( H=h )) k w k S ; 8 w 2 W: Finally summarizing, k B T B ~ w k S k B T Bw k S + k B T Bz k S C (1 + log H=h ) k w k S : F rom this and (4.21), the statemen t of the lemma follo ws. 4.3.3 Condition Num ber Estimate W e are no w ready to pro v e the nal result. It follo ws from the abstract estimate in Lemma 35 with the assumptions v eried b y Lemma 39 and Lemma 41. 89

PAGE 99

Theorem 42 Under the assumptions of section 4.3.1, the condition n um ber of the FETI method with the Diric hlet preconditioner satises = max ( PDPF ) min ( PDPF ) C (1 + log H h ) r with r = 3, and r = 2 in the special cases listed in Lemma 39. 4.4 Analysis of the Method b y P ark W e consider the same assumptions as for the original FETI for second order problems except for the c hoice of B W e consider B giv en b y projection matrix P L (3.43). The essen tial propert y of this matrix B used here is that B = B T B = B T since B is an orthogonal projection. W e presen t t w o lemmas that v erify assumptions of Theorem 38. Lemma 43 F or all 2 Ker G T = Ker B T and all w 2 W suc h that Bw 2 Ker G T = Ker B T there exists ~ w 2 W suc h that h ; B ~ w i = h ; Bw i ; k ~ w k S = k B T Bw k S Pr o of. Since B T = B and B 2 = B w e ma y c hoose ~ w = w Lemma 44 Let 2 Ker G T = Ker B T Then for all w 2 W w ? Ker S there is a ~ w 2 W suc h that B T B ~ w ? Ker S and h ; Bw i = h ; B ~ w i ; k B T B ~ w k 2 S C (1 + log H=h ) 2 k w k 2 S Pr o of. Let w 2 W and B ~ w = PBw Then h ; Bw i = h ; B ~ w i (4.22) Let z 2 Ker S Bz = Bw B ~ w Then, since P is an orthogonal projection, k Bz k k Bw k : 90

PAGE 100

No w, from the denition of B and the P oincar e inequalit y w e obtain h k Bw k 2 2 h k w k 2 CH k w k 2 S : Also, since z 2 Ker S it is constan t on eac h @ n i and using the connectivit y assumption, w e ha v e the follo wing b y Lemma 24 k Bz k 2 S C h H k Bz k 2 (1 + log H=h ) 2 : T ogether this yields k Bz k 2 S C (1 + log H=h ) 2 k w k 2 S : By the denition of B ( Bw ) i on @ n i [ @ n j is a linear com bination (with bounded coecien ts) of (a bounded n um ber of) w k from all @ n k adjacen t to @ n i [ @ n j F rom Lemma 24 and the P oincar e-F riedric hs inequalit y (Lemma 27), k Bw k S C (1 + log( H=h )) k w k S ; 8 w 2 W: Finally summarizing, k B ~ w k S k Bw k S + k B T Bz k S C (1 + log H=h ) k w k S : F rom this, (4.22), and B T B = B the result follo ws. Theorem 45 Under the assumption of Section 4.3.1, the condition n um ber of the P ark's v arian t of FETI method with the Diric hlet preconditioner (3.39) satises = max ( QDPF ) min ( QDPF ) C 1 + log H h 2 : Pr o of. Lemmas 43 and 44 v erify the assumptions (i) and (ii) of Theorem 38 with C 1 = 1 and C 2 = C (1 + log H=h ) 2 respectiv ely 91

PAGE 101

4.5 Con v ergence Estimates for Plate Bending The follo wing appr oximate p ar ametric variational principle form ulated in [42 ] will allo w us to considerate to the plate bending problem and the biharmonic equation. Assumption 46 ([41, 42]) W e consider elemen ts with displacemen t and rotation degrees of freedom at the v ertices only and assume that there exist constan ts c 1 > 0, c 2 suc h that if the plate thic kness t satises 0 < t h then for eac h elemen t T the local stiness matrix K T satises c 1 K HCT T K T c 2 K HCT T (4.23) where K HCT T is the HCT elemen t lev el stiness matrix of the biharmonic equation [15], with the rotations in terpreted as deriv ativ es of the transv ersal displacemen t in the HCT elemen t. That is, as the thic kness of the plate goes to zero, the stiness matrix of the elemen t should be spectrally equiv alen t to that of the HCT elemen t for the biharmonic equation. W e refer to [14 ] and Section 4.1.2 for more details on HCT elemen ts. Here w e just summarize that HCT elemen ts are in C 1 and they use cubic splines for v alues on elemen t sides, linear in terpolation for normal deriv ativ es on the sides, and piecewise polynomial extension in to the elemen t in terior. In [41 ], Assumption 46 is v eried for the particular case of the DKT elemen t [5]. Assumption 46 also holds for the follo wing general class of non-loc king P 1 Reissner-Mindlin elemen ts. Theorem 47 ([41 ]) Assume that the energy functional for an elemen t T is spectrally equiv alen t to Z T jr j 2 dx + 1 t 2 + h 2 Z T j r u j 2 dx; (4.24) 92

PAGE 102

where h = diam ( T ), u 2 P 1 ( T ) is the transv ersal displacemen t, and 2 ( P 1 ( T )) 2 is the rotation. Then (4.23) holds. Elemen ts with the energy functional of the form (4.24) include the DKT elemen t as restated in [60 ]. It should be noted that for the related Timoshenk o beam elemen t, the thin limit is exactly the discretization b y cubic splines of the biharmonic equation [35]. 4.5.1 Assumptions Before pro ving a bound on the condition n um ber of the generalized FETI method, w e need to in troduce some specic assumptions. W e refer to the model problem (2.10) and using the spectral equiv alence, that w e discussed in the previous section, w e consider the biharmonic boundary v alue problem in a v ariational form. Find u 2 H 2 0 (n) suc h that a ( u; v ) = f ( v ) ; 8 v 2 H 2 0 (n) ; where a ( u; v ) = Z n @ 11 u@ 11 v + @ 12 u@ 12 v + @ 22 u@ 22 v; 8 u; v 2 H 2 0 (n) Let the domain n be divided in to non-o v erlapping subdomains n i i = 1 ; :::; N s that are shape regular of diameter O ( H ) (Denition 21). Without loss of generalit y w e suppose that H < 1. F urthermore, w e require that, for all r ; s = 1 ; : : :; N s if @ n r \ @ n s is an edge, then its length is larger than a certain prescribed fraction of the length of @ n r W e note that this condition implies the connectivit y condition as that of the original FETI. It is possible to pro v e the result using the general condition, but, for the sak e of simplicit y w e pro v e it using the stronger assumption. 93

PAGE 103

W e assume that the problem is discretized using reduced HCT elemen ts. The general case of plate bending then follo ws from spectral equiv alence of the local elemen t stiness matrices follo wing Assumption 46. Let V HCT h (n) denote the corresponding nite elemen t space, and h denote the c haracteristic elemen t size. Eac h subdomain n i is assumed to be a union of some of the elemen ts. The degrees of freedom are v alues of the transv ersal displacemen t and its deriv ativ es (rotations) at the nodal poin ts of the discretization. W e dene B as follo ws, cf., Fig. 4.3. F or eac h node x on subdomain in terface and eac h pair ( r ; s ) suc h that x 2 n r \ @ n s and @ n r and @ n s share an edge (i.e, do not meet only at a crosspoin t), Bw includes one sub v ector of three elemen ts, ( Bw ) rs ( x ) = rs ( w r ( x ) w s ( x )) ; where w r ( x ) and w s ( x ) denote sub v ectors of w con taining the three degrees of freedom associated with node x and subdomain n r and n s respectiv ely and rs is c hosen to be either 1 or 1. W e poin t out that rs does not depend on x 2 @ n r \ @ n s i.e., coecien ts rs are uniform along the edge bet w een n r and n s Note that, unlik e in the case described in Section 4.3.1, this denition implies redundan t constrain ts at subdomain crosspoin ts. Let us mak e the denition of the matrix C from Section 3.2.3 more precise. Let m = dim and f k j ; j = 1 ; : : :; n g be the the complete list of indices of Lagrange m ultipliers corresponding to conditions (expressed through B ) enforcing con tin uit y of the solution (but not the con tin uit y of the deriv ativ es) at crosspoin ts. Then C = [ c ij ] is the m n matrix satisfying c ij = 1 if i = k j 94

PAGE 104

c ij = 0 otherwise : 4.5.2 Discrete Norm Bounds Using Lemma 33 for subdomain stiness matrices K i for eac h subdomain and summing o v er all subdomains, w e obtain c 1 N s X i =1 jr I HCT R i w j 2 1 2 ; 2 ;@ n i k w k 2 S c 2 N s X i =1 jr I HCT R i w j 2 1 2 ; 2 ;@ n i 8 w 2 W ; (4.25) where the positiv e constan ts c 1 ; c 2 are independen t of the c haracteristic mesh size h and the subdomain diameter H The constan ts c 1 ; c 2 ma y depend on the regularit y of the shape of subdomains. F or denitions of Sobolev spaces, see, for example [1]. Throughout this section, c; C; c 1 ; c 2 ; c 3 ; c 4 ; c 5 and c 6 denote positiv e constan ts independen t of H and h The next t w o lemmas con tain the principal tec hnical estimates. Lemma 48 F or all 2 V 0 and all w 2 W suc h that Bw 2 V there exists a ~ w 2 W suc h that h ; B ~ w i = h ; Bw i and k ~ w k 2 S C (1 + log H=h ) k B T Bw k 2 S : where = 1, and = 0 if 1 2 B B T = I whic h happens when there are no nodes shared b y more than t w o subdomains. Pr o of. Let us rst pro v e that, in the general case, w e obtain 1. Let w 2 W and Bw 2 V That is Z T B T Bw = 0 and C T Bw = 0. W e dene ~ w = B T ( B B T ) + Bw Then, B ~ w = Bw By the triangle inequalit y w e ma y write k ~ w k S k 1 2 B T Bw k S + k 1 2 B T ( I ( 1 2 B B T ) + ) Bw k S : (4.26) Denote z = 1 2 B T ( I ( 1 2 B B T ) + ) Bw F rom the denition of B z is zero at all nodes that belong to at most t w o subdomains. The remaining 95

PAGE 105

nodes lie on subdomain crosspoin ts. A t ev ery suc h node, I HCT R i z ( x ) is a linear com bination of the en tries of B T Bw that correspond to the same node x and the coecien ts of the linear com binations are bounded only in terms of the n um ber of subdomains to whic h the node belongs. In addition, since C T Bw = 0, the transv ersal displacemen t componen ts of z at crosspoin ts are zero. Using (4.25) and Lemma 24 for the subdomain crosspoin t v ertices, one subdomain at a time, w e obtain k z k 2 S C (1 + log H h ) N s X i =1 1 H kr I HCT R i B T Bw k 2 0 ; 2 ;@ n i + jr I HCT R i B T Bw j 2 1 2 ; 2 ;@ n i This together with the P oincar e inequalit y Lemma 34, and (4.26) yields the result. If 1 2 B B T = I w e simply c hoose ~ w = 1 2 B T Bw = w No w w e deriv e the con v erse bound. Lemma 49 F or all 2 V 0 and w 2 W; w ? Ker S there is a ~ w 2 W suc h that B ~ w 2 V h ; Bw i = h ; B ~ w i ; and k B T B ~ w k 2 S C 1 + log H h 2 k w k 2 S : Pr o of. Let 2 V 0 w 2 W; w ? Ker S and consider B ~ w = P Bw + PFPC Since 2 V 0 it is easily v eried that h ; Bw i = h ; B ~ w i V ector can be found from the condition C T B ~ w = 0, whic h can be rewritten as h PFPC; C ~ i = h P Bw; C ~ i 8 ~ : By the denition of F w e can rewrite this as h S + B T PC; B T PC ~ i = h w; B T PC ~ i 8 ~ : 96

PAGE 106

Since S + is positiv e semidenite, the equation yields suc h that h S + 1 = 2 B T PC; S + 1 = 2 B T PC i = h S 1 = 2 w; S +1 = 2 B T PC i and, from the Cauc h y inequalit y k S + 1 = 2 B T PC k k S 1 = 2 w k Therefore, k S + B T PC k S k w k S : (4.27) W e need to estimate k B T B ~ w k S Let B ij be the matrix constructed from B b y zeroing out all the ro ws that do not correspond to the in terface conditions bet w een n i and n j Then, B T B ~ w = N s X i;j =1 ;i
PAGE 107

Denote u = w + S + B T PC Then, B ~ w = P Bu The triangle inequalit y and (4.27) yield k u k S k w k S + k S + B T PC k S 2 k w k S : So no w w e only ha v e to pro v e that X k 1 H kr I HCT P Bu k 2 0 ; 2 ; k + jr I HCT P Bu j 2 1 2 ; 2 ; k c 3 k u k 2 S : W e rst use the triangle inequalit y to get jr I HCT P Bu j 1 2 ; 2 ; k jr I HCT Bu j 1 2 ; 2 ; k + jr I HCT ( I P ) Bu j 1 2 ; 2 ; k : F rom the denition of B eac h en try reduces to a linear com bination with coecien ts 1,constan t along ev ery edge. Since ( I P ) Bu 2 Im G I HCT ( I P ) Bu is a restriction of a linear function on ev ery edge, jr I HCT ( I P ) Bu j 1 2 ; 2 ; k = 0. F urthermore, X k jr I HCT Bu j 2 1 2 ; 2 ; k 2 N s X i =1 jr I HCT u j 2 1 2 ; 2 ;@ n i : (4.28) T o estimate the L 2 terms, w e again rst use the triangle inequalit y: kr I HCT P Bu k 0 ; 2 ; k 2 kr I HCT u k 0 ; 2 ; k + kr I HCT ( I P ) Bu k 0 ; 2 ; k : Since r I HCT ( I P ) Bu is a linear function on ev ery edge and P is an orthogonal projection, w e ha v e X k kr I HCT ( I P ) Bu k 2 0 ; 2 ; k c 4 h k ( I P ) Bu k 2 c 5 h k Bu k 2 2 c 5 h k u k 2 : F urthermore, Lemma 32 sho ws that h k u k 2 c 6 (1 + H 2 ) N s X i =1 kr I HCT u k 2 0 ; 2 ;@ n i Finally since u ? Ker S the P oincar e inequalit y (Lemma 34) and the equiv alence of the norms (4.25) conclude the estimate. 98

PAGE 108

4.5.3 Condition Num ber Estimate W e ha v e no w ev erything ready to pro v e the nal result. Theorem 50 Under the assumptions of Section 4.5.1, the condition n um ber of the generalized FETI method with the Diric hlet preconditioner (3.39) satises = max ( QDPF ) min ( QDPF ) C 1 + log H h r with r = 3, and r = 2 if there are no crosspoin ts bet w een more than t w o subdomains. Pr o of. Since V and V 0 can be replaced b y the factorspaces ~ V and ~ V 0 Lemmas 48 and 49 v erify the assumptions (i) and (ii), of Theorem 38, with C 1 = C (1 + log H=h ) = 0 or 1, and C 2 = C (1 + log H=h ) 2 respectiv ely 99

PAGE 109

Figure 4.3. Denition of B for plate bending t t t t t t t t t t t t t t t t t r r r r r r r r r r n 1 n 2 n 3 n 4 @ n 1 @ n 2 @ n 3 @ n 4 w 4 ( x A ) w 4 ( x B ) w 4 ( x C ) w 3 ( x C ) w 3 ( x D ) w 3 ( x E ) w 1 ( x A ) w 1 ( x B ) w 1 ( x C ) w 2 ( x C ) w 2 ( x D ) w 2 ( x E ) w 1 ( x F ) w 2 ( x F ) w 4 ( x G ) w 3 ( x G ) w 1 ( x A ) w 4 ( x A ) w 1 ( x B ) w 4 ( x B ) w 1 ( x C ) w 4 ( x C ) w 2 ( x C ) w 3 ( x C ) w 2 ( x D ) w 3 ( x D ) w 2 ( x E ) w 3 ( x E ) w 1 ( x F ) w 2 ( x F ) w 1 ( x C ) w 2 ( x C ) w 4 ( x C ) w 3 ( x C ) w 4 ( x G ) w 3 ( x G ) 100

PAGE 110

5. Computational Results W e refer to the paper b y F arhat et al. [32] for n umerical results for second order problems conrming our theoretical results. The paper also deals with parallel implemen tation issues and performance. Here w e summarize results for plate bending problems based on the paper b y Mandel, T ezaur and F arhat [53 ], due to F arhat. W e consider the plate bending problem on a unit square. The plate is discretized b y a uniform mesh of three-node triangular DKT plate elemen ts and subjected to a uniform pressure load. The thic kness of the plate is 10 3 the Y oung modulus E = 10 6 and the P oisson coecien t = 0 : 3. The FETI and generalized FETI methods are preconditioned b y the Diric hlet preconditioner (Section 3.3) [32, 34], and the follo wing stopping criterion is used k z k 1 k k f k = 10 3 (5.1) where z k 1 is the preconditioned residual in Algorithm 20. This condition is a good estimator of the global error as sho wn in [34]. W e construct sev eral meshes with dieren t mesh size h and sev eral mesh partitions with dieren t subdomain size H and report the performance of the original FETI and the generalized FETI. Three series of computational experimen ts are reported. The condition n um ber is denoted b y and the n um ber of iterations b y n it in the tables. The rst series of experimen ts sho ws the performance of the method 101

PAGE 111

for three dieren t n um bers of subdomains ( H = 1 = 2 ; 1 = 4 ; 1 = 8) and three differen t meshes corresponding to h = 1 = 10, h = 1 = 20, and h = 1 = 40. Results in T able 5.1) demonstrate that for a giv en H the condition n um ber of the original FETI gro ws fast with the mesh size h while that of the new FETI method is m uc h smaller and gro ws only w eakly with h F or large n um ber of subdomains ( N s = 64) the new FETI method con v erges about 7 times faster than the original one. In the second series of experimen ts (T able 5.2), the mesh size is xed to h = 1 = 120 (28800 elemen ts and 86400 degrees of freedom), and H is v aried bet w een H = 1 = 2 (4 subdomains) and H = 1 = 12 (144 subdomains). In this case, the condition n um bers of both FETI methods are sho wn to decrease with the n um ber of subdomains. This is an expected result because when h is xed and H is decreased, the size of the coarse problem increases for both algorithms. The generalized FETI sho ws only a small dependence on the n um ber of subdomains as predicted b y the theory F or large n um ber of subdomains, the new FETI method clearly outperforms the original FETI. Finally the subdomain problem size is xed to h=H = 1 = 15, and the n um ber of subdomains as w ell as the size of the global problem are increased. The performance results reported in T able 5.3 sho w that in this case too, the new FETI method outperforms signican tly the original one. Since solution time is ultimately the most importan t criterion for assessing performance, Both FETI methods w ere also benc hmark ed for the same plate bending problem described abo v e with 960000 degrees of freedom and 64 subdomains. The performance results obtained on a 64-processor IBM SP2 are summarized in T able 5.4. They sho w that ev en though the new FETI 102

PAGE 112

T able 5.1. Fixed n um ber of subdomains, series of rened meshes FETI Generalized FETI h H n it n it 2x2 subdomains, H= 1 2 1 10 18 2578 12 7.6 1 20 22 30101 15 12.6 1 40 26 409987 17 18.6 4x4 subdomains, H= 1 4 1 10 61 6795 21 11.5 1 20 86 84199 27 17 1 40 119 1038120 36 24.4 8x8 subdomains, H= 1 8 1 10 172 21707 25 13 1 20 247 275004 34 19.4 1 40 323 3920613 42 27.6 T able 5.2. Fixed global mesh 120x120, h= 1 120 series of rened mesh partitions Decomposition FETI Generalized FETI H h H n it n it 1 2 1 60 27 2079032 18 23.2 1 3 1 40 64 839240 29 22.4 1 4 1 30 104 391470 32 21 1 5 1 24 135 234504 33 19.9 1 6 1 20 164 160173 32 18.6 1 8 1 15 222 94285 31 16.6 1 10 1 12 255 63896 29 14.9 1 12 1 10 245 46921 27 13.6 103

PAGE 113

T able 5.3. Fixed local mesh (15x15, h H = 1 15 ), series of rened meshes and mesh partitions Decomposition FETI Generalized FETI H h n it n it 1 2 1 30 20 11088 13 10 1 3 1 45 49 19004 21 13.4 1 4 1 60 74 29041 25 14.6 1 5 1 75 109 40120 28 15.4 1 6 1 90 145 55068 29 15.9 1 8 1 120 222 94285 31 16.6 1 10 1 150 318 144556 32 16.9 104

PAGE 114

T able 5.4. P erformance results for 960000 degrees of freedom and 64 subdomains FETI Generalized FETI n it T otal time Time per iter. n it T otal time Time per iter. 314 265 s 0.8 s 45 105 s 1.1 s method consumes an amoun t of CPU time (55.5 s.) equiv alen t to that of 50 of its iterations to set up and preprocess the coarse problem (3.33), and ev en though eac h of its iterations is 1.3 times more expensiv e than an iteration of the original FETI method, the new FETI method is 2.5 times faster than the original one at solving the system of 960000 plate bending equations on a 64-processor IBM SP2. 105

PAGE 115

REFERENCES [1] R. A. Adams. Sob olev sp ac es Academic press, New Y ork, 1975. [2] V. I. Agoshk o v. P oincar e{Steklo v's operators and domain decomposition methods in nite dimensional spaces. In R. Glo winski, G. H. Golub, G. A. Meuran t, and J. P eriaux, editors, First International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations pages 73{112, Philadelphia, 1988. SIAM. [3] Iv o Babu sk a, Alan W. Craig, Jan Mandel, and Juhanni Pitk aran ta. Efcien t preconditioning for the p -v ersion nite elemen t method in t w o dimensions. SIAM J. Numer. A nal. 28:624{662, 1991. [4] N. S. Bakh v alo v and A. V. Kn y azev. Preconditioned iterativ e methods in a subspace for linear algebraic equations with large jumps in the coecien ts. In Domain De c omp osition Metho ds in Scientic and Engine ering Computing: Pr o c e e dings of the Seventh International Confer enc e on Domain De c omp osition v olume 180 of Contemp or ary Mathematics pages 157{162, Pro vidence, Rhode Island, 1994. American Mathematical Societ y [5] J. L. Batoz. An explicit form ulation for an ecien t triangular plate bending elemen t. Int. J. Numer. Meth. Engr g. 18:1077{1081, 1982. [6] P E. Bjrstad and J. Mandel. Spectra of sums of orthogonal projections and applications to parallel computing. BIT 31:76{88, 1991. [7] J. F. Bourgat, R. Glo winski, P LeT allec, and M. Vidrascu. V ariational form ulation and algorithm for trace operator in domain decomposition calculations. In T. F. Chan, R. Glo winski, J. P eriaux, and O. B. Widlund, editors, Domain De c omp osition Metho ds pages 3{16. SIAM, Philadelphia, 1989. [8] J. H. Bram ble, J. E. P asciak, and A. H. Sc hatz. An iterativ e method for elliptic problems on regions partitioned in to substructures. Math. Comp. 46:361{369, 1986. 106

PAGE 116

[9] J. H. Bram ble, J. E. P asciak, J. W ang, and J. Xu. Con v ergence estimates for product iterativ e methods with applications to domain decomposition. Math. Comp. 57:1{21, 1991. [10] James H. Bram ble, Joseph E. P asciak, and Alfred H. Sc hatz. The construction of preconditioners for elliptic problems b y substructuring, I. Math. Comp. 47(175):103{134, 1986. [11] James H. Bram ble, Joseph E. P asciak, and Alfred H. Sc hatz. The construction of preconditioners for elliptic problems b y substructuring, IV. Math. Comp. 53:1{24, 1989. [12] Sue Brenner. Tw o-lev el additiv e Sc h w arz preconditioners for plate bending problems. T o appear in Numer. Math. [13] T. F. Chan and T. P Mathew. The in terface probing tec hnique in domain decomposition. SIAM J. Matrix A nal. Applic. 13:212{238, 1992. [14] P G. Ciarlet. Basic error estimates for elliptic problems. In P .G. Ciarlet and J. L. Lions, editors, Handb o ok of Numeric al A nalysis v olume II, pages 17{352. North-Holland, Amsterdam, 1989. [15] R. W. Clough and J. L. T oc her. Finite elemen t stiness matrices for analysis of plate bending. In Pr o c. 1965 Conf. Matrix Metho ds Struct. Me ch., Wright-Patterson AFB, Ohio, AFFDL-TR-66-80 pages 515{546, 1966. [16] La wrence Co wsar. P ersonal comm unication, 1992. [17] Q. V. Dinh, R. Glo winski, J. He, V. Kw oc k, T. W. P an, and J. P eriaux. Lagrange m ultiplier approac h to ctitious domain methods: application to ruid dynamics and electro{magnetics. In D. E. Key es, T. F. Chan, G. A. Meuran t, J. S. Scroggs, and R. G. V oigt, editors, Fifth International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations pages 151{194, Philadelphia, 1992. SIAM. [18] Q. V. Dinh, R. Glo winski, and J. P eriaux. Solving elliptic problems b y domain decomposition methods with applications. In G. Birkho and A. Sc hoenstadt, editors, Elliptic Pr oblem Solvers II Academic Press, 1984. [19] M. R. Dorr. On the discretization of in terdomain coupling in elliptic boundary{v alue problems. In T. F. Chan, R. Glo winski, J. P eriaux, and 107

PAGE 117

O. B. Widlund, editors, Domain De c omp osition Metho ds pages 17{37. SIAM, Philadelphia, 1989. [20] M. Dryja. An additiv e Sc h w artz algorithm for t w o{and three{ dimensional nite elemen t elliptic problems. In T. F. Chan, R. Glo winski, J. P eriaux, and O. B. Widlund, editors, Domain De c omp osition Metho ds pages 168{ 172. SIAM, Philadelphia, 1989. [21] M. Dryja and O. B. Widlund. Domain decomposition algorithms with small o v erlap. SIAM J. Sci. Comput. 15:604{620, 1994. [22] Maksymilian Dryja. A method of domain decomposition for 3-D nite elemen t problems. In Roland Glo winski, Gene H. Golub, G erard A. Meuran t, and Jacques P eriaux, editors, First International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations Philadelphia, P A, 1988. SIAM. [23] Maksymilian Dryja, Barry F. Smith, and Olof B. Widlund. Sc h w arz analysis of iterativ e substructuring algorithms for elliptic problems in three dimensions. SIAM J. Numer. A nal. 31(6):1662{1694, Decem ber 1994. [24] Maksymilian Dryja and Olof B. Widlund. T o w ards a unied theory of domain decomposition algorithms for elliptic problems. In T on y Chan, Roland Glo winski, Jacques P eriaux, and Olof Widlund, editors, Thir d International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations, held in Houston, T exas, Mar ch 20-22, 1989 SIAM, Philadelphia, P A, 1990. [25] Maksymilian Dryja and Olof B. Widlund. Sc h w arz methods of NeumannNeumann t ype for three-dimensional elliptic nite elemen t problems. Comm. Pur e Appl. Math 48:121{155, 1995. [26] C. F arhat, P S. Chen, and J. Mandel. Scalable Lagrange m ultiplier based domain decomposition method for time-dependen t problems. Int. J. Numer. Meth. Engr g. 38:3831{3853, 1995. [27] C. F arhat and F.-X. Roux. A method of nite elemen t tearing and interconnecting and its parallel solution algorithm. Int. J. Numer. Meth. Engng. 32:1205{1227, 1991. [28] Charbel F arhat. Lagrange m ultiplier based divide and conquer nite elemen t algorithm. J. Comput. Sys. Engr g. 2:149{156, 1991. 108

PAGE 118

[29] Charbel F arhat, P o-Sh u Chen, Jan Mandel, and F rancois-Xa vier Roux. The t w o-lev el FETI method P art II: Extension to shell problems, parallel implemen tation and performance results. Comp. Meth. Appl. Mec h. Engrg, submitted, 1995. [30] Charbel F arhat, P o-Sh u Chen, F rancois Roux, and Jan Mandel. The t w o-lev el FETI method-part ii: Extension to shell problems, parallel implemen tation, and performance results. CU-CAS Report 96-10, Cen ter for Aerospace Structures, Univ ersit y of Colorado, 1996. [31] Charbel F arhat and Jan Mandel. The t w o-lev el FETI method for static and dynamic plate problems P art I: An optimal iterativ e solv er for biharmonic systems. T ec hnical Report UCB/CAS Report CU-CAS-95-23, Cen ter for Aerospace Structures, Univ ersit y of Colorado at Boulder, 1995. Comp. Meth. Appl. Mec h. Engrg, submitted. [32] Charbel F arhat, Jan Mandel, and F rancois-Xa vier Roux. Optimal conv ergence properties of the FETI domain decomposition method. Comput. Metho ds Appl. Me ch. Engr g. 115:367{388, 1994. [33] Charbel F arhat and F rancois-Xa vier Roux. An uncon v en tional domain decomposition method for an ecien t parallel solution of large-scale nite elemen t systems. SIAM J. Sci. Stat. Comput. 13:379{396, 1992. [34] Charbel F arhat and F rancois-Xa vier Roux. Implicit parallel processing in structural mec hanics. Comput. Me ch. A dvanc es 2:1{124, 1994. [35] Carlos A. F elippa. A surv ey of parametrized v ariational principles and application to computational mec hanics. Comput. Metho ds Appl. Me ch. Engr g 113:109{139, 1994. [36] R. Glo winski and M. F. Wheeler. Domain decomposition and mixed nite elemen t methods for elliptic problems. In R. Glo winski, G. H. Golub, G. A. Meuran t, and J. P eriaux, editors, First International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations pages 144{172, Philadelphia, 1988. SIAM. [37] Roland Glo winski and Mary F. Wheeler. Domain decomposition and mixed nite elemen t methods for elliptic problems. In Roland Glo winski, Gene H. Golub, G erard A. Meuran t, and Jacques P eriaux, editors, First International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations Philadelphia, P A, 1988. SIAM. 109

PAGE 119

[38] Gene H. Golub and Charles F. V an Loan. Matrix Computations Johns Hopkins Univ. Press, 1989. Second Edition. [39] Gene H. Golub and Charles F. V an Loan. Matrix Computations John Hopkins Univ ersit y Press, second edition, 1989. [40] Manoel R. Justino, K. C. P ark, and Carlos A. F elippa. An algebraically partitioned FETI method for parallel structural analysis: Implemen tation and n umerical performance ev aluation. International Journal of Numeric al Metho ds in Engine ering. 40:2739{2758, 1997. [41] P atric k Le T allec, Jan Mandel, and Marina Vidrascu. A NeumannNeumann domain decomposition algorithm for solving plate and shell problems. SIAM J. Numer. Anal., to appear. [42] P atric k Le T allec, Jan Mandel, and Marina Vidrascu. Balancing domain decomposition for plates. Contemp or ary Mathematics 180:515{524, 1994. Proceedings of the 7th In ternational Symposium on Domain Decomposition Methods, P enn State, No v em ber 1993. [43] P LeT allec, Y.-H. de Roec k, and M. Vidrascu. Domain{decomposition methods for large linearly elliptic three dimensional problems. J. Comput. Appl. Math. 34:93{117, 1991. [44] Pierre Louis Lions. On the Sc h w arz alternating method. I. In Roland Glo winski, Gene H. Golub, G erard A. Meuran t, and Jacques P eriaux, editors, First International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations Philadelphia, P A, 1988. SIAM. [45] J. Mandel. Balancing domain decomposition. Comm. Numer. Meth. Engr g. 9:233{241, 1993. [46] Jan Mandel. Balancing domain decomposition. Comm. in Numeric al Metho ds in Engr g. 9:233{241, 1993. [47] Jan Mandel. Hybrid domain decomposition with unstructured subdomains. In Y uri A. Kuznetso v, Jacques P eriaux, Alo Quarteroni, and Olof B. Widlund, editors, Domain De c omp osition Metho ds in Scienc e and Engine ering: The Sixth International Confer enc e on Domain De c omp osition v olume 157. AMS, 1994. Held in Como, Italy June 15{19,1992. [48] Jan Mandel and Marian Brezina. Balancing domain decomposition: Theory and computations in t w o and three dimensions. UCD/CCM Report 2, 110

PAGE 120

Cen ter for Computational Mathematics, Univ ersit y of Colorado at Den v er, 1993. [49] Jan Mandel and Marian Brezina. Balancing domain decomposition for problems with large jumps in coecien ts. Mathematics of Computation 65:1387{1401, 1996. [50] Jan Mandel and Radek T ezaur. Con v ergence of a substructuring method with Lagrange m ultipliers. Numerische Mathematik 73:473{487, 1996. [51] Jan Mandel, Radek T ezaur, and Charbel F arhat. Optimal Lagrange m ultiplier based domain decomposition method for plate bending problems. UCD/CCM Report 61, Cen ter for Computational Mathematics, Univ ersit y of Colorado at Den v er, 1995. [52] Jan Mandel, Radek T ezaur, and Charbel F arhat. Optimal Lagrange m ultiplier based domain decomposition method for plate bending problems. UCD/CCM Report 61, Cen ter for Computational Mathematics, Univ ersit y of Colorado at Den v er, 1995. [53] Jan Mandel, Radek T ezaur, and Charbel F arhat. A scalable substructuring method b y Lagrange m ultipliers for plate bending problems. UCD/CCM Report 84, Cen ter for Computational Mathematics, Univ ersit y of Colorado at Den v er, 1996. Submitted to SIAM J. Sci. Comput. [54] A. M. Matsokin and S. V. Nepomn y asc hikh. A Sc h w arz alternating method in subspaces. In. V uzov 10:61{66, 1985. Also in So viet Mathematics, 10 (1985), pp. 78{84. [55] S. V. Nepomn y asc hikh. Domain de c omp osition and the Schwarz metho d in subsp ac e for appr oximate solution of elliptic b oundary pr oblems PhD thesis, Computing Cen ter of the Siberian Branc h of the USSR Academ y of Sciences, No v osibirsk, Russia, 1986. [56] Jind ric h Ne cas. L es m etho des dir e ctes en th eorie des equations elliptiques Academia, Prague, 1967. [57] J. T. Oden, Abani P atra, and Y usheng F eng. P arallel domain decomposition solv er for adaptiv e hp nite elemen t methods. T ec hnical Report TICIAM 94-11, Univ ersit y of T exas at Austin, 1994. [58] K. C. P ark and Carlos A. F elippa. A v ariational framew ork for solution method dev elopmen ts in structural mec hanics. CU-CAS Report 96-21, 111

PAGE 121

Cen ter for Aerospace Structures, Univ ersit y of Colorado, 1996. T o appear in Journal of Applied Mathematics. [59] K. C. P ark, Manoel R. Justino, and Carlos A. F elippa. An algebraically partitioned FETI method for parallel structural analysis: Algorithm description. International Journal of Numeric al Metho ds in Engine ering. 40:2717{2737, 1997. [60] Juhanni Pitk aran ta. On a simple nite elemen t method for plate bending problems. In W olfgang Hac kbusc h and Kristian Witsc h, editors, Numeric al T e chniques in Continuum Me chanics n um ber 16 in Notes on Numerical Fluid Mec hanics, pages 84{101. View eg, Braunsc h w eig/Wiesbaden, 1987. Proceedings of 2nd GAMM-Seminar, Kiel, Jan uary 1986. [61] Y.-H. de Roec k and P LeT allec. Analysis and test of a local domain decomposition. In R. Glo winski, Y u. A. Kuznetso v, G. A. Meuran t, J. P eriaux, and O. B. Widlund, editors, F ourth International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations pages 112{128, Philadelphia, 1991. SIAM. [62] F.-X. Roux. Acceleration of the outer conjugate gradien t b y reorthogonalization for a domain decomposition method with Lagrange m ultiplier. In T. F. Chan, R. Glo winski, J. P eriaux, and O. B. Widlund, editors, Thir d International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations pages 314{321, Philadelphia, 1990. SIAM. [63] H. A. Sc h w arz. Uber einige Abbildungsaufgaben. Ges. Math. A bh. 11:65{ 83, 1869. [64] Radek T ezaur, P etr V an ek, and Marian Brezina. Tw o-lev el method for solids on unstructured meshes. UCD/CCM Report 73, Cen ter for Computational Mathematics, Univ ersit y of Colorado at Den v er, 1996. Submitted to SIAM J. Sci. Comput. [65] P V an ek, J. Mandel, and M. Brezina. Algebraic m ultigrid b y smoothed aggregation for second and fourth order elliptic problems. T o appear in Computing. [66] P etr V an ek and Jitk a K r zk o v a. Tw o{lev el method on unstructured meshes with con v ergence rate independen t of the coarse-space size. T ec hnical Report 35, Cen ter of Computational Mathematics/UCD, 1995. 112

PAGE 122

[67] P etr V an ek, Jan Mandel, Marian Brezina, and Radek T ezaur. Tw o-lev el substructuring tec hnique for problems with coecien t jumps within subdomains. In preparation, 1996. [68] P etr V an ek, Radek T ezaur, Marian Brezina, and Jitk a K r zk o v a. Tw o{ lev el method with coarse space size independen t con v ergence. In R. Glo winski, J. P eriaux, Z. Shi, and O. Widlund, editors, Domain Dec omp osition Metho ds in Scienc es and Engine ering John Wiley & Sons Ltd., New Y ork, N.Y., 1995. T o appear. [69] Olof B. Widlund. An extension theorem for nite elemen t spaces with three applications. In W olfgang Hac kbusc h and Kristian Witsc h, editors, Numeric al T e chniques in Continuum Me chanics pages 110{122, Braunsc h w eig/Wiesbaden, 1987. Notes on Numerical Fluid Mec hanics, v. 16, F riedr. View eg und Sohn. Proceedings of the Second GAMM-Seminar, Kiel, Jan uary 1986. [70] J. Xu. Iterativ e methods b y space decomposition and subspace correction: A unifying approac h. SIAM R eview 34:581{613, 1992. [71] Xuejun Zhang. Domain decomposition algorithms for the biharmonic Diric hlet problem. In T on y F. Chan, Da vid E. Key es, G erard A. Meuran t, Jerey S. Scroggs, and Robert G. V oigt, editors, Fifth International Symp osium on Domain De c omp osition Metho ds for Partial Dier ential Equations Philadelphia, P A, 1992. SIAM. 113