An efficient direct solver for the boundary concentrated FEM in 2D.
The boundary concentrated FEM, a variant of the hp-version of the finite element method, is proposed for the numerical treatment of elliptic boundary value problems. It is particularly suited for equations with smooth coefficients and non-smooth boundary conditions. In the two-dimensional case it is shown that the Cholesky factorization of the resulting stiffness matrix requires O(Nlog4N) units of storage and can be computed with O(Nlog8N) work, where N denotes the problem size. Numerical results confirm theoretical estimates.
Introduction
The recently introduced boundary concentrated nite element method of [8] is a numerical
method that is particularly suited for solving elliptic boundary value problems with the following
two properties: a) the coe-cients of the equations are analytic so that, by elliptic regularity,
the solution is globally, the solution has low Sobolev regularity due to, for example,
boundary conditions with low regularity or non-smooth geometries. The boundary concentrated
FEM exploits interior regularity in the framework of the hp-version of the nite element method
(hp-FEM) by using special types of meshes and polynomial degree distributions: meshes that
are strongly rened toward the boundary (see Figs. 5, 6 for typical examples) are employed in
order to cope with the limited regularity near the boundary; away from the boundary where
the solution is smooth, high approximation order is used on large elements. In fact, judiciously
linking the approximation order to the element size leads to optimal approximation results (see
Theorem 2.4 and Remark 2.6 for the precise notion of optimality).
In the present paper, we focus on the boundary concentrated FEM in two space dimensions and
present a scheme for the Cholesky factorization of the resulting stiness matrix that requires
O(N log 4 N) units of storage and O(N log 8 N) work; here, N is the problem size. The key to
this e-cient Cholesky factorization scheme is an algorithm that numbers the unknowns such
that the prole of the stiness matrix is very small (see Fig. 1 for a typical sparsity pattern of
the Cholesky factor). Numerical examples conrm our complexity estimates.
The boundary concentrated FEM can be used to realize a fast (i.e., with linear-logarithmic
complexity) application of discrete Poincare-Steklov and Steklov-Poincare operators as we will
discuss in Section 2.5. This use of the boundary concentrated FEM links it to the classical
boundary element method (BEM). Indeed, it may be regarded as a generalization of the BEM:
While the BEM is eectively restricted to equations with constant coe-cients, the boundary
concentrated FEM is applicable to equations with variable coe-cients yet retains the rate of
convergence of the BEM. Since the Cholesky factorization of the stiness matrix allows for an
exact, explicit data-sparse representation of boundary operators such as the Poincare-Steklov
operator with linear-logarithmic complexity, the boundary concentrated FEM provides a sparse
direct solver for the 2D BEM (because it directly computes the full set of Cauchy data on the
boundary corresponding to L-harmonic functions) that has almost linear complexity. It is also
a new alternative to modern matrix compression techniques now used in BEM.
We present a Cholesky factorization scheme for the boundary concentrated FEM in two dimen-
sions, that is, a direct method. We mention that iterative methods for solving the system of
linear equations arising in the boundary concentrated FEM are considered in [8]. Depending
on the boundary conditions considered, dierent preconditioners are required for an e-cient
iterative solution method. For example, while for Dirichlet problems the condition number
of the stiness matrix grows only polylogarithmically with the problem size, [8], Neumann
problems require more eective preconditioning. A suitable preconditioner, which has block-diagonal
structure, was proposed in [8]. We point out, however, that an application of this
preconditioner requires an inner iteration making our direct solver an attractive option.
We conclude this introduction by mentioning that, although we restrict our exposition to the
case of symmetric positive denite problems, our procedure can be extended to non-symmetric
problems by constructing the LU-decomposition rather than the Cholesky factorization.
Boundary Concentrated FEM
In this section, we present a brief survey of the boundary concentrated FEM. We refer to [8]
for a detailed description of this technique and complete proofs. For the sake of concreteness,
we discuss Dirichlet problems although other types of boundary conditions such as Neumann
or mixed boundary conditions can be treated analogously.
2.1 Problem class and abstract Galerkin FEM
For a polygonal Lipschitz
domain
we consider the Dirichlet problem
in
@
where the dierential operator L is given by
with uniformly (in x 2
symmetric positive denite matrix
i;j=1 . Moreover, in the
boundary concentrated FEM, we assume that A and the scalar functions a are analytic on
The operator
is the trace operator that restricts functions
on
to
the boundary @ We assume that the operator L generates an H
i.e.,
The boundary value problem (2.1) is understood in the usual, variational sense. That is, solving
is equivalent to the problem:
Find
with
Z
The standard FEM is obtained from the weak formulation (2.5) by replacing the space H
withanitedimensionalsubspace VN  H
8 For the Dirichlet problem (2.1), we introduce
the trace space
@
For an approximation N 2 YN to  we can then dene the FEM for (2.5) as follows:
Find uN 2 VN s.t.
Z
The coercivity assumption (2.4) ensures existence of the nite element approximation uN .
Furthermore, by Cea's Lemma there is a C > 0 independent of VN such that uN satises
C inf
In practice, the approximations N are obtained with the aid of the L 2 -projection operator
YN by setting
the function QN  is dened
by
In the next section, we specify the approximation spaces VN , the proper choice of which is
intimately linked to the regularity properties of the solution u of (2.1). The analyticity of
the data A, a 0 , and f implies by interior regularity that the solution u is analytic on
if
furthermore
for some - 2 (0; 1], then the blow-up of higher order derivatives near
the boundary can be characterized precisely in terms of so-called countably normed spaces
(see [8] for the details). This regularity allows us to prove an optimal error estimate for the
boundary concentrated FEM in Theorem 2.4 below.
2.2 Geometric meshes and linear degree vectors
For ease of exposition, we will restrict our attention to regular triangulations (i.e., no hanging
nodes) consisting of a-ne triangles. (We refer, for example, to [?] for the precise denition
of regular triangulations.) We emphasize, however, that an extension to quadrilateral and
curvilinear elements is possible. The triangulation of the
domain
consists of
elements K, each of which is the image FK
K) of the equilateral reference triangle
under the a-ne map FK . We furthermore assume that the triangulation T is
-shape-regular,
i.e.,
Here, hK denotes the diameter of the element K. Of particular importance to us will be
geometric meshes, which are strongly rened meshes near the boundary @
Denition 2.1 (geometric mesh) A
-shape-regular (cf. (2.10)) mesh T is called a geometric
mesh with boundary mesh size h if there exist c 1 , c 2 > 0 such that for all K
1. if K \
@
2. if K \
@
dist (x; @  hK  c 2 sup
dist (x; @
Typical examples of geometric meshes are depicted in Figs. 5, 6. Note that the restriction to
the boundary
@
of a geometric mesh is a quasi-uniform mesh, which justies speaking of a
\boundary mesh size h".
In order to dene hp-FEM spaces on a mesh T , we associate a polynomial degree p K 2 N with
each element K, collect these p K in the polynomial degree vector p := (p K ) K2T and set
where for p 2 N we introduce the space of all polynomials of degree p as
The linear degree vector is a particularly useful polynomial degree distribution in conjunction
with geometric meshes:
Denition 2.2 (linear degree vector) Let T be a geometric mesh with boundary mesh size
h in the sense of Denition 2.1. A polynomial degree vector is said to be a linear
degree vector with slope  > 0 if
An important observation about geometric meshes and linear degree vectors is that the dimension
dimS
of the space S
proportional to the number of points
on the boundary:
Proposition 2.3 ([8]) Let T be a geometric mesh with boundary mesh size h. Let p be a
linear degree vector with slope  > 0 on T . Then there exists C > 0 depending only on the
shape-regularity constant
and the constants c 1 , c 2 ,  of Denitions 2.1, 2.2 such that
dimS
K2T
K2T
2.3 Error and complexity estimates
We formulate an approximation result for the hp-FEM on geometric meshes applied to (2.1):
Theorem 2.4 ([8]) Let u be the solution to (2.1) with coe-cients A, a 0 , and right-hand side f
analytic
on
. Assume additionally that u 2 H
for some - 2 (0; 1). Let T be a geometric
mesh with boundary mesh size h and let p be a linear degree vector on T with slope  > 0.
Then the FE solution uN given by (2.7) satises
C
The constants C, b > 0 depend only on the shape-regularity constant
, the constants c 1 , c 2
appearing in Denition 2.1, the data A, c, f
, and -, kuk H
. For  su-ciently large the
boundary concentrated FEM achieves the optimal rate of convergence
number of boundary points:
Remark 2.5 Theorem 2.4 is formulated for the Dirichlet problem (2.1). Analogous approximation
results hold for Neumann or mixed boundary conditions as well.
Remark 2.6 Theorem 2.4 asserts a rate O(n - ) for the boundary concentrated FEM, where
This rate is optimal in the following sense: Setting for - 2 (0; 1)
on
we can introduce the n-width
En
sup
v2En
where the rst inmum is taken over all subspaces E n  H
of dimension n. It can then be
shown that Cn -  d n for some C > 0 independent of n.
2.4 Shape functions and stiness matrix
In order to convert the variational formulation (2.7) into a system of linear equations, a basis
of the nite element space S
has to be chosen. Several choices of basis
functions (\shape are standard in hp-FEM, [2, 12, 7]. Their common feature is
that the shape functions can be associated with the topological entities \vertices," \edges,"
and \elements" of the triangulation T . This motivates us to introduce the following notion of
\standard" bases:
Denition 2.7 A basis B of S
P(
said to be a standard hp-FEM basis if each shape
function exactly one of the following three categories:
1. vertex shape functions: ' is a vertex shape function associated with vertex V if supp '
consists of all elements that have V as a vertex;
2. edge shape functions: ' is an edge shape function associated with edge e of T if supp '
consists of the (at most two) elements whose edge includes e;
3. internal shape functions: ' is an internal shape function associated with element K if
For a standard basis in this sense, we assign spatial points, called nodes, to degrees of freedom
as follows:
1. we assign to the shape function associated with vertex V the point
2. we assign to the side shape functions associated with edge e the midpoint of e;
3. we assign to the internal shape functions associated with element K the barycenter of K.
One example of a standard basis in the sense of Denition 2.7 is obtained by assembling (see,
e.g., [2, 12, 7]) the so-called hierarchical shape functions:
Example 2.8 We construct a basis of the space S
in two steps: In the rst step, we
dene shape functions on the reference element ^
K. In the second step, we dene the basis of
by an assembling process.
1. step: Dene one-dimensional shape functions  i on the reference interval ( 1; 1) by
where the functions L i are the standard Legendre polynomials and the scaling factors c i are
given by c
Denote by v i , the three vertices of ^
K and by i , the three edges (we
assume polynomial degrees that we associate
with the edges i and let p 2 N be the polynomial degree of the internal shape functions. We
then dene vertex shape functions V, side shape functions
functions I as follows:
V := the usual linear nodal shape functions n i with
y
I
The side shape functions S 2 , S 3 are obtained similarly with p 1 replaced with
suitable coordinate transformation. Note that internal shape functions vanish on @ ^
K and that
edge shape functions vanish on two edges.
2. step: Shape functions as dened on the reference element ^
K are now assembled to yield a
basis of S
First, the standard piecewise linear hat functions are obtained by simply
assembling the shape functions V of each element. The internal shape functions are simply
taken as
else
I is the internal shape function on the reference element ^
K dened above. It
remains to assemble the side shape functions. To that end, we associate with each edge e a
polynomial degree p e := min fp K j e is edge of Kg. Let e be an edge shared by two elements
K 0 . For simplicity of notation assume that the element maps FK , FK 0 are such that FK( 1
and that additionally FK (x; (we refer to [2, 7] for
details on how to treat the general case). We then dene p e 1 edge shape functions ' i;e
associated with edge e by setting
else
An analogous formula holds for edges e with e  @
Once a basis of
chosen, the hp-FEM (2.7) can be formulated as seeking the
solution U of a system of linear equations
We mention in passing that computing the stiness matrix A and the load vector F to su-cient
accuracy can be accomplished with work O(dimVN ), [8].
2.5 Sparse Factorization of the Schur complement
We discuss how the Cholesky factorization of the stiness matrix A leads to the explicit sparse
representation of discrete Poincare-Steklov operators. Let T be the Poincare{Steklov operator
(Dirichlet-Neumann map)
where
1 u is the co-normal derivative of the solution u to the equation
condition
and the trace space YN :=
VN , the discrete
approximation
N is dened as follows: For  2 YN , the approximation TN  2 Y 0
is given by
where ev 2 VN is an arbitrary extension of v satisfying
An analysis of the error T TN was presented in [8]:
Theorem 2.9
Let
be a polygon. Then the following two statements hold:
1. There exists - 0 > 1=2 such that the Poincare-Steklov operator T maps continuously from
2. Under the hypotheses of Theorem 2.4 (with - 2 (0; 1) as in the statement of Theorem 2.4)
there holds for arbitrary - 2 [0; -] \ [0; - 0 )
1=2;@
If a standard basis in the sense of Denition 2.7 of the space S
chosen, then the shape
functions can be split into \interior" and \boundary" shape functions. A shape functions is
said to be \interior" if its node (see Denition 2.7) lies in
it is a \boundary" shape function
if its node lies on @ To this partitioning of basis functions corresponds a block partitioning of
the stiness matrix AN 2 R dim VN dimVN for the unconstrained space VN of the following form:
A A I
A I A II
The subscript I indicates the interior shape functions and marks the boundary shape func-
tions. If we choose Y 0
then the matrix representation of the operator TN is given by the
Schur complement
TN := A A I A 1
II A I :
Inserting the Cholesky factorization LL leads to the desired direct FE method for the
Because the Cholesky factor L can be computed with linear-logarithmic complexity (see Section
provides an e-cient representation of the Poincare-Steklov operator. We
nally mention that our factorization scheme for the Schur complement carries over verbatim
to the case of the Neumann-Dirichlet map and also to the case of mixed boundary conditions.
3 Cholesky factorization of the stiness matrix
This section is devoted to the main result of the paper, the development of an e-cient Cholesky
factorization scheme for the stiness matrix arising in the 2D-boundary concentrated FEM.
The key issue is the appropriate numbering of the degrees of freedom. First, we illustrate this
numbering scheme for the case of geometric meshes and constant polynomial degree
(Sections 3.2, 3.3). We start with this simpler case because our numbering scheme for the
degrees of freedom (Algorithm 3.10) is based on a binary space partitioning; in the case
the degrees of freedom can immediately be associated with points in space, namely, the vertices
of the mesh. The general case of linear degree vectors is considered in Section 3.4.
For the case we recall that by Denition 2.7, the vertices of the mesh are called nodes.
We denote by V the set of nodes and say that a node V 0 is a neighbor of a node V if there exists
an element K 2 T such that V and V 0 are vertices of K. It will prove useful to introduce the
set of neighbors of a node V 2 V as
is a neighbor of V
because the sparsity pattern of the stiness matrix can be characterized with the aid of the
sets N .
3.1 Nested dissection in Direct Solvers
Let A 2 R NN be a symmetric positive denite matrix. We denote by L its Cholesky factor,
i.e., the lower triangular matrix L with LL A. For sparse, symmetric positive matrices
A 2 R NN it is customary, [5], to introduce the i-th bandwidth  i and the j-th frontwidth ! j of
A as
It can be shown (see Proposition 3.1(i)) that in each row i, only the entries L ij with i  i  j  i
are non-zero. The j-th frontwidth ! j measures the number of non-zero subdiagonal entries in
column j of L, i.e.,
The frontwidth ! is given by
The cost of computing the Cholesky factor L can then be quantied in terms of the numbers
Proposition 3.1 Let A 2 R
NN be symmetric, positive denite. Then
(i) the storage requirement for L is
(ii) the number of
oating point operations to compute L is
N square roots for the diagonal entries L ii ,
1P N 1
multiplications.
In particular, the storage requirement nnz and the number of
oating point operations W can
be bounded by
Proof. See, for example, [5, Chapters 2, 4].
In view of the estimate (3.5), various algorithms have been devised to number the unknowns so
as to minimize the frontwidth !; the best-known examples include the Reverse Cuthill-McKee
algorithm, [4], the algorithm of Gibbs-Poole-Stockmeyer, [6], and nested dissection. For stiness
matrices arising in the 2D-boundary concentrated FEM, we will present an algorithm based
on nested dissection in Section 3.3 that numbers the nodes such that
The basic nested dissection algorithm in FEM reads as follows:
Algorithm 3.2 (nested dissection)
nested dissection(V; N 0 )
input: Set of nodes V, starting number N 0
output: numbering of the nodes V starting with N 0
label the element of V with the number N 0
else f
1. partition the nodes V into three mutually disjoint sets V left , V right , V bdy such that
(a) jV left j  jV right j
(b) jV bdy j is \small"
(c) right [V bdy for all right
2. if V left nested dissection(V left
3. if V right nested dissection(V right
4. enumerate the elements of V bdy starting with the number
return
The key property is (1c). It ensures that the stiness matrix A has the following block structure:
A =@ A lef t;left 0 A >
bdy;left
bdy;right
A bdy;left A bdy;right A bdy;bdyA (3.6)
sparsity pattern of L

Figure

1: Left: mesh and initial geometric partitioning (thick line). Right: sparsity pattern of
Condition (1b) is imposed in view of the fact that the frontwidth !(A) of A can be bounded
by
where are the frontwidths of the submatrices A left , A right . Thus, if
the recursion guarantees that jV bdy j is small and the frontwidths of these submatrices are small,
then the numbering scheme is e-cient. Since nested dissection operates recursively on the
sets V left , V right , its eectiveness hinges on the availability of partitioning strategies for Step 1
of Algorithm 3.2 that yield small jV bdy j. The particular node distributions appearing in the
boundary concentrated FEM will allow us to devise such a scheme in Section 3.3.
3.2 Nested dissection: an Example
We show that for meshes that arise in the boundary concentrated FEM, it is possible to perform
the partitioning of Step 1 in Algorithm 3.2 such that the set V bdy is very small compared to V.
We illustrate this in the following example.
Example 3.3 Consider meshes that are rened toward a single edge as shown in Fig. 1. The
thick line partitions R 2 into two half-spaces H< , H> , and the nodes are partitioned as follows:
has a neighbor in H> g [ has a neighbor in H<
Note that due to choosing the partitioning line as the center line, we have
(Estimates of this type are rigorously established in Lemma 3.4 below.) We then proceed as
in Algorithm 3.2 by partitioning along a \center line" of the subsets of nodes (a more rigorous
realization of this procedure is Algorithm 3.9 below that is based on binary space partitioning).
We note that the subsets V left , V right have a similar structure as the original set V; thus, they
are partitioned satisfying an estimate analogous to (3.8). Using this partitioning scheme in
Algorithm 3.2 leads to very small frontwidths: For a mesh with nodes, a frontwidth
obtained (see Fig. 1 for the actual sparsity pattern). We analyze this example in
more detail in Examples 4.1, 4.2 below.
The bounds in (3.8) are \geometrically clear." A more rigorous proof is established in the
following lemma.
Lemma 3.4 Let T be a geometric mesh with boundary mesh size h on a
domain
. Fix b 2
@
and choose a partitioning vector t 6= 0 such that the following cone condition is satised (cf.
Fig.
ni > -jx bj j~njg \ B  (b)
Dene the half-spaces
and set
is a node of T and
has a neighbor in H> g [ has a neighbor in H<
Then
log jV  j;
where
depends only on -, , and the constants describing the geometric mesh T . In
particular,
is independent of the point b.
Proof. Let l be the line passing through the point b with direction ~ n,
0g. Next, dene the function
dist
The key property of d is that
Denoting by K bdy the set of all elements that intersect the line l, K bdy
we can bound
K2K bdy
Z
Z
In order to proceed, we need two assertions:
Assertion 1: For - 2 (0; 1) given by the cone condition (3.9) there holds
dist (x; @  dist (x; b) p
dist (x; @ 8x\ l: (3.12)
PSfrag replacements
@
H<
H>

Figure

2: Notation for partitioning at a boundary point b.
The rst estimate of (3.12) is obvious. For the second one, geometric considerations show for
x\ l
dist
dist
dist
where C 0 is the lateral part of @C t . This proves (3.12).
Assertion 2: There exists C > 0 depending only on the parameters describing the geometric
mesh and the parameters of the cone condition such that
Again, the rst bound is obvious. For the second bound, let x 2 K for some K 2 K bdy and
choose xK 2 K \ l. Then
dist dist
dist
where we used (3.12). Next, we use the properties of geometric meshes and (3.10) to get
dist (x; b)g  max fh; hK g +p
dist
Inserting this bound in (3.11) gives
Z
dx  C
Z
r dr  Cj log hj:
Since jV  j  h 1 (cf. [8, Prop. 2.7]), we have proved the rst estimate of the lemma.
For the second estimate of the lemma, we note that the boundary parts < :=
@
\H< \B  (b)
and >
@
\H> \B  (b) have positive lengths. Thus we have jV  \ < j  h 1 and jV  \ > j  h 1
and a fortiori jV left j  jV  j  jV right j; the proof of the lemma is now complete.
The reason for the eectiveness of the partitioning strategy in Example 3.3 is that at each stage
of the recursion, Algorithm 3.2 splits the set of nodes into two sets of (roughly) equal size, and
a set of boundary nodes V bdy that is very small. The property (3.8), proved in Lemma 3.4,
motivates the following denition:
Denition 3.5 ((
q)-balanced partitioning) The nested dissection Algorithm 3.2 is said
to be (
q)-balanced for a set V if at each stage i of the recursion there holdsjV (i)
right j
left
Here, the superscripts i indicate the level of the recursion.
For a (
q)-balanced nested dissection algorithm, we can then show that the frontwidth grows
only moderately with the problem size:
Proposition 3.6 If the nested dissection Algorithm 3.2 is (
q)-balanced for V, then the numbering
generated by Algorithm 3.2 leads to a frontwidth !(A) of the stiness matrix A with

Proof. The assumption that the algorithm is (
q)-balanced implies easily
right j
Thus, jV (i) j
jVj, and the depth of the recursion is at most
since the recursion stops if jV (i) 1. The bound (3.7) then implies

where we used the denition of C
of the statement of the proposition.
In Example 3.3, we studied the model situation of meshes rened toward a straight edge.
In view of Lemma 3.4, we expect the partitioning strategy of be (
1)-balanced. Hence, we
expect the frontwidth to be of the order O(log 2 jVj). In Examples 4.1, 4.2, we will conrm this
numerically.
3.3 Node Numbering for geometric meshes: the case
We now present a partitioning strategy that allows Algorithm 3.2 to be (
1)-balanced for node
sets that arise in the boundary concentrated FEM. The partitioning rests on the binary space
partitioning (BSP), [3], which is reproduced here for convenience's sake:
Algorithm 3.7 (BSP)
input: Set of points X , partitioning vector t
output: partitioning of X into X< , X> , X= with jX < j  jX > j and
1. determine the median m of the set fhx; ti
2. X< := fx
return
Remark 3.8 Since the median of a set can be determined in optimal (i.e., linear in the number
of elements) complexity (see, e.g., [1, 9]), Algorithm 3.7 can be realized in optimal complexity.
The next algorithm formalizes our procedure of the example in Section 3.2.
Algorithm 3.9 (subdomain numbering)
nodes V, vector t, starting number N 0
%output: numbering of nodes V
label the element of V with the number N 0
else f
1.
2.
right
3. if V left
4. if V right
5. enumerate the elements of V bdy starting with the number
return
Algorithm 3.9 allows us to number e-ciently nodes of a mesh that is rened toward a line as
in Fig. 1. Our nal algorithm splits the domain into subdomains, each of which can be treated
e-ciently by Algorithm 3.9.
Algorithm 3.10 (node numbering)
1. split the
domain
into
subdomains
choose vectors t i
2.
3.
4. N := 1
5. for do f
call numbering
6. number the nodes of V bdy starting at N
The
subdomains
i and the partitioning vectors t i should be chosen such that
(a) jV bdy
(b) the partitioning in the subsequent calls of
1)-balanced in the
sense of Denition 3.5.
To obtain guidelines for the selection of
subdomains
i and partitioning vectors t i , it is valuable
to study examples where Condition (a) or Condition (b) are not satised. This is the purpose
of the next example.
Example 3.11 The left and center pictures in Fig. 3 illustrate situations in which Conditions
(a), (b) are violated: In the left picture of Fig. 3, the common boundary
@
@
j is
tangential to
@
at b, and thus we cannot expect jV bdy (cf. the cone condition
(3.9)). In the center picture of Fig. 3, the partitioning vector t i is parallel to the outer normal
vector n(x) at the boundary point x 2 @ This prevents the partitioning from being
1)-balanced, since at some stage of the recursion, Condition (3.15) will be violated (note
again the cone condition (3.9)). We refer to Example 4.5 below, where this kind of failure is
demonstrated numerically. Finally, we point out that in the center picture of Fig. 3 the vector
for the
subdomain
@
>From the two cases of failure in Example 3.11, we draw the following guidelines:
1. The subdomains should be such that
@
@
j is non-tangential to @
2. For each
subdomain
i , the partitioning vector t i should be chosen such that the cone
condition holds uniformly in b 2
@
are independent of b.
A partition chosen according to these rules is depicted in the right part of Fig. 3.
Remark 3.12 The guideline for choice t i such that the cone condition (3.9) is satised at each
point
@
\@
guarantees that in the partitioning, jV (i)
O(log jVj) at each stage i of the
recursion (see Lemma 3.4). Thus, the partitioning is (
1)-balanced if we can ensure (3.14),
that is, that V left and V right are comparable in size. Note that this could be monitored during
run time.
Remark 3.13 In all steps of the recursion in Algorithm 3.9, we use a xed partitioning vector
. This is done for simplicity of exposition. In principle, it could be chosen dierently in each
step of the recursion depending on the actual set to be partitioned. Since the partitioning
strategy should be (
1)-balanced, one could monitor this property during run time and adjust
the vector t as necessary.
We conclude this section with a work estimate for the case
PSfrag replacements
@

@
@
PSfrag
replacements
@

@
@
@

x
PSfrag
replacements
@

@
@
@

x

Figure

3: Choosing subdomains and partitioning vectors. Left and center: cases of failure.
Right: possible partitioning with arrows indicating a good choice of vector t i ; in the three
subdomains without an arrow, t i may be chosen arbitrary.
Proposition 3.14 Let V be the set of nodes corresponding to a geometric mesh on a domain
. Assume that the
subdomains
Algorithm 3.10 are chosen
such that a) jV bdy the partitioning in each call
1)-balanced. Then the frontwidth !(A) of the stiness matrix on geometric meshes with
bounded by
where the constant C is independent of jVj. The storage requirements nnz and work W for the
Cholesky factorization are bounded by
nnz  CjVj log 2 jVj; W  jVj log 4 jVj:
Proof. The hypothesis that Algorithm 3.10 is based on a (
1)-balanced partitioning together
with Proposition 3.6 implies jVj). The estimates concerning storage requirement
and work then follow from Proposition 3.1.
Remark 3.15 On each level the nodes of V (i)
bdy are numbered arbitrarily. Suitable numbering
strategies of these sets could further improve the frontwidth !(A).
3.4 Node Numbering: geometric meshes and linear degree vectors
We now consider the case of geometric meshes and linear degree vectors. We proceed as in
Section 3.3 for the case identifying degrees of freedom with points in space. We use
the notion of nodes introduced in Denition 2.7 and denote by V the set of all nodes. We count
nodes according to their multiplicity, that is, the number of shape functions corresponding to
that node. This procedure is justied by the fact that shape functions associated with the same
node have the same support and therefore the same neighbors. As in the case
that node V 0 is the neighbor of a node V , if the intersection of the supports of the associated
shape functions has positive measure. The set of neighbors of a node is dened as in (3.1).
Remark 3.16 Nodes are counted according to their multiplicity. If a one-to-one correspondence
between points in space and degrees of freedom is desired, one could choose distinct
nodes on an edge (e.g., uniformly distributed) to be assigned to the shape functions associated
with that edge; likewise distinct nodes in an element could be selected to be assigned to shape
functions associated with that element. The performance of the algorithms below will be very
similar.
To this set of nodes and this notion of neighbors, we can apply Algorithm 3.10. In order to
estimate the resulting frontwidth, we need the analog of Lemma 3.4.
Lemma 3.17 Let T be a geometric mesh on a
domain
, p be a linear degree vector, and
assume the cone condition (3.9). Dene the half-spaces
and set
is a node and
has a neighbor in H> g [ has a neighbor in H<
Then
log 3 jV  j;
where
depends only on -, , and the constants describing the geometric mesh T and the
linear degree vector p.
Proof. The proof of this lemma is very similar to that of Lemma 3.4. For the bound on jV bdy j
we have to estimate (using the notation of the proof of Lemma 3.4)
K2K bdy
The desired bound then follows as in the proof of Lemma 3.4 if we observe that
In view of the appearance of the exponent 3 in Lemma 3.17, we expect Algorithm 3.10 to be
3)-balanced. In this case, we can obtain the following result for the performance of the
numbering obtained by Algorithm 3.10:
Theorem 3.18 Let T be a geometric mesh and p be a linear degree vector. Set
.
Assume that
subdomains
i and partitioning vectors t i , Algorithm 3.10 are
chosen such that a) jV bdy the partitioning in each call
is (
3)-balanced. Then the frontwidth !(A) of the stiness matrix is bounded by
The storage requirements nnz and work W for the Cholesky factorization are bounded by
nnz  CN log 4 N; W  N log 8 N:
Remark 3.19 In view of Remark 3.8, Algorithm 3.10 requires O(N log N) work (i.e., optimal
complexity) to compute the numbering.
Remark 3.20 We assumed that the mesh consists of triangles only. However, Algorithm 3.10
can be applied to meshes containing quadrilaterals and curved elements. Theorem 3.18 holds
verbatim in these cases as well.
frontwidth for refinement towards a single edge
log(N)
frontwidth
frontwidth, t=(1,1)
log 2 (N)

Figure

4: Examples 4.1, 4.2: in
uence of partitioning vector in BSP on frontwidth.
In this section, we conrm that the numbering obtained by Algorithm 3.10 allows for computing
the Cholesky factorization with O(N log q N ), q 2 f4; 8g, work. We restrict ourselves to the case
that is, we illustrate Proposition 3.14. In all examples, the nodes on the
boundary correspond to unknowns, i.e., we consider Neumann problems. In all examples, we
use Algorithms 3.10, 3.9 to obtain the numbering of the nodes.
In our computational experiments, the meshes are generated with the code Triangle of
J.R. Shewchuck, [11]. Triangle is a realization of Ruppert's Algorithm, [10], which creates
triangulations with a guaranteed minimum angle of 20 - . Our reason for working with
this particular triangulation algorithm is that it automatically produces meshes
the desired property diamK  dist (K; @
if its input consists of quasi-uniformly distributed
boundary nodes only (the meshes in Figs. 5, 6, for example, are obtained with Triangle from
200 uniformly distributed points on the boundary).
In all tables and gures, N stands for the number of nodes of the mesh generated by Triangle,
! is the frontwidth of the stiness matrix, and nnz the storage requirement for the Cholesky
ops is the number of multiplications, and t chol the CPU-time required to perform the
We implemented the Cholesky factorization in \inner product form,"
where L is computed columnwise and the sparsity pattern of L is exploited.
The basic building block of our procedure is Algorithm 3.9. Our rst example, therefore,
analyzes in detail the situation already discussed in Example 3.3.
Example 4.1 Let be the unit square. For n 2 N the initial input for Triangle
are the points f(i=n; ng [ f(1; 0:3); (1; 1); (0:7; 1); (0; 1); (0; 0:7)g (see Fig. 4 for
Triangle's output for 50). The node numbering is then obtained by applying Algorithm
3.9 with the partitioning vector . The results are collected in Table 1 and
Fig. 4. In view of Proposition 3.6 we expect the frontwidth ! to be O(log 2 N ). The results of

Table

are plotted in Fig. 4, and the observed growth is indeed very close to O(log 2 N ). In
ops t chol [sec] ! nnz

Table

1: Examples 4.1, 4.2: points on edge,
view of Proposition 3.1,
ops  t chol
these estimates.
Example 4.2 The choice in Example 4.1 of the partitioning vector
well-suited to the case of renement toward a straight edge. In view of Lemma 3.4, we expect
Algorithm 3.9 to be still (
1)-balanced for partitioning vectors t that are not parallel to the
normal vector. In order to illustrate that \non-optimal" choices of partitioning vectors t still
lead to (
1)-balanced nested dissection, we consider the same meshes as in Example 4.1 but
employ Algorithm 3.9 with the vector . The results are presented in the right part of

Table

1 and in Fig. 4. We observe that this choice of partitioning strategy leads to very similar
results as in Example 4.1, showing robustness of our algorithm with respect to the choice of
partitioning vector t.
Example 4.3 In this example, the domain is the unit
square
input are n uniformly distributed nodes on the boundary
@
(see Fig. 5 for Triangle's output
200). The node numbering is achieved with Algorithm 3.10 for
subdomains
corresponding vectors t i given
by:
The numerical results are collected in Table 2. We expect the frontwidth to be O(log 2 N ),
which is visible in Fig. 5.
refinement towards bdy of square and clover leaf
log(N)
frontwidth
square
clover leaf
log 2 (N)

Figure

5: Examples 4.3, 4.4: meshes with 200 points on boundary (left) and frontwidth vs.
log N (right).
Example 4.4 In this example, we replace the square of Example 4.3 with
(See Fig. 5 for Triangle's output for boundary points.) The partitioning into four
subdomains and the choice of the partitioning vectors t i is given by (4.1). The numerical results
are collected in Table 3. The expected relation visible in Fig. 5.
Example 4.5 The key feature of the choice of the partitioning vectors t i in Examples 4.3, 4.4
is that, for each i,
sup
x 2
@
@
is the normal at a boundary point x). (4.2)
This condition was identied in Example 3.11 as necessary for the binary space partitioning
strategy with (xed) vector t i to be (
1)-balanced. In this last example, we illustrate that
(4.2) is indeed necessary. To that end, we replace the square of Example 4.3 with
The
subdomains
i and the partitioning vectors t i are chosen as in Example 4.3 and given by
(4.1). A calculation shows that condition (4.2) is satised for c 2 (0; 2=3) and is violated for
subdomains,
4 in Fig. 6). For c close to 2=3 we therefore expect our binary
space partitioning to perform poorly in the sense that the sets V (i)
bdy become large, thus resulting
in large frontwidths. This is illustrated in Table 4, where the frontwidths for dierent values
of c are shown in dependence on n, the number of points on the boundary. Fig. 6 shows that
in the limit 2=3, the frontwidth ! does not grow polylogarithmically in N . In Table 4,
we reported the number N of nodes for the case however, that the
meshes produced by Triangle for the three cases are very similar.

Acknowledgments

: The authors would like to thank Profs. W. Hackbusch and L.N. Trefethen
for valuable comments on the paper.
ops t chol [sec]

Table

2: Example 4.3:
ops t chol [sec]

Table

3: Example 4.4:
-0.4
c=2/3
c=0.3
effect of subdomain choice on frontwidth; domain: dumbbell
log(N)
frontwidth
c=2/3
log 2 (N)

Figure

Example 4.5: mesh with 200 points on boundary (top left), the four subdomains for
(bottom left) and frontwidth vs. log N (right).
1024 2022 107 118 118
2048 4079 128 159 155
131072 264600 370 577 768
262144 529340 423 723 1043

Table

4: Example 4.5:



--R

Time bounds for selection.
A general and exible fortran 90 hp-FE code
On visible surface generation by a priori tree structures.
Computer implementation of the
Computer solution of large sparse positive de
An algorithm for reducing the bandwidth and pro
Element Methods for CFD.
Boundary concentrated
The Art of Computer Programming
A Delaunay re
2d mesh generator

--TR
A Delaunay refinement algorithm for quality 2-dimensional mesh generation
The art of computer programming, volume 3
Computer Solution of Large Sparse Positive Definite
Finite Element Method for Elliptic Problems
Computer implementation of the finite element method
hp-finite element methods;direct solvers;schur complement;meshes refined towards boundary