In a similar way, the variational method can mimic patching if the variational principle

is slightly modi¬ed. Again, however, the modi¬cations do not raise or lower the order

of the method. There is no compelling reason for choosing patching over the variational

approach or vice versa; it is a matter of preference and taste.

22.8 Matrix Inversion

As reviewed by Canuto et al. (1988), linear algebra methods used with spectral elements

include:

(i) Alternating Schwarz method

(ii) Static condensation

(iii) In¬‚uence matrix method

(iv) Conjugate gradient and

(v) Multigrid.

The Schwarz algorithm de¬nes overlapping elements and then iterates to match the solu-

tions in different subdomains. Canuto & Funaro (1988) report great success and Schwarz

iteration has become popular with ¬nite difference methods, too.

“Static condensation” is a two-stage application of Gaussian elimination. On a given

subdomain, the numerical solution has both internal and external degrees of freedom. The

external degrees of freedom are those which can be computed only simultaneously with

the corresponding quantities for all the other elements. For the second order ordinary

differential equation of Sec. 4, for example, the general solution can be written as

x ∈ [dj’1 , dj ] (22.28)

uj (x) = pj (x) + Uj’1 hLj (x) + Uj hRj (x)

where Uj ≡ u(xj ) and the particular solution pj (x) and the homogeneous solutions, hLj (x)

& hRj (x), satisfy (22.12) and (22.13). The internal degrees of freedom are the polynomial

coef¬cients of pj (x), hLj (x) & hRj (x). These coef¬cients are “internal” because one can

compute these three functions on the j-th element without using any information outside

the j-th element. This is possible because these three functions are arbitrarily assigned

boundary values of either 1 or 0.

The external degrees of freedom are the values of uj (x) on the boundaries of the ele-

ment, Uj . These can be determined only by matching the elemental solutions across the

global domain to form a single continuous solution.

The ¬rst stage of static condensation is to eliminate the internal degrees of freedom. The

second stage is to solve a reduced matrix system (i. e. (22.16)) for the external degrees of

freedom. The dimension of this “condensed” system is small in comparison to N , the total

number of degrees of freedom, resulting in huge savings.

The in¬‚uence matrix method, which is the theme of the next section, is very similar to

static condensation. Both are direct methods that perform as much of the work as possible

within each element before solving a matrix problem of moderate size to link the piecewise

polynomials together.

The conjugate gradient and multigrid algorithms are iterations. However, both share

a common philosophy with the direct methods: exploit the weak inter-element coupling

through heavy use of intra-element operations. As emphasized in Sec. 4, this strategy is

fundamental to all spectral element procedures.

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

488

Both multigrid and a conjugate gradient-like algorithm have already been discussed

in an earlier chapter. However, Rønquist & Patera (1987b) and Maday & Munoz (1988)

have shown that the incomplete-factorizations employed with global Chebyshev multigrid

are unnecessary for Legendre spectral elements. Instead, a simple pointwise Richardson™s

iteration (which these authors label with Jacobi™s name) is extremely effective, at least for

simple problems.

It is usual to abandon iteration when multigrid has descended to a grid with only a

small number of points and solve the matrix problem on the coarsest grid via Gaussian

elimination. In Rønquist & Patera (1987b), the multigrid/Richardson™s iteration polyalgo-

rithm is applied only within each element. The direct-solution-on-the-coarsest-grid is the

step at which the various elemental solutions are matched together. Thus, the non-iterative

part of multigrid is identical with the second stage of static condensation or the in¬‚uence

matrix algorithm described below.

22.9 The In¬‚uence Matrix Method

If we knew the values of u(x, y) on the walls of an element, then the problem on that

subdomain would be completely self-contained. If we had a subroutine that solved the

differential equation using a global spectral series, we could call this subprogram once and

be done.

The dif¬culty is that when one or more of the walls of an element is an internal bound-

ary, the values of u(x, y) on that wall are unknown. Nevertheless, if we discretize the bound-

ary, we can calculate the in¬‚uence of each boundary point on the rest of the solution. This

calculation of “in¬‚uence” requires nothing more nor less than calling the subroutine to

solve the Dirichlet problem on the j-th element.

Let the elliptic equation be

(22.29)

Lu = f

For simplicity, assume homogeneous Dirichlet conditions. (If the boundary conditions

are inhomogeneous, we may always modify u and f to recast the problem so that the

boundary conditions are homogeneous as described in Chapter 6, Secs. 5 and 6). The “dis-

crete boundary Green™s functions on the j-th element”, which we shall call the “elemental

Green™s functions” for short, are de¬ned as the pseudospectral solutions of

(x, y) ∈ „¦j

on (22.30a)

Lγi = 0

(xk , yk ) ∈ ‚„¦j (22.30b)

γi (xk , yk ) = δik

In words, each function γi (x, y) vanishes at all but one of the points on the boundary of

the element. Strictly speaking, we should add subscripts “j” to both the boundary Green™s

functions and to the boundary point labels to specify the j-th element, but we have sup-

pressed these for simplicity. We shall append the element index below wherever needed to

avoid confusion.

In addition, it is usually necessary to cover the domain with curved quadrilaterals and

then map each into the unit square as described in the next section. After such a trans-

formation, L, f , „¦j , and ‚„¦j should have tildes to denote that they have been altered by

the mapping. For simplicity, such tildes have been omitted. The mapping is important,

however, and forms the theme of the next section.

To enforce continuity of u(x) at the inter-element boundaries, it is helpful to assemble

the elemental Green™s functions into global Green™s functions that display the total in¬‚u-

ence of a given boundary point. Let K denote an index that runs over all interior boundary

22.9. THE INFLUENCE MATRIX METHOD 489

1

1 2

3

2 3

4

Figure 22.2: Schematic of an L-shaped domain which has been divided into three square

elements. The boundaries between elements are shown by dotted lines. The grid points

on the interior of these intra-element boundaries are shown by numbers in circles. We need

four “boundary elemental Green™s functions”. The Green™s function G2 (x, y) is equal to

one at boundary grid point 2 and is zero at each of the other three circled points. G2 is

also zero at all grid points on the boundaries of the L-shaped domain (not marked). The

cross-hatching shows the two elements where G2 (x, y) is non-zero.

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

490

points and let J(K) and J (K) denote the elements whose common wall contains the point

labelled by K. Fig. 22.2 schematically illustrates this system for an L-shaped domain par-

titioned into three elements. The index K runs from 1 to 4. The point with K = 2, for

example, lies on the boundary between elements 1 and 2, so that J(2) = 1 and J (2) = 2.

Further, let i(K) and i (K) denote the label of the K-th point in the local number system

on the elements j = J(K) and j = J (K). Then, taking j(x, y) to denote the number of the

element which contains the point (x, y),

±

γi(K),J(K) j = J(K)

GK (x, y) ≡ (22.31)

0 j = J(K), J (K)

γi (K),J (K) j = J (K)

At corner points where four domains come together, the de¬nition of GK must be general-

ized to include four non-zero pieces.

We also must compute a single particular integral pj (x, y) for each element; these solve

x, y ∈ „¦j on ‚„¦j (22.32)

L pj = f & pj = 0

Then

Kmax M

(22.33)

u(x, y) = uK GK (x, y) + pj (x, y)

j=1

K=1

where Kmax is the total number of interior boundary points (=4 in Fig. 22.2), M is the num-

ber of subdomains, and uK denotes the value of u(x, y) at the K-th boundary point. Eq.

(22.33) is the two-dimensional generalization of (22.28); hLj (x) and hRj (x) are the elemen-

tal boundary Green™s functions for that one-dimensional example.

The usual patching conditions are that both u and its derivatives should be continuous

at all boundary points. (Strictly speaking, one wants continuity along the whole interval of

each inter-element boundary, but with any discretization, the best one can do is to impose

such conditions only at the grid points.) The GK (x, y) are constructed to be continuous at

the boundary points. The particular integrals pj (x, y) are all zero at the boundary points

and therefore are continuous also. Thus, the only non-trivial conditions are those on the

derivatives.

Imposing continuity of the normal derivative, that is, the derivative in the direction per-

pendicular to the wall, gives us one condition at each of the Kmax interior boundary points.

These continuity conditions give a matrix problem of dimension Kmax to determine the free

parameters in (22.33), which are the uK . Note that because Kmax N , the dimension of

this element-coupling matrix system is relatively small.

Naturally, one wants continuity of the tangential derivative, too. However, continuity

of u implies that the numerical solutions in two neighboring elements agree at each of

the boundary points along their common wall. It follows that if we expand ui and uj as

one-dimensional functions of the coordinate tangential to the wall, their interpolating poly-

nomials will be identical. It follows that the derivatives tangential to the wall are forced

to match because ui and uj match along the wall. (For simplicity, this argument assumed

that the number of boundary points is the same for each element, but this assumption may

be relaxed without altering the conclusion: the jump in the tangential derivative across the

wall is either zero or spectrally small.)

At a corner point, this argument fails. It seems that one would need to explicitly im-

pose continuity of the derivatives with respect to both x and y, which would give more

conditions than unknowns uK , creating an overdetermined matrix system. It turns out,

22.10. TWO-DIMENSIONAL MAPPINGS & SECTORIAL ELEMENTS 491

however, that it is suf¬cient to impose continuity of either derivative at a corner point; the

jump in the other derivative will be spectrally small (Canuto et al., 1988, pg. 454).

Once the uK have been found, the R. H. S. of (22.33) is completely known. The matrix

systems which must be solved to determine the solution within each element are rather

small: the dimension is Nj , the number of grid points that lie in the interior (and not on

the boundaries) of the j-th element. In most practical applications, Nj ≈ 50“100 (Rønquist

& Patera, 1987b).

22.10 Two-Dimensional Mappings & Sectorial Elements

Spectral elements, like ¬nite elements, may accomodate very irregular or complicated ge-

ometry through mapping. Each element is square in the computational coef¬cients (r, s),

but may be a highly distorted quadrilateral in physical space where the coordinates are

(x, y).

The mapping mechanics are a straightforward (if messy) extension of the one-dimensional

transformations of Chapter 16. One writes

r, s ∈ [’1, 1] (22.34)

x = F (r, s) & y = G(r, s)

and then applies the chain rule to transform derivatives with respect to the physical coor-

dinates x and y into derivatives with respect to the computational coordinates (Table E.9,

Appendix E).

The simplest map is that which transforms the rectangle [ax , bx ] — [ay , by ] in physical

space into the unit square:

1 1

[(bx + ax ) + (bx ’ ax ) r] y = [(by + ay ) + (by ’ ay ) s] (22.35)

x= &

2 2

However, one is free to use analytical maps of almost any form desired.

When the geometry is so complex so that a closed-form mapping for each element can-

not be found, the usual strategy is to employ a so-called “isoparametric mapping”. The

mapping functions F (r, s) and G(r, s) are approximated by double Chebyshev series of

the same order as will be used for the solution u(x, y), i. e.

N N

(22.36)

x= Xij Ci (r) Cj (s)

i=0 j=0

where the Ci are the usual cardinal functions on the “Lobatto” grid (Appendix F, Sec. 6),

plus a similar expansion for y(r, s).

The earliest strategies for coping with curving boundaries were quite crude. A circle

might be approximated by an octagon, or a curving boundary replaced by a zigzag of line

segments that alternately paralleled the x and y axis so that the computational domain was

a union of rectangles. As higher order ¬nite elements were developed, such crudeness

was no longer adequate: the poor approximation to the boundary would erase the high

accuracy inherent in the basis set. At the other extreme, using an eighth order polynomial

for the map functions while employing only a fourth order basis for u(x, y) is also stupid;

the numerical solution cannot resolve the ¬ne structure in the mapping functions. The

isoparametric strategy is the common-sense tactic of approximating all parts of the problem

” solution & mapping functions ” by polynomials of the same order.

Generating numerical grids is a complex topic in and of itself; (Thompson, Warsi, and

Mastin, 1985). However, when the geometry is not complex, or when the number of ele-

ments is suf¬ciently large, it is not dif¬cult to devise a systematic strategy for mapping the

whole domain, one curvilinear quadrilateral at a time.

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

492

The “Global Element School” (Delves & Hall, 1979, Delves & Phillips, 1980, Kermode,

McKerrell & Delves, 1985, Delves et al., 1986, and McKerrell, 1988) have creatively used

mapping to neutralize corner singularities. The basic strategy, which is borrowed from

Wait, is to map a sectorial or triangular element in physical space into the unit square in the

computational coordinates r and s using a map which is designed to remove a speci¬c type

of singularity.

To explain this strategy, it is convenient to use the polar coordinates (in physical space)

y

ρ≡ x2 + y 2 (22.37)

& θ = arctan

x

where ρ = 0 is the location of the boundary corner where u(x, y) is singular. To apply the

method, one must be able to characterize the branch point using local analysis. Assume,

for example, that the singularity is of the form

u(x, y) = [constant] ρw f (θ) + less singular functions (22.38)

for some non-integral constant w.

When the element in (x, y) is the sector

ρ ∈ [0, ±] θ ∈ [0, β] (22.39)

&

where ± is the radius of the sector and β is the sector™s angular width, then the mapping

which (i) transforms the sector into the unit square and (ii) neutralizes a singularity of the