. 100
( 136 .)


patching is very close to the variational formalism without being identical.
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
(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.

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

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

1 2
2 3

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.

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

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

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
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,

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.
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.

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)
ρ≡ x2 + y 2 (22.37)
& θ = arctan
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


. 100
( 136 .)