A DOMAIN DECOMPOSITION ALGORITHM FOR ELLIPTIC PROBLEMS IN THREE DIMENSIONS

BARRY F. SMITH^{∗}

Mathematics and Computer Science Division Preprint MCS-P185-1090 Argonne National Laboratory November 1990

Summary. Most domain decomposition algorithms have been developed for problems in two dimensions. One reason for this is the difficulty in devising a satisfactory, easy-to-implement, robust method of providing global communication of information for problems in three dimensions. Several methods that work well in two dimensions do not perform satisfactorily in three dimensions.

A new iterative substructuring algorithm for three dimensions is proposed. It is shown that the condition number of the resulting preconditioned problem is bounded independently of the number of subdomains and that the growth is quadratic in the logarithm of the number of degrees of freedom associated with a subdomain. The condition number is also bounded independently of the jumps in the coefficients of the differential equation between subdomains. The new algorithm also has more potential parallelism than the iterative substructuring methods previously proposed for problems in three dimensions.

Subject classifications: AMS(MOS): 65F10, 65N30

1. Introduction. In all domain decomposition algorithms with more than just a few subdomains, there is a crucial need for a mechanism to provide for global communication of information in each iteration. The convergence rate will decay rapidly with an increasing number of subdomains if communication is purely local, e.g. between neighboring subdomains only.

In two dimensions, global communication can be provided for by solving a problem that is essentially the discretization of the operator on the coarse grid defined by the subdomains. Such a method will be called a vertex-based algorithm. The coarse grid solution provides values at the vertices of the subdomains; values at all other nodes, for the coarse problem contribution to the solution, are obtained by interpolation. In three dimensions, this interpolation can produce a component of large energy, resulting in a

∗ Argonne National Laboratory, 9700 South Cass Ave, Argonne, IL 60439; formerly of the Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, NY 10012. This work was supported in part by the National Science Foundation under grant NSF-CCR-8903003 and by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S. Department of Energy, under Contract W-31-109-Eng-38.

poorly conditioned problem.

Considerations such as these motivated Bramble, Pasciak, and Schatz [3] to work with coarse problems that involve all the nodes of the wirebasket. The wirebasket consists of all nodes that belong to the closure of more than two subdomains; see Fig.

1. It is then possible to design an interpolation procedure that is much less harmful. Other wirebasket-based algorithms have been suggested by Dryja [5] and Mandel [7], [8]. A problem with these algorithms is that they require the solution of the coarse wirebasket problem before the solution of the local problems can commence. With the vertex-based algorithms, all the problems can be solved simultaneously.

We propose and analyze a new wirebasket-based method, which has the same superior convergence properties as the previous algorithms, but with the added feature that the local problems and the coarse problem can be solved simultaneously.

The new algorithm, like those of Bramble, Pasciak, and Schatz [3] and Mandel [7], [8], is an iterative substructuring algorithm. The preconditioner is constructed in such a manner that the condition number is independent of the number of substructures and, very important, is also independent of the jumps in the coefficients of the differential equation between subdomains. We note that the additive Schwarz algorithms considered by Dryja and Widlund [6], and Smith [9], [10], which use overlapping subdomains, converge very quickly for problems with slowly varying coefficients. However, their behavior for problems with large jumps in the coefficients is not fully understood.

This paper is organized as follows. In the next section, we introduce the elliptic and finite element problems. Section 3 contains a discussion of a vertex-based method that works well in two dimensions. Our new algorithm is introduced in Section 4 and is briefly contrasted with those of Bramble, Pasciak, and Schatz [3] and Mandel [7], [8]. Finally, in Section 5, we state and prove our theoretical results.

2. The Finite Element Problem. We consider a scalar, second-order, self adjoint, coercive, bilinear form a_{Ω}(u, v) on Ω ⊂ R^{3} and, for simplicity, impose a homogeneous Dirichlet condition on Γ_{0} ⊂ ∂Ω and Neumann boundary conditions on ∂Ω \ Γ_{0}. We assume the underlying elliptic operator has no zero order terms. The variational problem is then to find u ∈ H_{Γ_{0}}^{1}(Ω) such that,

a(u, v) = (f, v), ∀ v ∈ H_{Γ_{0}}^{1}(Ω).

We introduce two levels of triangulations into substructures Ω_{i} of diameter O(H) and into elements of diameter O(h). We assume shape regularity and that the substructures and elements satisfy the usual rules of finite element triangulations; see Ciarlet [4]. V^{H}(Ω) and V^{h}(Ω) are the spaces of continuous, piecewise linear functions, on the two triangulations, which vanish on Γ_{0}. Our work can be extended to other finite element models.

The discrete problem is then to find u^{h} ∈ V^{h}(Ω) such that

a(u^{h}, v^{h}) = (f, v^{h}), ∀ v^{h} ∈ V^{h}(Ω). (1)

If we expand u^{h} in the standard nodal basis, u^{h} = ∑_{k} x_{k}ϕ_{k}, the variational problem

[figure: for each node on W.; z_{W}: vector with a component of one; It is of the same length as z_{W} .; with ∂F^{j} are one; all others are zero.; z_{∂F^{j}}: vector whose components associated; z_{F^{j}}: vector whose components associated; with F^{j} are one; all others are zero.; nodes of F^{j}.; lies on F^{j} and is of value one on the; θ_{F^{j}}: finite element function whose support; lies on F^{j}.; u_{F^{j}}: finite element function whose support; Nodes on the wirebasket: W; Nodes on a face: F^{j}; Nodes on the boundary of a face: ∂F^{j}]

Fig. 1. Nodel subsets associated with a substructure Ω

(1) can be written as the linear system

Kx = f.

The elements of the stiffness matrix K are given by

K_{ij} = a(ϕ_{i}, ϕ_{j})

and those of the right-hand side f by

f_{i} = (f, ϕ_{i}).

The local contribution to the stiffness matrix and right-hand side can be formed one subdomain at a time. The stiffness matrix is then obtained by subassembly of these parts. If we order the nodes interior to the subdomains first followed by the nodes on the subdomain’s boundaries, we can write the linear system symbolically as

∑_{i} ( \begin{matrix}K_{II}^{(i)} & K_{IB}^{(i)} \\ K_{IB}^{(i)^{T}} & K_{BB}^{(i)} \end{matrix} ) ( \begin{matrix}x_{I}^{(i)} \\ x_{B}^{(i)} \end{matrix} ) = ∑_{i} ( \begin{matrix}f_{I}^{(i)} \\ f_{B}^{(i)} \end{matrix} ) .

In the first step of most iterative substructuring algorithms, the unknowns in the interior of the subdomains are eliminated. Elimination is done by calculating the Schur complement with respect to the variables associated with the boundaries of the substructures. The resulting linear system is

∑_{i} ( \begin{matrix}K_{II}^{(i)} & K_{IB}^{(i)} \\ 0 & S_{BB}^{(i)} \end{matrix} ) ( \begin{matrix}x_{I}^{(i)} \\ x_{B}^{(i)} \end{matrix} ) = ∑_{i} ( \begin{matrix}f_{I}^{(i)} \\ f_{B}^{(i)} − K_{IB}^{(i)^{T}} K_{II}^{(i)^{−1}}f_{B}^{(i)} \end{matrix} ) ,

where

S^{(i)} = S_{BB}^{(i)} = K_{BB}^{(i)} − K_{IB}^{(i)^{T}} K_{II}^{(i)^{−1}}K_{IB}^{(i)}.

We write the reduced system as

Sx = f. (2)

The matrix S is obtained, from the S^{(i)}, by subassembly.

An iterative substructuring algorithm is obtained by solving this reduced, Schur complement system using a preconditioned iterative method. Many of these algorithms do not require that the Schur complement be formed, since it needs to be available only in terms of a matrix-vector product. This product can be calculated, using the definition of S^{(i)}, at the expense of solving a Dirichlet problem on each subdomain. The new algorithm can be formulated in such a way as not to require the explicit calculation of S.

3. The Vertex-based Method. The vertex-based method can be understood as a block diagonal preconditioner in a new basis. We give the description for problems in three dimensions; the same approach also works for two dimensions. We are solving the reduced, Schur complement problem, i.e. equation (2), and therefore the interior variables x_{I} play no role in the description. Let Ṽ^{h} be the subspace of functions in V^{h} that are discrete harmonic in the interiors of the subdomains. Such functions satisfy the relation K_{II}x_{I} + K_{IB}x_{B} = 0 and are completely defined by their values, x_{B}, on the interfaces between the substructures. We introduce basis functions ϕ̃^{k} of Ṽ^{h} which take on the same values as ϕ^{k} on the interfaces between substructures and are extended as discrete harmonic in the interiors of the subdomains.

We first consider the substructures, which have diameters on the order of H, as large elements and consider the corresponding nodal functions, {ψ^{j}}, as a basis for the space V^{H}. These basis functions are discrete harmonic, and hence V^{H} is a subspace of Ṽ^{h}. We decompose the space Ṽ^{h} into

Ṽ^{h} = V^{H} ⊕ (Ṽ^{h} \ V^{H}).

Any function in Ṽ^{h} can then be represented in the basis {ϕ̃^{k}} but also in the new partial hierarchical basis. We group the nodal basis functions into two groups, those associated with a vertex of a substructure, {ϕ̃_{V}^{k} }, and all the rest,{ϕ̃_{E}^{k}}, and represent u^{h} by

u^{h} = ∑_{j} y_{V_{j}}ψ^{j} + ∑_{k} y_{E_{k}} ϕ̃_{E}^{k},

and

u^{h} = ∑_{j} x_{V_{j}} ϕ̃_{V}^{j} + ∑_{k} x_{E_{k}} ϕ̃_{E}^{k}.

For a single substructure the mapping between the two sets of coefficients is given by

( \begin{matrix}y_{E}^{(i)} \\ y_{V}^{(i)} \end{matrix} ) = ( \begin{matrix}I & −R^{(i)^{T}} \\ 0 & I \end{matrix} ) ( \begin{matrix}x_{E}^{(i)} \\ x_{V}^{(i)} \end{matrix} ) , ( \begin{matrix}x_{E}^{(i)} \\ x_{V}^{(i)} \end{matrix} ) = ( \begin{matrix}I & R^{(i)^{T}} \\ 0 & I \end{matrix} ) ( \begin{matrix}y_{E}^{(i)} \\ y_{V}^{(i)} \end{matrix} ) .

R^{(i)^{T}} simply represents linear interpolation from the space Ṽ^{h} to V^{H}. The values of the interpolant on an edge or face depend only on the values at the subdomain vertices that belong to the boundary of that edge or face. Therefore, we can express the global basis change as

( \begin{matrix}y_{E} \\ y_{V} \end{matrix} ) = ( \begin{matrix}I & −R^{T} \\ 0 & I \end{matrix} ) ( \begin{matrix}x_{E} \\ x_{V} \end{matrix} ) , ( \begin{matrix}x_{E} \\ x_{V} \end{matrix} ) = ( \begin{matrix}I & R^{T} \\ 0 & I \end{matrix} ) ( \begin{matrix}y_{E} \\ y_{V} \end{matrix} ) .

The elements of R^{T} are not obtained by subassembly, but rather as the common value of the corresponding elements in the appropriate matrix R^{(i)^{T}} .

The preconditioner for the vertex-based methods is obtained by first making this change to a partial hierarchical basis, and then using a simple block diagonal splitting. Let S^{(i)} be the contribution to the Schur complement from the ith subdomain,

S^{(i)} = ( \begin{matrix}S_{EE}^{(i)} & S_{EV}^{(i)} \\ S_{EV}^{(i)^{T}} & S_{V V}^{(i)} \end{matrix} ) .

We make a partial change to hierarchical basis, by post-multiplying by

Q^{(i)} = ( \begin{matrix}I & R^{(i)^{T}} \\ 0 & I \end{matrix} )

and by pre-multiplying by Q^{(i)^{T}} . This gives

Q^{(i)^{T}} S^{(i)}Q^{(i)} = ( \begin{matrix}S_{EE}^{(i)} & Non-zero \\ Non-zero & S̃_{V V}^{(i)} \end{matrix} ) .

We now replace S_{EE}^{(i)} with its block diagonal part with a block for each edge and each face, and we also drop the coupling between the vertices and the other nodes. This results in

( \begin{matrix}Ŝ_{EE}^{(i)} & 0 \\ 0 & S̃_{V V}^{(i)} \end{matrix} ) .

Other iterative substructuring algorithms can be obtained by replacing the diagonal blocks with spectrally equivalent matrices which might be more computationally efficient to work with.

Finally, we return to the usual nodal basis by post-multiplying by Q^{(i)^{−1}} and premultiplying by Q^{(i)^{−T}} .

Ŝ^{(i)} = ( \begin{matrix}I & 0 \\ −R^{(i)} & I \end{matrix} ) ( \begin{matrix}Ŝ_{EE}^{(i)} & 0 \\ 0 & S̃_{V V}^{(i)} \end{matrix} ) ( \begin{matrix}I & −R^{(i)^{T}} \\ 0 & I \end{matrix} ) . (3)

We construct the global preconditioner by subassembly,

Ŝ = ( \begin{matrix}I & 0 \\ −R & I \end{matrix} ) ( \begin{matrix}Ŝ_{EE} & 0 \\ 0 & S̃_{V V} \end{matrix} ) ( \begin{matrix}I & −R^{T} \\ 0 & I \end{matrix} ) . (4)

The preconditioner is easily inverted to give

Ŝ^{−1} = ( \begin{matrix}R^{T} \\ I \end{matrix} ) S̃_{V V}^{−1} ( R I ) + ( \begin{matrix}I \\ 0 \end{matrix} ) Ŝ_{F F}^{−1} ( I 0 ) .

Note that there are independent problems associated with each edge, each face, and the coarse problem.

We return briefly to the local representation of the preconditioner. We note that the local contribution to the preconditioner, Ŝ^{(i)}, has the same null space as S^{(i)}. This leads to the very useful observation that it is possible to bound the condition number of the preconditioned problem by bounds obtained on individual substructures. This important idea is used in Bramble, Pasciak, and Schatz [2], [3] and presented very clearly in Mandel [7], [8]. This also makes it possible to show that the condition number is independent of the jumps in the coefficients between subdomains. Thus, if

c_{i}Ŝ^{(i)} ≤ S^{(i)} ≤ C_{i}Ŝ^{(i)}, ∀i,

then

min_{i} c_{i}Ŝ ≤ S ≤ max_{i} C_{i}Ŝ,

or equivalently

κ(Ŝ^{−1}S) ≤ \frac{max_{i} C_{i}}{min_{i} c_{i}} .

We note that this analysis will work if and only if

null(Ŝ^{(i)}) = null(S^{(i)}), ∀i.

The observation above appears as a theorem in Mandel [7].

For scalar, second-order, elliptic problems without zero-order terms, the subdomains that have no given Dirichlet boundary data have a null space of the constant functions. Those with any given Dirichlet boundary data do not have a null space. Using the above argument, we can conclude that the vertex-based method has a condition number that is bounded independently of the number of subdomains.

Dryja and Widlund [6] have shown that, in two dimensions, the vertex-based method satisfies

κ ≤ C(1 + log(H/h))^{2}.

In three dimensions the bound is given by

κ ≤ C(H/h)(1 + log(H/h)).

Throughout this paper, c and C are constants, independent of h and H. These bounds are sharp. The result for two dimensions is quite satisfactory, while that for three dimensions is disappointing. In the next section we discuss why there is so much difference and show how this understanding can lead to better algorithms in three dimensions.

4. The Wirebasket-based Algorithms. The following Sobolev type inequality (see Bramble [1]) holds for finite element functions in two dimensions:

||u^{h} − α||_{L^{∞}(Ω_{i})}^{2} ≤ C(1 + log(H/h))|u^{h}|_{H^{1}(Ω_{i})}^{2}.

Here α is any convex combination of the values of u^{h} in Ω_{i}. Using this inequality, we can show that the energy of the coarse mesh interpolant of u^{h} can exceed that of u^{h} by at most a factor of C(1 + log(H/h)).

However, in three dimensions we have a much weaker bound

||u^{h} − α||_{L^{∞}(Ω_{i})}^{2} ≤ C(H/h)|u^{h}|_{H^{1}(Ω_{i})}^{2}.

Therefore, in this case, interpolating the value of u^{h} from a vertex can result in a O(H/h) change in the energy. See Smith [10] for a detailed discussion of this point.

Bramble, Pasciak, and Schatz [3] observed that, for three dimensions, there exists an alternative bound. Let α be any convex combination of values of u^{h} on substructure Ω_{i}. Then

h||u^{h} − α||_{l^{2}(W_{i})}^{2} ≤ C(1 + log(H/h))|u^{h}|_{H^{1}(Ω_{i})}^{2}. (5)

Here W_{i} represents the wirebasket on substructure Ω_{i}. The wirebasket-based schemes of Bramble, Pasciak, and Schatz [3] and Dryja [5] are based on calculating the average value of u^{h} on Ω_{i} and then interpolating these values onto the faces of the individual substructures.

The use of the average value of u^{h} over the substructure makes it impossible to solve the coarse problem in parallel with the local problems. We can see this difficulty by noting that the local contribution to the preconditioner can be written in the form of (3). However, the action of the two transformation matrices R^{(i)^{T}} generally provides different values on the face, since each R^{(i)^{T}} also depends on values that do not belong to both substructures. We therefore cannot express the global form of the preconditioner as in (4), and the coarse problem and face problems cannot be solved in parallel.

4.1. The New Method. To introduce the new algorithm, we observe that the estimate (5) remains valid if we replace the entire wirebasket by a line segment of length cH. In the new algorithm we will interpolate averages for the parts of the wirebasket adjacent to each face; see Fig. 1.

We construct the preconditioner locally, using a method similar to that of Section 3, one substructure at a time, and obtain the preconditioner by subassembly. We first order the nodes on the faces and then those on the wirebasket. The local contribution of substructure Ω_{i} to the Schur complement is given by

S^{(i)} = ( \begin{matrix}S_{F F}^{(i)} & S_{F W}^{(i)} \\ S_{F W}^{(i)^{T}} & S_{W W}^{(i)} \end{matrix} ) .

Let T^{(i)^{T}} map the weighted average of the values of the boundary nodes of each face (the adjacent vertices and edges) to the nodes on the corresponding face; see Fig. 1. Let Ŝ_{F F}^{(i)} be the block diagonal part of S_{F F}^{(i)} with a block for each face. The local contribution to the preconditioner is then

Ŝ^{(i)} = ( \begin{matrix}I & 0 \\ −T^{(i)} & I \end{matrix} ) ( \begin{matrix}Ŝ_{F F}^{(i)} & 0 \\ 0 & Ĝ^{(i)} \end{matrix} ) ( \begin{matrix}I & −T^{(i)^{T}} \\ 0 & I \end{matrix} ) .

Note that the structure is similar to that of the preconditioner introduced in Section 3. We define Ĝ^{(i)} by

x_{W}^{(i)^{T}} Ĝ^{(i)}x_{W}^{(i)} = min_{w̄^{(i)}}(x_{W}^{(i)} − w̄^{(i)}z_{W}^{(i)})^{T} G^{(i)}(x_{W}^{(i)} − w̄^{(i)}z_{W}^{(i)}). (6)

G^{(i)} is given by δ(H/h)I or δ(H/h) times a constant (block) diagonal matrix and z_{W}^{(i)} is a vector of all ones of the same dimension as x_{W}^{(i)}; see Fig. 1. The optimizing w̄^{(i)} corresponds to the weighted average of the values of x on the wirebasket of substructure Ω_{i}. The Ĝ^{(i)} is constructed in this manner to force Ĝ^{(i)}, and hence Ŝ^{(i)}, to have a null space consisting of the constants. Since, for interior subdomains, this is also the null space for S^{(i)}, we immediately obtain that the condition number is independent of the number of subdomains (see Section 3). For subdomains with any prescribed Dirichlet boundary data, the construction of the preconditioner is identical, except that the value of u^{h} at all the Dirichlet boundary nodes is constrained to be zero. In this case, the S^{(i)} and Ŝ^{(i)} do not have null spaces.

The δ(H/h), a scalar function of (H/h), is chosen to optimize the scaling between the ‘coarse’ problem and the ‘face’ problems. Our results show that when δ(H/h) = (1 + log(H/h)), the condition number of the preconditioned problem is bounded by C(1 + log(H/h))^{2} in three dimensions. If δ(H/h) is chosen to be a constant, then the bound is C(1 + log(H/h))^{3}.

The global preconditioner is obtained by subassembly:

Ŝ = ( \begin{matrix}I & 0 \\ −T & I \end{matrix} ) ( \begin{matrix}Ŝ_{F F} & 0 \\ 0 & Ĝ \end{matrix} ) ( \begin{matrix}I & −T^{T} \\ 0 & I \end{matrix} ) .

Here Ŝ_{F F} and Ĝ are obtained by subassembly from Ŝ_{F F}^{(i)} and Ĝ^{(i)}, respectively. As in Section 3, we have used the fact that the actions of the T^{(i)^{T}} for two adjacent subdomains along a shared face are identical. The preconditioner is easily inverted to give

Ŝ^{−1} = ( \begin{matrix}T^{T} \\ I \end{matrix} ) Ĝ^{−1} ( T I ) + ( \begin{matrix}I \\ 0 \end{matrix} ) Ŝ_{F F}^{−1} ( I 0 ) .

Note that the resulting operator has independent parts associated with each face and the coarse problem. This is not true of the methods considered by Bramble, Pasciak, and Schatz [3] and Mandel [7], [8].

To derive an explicit formula for T^{(i)}, we need to introduce some more notation. Let z_{F^{j}}^{(i)} be a vector of the same length as x_{F}^{(i)} with zeros at all nodes except those associated with the face F^{j} where the coefficients are one. Let z_{∂F^{j}}^{(i)} be a vector of the same length as x_{W}^{(i)}, which is zero at all nodes except those on the boundary of the face F^{j} where the coefficients are one (see Fig. 1). Then

T^{(i)^{T}} = ∑_{j} z_{F^{j}}^{(i)}(z_{∂F^{j}}^{(i)^{T}} G^{(i)}z_{∂F^{j}}^{(i)})^{−1}z_{∂F^{j}}^{(i)^{T}} G^{(i)}.

4.2. The Previous Wirebasket Method. The original iterative substructuring wirebasket-based algorithm was proposed by Bramble, Pasciak, and Schatz [3]. Dryja [5] adapted the same type of wirebasket-based coarse problem to a Neumann-Dirichlet algorithm. In [7] and [8] Mandel introduced a more general presentation and interpretation of the Bramble, Pasciak, and Schatz algorithm, and we adopt this approach here. Motivated by the desire to obtain condition numbers independent of the number of substructures, Mandel proposed the following preconditioner. Let z^{(i)} be the vector of all ones of the same length as x^{(i)}. The local contribution to the preconditioner is written as

x^{(i)^{T}} Ŝ^{(i)}x^{(i)} = min_{w̄^{(i)}}(x^{(i)} − w̄^{(i)}z^{(i)})^{T} ( \begin{matrix}Ŝ_{F F}^{(i)} & 0 \\ 0 & G^{(i)} \end{matrix} ) (x^{(i)} − w̄^{(i)}z^{(i)}).

Here G^{(i)} is given by

I  or  diag(S_{W W}^{(i)} )

or a block diagonal part of S_{W W}^{(i)} . We no longer use a scale factor δ(H/h).

Remark: In Bramble, Pasciak, and Schatz [3] the problems associated with the faces are solved using a fast, approximate technique, which does not require the explicit calculation of S_{F F}^{(i)} . If we restrict our substructures, as in [3], to be brick-shaped, a similar approach can be applied for this new algorithm as well. We would, in the preconditioner, replace each block associated with a face with the l_{0}^{1/2} operator considered in Bramble, Pasciak, and Schatz [3].

4.3. Solving the Coarse Problem. For all of these methods based on averaging, we need to solve a quadratic problem of the form

min_{x} ∑_{i} \frac{1}{2}(x^{(i)} − w̄^{(i)}z^{(i)})^{T} B^{(i)}(x^{(i)} − w̄^{(i)}z^{(i)}) − x^{T} f,

with w̄^{(i)} defined by

arg min_{v̄}(x^{(i)} − v̄z^{(i)})^{T} B^{(i)}(x^{(i)} − v̄z^{(i)}).

For the new method this is the system associated with the matrix Ĝ. For the other wirebasket-based methods, all the variables associated with the preconditioner are involved.

We use the solution technique developed by Mandel [7],[8]. Bramble, Pasciak, and Schatz [3] used different tools and provided a technique that appears less generally applicable.

We write the problem as

min_{x} ∑_{i} min_{w̄^{(i)}} \frac{1}{2}(x^{(i)} − w̄^{(i)}z^{(i)})^{T} B^{(i)}(x^{(i)} − w̄^{(i)}z^{(i)}) − x^{T} f.

We then take derivatives with respect to w̄^{(i)} and x. We obtain the linear system

z^{(i)^{T}} B^{(i)}(x^{(i)} − z^{(i)}w̄^{(i)}) = 0, ∀i,

Bx − ∑_{i} B^{(i)}z^{(i)}w̄^{(i)} = f. (7)

We then eliminate x and get the following system for the w̄^{(i)} :

(z^{(i)^{T}} B^{(i)}z^{(i)})w̄^{(i)} − z^{(i)^{T}} B^{(i)}B^{−1} ∑_{j} B^{(j)}z^{(j)}w̄^{(j)} = z^{(i)^{T}} B^{(i)}B^{−1}f.

Once the w̄^{(i)} are known, x can be found by solving (7).

5. Proofs of the Near Optimality of the Methods. In this section we prove the convergence properties of the new method and Mandel’s method. We first give a series of lemmas needed, then state and prove our results. Many of the techniques are similar to those of Bramble, Pasciak, and Schatz [3] and Dryja [5]. However, our algorithm is quite different, as is some of the analysis.

We first use the approach of Section 3 to reduce the calculation of the condition number to that of a single subdomain.

Lemma 5.1. The vertex-based method, the new method and the other wirebasketbased method satisfy

κ(Ŝ^{−1}S) ≤ \frac{max_{i} λ_{max}(Ŝ^{(i)^{+}}S^{(i)})}{min_{i} λ_{min}(Ŝ^{(i)^{+}}S^{(i)})}.

The starting point for the analysis of a large body of domain decomposition algorithms is the lemma.

Lemma 5.2. In two dimensions

||u^{h}||_{L^{∞}(Ω̄)}^{2} ≤ C(1 + log(H/h))||u^{h}||_{H^{1}(Ω)}^{2}.

Proof. See, e.g., Bramble, Pasciak, and Schatz [2] or Dryja [5]. □ In the rest of this section Ω ⊂ R^{3} is a substructure of diameter H, and Γ is its

boundary. We use the weighted H^{1} norm

||u||_{H^{1}(Ω)}^{2} = |u|_{H^{1}(Ω)}^{2} + \frac{1}{H^{2}} ||u||_{L^{2}(Ω)}^{2},

and the weighted H^{1/2} norm

||u||_{H^{1/2}(Γ)}^{2} = |u|_{H^{1/2}(Γ)}^{2} + \frac{1}{H} ||u||_{L^{2}(Γ)}^{2}.

The discrete l^{2} norm along an interval I is defined simply by

||u^{h}||_{l^{2}}^{2} = ∑_{j} u^{h}(x_{j})^{2},

where the x_{j} are the nodes along the interval I.

The following lemma is a consequence of Lemma 5.2 and is obtained by applying it to two-dimensional slices one at a time. Proofs of the first and second parts are given in Bramble, Pasciak, and Schatz [3] and Dryja [5], respectively.

Lemma 5.3. Let I be a line segment in Ω, of length cH. Then

h||u^{h}||_{l^{2}(I)}^{2} ≤ C(1 + log(H/h))||u^{h}||_{H^{1}(Ω)}^{2}.

If, in addition, ū^{h} is the average of u^{h} on I, then

h||u^{h} − ū^{h}||_{l^{2}(I)}^{2} ≤ C(1 + log(H/h))|u^{h}|_{H^{1}(Ω)}^{2}.

Lemma 5.4. Let ū^{h} be the average of u^{h} over an interval I, of length cH. Then

(ū^{h})^{2} ≤ C \frac{(1 + log(H/h))}{H}||u^{h}||_{H^{1}(Ω)}^{2}.

In addition,

(ū^{h})^{2} ≤ C \frac{h}{H} ||u^{h}||_{l^{2}(I)}^{2}.

Proof.

(ū^{h})^{2} ≤ C \frac{h^{2}}{H^{2}} (∑_{i∈I} u^{h}(x_{i}))^{2}
≤ C \frac{h^{2}}{H^{2}} h^{−1}H ∑_{i∈I} u^{h}(x_{i})^{2}
= C \frac{h}{H} ||u^{h}||_{l^{2}(I)}^{2}
≤ C \frac{(1 + log(H/h))}{H}||u^{h}||_{H^{1}(Ω)}^{2}. □

Let Γ be the boundary of Ω, and let F^{i} ⊂ Γ be a face. Let θ_{F^{i}} be the finite element function that is one at the interior nodes of F^{i} and vanishes at all of the other nodes on Γ.

Lemma 5.5. For θ_{F^{i}} as defined above,

|θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ C(1 + log(H/h))H.

Proof. This is an intermediate result in the proof of Lemma 5.6. □ Lemma 5.6. Let u_{F^{i}}^{h} be that finite element function that is equal to u^{h} on the

interior nodes of F^{i} but zero on ∂F^{i}. Then

|u_{F^{i}}^{h}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ C(1 + log(H/h))^{2}||u^{h}||_{H^{1/2}(Γ)}^{2}.

Proof. This is Lemma 4.3 in Bramble, Pasciak, and Schatz [3]. □ Lemma 5.7. Let I be a line segment in Ω, of length cH, and let u_{W}^{h} be the finite

element function that equals u^{h} on I but that vanishes at all nodes not on I. Then

|u_{W}^{h} |_{H^{1/2}(Γ)}^{2} ≤ Ch||u^{h}||_{l^{2}(I)}^{2}.

Proof. This is Lemma 4.1 in Bramble, Pasciak, and Schatz [3]. □ The following extension theorem is given in Bramble, Pasciak, and Schatz [3]. Lemma 5.8. Let u^{h} be a discrete harmonic, piecewise linear, finite element func-

tion in Ω. Then there exists a C independent of u^{h},h and H such that

|u^{h}|_{H^{1}(Ω)} ≤ C|u^{h}|_{H^{1/2}(Γ)}.

Lemma 5.9. Let ū^{h} be the average of u^{h} on an interval I, of length cH. Then

|u_{F^{i}}^{h} − ū^{h}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ C(1 + log(H/h))^{2}|u^{h}|_{H^{1/2}(Γ)}^{2}.

Proof. We use Lemma 5.6, Lemma 5.4, and Lemma 5.5 to obtain

|u_{F^{i}}^{h} − ū^{h}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ 2|u_{F^{i}}^{h}|_{H_{00}^{1/2}(F^{i})}^{2} + 2|ū^{h}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2}
≤ C(1 + log(H/h))^{2}||u^{h}||_{H^{1}(Ω)}^{2}.

To conclude the argument, we note that the left-hand side is unchanged if a constant is added to u^{h}. We therefore shift u^{h} by a constant so that its average is zero. We can then replace the norm on the right-hand side with a semi-norm using a Poincaré inequality. We return to the H^{1/2} semi-norm using Lemma 5.8. □

The main result for the new algorithm is given in the following theorem.

Theorem 5.1. For the new wirebasket method using G^{(i)} = (1 + log(H/h))I, the condition number is bounded by

κ(Ŝ^{−1}S) ≤ C(1 + log(H/h))^{2}.

Moreover, if G^{(i)} = I, then

κ(Ŝ^{−1}S) ≤ C(1 + log(H/h))^{3}.

Proof. We prove the first result; the second result follows easily. We will establish

\frac{x^{(i)^{T}} Ŝ^{(i)}x^{(i)}}{c(1 + log(H/h))^{2}} ≤ x^{(i)^{T}} S^{(i)}x^{(i)} ≤ Cx^{(i)^{T}} Ŝ^{(i)}x^{(i)}. (8)

The proofs for interior subdomains and boundary subdomains are nearly identical. The only difference is that, for the boundary subdomains, the u^{h} are constrained to be zero on all nodes with prescribed Dirichlet data and hence lie in a subspace of Ṽ^{h}. We prove the result for all u^{h} ∈ Ṽ^{h}.

We need the following definitions. Let w̄ be the average value of u^{h} on the nodes of the wirebasket and w̄^{∂F^{i}} be the average value of u^{h} on ∂F^{i}. Recall that u_{F^{i}}^{h} is the finite element function that is equal to u^{h} at the interior nodes of face F^{i}, but is zero on ∂F^{i}. W = ∪∂F^{i} is the union of the line segments that make up the wirebasket, and u_{W}^{h} is the finite element function that is equal to u^{h} at the nodes of W, but is zero at all other nodes of Γ. Using this notation the lower bound of (8) is equivalent to

δ(H/h)h||u^{h} − w̄||_{l^{2}(W)}^{2} + ∑_{i} |u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ c(1 + log(H/h))^{2}|u^{h}|_{H^{1/2}(Γ)}^{2}.

We bound each of the left-hand terms separately. We first note that

δ(H/h)h||u^{h} − w̄||_{l^{2}(W)}^{2} ≤ cδ(H/h)(1 + log(H/h))|u^{h}|_{H^{1}(Ω)}^{2}
≤ c(1 + log(H/h))^{2}|u^{h}|_{H^{1}(Ω)}^{2}
≤ c(1 + log(H/h))^{2}|u^{h}|_{H^{1/2}(Γ)}^{2}.

This follows from Lemma 5.3 and Lemma 5.8. We bound the other terms using Lemma 5.9,

|u_{F}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ c(1 + log(H/h))^{2}|u^{h}|_{H^{1/2}(Γ)}^{2}.

We must now prove the upper bound,

x^{(i)^{T}} S^{(i)}x^{(i)} ≤ Cx^{(i)^{T}} Ŝ^{(i)}x^{(i)}.

We shift u^{h}, by a constant, so that the average of u^{h} on the wirebasket is zero, i.e., w̄ = 0. The upper bound is then equivalent to

|u^{h}|_{H^{1/2}(Γ)}^{2} ≤ C(δ(H/h)h||u^{h}||_{l^{2}(W)}^{2} + ∑_{i} |u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2}).

We establish this bound by

|u^{h}|_{H^{1/2}(Γ)}^{2} = | ∑_{i}(u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}) + u_{W}^{h} + ∑_{i} w̄^{∂F^{i}}θ_{F^{i}}|_{H^{1/2}(Γ)}^{2}
≤ C(∑_{i} |u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H^{1/2}(Γ)}^{2} + |u_{W}^{h} |_{H^{1/2}(Γ)}^{2} + ∑_{i} |w̄^{∂F^{i}}θ_{F^{i}}|_{H^{1/2}(Γ)}^{2})
≤ C(∑_{i} |u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} + h||u^{h}||_{l^{2}(W)}^{2} + (1 + log(H/h))h||u^{h}||_{l^{2}(W)}^{2})
≤ C(∑_{i} |u_{F^{i}}^{h} − w̄^{∂F^{i}}θ_{F^{i}}|_{H_{00}^{1/2}(F^{i})}^{2} + δ(H/h)h||u^{h}||_{l^{2}(W)}^{2}).

The bound follows from Lemma 5.7, Lemma 5.4, and Lemma 5.5. □

This trick of scaling a portion of the preconditioner with δ(H/h) has been previously used by Dryja [5] for a version of the Neumann-Dirichlet algorithm.

We conclude by proving the result for the previous wirebasket method. Note that this is essentially the result in Bramble, Pasciak, and Schatz [3].

Theorem 5.2. For Mandel’s version of the wirebasket algorithm, using G^{(i)} = I, the condition number is bounded by

κ(Ŝ^{−1}S) ≤ C(1 + log(H/h))^{2}.

Proof. We will prove that

\frac{x^{(i)^{T}} Ŝ^{(i)}x^{(i)}}{c(1 + log(H/h))^{2}} ≤ x^{(i)^{T}} S^{(i)}x^{(i)} ≤ Cx^{(i)^{T}} Ŝ^{(i)}x^{(i)}.

We shift u^{h} by a constant so that its weighted average is zero. If we use the definitions given in the previous theorem, the lower bound follows from

h||u^{h}||_{l^{2}(W)}^{2} + ∑_{i} |u_{F^{i}}^{h}|_{H_{00}^{1/2}(F^{i})}^{2} ≤ c(1 + log(H/h))^{2}|u^{h}|_{H^{1/2}(Γ)}^{2}.

This inequality follows from Lemma 5.3, Lemma 5.6, Lemma 5.8, and a quotient space argument.

Lemma 5.7 implies the upper bound:

|u^{h}|_{H^{1/2}(Γ)}^{2} = | ∑_{i} u_{F^{i}}^{h} + u_{W}^{h} |_{H^{1/2}(Γ)}^{2}
≤ C(∑_{i} |u_{F^{i}}^{h}|_{H^{1/2}(Γ)}^{2} + |u_{W}^{h} |_{H^{1/2}(Γ)}^{2})
≤ C(∑_{i} |u_{F^{i}}^{h}|_{H_{00}^{1/2}(F^{i})}^{2} + h||u^{h}||_{l^{2}(W)}^{2}). □

Remark: The bound on the condition number for both wirebasket-based algorithms is essentially the same. However, the analysis gives no indication how the condition numbers for the two methods compare. This can only be decided by numerical computations.

Acknowledgments. I would like to thank my adviser Olof Widlund for the assistance he has given me in revising this paper. This work is based on Chapter 5, Section 2, of my Ph.D. thesis [10].

REFERENCES

[1] J. H. Bramble, A second order finite difference analogue of the first biharmonic boundary value problem, Numer. Math., 9 (1966), pp. 236–249.

[2] J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring, I, Math. Comp., 47 (1986), pp. 103–134.

[3] , The construction of preconditioners for elliptic problems by substructuring, IV, Math. Comp., 53 (1989), pp. 1–24.

[4] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.

[5] M. Dryja, A method of domain decomposition for 3-D finite element problems, in First International Symposium on Domain Decomposition Methods for Partial Differential Equations, R. Glowinski, G. H. Golub, G. A. Meurant, and J. Périaux, eds., Society for Industrial and Applied Mathematics, Philadelphia, PA, 1988, pp. 43–61.

[6] M. Dryja and O. Widlund, Some domain decomposition algorithms for elliptic problems, in Iterative Methods for Large Linear Systems, Academic Press, San Diego, CA, 1989, pp. 273–291.

[7] J. Mandel, Iterative solvers by substructuring for the p-version finite elment method, Comput. Methods Appl. Mech. Engrg., 80 (1990), pp. 117-128.

[8] , Two-level domain decomposition preconditioners for the p-version finite element method in three dimensions, Int. J. Numer. Methods Engrg., 29 (1990), pp. 1095-1108.

[9] B. F. Smith, An optimal domain decomposition preconditioner for the finite element solution of linear elasticity problems, Tech. Report 482, Department of Computer Science, Courant Institute, New York, NY, 1989; to appear in SIAM J. Sci. Stat. Comput.

[10] , Domain Decomposition Algorithms for the Partial Differential Equations of Linear Elasticity, Ph.D. Thesis, Courant Institute, New York, NY, 1990.