ERROR ANALYSIS OF THE FINITE VOLUME ELEMENT
METHOD FOR ELLIPTIC AND PARABOLIC
PARTIAL DIFFERENTIAL EQUATIONS
by
Rick Vance Trujillo
B.S., Metropolitan State College of Denver, 1988
M.S., University of Colorado at Denver, 1991
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1996
This thesis for the Doctor of Philosophy
degree by
Rick Vance Trujillo
has been approved
by
Thomas F. Russell
Leopoldo P. Franca
Chaoqun Liu
Thomas A. Manteuffel
Stephen F. McCormick
Date
Trujillo, Rick Vance (Ph.D., Applied Mathematics)
Error Analysis of the Finite Volume Element Method
for Elliptic and Parabolic Partial Differential Equations
Thesis directed by Professor Thomas F. Russell
ABSTRACT
An a priori error analysis of the finite volume element method, a
locally conservative, PetrovGalerkin, finite element method for the numerical
solution of elliptic and parabolic partial differential equations arising in fluid
dynamics applications, is presented. Existing error estimates apply to dis
cretizations of steady diffusion equations by linear finite elements in two spatial
dimensions. These results are extended to steady advectionreactiondiffusion
equations and are generalized to polynomial finite elements of arbitrary order
in three spatial dimensions and to the full range of admissible regularities for
the exact solution. Optimalorder error estimates for h, p, and hp versions of
the method with uniform refinement are derived in a discrete H1 norm, un
der minimal regularity assumptions for the exact solution, the finite element
triangulation, and the finite volume construction. With additional symmetry
assumptions for the finite volumes, multidimensional H1 superconvergence
results are obtained for linear finite elements. Using an elliptic projection ar
gument, we outline the analysis of continuoustime and backward Euler and
CrankNicolson discretetime methods for general parabolic equations.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Thomas F. Russell
IV
DEDICATION
To my Family, with love:
To my Father, who finds me when I am lost
and welcomes me home with open arms;
To my elder brother Yeshua, who selflessly pours out his lifeblood
to save me and bears my burdens upon his shoulders;
To my dear Comforter, who always illuminates the way of
truth and life in the darkest nights;
To my Mother, who graciously intercedes on my behalf
and keeps me close to her heart;
To my Guardian, who subtly guides and protects me
all the days of my life;
To my Brothers and Sisters, who encourage and admonish me
with the example of their lives.
Please never let me be separated from thee,
my precious loving Family!
CONTENTS
Chapter
1. The Finite Volume Element Method................................. 1
1.1 Introduction................................................. 1
1.1.1 Contents of the Thesis................................. 3
1.2 Background................................................... 7
1.3 Numerical Methods ........................................... 9
1.4 The Finite Volume Method (FV)............................... 12
1.5 The Finite Element Method (FE) ............................. 14
1.6 The Finite Volume Element Method (FVE)...................... 18
1.7 Historical Development...................................... 22
1.8 Summary..................................................... 25
2. Theoretical Preliminaries..................................... 26
2.1 Triangulation............................................... 26
2.1.1 FE Triangulation...................................... 27
2.1.2 FVE Primal Mesh ...................................... 28
2.2 Volumization................................................ 32
2.2.1 Volume Construction................................... 32
2.2.2 Volume Symmetry....................................... 38
2.2.3 Volume Regularity..................................... 38
2.3 Versions h, p, and hp...................................... 41
2.4 FVE vs. FE ................................................. 43
vi
2.5 Sobolev Spaces................................................ 44
2.5.1 IntegerOrder Spaces.................................... 45
2.5.2 FractionalOrder Spaces................................. 46
2.5.3 Sobolev Embedding ...................................... 46
2.5.4 Trace Theorem........................................... 47
2.5.5 Poincare Inequality..................................... 47
2.6 Finite Element Analysis....................................... 49
2.6.1 Boundedness ............................................ 50
2.6.2 Approximation Theory.................................... 51
2.6.3 Ellipticity............................................. 52
2.6.4 Error Analysis ......................................... 52
2.6.4.1 Error Splitting ................................ 52
2.6.4.2 Zero Property................................... 53
2.6.4.3 Error Estimate.................................. 53
2.7 Summary....................................................... 54
3. Analysis for Diffusion Equations.................................. 55
3.1 Preliminaries................................................. 56
3.2 Motivation.................................................... 59
3.3 Discrete Norms................................................ 61
3.3.1 Discrete L2 Norm........................................ 62
3.3.2 Discrete H1 Norm ....................................... 63
3.3.3 Discrete Poincare Inequality............................ 64
3.4 Boundedness................................................... 65
3.5 BrambleHilbert Lemma......................................... 68
vii
3.5.1 Applications............................................ 69
3.6 Approximation Theory ....................................... 72
3.6.1 Optimalorder Approximation............................. 72
3.6.2 Superconvergent Approximation........................... 77
3.6.2.1 Extension to Three Dimensions................... 81
3.6.2.2 Degenerate Volume Paradigm.................... 84
3.7 Ellipticity................................................... 87
3.8 Error Analysis................................................ 91
3.8.1 Error Splitting ........................................ 91
3.8.2 Zero Property........................................... 92
3.8.3 Error Estimate.......................................... 92
3.9 Summary ...................................................... 94
4. Analysis for General Elliptic Equations............................ 97
4.1 Preliminaries................................................. 98
4.2 Boundedness.................................................. 100
4.2.1 Advection.............................................. 100
4.2.2 Reaction............................................... 101
4.3 Approximation Theory ........................................ 103
4.3.1 Advection.............................................. 103
4.3.2 Reaction............................................... 104
4.4 Ellipticity.................................................. 108
4.4.1 Reaction............................................... 108
4.4.2 Advection.............................................. 110
4.5 Error Analysis............................................... 113
viii
4.5.1 ReactionDiffusion.................................... 113
4.5.2 AdvectionReactionDiffusion.......................... 114
4.5.3 General Error Estimate ............................... 115
4.6 Elliptic Projection......................................... 116
4.7 Summary .................................................... 118
5. Analysis for Parabolic Equations................................ 126
5.1 Numerical Methods .......................................... 128
5.1.1 ContinuousTime....................................... 128
5.1.2 DiscreteTime......................................... 129
5.2 Components of Analysis...................................... 131
5.2.1 Boundedness .......................................... 131
5.2.2 Approximation Theory.................................. 132
5.2.2.1 Spatial....................................... 133
5.2.2.2 Temporal...................................... 134
5.2.3 Ellipticity........................................... 134
5.2.3.1 Spatial....................................... 135
5.2.3.2 Temporal...................................... 135
5.3 Error Analysis.............................................. 137
5.3.1 Error Splitting ...................................... 137
5.3.2 Zero Property......................................... 138
5.3.3 Energy Argument....................................... 138
5.4 Summary .................................................... 141
5.5 Future Work................................................. 149
References.......................................................... 152
IX
FIGURES
Figure
2.1 An example of the triangulation .............................. 29
2.2 7~h!p restricted to a single element of Th.................... 30
2.3 Construction of the subvolume Vi(Tp)......................... 34
2.4 A circumcenter finite volume U................................ 36
2.5 A general finite volume U..................................... 37
2.6 Construction of a pair of elements in Th/p.................... 40
3.1 Formation of the macroelement Ea = T} U Tf................. 79
3.2 An example of Vhl2 degenerating to Vh......................... 85
x
ACKNOWLEDGMENTS
I gratefully acknowledge the financial support of the CUDenver Math
Department through many teaching assistantships, of Pat Roache through Eco
dynamics Research Associates, Inc. and Sandia National Laboratories of Albu
querque, NM, and of Tom Russell through many research assistantships: most
recently through National Science Foundation Grant No. DMS9312752 and
U.S. Army Research Office Grant No. 34593GS.
I would like to acknowledge some of the people who encouraged and
admonished me during research seminar presentations: Teri Barth, Dave Dean,
Rick Healy, Randy LeVeque, Tom Manteuffel, Steve McCormick, Suely Oliveira,
Jim Otto, Rossen Parashkevov, Pat Roache, Tom Russell, and Hong Wang. I
would also like to acknowledge some of the people who encouraged and ad
monished me in my true vocation as a teacher: Bill Briggs, Tom Kelley, and
my students! And a special thanks to Debbie Wangerin, who is always helping
behind the scenes.
1. The Finite Volume Element Method
1.1 Introduction
This thesis presents an a priori error analysis of the finite volume
element method (FVE) for elliptic and parabolic partial differential equations.
As introduced by Baliga and Patankar [4] and elaborated by McCormick [39]
and in the references cited therein, we can view FVE as a combination of the
standard finite volume method (FV), also known as cellcentered finite differ
ences (see Mitchell and Griffiths [36] for details), and the standard Galerkin
finite element method (FE) (see Ciarlet [11] for details). In FVE, the governing
partial differential equation is posed as a local (integral) conservation law on
a set of control or finite volumesregions of local conservationthat parti
tions the problem domain, as in FV. The unknown FVE solution and known
(boundary and source term) data are systematically discretized by a C piece
wise polynomial trial space of degree p > 1 that is defined on a finite element
triangulation with maximum mesh diameter h that partitions the problem do
main, as in FE. Thus, the finite volume element method can be viewed as
either a systematic generalization of the standard finite volume method or as a
locally conservative, as well as globally conservative, (PetrovGalerkin) variant
of the standard finite element method. In analogy to the h, p, and hp version
FE of Babuska (see [2, 3, 27, 51] and the references cited therein), we consider
h, p, and hp version FVE where increased accuracy is obtained by refining
the computational mesh parametrized by h and/or by increasing the order p
1
of the FE polynomial space. Other methods that are either closely related or
identical to the hversion FVE include: control volume finite element methods
[4], generalized box methods [5], and the method of support operators [48],
Although a large body of theoretical infrastructure and of results for
elliptic and parabolic equations exists for both FV and FE, analysis for the h
version FVE is limited to the foundation laid by Cai and McCormick [8], Cai,
Mandel, and McCormick [9], and culminating in Cai [10] for the numerical
solution of steady diffusion equations by linear finite elements in two spatial
dimensions. Starting from this foundation, we extend and generalize FVE error
analysis to steady and transient advectionreactiondiffusion equations and to
polynomial finite elements of arbitrary order in three spatial dimensions. To
construct an error analysis for FVE, we are guided by the large body of FE
theory for elliptic and parabolic problems and their solution by arbitrary order
finite elementsas contained in Ciarlet [11], Douglas [12], Douglas and Dupont
[13, 14], Wheeler [55], and in references cited thereinand by modern FV error
analysis for the lowestorder methods based on FEstyle arguments for linear
finite elements: e.g., Bank and Rose [5], Ewing et al. [17, 18], Hackbusch [21],
Heinrich [23], Herbin [24], Kreiss et al. [31], Lazarov et al. [32, 33], Manteuffel
and White [35], Morton and Siili [38], Samarskii et al. [46], Siili [50], and Weiser
and Wheeler [54],
By adapting FV and FE arguments to fit the context of FVE and
developing some new arguments unique to FVE, we derive optimalorder error
estimates under the full range of admissible regularities for the exact solu
tion for h, p, and hp version FVE with uniform refinement in a discrete H1
2
norm for Lagrange finite elements of arbitrary orderunder minimal regular
ity assumptions for the exact solution, the finite element triangulation, and
the finite volume construction. With additional uniformity assumptions for
the finite volumes, Hl superconvergence results are obtained for linear finite
elements. Using an elliptic projection argument, we outline the analysis of
continuoustime and backward Euler and CrankNicolson discretetime meth
ods for general parabolic equations.
1.1.1 Contents of the Thesis
The remainder of the thesis is organized as follows. In this chapter, we
give the necessary introduction for the finite volume element method (FVE)
for elliptic equations before we proceed with an analysis of the method. In
Section 1.2, we outline how FVE arises from the approximation of certain
integral conservation laws that lead to elliptic partial differential equations.
In Section 1.3, we detail a three step procedure describing the construction
of numerical methods for elliptic equations. In Sections 1.4, 1.5, and 1.6, we
compare and contrast three specific numerical methods that arise from our
three step procedure: the finite volume method (FV) in Section 1.4, the finite
element method (FE) in Section 1.5, and most importantly, the finite volume
element method (FVE) in Section 1.6. Finally, a brief outline of the historical
development of computational fluid dynamics and roles played by the three
methods is given in Section 1.7.
In Chapter 2, we present some theoretical preliminaries. In Sections
2.1 and 2.2, we define and discuss the computational meshes (i.e., the FE
triangulation Th, the FVE primal mesh T^, and the FVE dual mesh Vh!p)
3
that are fundamental to the implementation and analysis of FVE. In Sections
2.3 and 2.4, we define h, p, and hp version FE and FVE and discuss their
relative merits. In Section 2.5, we introduce standard definitions and results
for Sobolev spaces that are fundamental to FE and FVE analysis. In Section
2.6, we briefly outline an FE analysis for the steady diffusion equation to set a
reference point for an FVE analysis in the next two chapters.
In Chapter 3, we present an h, p, and hp version FVE analysis for
diffusion equations. In Section 3.1, we provide some preliminary definitions and
notation that facilitates the transition from FE to FVE analysis. In Section
3.2, we motivate our FVE analysis and its relation to previous FE work. In
Section 3.3, we define some discrete norm counterparts to the continuous L2
and H1 Sobolev norms of Section 2.5 and demonstrate a discrete Poincare
inequality. As in the FE analysis of Section 2.6, these tools are crucial to FVE
analysis. In Section 3.4, we determine an upper bound or continuity result
in terms of an FVE diffusion functional and our discrete H1 seminorm for
the FVE bilinear form that is used to discretize the model diffusion equation
discussed in Section 2.6. In Section 3.5, we state and outline applications of
the linear BrambleHilbert lemma which is at the heart of our FVE analysis.
In Section 3.6, we use the BrambleHilbert lemma to derive error estimates
for the FVE approximation of elements from fractionalorder Sobolev spaces
by elements from the piecewise polynomial trial spaces employed by the h,
p, and hp version FVE methods: in terms of our FVE diffusion functional,
i71(0)equivalent optimalorder approximation results are derived for pthorder
trial spaces in Section 3.6.1; robust and practical superconvergence results for
4
linear trial spaces on a symmetric volumization are derived in two and three
dimensions in Section 3.6.2. In Section 3.7, we derive an asymptotic discrete
ellipticity result for the FVE diffusion bilinear form in terms of our discrete
Hl seminormaided by the approximation theory from the previous section.
Finally in Section 3.8, we construct an FVE analysis for diffusion equations
with the components developed in the preceding sections: optimalorder and
superconvergence results are obtained for h, p, and hp version FVE in our
discrete H1 seminorm. The key results of this chapter are summarized in
Section 3.9.
In Chapter 4, we note modifications for an h, p, and hp version FVE
analysis for general elliptic equations and elliptic projections. After presenting
some preliminary theoretical infrastructure in Section 4.1, we develop upper
bounds in terms of FVE advection and reaction functionals and our discrete
L2 and H1 norms for the FVE bilinear forms corresponding to the treatment
of advection and reaction terms in a general elliptic equation in Section 4.2.
In Section 4.3, we develop T2equivalent optimalorder approximation theory
results for our FVE advection and reaction functionals. Our goal here is to find
estimates for advection and reaction that do not interfere or add any additional
volume symmetry constraints to the diffusion optimalorder and supercon
vergence results of the previous chapter. In Section 4.4, we derive asymptotic
discrete ellipticity results for the FVE reaction and advection bilinear forms.
Finally in Section 4.5, we construct an FVE analysis for reactiondiffusion and
advectionreactiondiffusion equations with the components developed in the
preceding sections: optimalorder and superconvergence results are obtained
5
for h, p, and hp version FVE in our discrete H1 seminorm. In Section 4.6,
we extend the theory developed so far to analyze the elliptic map which is the
key theoretical tool in our analysis of parabolic equations. The key results of
this chapter are summarized in Section 4.7.
In Chapter 5, we outline an FVE analysis for parabolic equations.
First, we define continuoustime and discretetime numerical methods for gen
eral parabolic equations in Section 5.1: the continuoustime method is an FVE
variant of the method of lines; the discretetime methods are FVE variants of
Backward Euler and CrankNicholson time marching methods. In Sections
5.2 and 5.3, the components of FVE parabolic analysis are developed and as
sembled to yield optimalorder and superconvergence results for h, p, and hp
version continuoustime and discretetime FVE. The key results of this chapter
are summarized in Section 5.4 and future work based on the thesis is outlined
in Section 5.5.
The main results of the thesis are summarized in Theorems 3.3, 4.3,
4.4, 5.4, and 5.5.
6
1.2 Background
The finite volume element method for elliptic equations is based on
the fact that these equations arise from integral conservation laws [39]. For
definiteness and simplicity, we work in two spatial dimensions (modifications to
our analysis for one or three dimensions will be stated as needed) and consider
the following model elliptic equation for an unknown distribution u:
Cu = V (an DVu) + r u = f in 0, (1.1)
where a = (ai,a2), D, and r are advection, diffusion and reaction coefficients,
/ is a source term, and 0 C K2 is the spatial domain. In Chapter 5, we will
consider a model parabolic equation based on (1.1) augmented by the mass
storage term mdtu. For simplicity, we assume here that the coefficients are
smooth functions of space alone to avoid difficulties caused by discontinuities
and nonlinearities in our analysis, though these cases can be handled (see
Section 4.7).
For (1.1), we consider Dirichlet or flux boundary conditions:
u = g or (a u DVu) n = g on 90, (1.2)
where 90 is the boundary of 0, g is boundary data on 90, and n is the
outward unit normal on 90. To guarantee an unique solution in the case of
flux boundary conditions, the source term and boundary data must also satisfy
a compatibility condition:
f /dx+ f g dS = 0. (1.3)
Jil J dÂ£l
When coupled with boundary conditions (1.2), the elliptic equation (1.1) is ca
pable of accurately modeling a wide range of physical phenomena in disparate
7
application areas. Therefore, numerical methods for the elliptic problem de
fined by equation (1.1) subject to (1.2), like FVE, and their corresponding
analysis are vitally important.
Although we stated at the outset that we are investigating numerical
solutions of the elliptic partial differential equation (1.1), we are more precisely
studying local integral conservation laws on subdomains V C 0 that lead to
equations such as (1.1):
/ (a u D'Vu) ndS + f r udx. = / / dx, VhCj], (1.4)
JdV JV JV
where dV is the boundary of V, n is the outward unit normal on the boundary
dV, and dS is a spatial boundary measure. The first integral on the lefthand
side of equation (1.4) is the net advectivediffusive outflow or outflux of u
across dV. The second integral on the lefthand side of equation (1.4) is the
net internal reaction within V which can be thought of as a type of nonlinear
source term. And the single integral on the righthand side of equation (1.4)
is the net contribution due to external sources imposed within V. The local
(when V C 0) and global (when V = 0) conservation law (1.4) states that
the advectivediffusive flux across volume boundary dV is counterbalanced by
reactions and sources within the volume V. In Chapter 5, we will introduce
time dependence and augment (1.4) with the mass storage integral fv m dtu dx.
Equation (1.4) is called the primitive form of the elliptic equation
(1.1) because it contains the essence of the physical model within the least
restrictive mathematical model. Finally, it is equation (1.4)and not equation
(1.1) that we hope to mimic discretely in our numerical methods.
8
1.3 Numerical Methods
FVE is a hybrid method composed from standard finite volume (FV)
and finite element (FE) parent methods and is best understood in terms of
these methods. Therefore, we will detail a general framework for developing
numerical methods that includes FV, FE, and FVE so that the similarities and
differences between the methods may be clearly delineated and understood.
When we develop numerical methods, we take the middle ground
between the strong form equation (1.1) and the conservation law equation
(1.4) and work with the associated variational or weak form equation:
/ wCud^= / /rcdx, VrcGW, (1.5)
Jn Jn
where W is a test space of functions defined on the spatial domain 0. Setting
the stage for later developments, we let (, ) denote the standard L2{0) inner
product and rewrite equation (1.5) in short hand notation:
(Â£u, w) = (/, w), Vw
Now the procedure by which we produce numerical methods is easily
enumerated in a three step discretization process:
(1) Discretize the spatial domain 0 in order to
(a) define a finite number of subdomains in 0 on which to pose equa
tion (1.5) and to
(b) form a computational mesh whose finite number of degreeof
freedom (DOF) points will define the numerical approximation
to the distribution u of equation (1.5);
(2) Discretize the test space W of equation (1.5):
9
(a) replace W by a finite dimensional subspace of polynomials test
functions defined on the spatial subdomains of Step 1 so that
(b) equation (1.5) is localized on each subdomain of Step 1;
(3) Discretize the unknown distribution u and data / of equation (1.5) by
(a) finite difference or finite element approximations with respect to
the DOF points of the computational mesh of Step 1 in order to
(b) pose a discrete version of equation (1.5) on each subdomain of
Step 1.
In Step 1, we discretize the spatial domain 0 in two different but
interconnected ways: first, we discretize 0 by
standard finite element triangulations of 0: with triangular or rect
angular elements denoted by T; Th^p, the pfold refinement of Th that
connects DOF points for polynomials of order p defined each element T
in Th. Hence, Th and 'Thlp are basic or primal computational meshes.
The details of the triangulation constructions are given in Section 2.1;
second, we discretize 0 by
standard finite volume volumizations of 0, denoted by Vh or V^,
with convex or starshaped volumes, denoted by V, to form dual com
putational meshes corresponding to the triangulations Th and 'Thlp,
respectively. The details of the volumization constructions are given
in Section 2.2.
The finite number of degreesoffreedom located at the vertices the 'Thlp tri
angulation are used represent the numerical approximation to the unknown
distribution u of equation (1.5); the set of DOF points on the computational
10
mesh is denoted by {x8}8e/, where / is a DOF point indexing set. Each DOF
point X; of the computational mesh is surrounded by a finite volume, denoted
by Vi, so that there is a onetoone correspondence between DOF points and
finite volumesthe specifics of the finite volume construction are given in Sec
tion 2.2. Therefore, the set of finite volumes in the volumizations Vh or yh/p
can be denoted and enumerated by According to standard finite vol
ume methodology, the value of the unknown u at DOF point x4 represents the
volume average of u on volume V{.
In Step 2, we discretize the test space by replacing W by one of the
following choices:
Sp(Th), the space of piecewise polynomials of order p > 1 defined on
each triangle or rectangle T of the triangulation 'Th;
iS,(V/l) or S(Vh/p), the space of piecewise constants defined on the
volumizations Vh or V^, respectively.
In Step 3, we discretize the unknown u, derivatives of u, and the
source term / by replacing them with
finite element approximations uh, Vuh, and fh, where uh and fh are
representations of u and / in terms of a finite dimensional trial space
chosen to be Sp(Th).
finite difference approximations to u, Vu, and / formed by pthorder
interpolation and/or extrapolation between vertices of the triangu
lation.
Now we can easily describe FV, FE, and FVE by tracing the development of
each numerical method through each of the three steps discussed above.
11
1.4 The Finite Volume Method (FV)
In the finite volume method (FV), we create the computational mesh
Th and its dual mesh Vh in Step 1 of Section 1.3. In Step 2 of Section 1.3, we
choose iS,(V/l) as the test space; the set of basis elements for iS,(V/l) is the set
of characteristic functions {Xi}iei (recall that / is a DOF point indexing set)
with Xi defined by
V'(x)
1, if x G Vi,
(1.7)
0, otherwise.
Referring back to equation (1.5), we choose the Xi basis functions as
test functions and apply Greens Theorem (i.e., perform integration by parts)
to the advection and diffusion term. As a result, we obtain the FV weak form:
/ (an DVu) ndS + / rudx = / /fix, VfiT VA, (1.8)
JdVi JVi JVi
where dVi is the boundary of V, n is the outward unit normal on the bound
ary dV, and dS is a spatial boundary measure. Referring to the shorthand
notation of equation (1.6), we perform integration by parts on (Cu,Xi) t ob
tain the lefthand side integrals of equation (1.8). To set the stage for later
developments, we rewrite equation (1.8) in the shorthand notation:
B(u,x,) = (f,x>), VX, (1.9)
where the the bilinear form B(, ) on the lefthand side of equation (1.9) con
sists of all of the integrals on the lefthand side of (1.8).
In Step 3 of Section 1.3, FV uses a somewhat ad hoc finite difference
methodology that usually corresponds to a reduced or lumped quadrature
strategy for the integrals of equation (1.8). That is, expressions for u, Vrt, and
12
/ are obtained by interpolation and/or extrapolation between vertex values of
u and / in the triangulation (see Mitchell and Griffiths [36] or Roache [43]
for details).
At this point we will defer a discussion about the specifics of obtain
ing a finite volume numerical solution and the implementation of boundary
conditions until we address those issues in conjunction with the finite volume
element method.
The principal advantage of the finite volume approach is that its weak
form (1.8) closely mimics the integral conservation law (1.4). The method
will be locally conservative on each Vj, as well as globally conservative on
0 which is a crucial measure of the physical correctness of the numerical
solution. Of course, the conservative property does not imply the mathematical
correctness of the numerical solution. Usually, we fold that the crude finite
difference approximations to the integrands of (1.8) limit the effectiveness of
FV to relatively simple mesh structures (e.g., Cartesianproduct meshes) and
simple treatments of boundary conditions.
13
1.5 The Finite Element Method (FE)
To alleviate some of the difficulties due to FV, one usually turns to
the finite element method (FE). For FE, we create the computational mesh
with triangular or rectangular elements and its pfold refinement 'Thlp in Step
1 of Section 1.3. In Step 2 of Section 1.3, we choose Sp(l~h) as the test space
and denote this choice as Wh. The set of nodal basis elements for this space
Wh is the set {;};e/:
where 8ij is the Kronecker delta, and ;(x) is found by pthorder interpolation
between the DOF points at the vertices of 'Thlp. The support of ; is denoted
by l~h and it consists of all elements of the triangulation that have X; as a
vertex.
Referring back to equation (1.5), we choose the ; basis functions as
test functions and perform integration by parts on the advection and diffusion
terms. As a result, we obtain the FE weak form given by (1.11) on subdomains
away from the domain boundary <90:
[ DVu V4>i u a V4>i + r u fa dx = [ f fa dx, V ; Â£ Wh, (1.11)
JTth JTth
where %h ^ the support of 4. To describe modifications on the domain bound
ary dO, we define <90; as the intersection of the domain boundary with 7~hthe
closure of the ith test function support: <90; = dfinl~h. When <90; has positive
measure, the FE weak form (1.11) includes the boundary flux integral
/ (art DVu) n ; dS.
Jdti,
(1.12)
14
Regarding boundary modifications, let us pause for a few comments
regarding implementation of the boundary conditions specified in (1.2) with
respect to the boundary integral of (1.12)additional details can be found
the texts of Ciarlet [11] or Johnson [28]. If a Dirichlet boundary condition
is imposed, the boundary flux integrand of (1.12) is an unknown and can be
solved for separately if desired. The information from the Dirichlet boundary
condition is imposed directly upon the trial space discretization of u on <90 and
information from boundary nodes on <90 are passed to the righthand side of
our discrete equations (cf. equation (1.16) below) in an obvious manner that
is discussed in the texts cited above. However, if a flux boundary condition is
imposed, the boundary flux is known data and boundary flux integral of (1.12)
is moved to the righthand side of equation (1.11). In the discussion below,
well assume that a homogeneous flux boundary condition has been imposed
in order to simplify the presentation.
Referring to the shorthand notation of equation (1.6), we perform
integration by parts on (Cu, ;) to obtain the lefthand side integrals of equa
tion (1.11) and the boundary integral of (1.12). To set the stage for later
developments, we rewrite equation (1.11) in shorthand notation:
A(u,i) = (f,i), (1.13)
where the the bilinear form A(,) on the lefthand side of equation (1.13)
consists of all of the integrals on the lefthand side of (1.11).
In Step 3 of Section 1.3, the Galerkin FE uses Wh as a trial space, as
well as a test space, to represent the unknown u and source / of equation (1.11),
but other higher order polynomial spaces defined on the mesh can be used as
15
a trial space and the result is termed PetrovGalerkin FE. In the remainder of
this dissertation, well refer to the Galerkin FE as the finite element method
since it is the FE of choice in the bulk of applications. To describe our finite
element method, we define uh G Wh as the trial space representation of the
numerical solution:
h = ZuJ^(xh (1141
jei
where uh3 = uh(x.j) so that uh has a nodal representation. Similarly, we define
fh G Wh as our trial space representation of the source term:
jei
where fj = /(xj) so that fh = IIphf is the pthorder interpolant of / into
Sp(Th).
Referring to equation (1.13), the finite element method solution of
equation (1.1) is described as follows: find uh G Wh such that
A(uh, 4>i) = (fh, (f>i), V (j>i G yvh. (1.16)
Equivalently, the simultaneous solution to the finite dimensional problem given
by equation (1.16) is given by the corresponding matrix problem:
A u = fpE, (117)
where the finite element matrix A of (1.17) has entries [A]8J defined by
[A]tJ = A(^, fr), (1.18)
and the vector u of (1.17) contains the nodal values of the finite element solu
tion and has entries [u]j given by
R = y. (i.i9)
16
As a result of the definitions of uh in (1.14), A in (1.18), and u in (1.19), we
see that the ith entry of the vector formed by the matrix/vector inner product
on the lefthand side of equation (1.17), [Au],, is equivalent to the lefthand
side of equation (1.16) for the ith test function:
[An], = A(^,,) (120)
jei
Similarly, we define the ith entry of the finite element righthand side vector
of (1.17)[fpE]ias the righthand side of equation (1.16) for the ith test
function: recalling the definition of fh in (1.15), we have
[fFE], = {fhAi) = Â£ (121)
jei
The principal advantage of the finite element approach is its mathe
matical generality: FE triangulations can be fitted to very general geometries
and boundary conditions; the trial and test spaces defined on those triangula
tions provide accurate and robust approximations to the unknown u and data
/. However, the mathematical generality of FE can be a weakness as well as a
strength: the computational overhead required to set up the FE matrix system
of equations can be expensive compared to the set up costs for FV on general
computational meshes. Furthermore, this overhead is relatively fixed: the set
up cost for FE is largely problem independent for a general meshrelatively
simple diffusion problems require nearly as much work as complex advection
reactiondiffusion problems. From the physical standpoint, the boundary inte
gral (1.12) implies that FE will match the form of the governing conservation
law globally on 0 but not on any subdomain of 0; hence an important char
acteristic of the problemlocal conservationhas been lost.
17
1.6 The Finite Volume Element Method (FVE)
To realize the local conservation property of the finite volume method
within a finite element framework, we turn to the finite volume element method
(FVE)see Baliga and Patankar [4] and McCormick [39] for fundamentals. In
Steps 1 and 2 of Section 1.3, we follow the a combined FE and FV method
ology: in Step 1, we construct a FE triangulation with triangular or rect
angular elements, its pfold refinement T^, and a dual mesh Vh/p that is tied
to 'Thlp rather than Th as in FVthe importance of this distinction will be
discussed below; in Step 2, choose S(Vh/p) = Xh as a test space. After in
tegration by parts, FVE has the same weak form as FVcf. equations (1.8)
and (1.9)which preserves a finite dimensional analog of the local conservation
property of the integral conservation law (1.4). However in Step 3 of Section
1.3, we break away from the usual FV discretization of the unknown u and the
data / and turn to the usual FE discretization: replace u and / with repre
sentations uh and fh in the trial space Sp(l~h) = Wh as in equations (1.14)
and (1.15). We can do this because the volumization Vh/p establishes a oneto
one correspondence between the vertices of Thlp which are the DOF points for
Sp(Th). Hence, the use of the FE trial space is equivalent to systematic and
generalized FVstyle interpolations or extrapolations between nodes of Thlp.
Regarding boundary conditions, the general remarks for the finite
element method also apply to the finite volume element method. As in the FE
development of the previous section, well assume a homogeneous flux boundary
to simplify our presentation. Further discussion of boundary conditions and
their impact on our FVE analysis will be deferred to later chapters.
18
Referring to equation (1.9), the finite volume element method solution
for equation (1.1)or perhaps more appropriately equation (1.4)is described
as follows: find uh G Wh such that
B(uh,Xi) = (fh,Xi), (1.22)
where fh satisihes an L2 projection relationship with /,
(/\.Y,) = (/,X'.), Vx, ,t'\ (1.23)
to ensure that the numerical source term conserves mass on each volume in
Vh/p. Equivalently, the simultaneous solution to the finite dimensional problem
given by equation (1.22) is given by the corresponding matrix problem:
Bu = fFVE, (124)
where the finite volume element matrix B of (1.24) has entries [B]8J defined by
[B }ij = B(4>j,xi), (125)
and as before the vector u of (1.24) contains the nodal values of the finite
volume element solution and has entries [u]j of the form given by (1.19). As a
result of the definitions of uh in (1.14), B in (1.25), and u in (1.19), we see that
the Ah entry of the vector formed by the matrix/vector inner product on the
lefthand side of equation (1.24) is equivalent to the lefthand side of equation
(1.22) for the Ah test function:
[Bu], = $>jÂ£(^,Xt) = B{uh,Xi) (126)
jei
19
Similarly, we define the fth entry of the finite volume element righthand side
vector of (1.24)fpvEas the righthand side of equation (1.22) for the fth
test function: recalling the definition of fh in (1.15), we have
[fFVE]* = (/\ X*) = (fyiXi) (127)
jei
The principal advantage of the finite volume element approach is that
its FV weak form (1.8) mimics the integral conservation law (1.4) while incor
porating the mathematical generality of FE for very general geometries and
boundary conditions. From the FV perspective, FVE represents a systematic
and general approach to generating finite difference stencils that is only slightly
more complex than the standard FV approach. While from the FE perspec
tive, FVE represents a locally conservative variant of the standard FE by the
choice of a simpler, lower order, volumeoriented test space.
From the standpoint of application, FVE favors its FV parent as it
agrees with FV in Steps 1 and 2 of its formation. In Step 3, FVE can build upon
FV by using a FE trial space rather than interpolate or extrapolate between
nodes of 'Th. From the standpoint of analysis, FVE favors its FE parent as
it can be viewed as a PetrovGalerkin variant of the standard finite element
method. In classical finite volume analysis (see Mitchell and Griffiths [36],
Kreiss et al. [31], Manteuffel and White [35], and the references listed therein)
FV is viewed as a special type of finite difference method and analyzed in terms
of finite difference truncation error in C norms. However, in modern finite
volume analysis (see [5, 17, 18, 21, 23, 32, 33, 38, 46, 50, 54]) FV is viewed as a
special type of finite element method and its error is analyzed in terms of finite
element approximation theory results: results for the error in approximating
20
elements of a fractionalorder Sobolev space with piecewise polynomials. Since
FVE is even closer to FE than FV, our analysis will follow that of FE even
more closely.
Once we have established some theoretical infrastructure in the next
chapter, well come back to the theoretical similarities between FVE and FE.
But for now, we need to close the book with regard to the basic description of
the methods by briefly discussing the historical development of the computa
tional fluid dynamics and the roles played by FV, FE, and FVE.
21
1.7 Historical Development
Historically (see Batchelor [6]), the local and global conservation law
(1.4) was the starting point for theoretical fluid dynamics in the 18th century;
later with the development of continuum mechanics, it was shown (under suit
able conditions) that (1.4) implied (1.1). In the current century (see Roache
[43], Shashkov and Steinberg [48]), Unite difference numerical methods based
on (1.4) have enjoyed great success in a wide range of application problems in
fluid dynamics.
Starting in the 19th century (see Guenther and Lee [20]) with the
development of the calculus of variations and functional analysis, attention
in mathematical and physical circles shifted from local conservation laws to
integralbased variational principles that (under suitable conditions) lead to
equations like (1.1) in fluid dynamics as well as in other branches of contin
uum mechanics. In the current century (see Johnson [28] and Oden [40]),
Unite element numerical methods (discussed in Section 2.6) based on integral
variational principles have succeeded in cases where Unite differences have en
countered difficulties and even failures.
Coupled with a sophisticated yet elegant analysis (see [11, 12, 13,
14, 55]) that is one the crowning jewels of modern applied mathematics, it ap
peared likely that Galerkin Unite element methods would largely supplant finite
differences by the end of this century. However, this has not happened in fluid
dynamics and is not likely to happen in the future. Somewhat superficially,
the issue of computational efficiency has been problematic: the computational
22
setup costs for finite elements are relatively expensive when compared to fi
nite differences on general computational grids that discretize domains with
complex geometries.
More fundamentally, the problem lies in the choice of the underlying
model for computation. Standard finite volume or cellcentered finite difference
methods (FV) rely on the local and global conservation law (1.4), which is
always physically valid even when extrapolation to (1.1) is not. Standard
Galerkin finite element methods (FE) rely on an integral formulation of (1.1)
arising from variational principles and are not locally conservative; however,
FE is globally conservative. Hence, FE violate the spirit and the letter of the
fundamental local conservation laws that govern fluid dynamics.
Furthermore, the difficulties that arise in FV are not due to its con
servation law model but are usually due to somewhat ad hoc discretizations
of the unknown solution, known data, and boundary conditions for problems
posed on difficult geometries. Conversely, the successes of FE are largely due to
its systematic discretizations based on piecewise polynomial spaces defined on
FE triangulations for general geometries and not due to its variational model.
At the same time, variational principles are fundamental (and con
servation laws are not) in related disciplines of continuum mechanics such as
structural mechanics and electrostatics which are based on a principle of vir
tual work. In those applications, FE is the preferred numerical method (see
Babuska and Osborn [2] and the references cited therein).
In FVE, we try to combine the best of FV and FE. First, it is the lo
cal and global conservation law (1.4)and not a FE weak form of (1.1)that
23
we mimic discretely in FVE. Second, FVE approximates (1.4) using meshes
described in the next two sections, by replacing u and / in (1.4) with finite
element approximations uh and fh which are based on a finite element trian
gulation Th that partitions 0, and by posing equation (1.4) on a finite subset
Vh/p of volumes that partitions 0.
FVE is a fullfledged finite element method that is based on a local
conservation model rather than a variational model; but at the same time,
FVE can recover more computationally efficient equations like those of FV
with the use of certain reduced or lumped quadratures (see [5]). In general,
FVE enjoys some of the computational advantages of FV when compared to
FE. In FE, piecewise polynomial trial and test spaces of degree p > 1 are
used to discretize the unknown u and to setup FE equations, respectively.
Therefore, quadrature rules that are exact for polynomials of degree 2p are
used to evaluate multiple integrals in FE. In FVE, piecewise polynomial trial
spaces of degree p and piecewise polynomial test spaces of degree 0 are used to
discretize the unknown u and to setup FVE equations, respectively. Therefore,
quadrature rules that are exact for polynomials of degree p are used to evaluate
multiple integrals in FVE. If we let d be the number of spatial dimensions in
a given problem domain, we see that FVE quadrature will be 0(2d) times
less expensive than FE quadrature for a given problem (see Stroud [49] for
a discussion of quadrature for multiple integrals). For practical engineering
problems in three dimensions, the computational setup costs for FVE are
nearly an order of magnitude less expensive than FE.
24
1.8 Summary
In this chapter, we have given the necessary introduction for the fi
nite volume element method (FVE) for elliptic equations as a prelude to an
analysis of the method. In Section 1.2, we outlined how FVE arose from the
approximation of certain integral conservation laws that led to elliptic partial
differential equations. In Section 1.3, we detailed a three step procedure for the
construction of numerical methods for elliptic equations. In Sections 1.4, 1.5,
and 1.6, we compared and contrasted three specific numerical methods that
arose from our three step procedure: the finite volume method (FV) in Section
1.4, the finite element method (FE) in Section 1.5, and most importantly, the
finite volume element method (FVE) in Section 1.6. Finally, a brief outline of
the historical development of computational fluid dynamics and roles played
by the three methods was given in Section 1.7.
25
2. Theoretical Preliminaries
In this chapter, we make the discussion started in the previous chapter
more precise by providing some theoretical preliminaries: construction of the
FVE triangulations and volumizations, definitions of h, p, and hp version FVE,
a comparison between FE and FVE versions, Sobolev spaces to characterize the
exact solution to the model elliptic equation (1.1), and an outline of standard
FE analysis.
2.1 Triangulation
The most basic component in the implementation and analysis of
the finite volume element method is the discretization of the domain 0 into
computational meshes. The finite element triangulation of 0 and the finite
volume element volumization Vhlp of 0 are two differentyet interconnected
meshes or discretizations of 0. The two discretizations are connected by a
third mesh or discretization of 0: T^, the pfold refinement of Th The
nodes or vertices of 'Thlp are the locations of the degrees of freedom (DOF)
for uh G Sp(Th): uh is a C function on 0 that is a polynomial of degree
p > 1 when restricted to each element T of 'Th. In FVE terminology, the
triangulation Thlp is the primal mesh and the volumization Vhlp is the dual
mesh. Next, we describe these three meshes and their relationships.
26
2.1.1 FE Triangulation
Let Th be a nonoverlapping triangulation of 0, the closure of 0,
into a finite number of elements T. That is, satisfies the following general
properties [11]:
(1) For each T Â£ 'Th, the set T is closed and its interior, int(T), is
nonempty and connected;
(2) For each T Â£ T^, the boundary dT is Lipschitzcontinuous;
(3) Th partitions 0: i.e., 0 = (JTeTh T
(4) The elements of do not overlap: for each distinct Ti,T2 Â£ l~h, we
have int(Ti) fl int(T2) = 0.
(5) Any face of any Ti Â£ is either a subset of the boundary dfl or a
face of another T2 Â£ 'Th.
To simplify the discussion, we assume the elements of the triangulation are
triangles (possibly curved to accommodate nonpolygonal domains). Other
types of elements could be considered: in particular, rectangular elements for
a rectangular domain 0.
For each T of T^, we define the mesh parameters hy,h, and pj as
follows: hj is the diameter of the circumscribing circle for T; h is the maximum
value of hj; and py is the diameter of the inscribed circle in T.
We assume that the triangulation is shape regular: there exists a
positive constant a such that
<<7, V T Â£ ; (2.1)
Pt
27
the family (hx) is bounded and 0 is its unique accumulation pointi.e., h
approaches zero. The assumption of regularity is used to simplify error esti
mates and to avoid degenerate triangulations (i.e., a minimum angle condition
is satisfied). See Figure 2.1 for an example of the triangulation Th.
2.1.2 FVE Primal Mesh
Once a nodal or base triangulation is defined, we can define ele
ments of Sp(Th). The degrees of freedom (DOF) for uh Â£ Sp(l~h) are located
in a regular fashion in T: on the vertices, along edges, and in the interior of
T. Just as the FE triangulation can be determined by connecting the ver
tices of T to their nearest neighbors, an alternative triangulation Th^p can be
determined by connecting the DOF in T to their nearest neighbors within T:
let Tp denote the triangular elements of See Figure 2.2 for an example
of the triangulation J~h/p restricted to a single element in Th
To be specific, 'Thlp is generated from by the newest vertex bisec
tion algorithm developed by Sewell [47] and applied by Mitchell [37] to define
mesh hierarchies and the corresponding hierarchical basis functions used in the
Multigrid iterative solution of pversion FE problems. This refinement strategy
is very similar to the FAC (Fast Adaptive Composite grid) algorithms that were
successfully applied to hversion FVE in the monograph of McCormick [39].
Although we will not discuss the iterative solution of pversion FVE problem,
we use this mesh construction to define a volume mesh in Section 2.2.1, and to
define p and hp version FVE in Section 2.3. At the same time, we wish to set
the stage for the efficient numerical implementation p and hp version FVE at
a later date.
28
Figure 2.1: An example of the triangulation 'Th.
29
Figure 2.2: 'Thlp restricted to a single element of Th
30
Note that each element in 'Thlp is contained in a single element of
Th. The primal mesh 'Thlp can be seen as the pfold rehnement of Th: if
diam(T) = h, then diam(Tp) = h/p. Therefore, we write
yh/p g yh (2.2)
as a shorthand notation for this mesh hierarchy or nested property of the
triangulations.
31
2.2 Volumization
The FVE dual mesh Vh/p or volumization of 0 partitions 0 into a
finite number of nonoverlapping elements Rthe index i refers to a oneto
one correspondence between volumes and DOF for uh Â£ Sp(Th). To define Vh!p
from the primal mesh 7~h!p ^ we make frequent reference to 'T^v\ the union of
all Tp Â£ Thlp that have the location of the fth DOF, Xq as a vertex. In general,
Vh/p satisfies the following properties:
(1) For each V Â£ Vh^p, the set V is closed and its interior, int(V), is
nonempty and connected;
(2) For each V Â£ Vh//p, the boundary dV is Lipschitzcontinuous;
(3) Vh/p partitions 0: i.e., 0 = UveVh/pV
(4) The elements of Vh!p do not overlap: For each distinct Vi,Vj Â£ Vh^p,
one has int(V) 0 int(Vj) = 0.
(5) Any face of any V Â£ Vh/p is either a subset of the boundary dil or a
face of another Vj Â£ Vh^p]
(6) There is a onetoone correspondence between primal mesh nodes and
dual mesh volumes: x4 Â£ Vi and Vi C ', V i Â£ /.
Properties 1 through 5 are satisfied by any triangulation. However, Property
6 distinguishes the volumization Vh!p from triangulation Th/p and hints at the
concrete relationship between the FVE primal and dual meshes that will be
detailed below.
2.2.1 Volume Construction
To make our discussion more concrete, we follow the general volume
construction algorithm of Bank and Rose [5]:
32
(1) Select a point P G Tp, V Tp G T^V]
(2) Connect P by straight line segments to edge midpoints of Tp for the
two edges of Tp adjacent to the vertex Xy VTp Â£
(3) For each Tp Â£ , define a subvolume, Vi(Tp), as the region bounded
by the line segments formed in Step 2 and line segments connecting the
edge midpoints of Tp with the vertex x4. Note that Vi(Tp) = Vj fI Tp;
see Figure 2.3.
Finally, define the volume Vj as the region enclosed by all subvolumes
formed in Step 3:
Vi = U
TpeTp/p
Let vh/'p denote the set of all subvolumes; then, vh/p is a refinement
of or nested in yh/p.
Vh/P g yh/P' (2.4)
More importantly, each subvolume Vi(Tp) of vh/p is contained in a single ele
ment Tp of ,phlp; therefore, vh/p is nested in Th/p and in Th\
Vh/P g rh/P ^ jh' (2.5)
In our FVE analysis in Chapter 4, vh/p will be the focus of a local error analysis
for interpolation defined on elements of ph!p and then restricted to elements
of vhlp on an elementbyelement basis. Due to the nesting of the subvolumes
in the triangulations, we obtain error estimates on Th by summing over results
on vh!p.
33
Figure 2.3: Construction of the subvolume Vi(Tp).
34
The choice of P in Step 1 is crucial in volume construction. In practice
(see Cai [10]), we typically use one of the following volume types (characterized
by P):
(1) Circumcenter volume: P is the center of the circle circumscribed about
Tp; equivalently, P is the point of intersection for the perpendicular
bisectors of the edges of Tv. To ensure that P Â£ I), we require that
no interior angle of Tp exceeds tt/2 radiansno obtuse triangles are
permitted.
(2) Centroid volume: P is the center of mass or centroid of Tp; equivalently,
P is the point of intersection of the medians of Tp. A median of Tp
is a line segment connecting a vertex of Tp with the midpoint of its
opposite side.
(3) Incenter volume: P is the center of the circle inscribed in Tp; equiva
lently, P is the point of intersection of the angle bisectors of Tp.
(4) Orthocenter volume: P is the point of intersection of the altitudes of
Tp. An altitude of Tp is the line segment from a vertex of Tp that is
perpendicular to the line containing the opposite side. To ensure that
P G Tp, we require that no interior angle of Tp exceeds tt/2 radiansno
obtuse triangles are permitted.
In applications, advantages of circumcenter volumes are that they are always
convex and geometrically simple (see Figure 2.4 for an example), while other
types of volumes are usually nonconvex and geometrically complex (see Figure
2.5 for an example and Bank and Rose [5] for details).
35
Figure 2.4: A circumcenter finite volume Vfi
36
Figure 2.5: A general finite volume Vfi
37
2.2.2 Volume Symmetry
In the analysis to follow, superconvergence results for linear finite el
ements can be demonstrated for symmetric circumcenter volumes. To define
this symmetry precisely, we need additional notation: W7 is the edge or line
segment connecting nodes x4 and Xj of the triangulation 77 is the in
terface or volume boundary between volumes Vj and Vj of the volumization
Vh/ple. 7ij = V pi Vj. Following Cai et al. [9] and Cai [10], the volumization
yh/p
is symmetric to the triangulation 'Thlp if the following two symmetries
hold for all volumes in yh/p.
(Wsymmetry) 77 is a perpendicular bisector of W7;
(7symmetry) is a perpendicular bisector of 77.
We remark that for rectangular elements, only the first condition of
Wsymmetry is required for yh!p to be symmetric to 'Thlp (cf. Siili [50]).
2.2.3 Volume Regularity
Of more general significance is the notion of volume regularity,
which is assumed in all of our results (cf. Cai [10]). Form an auxiliary tri
angulation Th/p of triangular elements by connecting the endpoints of 77 with
the endpoints of W7 for every volume interface in yh!p. That is, we define a
macroelement C7 as the region formed by joining the endpoints of 77 and
Xij (See Figure 2.6 for an illustration). Then C7 is the union of two triangular
elements (with W7 as a common edge): i.e, C7 = T)} U T)}, where each element
Tj,k G {1,2}, is contained in a single element of 'Thlp. Therefore, 'Thlp is a
refinement of 'Thlp and diam(T^) = 0(h/p) and we have the following nesting
38
of triangulations:
fh/v c rh,p c rh. (2.6)
In our FVE analysis in Chapters 3 and 4, Th/p will be the focus of a local error
analysis for interpolation defined on elements of 'Thlp and then restricted to
elements of Th/p on an elementbyelement basis. Due to the nesting of the
triangulations, we obtain error estimates on Th by summing over results on
'fhjv.
For circumcenter volumes (the focus of superconvergent error esti
mates in Chapter 3), if this auxiliary triangulation Th/p is shape regular ac
cording to the definition in Section 2.1.1 (i.e., inequality (2.1) is satisfied) and
no interior angle of Tp Â£ J~h/p exceeds tt/2, then the volumization Vhlp is reg
ular. For the other volume types of Section 2.2.1, shape regularity of Thlp
implies regularity of yh!p. See Figure 2.6.
39
Figure 2.6: Construction of a pair of elements in Th!p.
40
2.3 Versions h, p, and hp
Motivated by the success of the h, p, and hp version FE of Babuska
in structural mechanics problems, we define and analyze h, p, and hp FVE
for fluid dynamics problems under uniform refinements of h and p in direct
analogy to the FE development outlined in [2, 3, 27, 51].
In the hversion FVE, we fix the order p of the FE polynomial trial
space and uniformly refine the computational mesh 'Th. The monograph of
McCormick details mesh refinement and the efficient numerical solution of h
version FVE discretizations. For p = 1, the paper of Cai and McCormick [8]
analyzes local refinement of the hversion FVE for a steady diffusion equation;
similarly, Ewing et al. [17, 18] and Lazarov et al. [32] analyze local refinement
for hversion FV for general elliptic and parabolic problems.
In the pversion FVE, we fix the computational mesh and uni
formly increase the order p of our FE trial space, which in turn refines the FVE
meshes Th!p and yh!p. In the hp version FVE, we refine uniformly in h and p.
As in p and hp version FE (see Szabo and Babuska [51] for background), hier
archical basis functions could be used which induce a mesh hierarchy for Th!p
and Vh/p in the p and hp version FVE as discussed in Mitchell [37]. However,
our focus here is on an a priori error analysis, so it is sufficient that we have a
set of basis functions (not necessarily hierarchical) on each element of that
is defined in terms of the vertex DOF on Thlp. So for simplicity, we can assume
that a standard (nonhierarchical) nodal basis is being used. However, in an a
posteriori error analysis and in realworld implementations, hierarchical basis
functions are essential to an efficient iterative solution (see Hu and Katz [25],
41
Mandel [34], Mitchell [37], and Pavarino [42] for background). To this end, we
note that our definitions of 7~h^p and Vh/p in the previous two sections were
stated with efficient implementation in mind and fully support basis hierarchy
according to the newest vertex bisection algorithm discussed in [37].
Although our error analysis in the remainder of the text is inde
pendent of any particular solution algorithm for FVE discretizations, we are
motivated and strongly influenced by the success of Multigrid and Multilevel
iterative algorithms (MG) in numerical partial differential equations (see Mc
Cormick [39] and the references cited therein for background). Out of the stan
dard hversion FV/FD, FE, and FVE discretizations, McCormick [39] found
that FVE was the most amenable to efficient MG solutions. Since the success
of MG and Domain Decomposition solution algorithms [25, 34, 37, 42] for p
and hp version FE discretizations is either directly based on or closely related
to the ideas presented in [39], we can reasonably expect similar success for p
and hp version FVE. In fact, it was the success of these methods for FE that
prompted the analysis of p and hp version FVE in this thesis. If the perfor
mance of FVE vs. FE for p and hp methods parallels the MG performance of
the hversion methods, FVE may be superior. Some anecdotal support for this
position was stated in Section 1.2: an operation count for the costs of equa
tion setup in FVE vs. FE (i.e., numerical quadrature of multiple integrals)
revealed that FVE was nearly an order of magnitude less expensive for three
dimensional problems. Further support for this position is given in the next
section.
42
2.4 FVE vs. FE
In Sections 1.1 and 1.2, we outlined the origin of fluid dynamics and
its basis on local conservation laws. Also, we motivated the need for FVE
which preserves local conservation on the discrete level. In addition, we saw
that FVE systematically discretizes the unknown solution, known source term
data and boundary conditions on domains with complex geometries via the
FE triangulation Th, the FVE triangulation 'Thlp, and the FVE volumizations
Vhlp. From this perspective, we can view FVE as a locally conservative Petrov
Galerkin FE. This prompted us to define h, p, and hp version FVE in direct
analogy to h, p, and hp FE in the previous section. In addition to local
conservation and the variety of versions, FVE has numerous other practical
advantages over FE that are outlined below.
In Section 1.2, we saw that the computational setup costs for FVE
were less expensive than the costs for Galerkin FE. If d is the spatial dimension
of the problem domain, the quadrature rules used to evaluate FVE integrals
and setup a system of linear equations were seen to be 0(2d) less expensive
than corresponding FE quadrature rules. That is, the quadrature operation
count for reaction and source term integrals is 0(p)d for FVE and 0(2p)d for
FE, since the order of polynomial trial and test space is p and 0 for FVE and
p and p for FE.
Furthermore, for reaction and source terms, the FVE integrals from
(1.4) are defined on volumes from Vh!p with diameters of 0(h/p) (see Section
2.2), whereas FE integrals are defined on elements from Th with diameters
of 0(h) (see Section 2.1); therefore, the FVE quadratures are performed on
43
regions that are 0(p)d smaller than their FE counterparts. This feature makes
the FVE equations more local in character and the FVE mass matrix arising
from the reaction term more strongly diagonally dominant than its FE coun
terpart; from the standpoint of numerical linear algebra and implementation,
this means that it should be easier to precondition and iteratively solve FVE
systems of linear equations than their FE counterparts.
In addition, FVE diffusive and advective fluxes are defined as sur
face integrals in (1.4) which are simpler (i.e., of lower dimension) than the
corresponding volume integrals for FE: i.e., the quadrature operation count
is 0(p)d~1 for FVE and 0(2p)d for FEa increased savings of 0(p2d). Cou
pled with the documented [39] success of Multilevel solution algorithms for
hversion FVE vs. FV discretizations and the documented [25, 35, 37, 42] suc
cess of similar Multilevel and Domain Decomposition solution algorithms for p
and hp version FE discretizations discussed in the previous section, we expect
FVE to outperform FE in standard measures of computational efficiency.
Furthermore, from the standpoint of the a priori error analysis pre
sented in the later chapters, we fold that the accuracy of FVE approximations
to the exact solution of general elliptic and parabolic problems will either match
or surpass FE asymptotic convergence rates for h, p, and hp versions of the
respective methods.
2.5 Sobolev Spaces
In this section, we present the Sobolev space infrastructure for our
FVE analysis. In the basic definitions and results, we have followed the well
known and comprehensive text of Ciarlet [11] unless otherwise noted. For a
44
more detailed and advanced discussion of Sobolev spaces, we refer the reader
to the authoritative text of Adams [1] for mathematical theory and to the text
of Kesavan [30] for applications to partial differential equations.
Throughout this section, 0 C $ld,d Â£ {1,2,3}, denotes a domain:
a bounded open connected set with Lipschitzcontinuous boundary T. We
also use the multiindex notation for denoting partial derivatives: let a multi
index a be denoted by a vector with nonnegative integer components, a =
(!,. . ay) Â£ Nd, equipped with the L1 vector norm, a = J2i=i i Taking
d = 2 for definiteness and letting x and y denote the spatial coordinates of $ft2,
we define dau as
dau = d^d^2u, (2.7)
so that dau represents the composition of the oqth and the a2th partial deriva
tives with respect to x and y.
2.5.1 IntegerOrder Spaces
For any integer k > 0, and number q Â£ [1, oo], the integerorder
Sobolev space Wk,q(Q) consists of those functions u Â£ Lq(fl) for which all
partial derivatives dau (in sense of distributions) with a < k belong to the
space Lq( fl). The space Wk,q(Q) is equipped with the norm:
 k,q,Â£l <
( Jn \dau\qd^y/q} g Â£ [1, oo),
\a\
maxjess supd"u(a:), q = oo.
(2.8)
The seminorm  \k,q,a is obtained by taking a = k in (2.8).
The space Wk,q(Q) is a Banach space and a Hilbert space when q = 2.
45
In the latter case, we denote Wk,2(Q) by Hk(0) and  11*,2,n by  fc,n If
u G Hk(0) for some k > 0, we say that u is kregular. In our analysis, we will
be using the Hk(0) spaces almost exclusively and will cite additional Sobolev
space results (e.g., fractionalorder spaces, trace theorem, embedding, etc.) in
terms of this q = 2 case to simplify notation. These results have generalizations
to q ^ 2 (see Adams [1] or Kesavan [30]), but they are not required in our
forthcoming work.
2.5.2 FractionalOrder Spaces
In addition to the usual integerorder Sobolev spaces, we also employ
fractionalorder Sobolev spaces Hs(0), where S > 0 is a real number, as in
Adams [1] and Kesavan [30]. In this discussion, S is decomposed into integer
and fractional parts: i.e., S = k + 0, where k > 0 is the integer part of S and
6 G [0,1) is the fractional part of S. We define Hs(0) in terms of Hk(fl) and
He{ 0):
Hs{0) = [u G Hk{0) : dau G tf0(O), V a = k}, (2.9)
H(Sl) = {,, i2(f!) : l(^)y~(Jh>,1 i2(f! x f!)}. (2.10)
If u G Hs(0) for some S' > 0, we say that u is Sregular.
2.5.3 Sobolev Embedding
In our analysis, we approximate elements from fractionalorder Sobolev
spaces by C piecewisepolynomial (Lagrange) interpolants. To ensure that in
terpolation is welldefined, we will only use the simplest embedding:
Hs(0) C^fi), if S' > d/2, (2.11)
46
where means that Hs(0) is contained in (7(0) with a compact injection.
Referring to the embedding (2.11), we see that if u is (1 + e)regular
for any e > 0, then u is represented by a (7 function in $ft2. Therefore, u
has point values and we can easily define (7 Lagrange interpolation for u.
Similarly, 3/2 + e regularity is required in three dimensions.
2.5.4 Trace Theorem
We will often employ the continuous and linear trace map of an el
ement of a Sobolev space from an interior domain 0 C d > 2, to a sub
dimensional boundary F C 1 that is sufficiently smooth (see Kesavan [30]
for details). The trace map, tr(), maps u Â£ iL0+1/2(fi) to tr(u) Â£ H&(T),
where 0 is a positive, noninteger, real number (i.e., 0 Â£ $R+\iV). When no
confusion should arise, we simply write tr(u) as u. The following trace theo
rem demonstrates the boundedness of the trace map: there exists a positive
constant C = (7(0) such that
IMkr < C' ^+i/2^, 7 K+\N. (2.12)
As a result of (2.12), we see that u Â£ Hs+1(Fl) for S > 1/2 ensures that u Â£
Hl(F) and that the diffusive fluxes in (1.4) are bounded. With this regularity,
(7 interpolation is easily defined (up to three dimensions) as discussed in the
embedding result above.
2.5.5 Poincare Inequality
The Sobolev space Hq(0) is the closure of the space DiFt)the space
of infinitely differentiable functions with compact support in 0in the space
Hl(Fl). More intuitively, Hq(FI) = {u Â£ Hl(0) : tr(u) = 0 on T}. As a
47
result of the boundedness of 0, we have the Poincare inequality, there exists
a positive constant 6Y(0) such that
Therefore,
pi2 is a norm over Hq(0) equivalent to  i,m
(2.13)
48
2.6 Finite Element Analysis
Here we present the rudiments of FE analysis to establish a reference
point for analogous developments to follow for FVE. For simplicity, consider a
steady diffusion equation with homogeneous Dirichlet boundary condition:
V(TV) = /, ;rÂ£0, (2.14)
u\dn = 0, (2,15)
where D is a continuously differentiable, bounded, nondegenerate (0 < Dm <
D < Dm < oo) diffusion coefficient and 0 C iff2 is a domain.
A weak solution u Â£ W = Hq(D) to (2.14), corresponding to / Â£
L2{0), satisfies:
A(u}w) = (/, w), V w Â£ W, (2,16)
where A(u,w) = [ DVu Vw dx, Jsi (2,17)
if,w) = / fwdx. Ja (2.18)
In our analysis, we assume additional regularity for u: i.e., u lies in the
(fractionalorder) Sobolev space W+ = Hs+1(0) fl W, where S Â£ (1/2, p] and
p > 1. The interval (l,p+ 1] for S' + 1 is the full range of admissible regularities
for u in FE analysis for a general unspecified problem domain. (Technically,
the lower bound on S can reduced to 0 in two dimensions; however, 1/2 is re
quired for our results to hold in three dimensions.) But for any specific problem
domain, the value of S' + 1 in that range is precisely determined: S = S(d0)
is a function of the smoothness of the boundary <90 and the geometry of any
corners in nonsmooth domains (see Grisvard [19] for details). For example, if
49
0 is an Lshaped domain, S + 1 = 5/3 eVe > 0. To justify the upper end
of the regularity scale (S + 1 > 3) associated with highorder polynomial trial
spaces (p > 2), we assume that curved elements are being used to approximate
the boundaries of smooth domains (see Ciarlet [11] for details).
In the FE analysis presented here and in the FVE analysis of the
following chapters, we invoke fractionalorder spaces to demonstrate error es
timates across the fullest range of admissible regularities to encompass any
particular problem domain.
Also, the numerical solution uh of (2.14) (see (2.19) below) will lie in
Wh = Sl(Th) C Wthe space of C functions that are polynomials of order
< p on each element T of the triangulation and vanish on dfl.
In direct analogy with (2.16), the (finite element) numerical solution
uh Â£ Wh to (2.14), corresponding to fh Â£ Sp(Th), satisfies:
A{uh, w) = (/\ w), V re Â£ Wh. (2.19)
Here, fh is the T2(0) projection of / into Sp(Th):
(ffh,w) = 0, VweSp(Th). (2.20)
The key components of finite element error analysis are boundedness,
approximation theory, and ellipticity results; we introduce and discuss each
below.
2.6.1 Boundedness
For the bilinear form A(, ) of (2.17)hereafter called the Aform
consider the following upper bound (or continuity) condition: there exists a
50
positive constant M such that
A(u, w) < M ^i,ntyi,n, V G tV. (2.21)
Such a condition is usually assumed to obtain convergence of the numerical
method as seen below.
2.6.2 Approximation Theory
At the heart of FE analysis are error estimates for approximation of
elements from Sobolev spaces by elements of piecewise polynomial FE trial
spaces. For the Hk(0) norms (k Â£ {0,1}) that we shall employ, we cite stan
dard approximation theory results for polynomial interpolation (see Babuska
et al.[2, 3], Jensen and Suri [27], and Dupont and Scott [16]): there exists a
constant C such that
\\vP(u)\\k,n < C(h/p)s+1k\\u\\s+1,n, S Â£ (1/2,p], (2.22)
where r]p(u) = Tlku u is the interpolation error for u in the polynomial
space Sp(Th). As discussed above for an Lshaped domain, (2.22) holds for
S' = 2/3 e. In three dimensions, the lower bound on S in (2.22) is necessitated
by the Sobolev embedding result (2.11), 3/2 + e regularity is required to embed
into C(fi) and easily define interpolation.
In accordance with the definitions of Section 2.3, the term (h/p)s
in (2.22) is understood as C(p)hs for hversion methods and as C(h)p~s for
pversion methods. In the remainder of the text, C will represent a (generic)
positive constant.
51
2.6.3 Ellipticity
Consider the following ellipticity (or lower bound) condition: there
exists a positive constant (3 such that
A(u,u)>P\u\la, VuGW. (2.23)
This condition is perhaps the most important technical component of FE error
analysis: it implies that the Aform is positive definite. Once ellipticity is
established, FE analysis reduces to the approximation theory of the previous
subsection.
2.6.4 Error Analysis
In the standard FE error analysis to follow, we apply the results
(2.23) and (2.21)(2.22); for (2.14), (2.23) and (2.21) hold with (3 = Dm and
M = Dm, respectively.
2.6.4.1 Error Splitting
Let e = u uh Â£ W be the error in the FE approximation to (2.14).
For theoretical purposes, define ( = 11))e as an interpolation of the error from W
into Wh, so that ( = H^u uh Â£ Wh. We refer to ( as the representation error
in Wh, i.e., the difference between interpolated and computed approximations
in Wh of the exact solution u. Then, we employ the error splitting e = ( rjp,
in which the FE error is split into its representation error and interpolation
error components: u + ?7Pi,n. With the estimate of \riP\i}a in
(2.22), estimating e1is reduced to estimating Ci,m Below, we develop an
estimate for Ci,n in terms of
52
2.6.4.2 Zero Property
For all w Â£ Wh, the exact solution satisfies (2.16), the numerical
solution satisfies (2.19), and the source term satisfies (2.20). Subtracting, we
have the zero property or orthogonality of the error,
A(uuh,w) = 0, VweWh, (2.24)
which in our context is more appropriately expressed as
A((, w) = A{r]p, w), V w Â£ Wh. (2.25)
2.6.4.3 Error Estimate
Since the representation error, (, is in Wh, we take w = ( in (2.25)
to obtain
p\c\U < a(c,o,
= Mnp,0,
< M^linCi,n. (2.26)
The three lines of (2.26) use (2.23), (2.25), and (2.21), respectively. By (2.22),
(2.26) implies the following H1 optimalorder error estimate for the finite ele
ment method:
\uuh\^ < C{h/p)s\\u\\s+i^ S Â£ (1/2,p\. (2.27)
In the FVE analysis to follow, we obtain error estimates that match or surpass
the pthorder estimate (2.27).
53
2.7 Summary
In this chapter, we presented some theoretical preliminaries. In Sec
tions 2.1 and 2.2, we defined and discussed the computational meshes (i.e., the
FE triangulation Th, the FVE primal mesh 'Thlp, and the FVE dual mesh Vh!p)
that are fundamental to the implementation and analysis of FVE. In Sections
2.3 and 2.4, we defined h, p, and hp version FE and FVE and discussed their
relative merits. In Section 2.5, we introduced standard definitions and results
for Sobolev spaces that are fundamental to FE and FVE analysis. In Section
2.6, we briefly outlined an FE analysis for the steady diffusion equation to set
a reference point for an FVE analysis in the next two chapters.
54
3. Analysis for Diffusion Equations
In this chapter, we present an h, p, and hp version FVE analysis for
diffusion equations. In Section 3.1, we provide some preliminary definitions and
notation that facilitates the transition from FE to FVE analysis. In Section
3.2, we motivate our FVE analysis and its relation to previous FE work. In
Section 3.3, we define some discrete norm counterparts to the continuous L2
and H1 Sobolev norms of Section 2.5 and demonstrate a discrete Poincare
inequality. As in the FE analysis of Section 2.6, these tools are crucial to FVE
analysis. In Section 3.4, we determine an upper bound or continuity result
in terms of an FVE diffusion functional and our discrete H1 seminorm for
the FVE bilinear form that is used to discretize the model diffusion equation
discussed in Section 2.6. In Section 3.5, we state and outline applications of
the linear BrambleHilbert lemma which is at the heart of our FVE analysis.
In Section 3.6, we use the BrambleHilbert lemma to derive error estimates
for the FVE approximation of elements from fractionalorder Sobolev spaces
by elements from the piecewise polynomial trial spaces employed by the h,
p, and hp version FVE methods: in terms of our FVE diffusion functional,
i71(0)equivalent optimalorder approximation results are derived for pthorder
trial spaces in Section 3.6.1; robust and practical superconvergence results for
linear trial spaces on a symmetric volumization are derived in two and three
dimensions in Section 3.6.2. In Section 3.7, we derive an asymptotic discrete
ellipticity result for the FVE diffusion bilinear form in terms of our discrete
55
H1 seminormaided by the approximation theory from the previous section.
Finally in Section 3.8, we construct an FVE analysis for diffusion equations
with the components developed in the preceding sections: optimalorder and
superconvergence results are obtained for h, p, and hp version FVE in our
discrete H1 seminorm.
3.1 Preliminaries
We employ the same continuous and discrete space notations, defini
tions, and structures as in the previous section for the FVE continuous and
discrete trial spaces that represent the solution to (2.14): u Â£ W = Hq(0)
and uh G Wh = Sq^T*1), for p > 1. Similarly in FVE analysis, we as
sume additional regularity for u: u lies in the (fractionalorder) Sobolev space
W+ = Hs+1(0) fl W, where S Â£ (1/2,p]. The 1/2 lower limit on S en
sures: a bounded diffusive flux, due to the Sobolev trace theorem (2.12) as
discussed below; and welldefined interpolations in up to three spatial dimen
sions, due to the Sobolev embedding (2.11). As stated at the start of Section
2.6, (3/2, p + 1] is the full range of admissible regularities in our analysis for
an unspecified problem domain. For any specified 0, the value of S in that
range is a function of the smoothness of dQ [19]: e.g., on an Lshaped domain,
S = 2/3 eVe > 0. In the remainder of the text, e will represent a (generic)
infinitesimal positive constant. Also, to justify the upper end of the regularity
scale (S + 1 > 3) associated with highorder polynomial trial spaces (p > 2),
we assume that curved elements are being used to approximate the boundaries
of smooth domains (see Ciarlet [11] for details).
56
In addition, let Xh denote the FVE test space Pq (Vh^p)the poly
nomials that are constant on each volume V of Vh!p and vanish on <90. Thus,
test functions corresponding to volumes that intersect <90 are identically zero.
Also, let xv denote the characteristic function associated with volume V C 0
and let V denote the collection of all such V.
Recalling the integral equation (1.4), a weak solution u of (2.14) sat
isfies the following problem: given / Â£ T2(0), find u Â£ W+ such that
B(u,Xv) = (f,Xv), VVÂ£V, (3.1)
where
B(u,xv) = [ D'S/u n dS, (32)
J dV
{f,Xv) = jvf dx.. (3.3)
FVE poses the weak form (3.1) on the volumization Vh/pa finite
subset of V with nonoverlapping volumesand employs standard FE repre
sentations of u and /. Thus, FVE replaces (3.1) with: given fh Â£ Sp('Th), find
uh Â£ Wh such that
B(u\ xv) = (/*,;tv), VV6V. (3.4)
Here, fh satishes a nonstandard L2(fl) projection relationship with /:
(ffh,xv) = 0, VVÂ£V^, (3.5)
which ensures that the numerical source term conserves mass on each volume.
Since Vh/p partitions ft and the volumes in yh/p
are related to the
finite element triangulation Th through Th/p as discussed in Section 2.2, (3.4)
57
is more amenable than (3.1) to an FEstyle analysis. Below, we point out
similarities and differences between FVE and FE from the perspective of an
FEstyle framework for analysis.
58
3.2 Motivation
To motivate our FVE analysis and its relation to previous FE work,
reconsider (2.14)(2.15) and the FE analysis of Section 2.6. We saw that the
crucial ellipticity result (2.23) for the Aform and boundedness result (2.21)
are stated in terms of the continuum space W, but the results also held for
the discrete space Wh (i.e., Wh C W). Superficially, it could appear that
the focus was on the continuum problem (2.16) rather than on the FE discrete
problem (2.19). However, in Section 2.6.4, the error analysis emphasized results
for the discrete FE problem (2.19) and the associated discrete spaces. Recall
that the exact weak solution u was invoked only in the zero property (2.25);
furthermore, recall that the key observation in (2.25) was that u satisfied the
FE discrete equation (2.19). Therefore, to begin an FEstyle analysis for FVE,
we need to generate results analogous to those of Sections 2.6.1, 2.6.2 and 2.6.3
for the FVE discrete problem (3.4) rather than for (3.1).
The continuum problems (2.16) and (3.1) are different in physical
character (minimization vs. conservation) as well as in mathematical content
and these differences discourage a single, unifying framework for their analysis.
For the discrete problems (2.19) and (3.4), the situation is more favorable for
a unified analysis. We now manipulate (3.4) to be more reminiscent of (2.19).
If w Â£ Xh, w has a volumewise constant representation:
w = Y,wiXi, (3.6)
iei
where yy = \vt and Vh/p = {R : i Â£ /}. In terms of w, the FVE solution
uh Â£ Wh of (3.4) also satisfies:
B(u\ w) = (fh,w), V re Â£ Xh, (3.7)
59
where
B(uh, w) = J2 B(uh, Wi Xi), (3.8)
iei
(/\m)=E(AiXi) (39)
iei
That is, if uh satishes the each of the local equations of (3.4), then, by (3.7), uh
also satishes any linear combination of those equations. Notice that (3.8) and
(3.9) are now global in extent, like their FE counterparts (2.17) and (2.18),
and (3.7) with w = 1 represents a global conservation law on the union or
aggregate of all volumes in Vh/p. Like the FE Aform, we refer to (3.8) as the
Fhform.
Similarly, the L2(0) projection relationship (3.5) for fh Â£ Sp(l~h) is
transformed into
(// = 0, V rc Â£ Xh. (3.10)
Now we establish some additional theoretical infrastructure to sup
port a FEstyle analysis of FVE based on (3.7).
60
3.3 Discrete Norms
In this section, we define some discrete norm counterparts to the
continuous integerordered Sobolev norms of Section 2.5 and demonstrate a
discrete Poincare inequality. As in FE analysis of Section 2.6, these tools are
crucial to FVE analysis.
In order to distinguish between continuous and discrete norms in our
definitions, we replace the physical domain 0 that was employed in the defini
tion of the continuous Sobolev norms with the set of computational mesh points
io (i.e., the vertices of the triangulation 'Th/p') contained in 0 in our definition
of discrete Sobolev norms. Recall that the continuous norms were defined in
terms of integrated functions over 0; however, our discrete norms will be de
fined in terms of function point values in u. Therefore, only Sobolev spaces
that can be embedded into (70(Q), like Hs(0) for S > 1/2, are admissible in a
discrete norm analysis.
In particular, we will restrict our attention to certain polynomial sub
spaces of C0(O) when defining certain discrete Sobolev norms and seminorms:
S(yhlp) for the L2(u) norm and Sl(Th^p) for the Hl(u) seminorm. This
choice of polynomial spaces will emphasize and exploit the relationships be
tween the FVE dual and primal meshes, which are at the heart of any FVE
implementation, in our FVE analysis. Thus, discrete norms will provide for
a unified analysis of both practical and theoretical questions in our study of
FVE: attention will be shifted from estimating the error over the whole domain
0 to estimating the error over only the subset of points lu where our numerical
method is approximating the exact solution to the original partial differential
61
equation (1.1). Starting with Section 3.6, we find that FVE discrete norm
analysis reveals nodal superconvergence in multiple dimensions as well as the
usual optimalorder convergence rates and is superior to the FE continuous
norm analysis of Section 2.6.
3.3.1 Discrete L2 Norm
First, we define a discrete variant of the standard continuous L2{0)
normreferred to as a L2(u) norm.
Definition 3.1 The L2(u) norm is given by
Mo,*= (E^2IK)1/2, (3.11)
iei
where
io = {x8 : i G 1} is the set of primal mesh points;
I is the primal mesh point indexing set;
Ui is the value of u at the primal mesh point xy
Vi G Vh/p is the volume that corresponds to xy
\Vi\ is the standard Lebesgue volume measure of V{.
Clearly, the L2(u) norm is equivalent (cf. Bank and Rose [5]) to the
T2(Q) norm restricted to S{Vh^p)the space of piecewise constant functions
on the FVE dual mesh Vh/p:
Mo,w = no/p^o,n, (3.12)
where II hJPu denotes interpolation into S(Vh^p). Then, L2(uj) is the space of
functions (generating primal mesh point sequences) that are bounded in the
L2(uj) norm.
62
3.3.2 Discrete H1 Norm
As in the FE analysis of Section 2.6, the H1 seminorm plays a pivotal
role in the course of FVE analysis. By analogy, we define a discrete variant
of the standard continuous i71(0) seminormreferred to as an Hl(u) semi
norm.
Definition 3.2 The Hl(u) seminorm is given by
lL=( Â£ (.))2M)I/2, (3.13)
{i,j }ew
where
W = {{*, j} : i Â£ I, j Â£ Wi with \^i3\ > 0} is a set of unordered index
pairs corresponding to the edges of each element Tp Â£ 7~h/'p]
Wi = {j : Xj Â£ Af(f), j Â£ /} is a set of indices for x4s nearest neighbor
mesh points;
. Af(i) is a set of x4s nearest neighbor mesh points in cj;
7ij is the standard Lebesgue surface measure of the volume interface
7ij = Vi fl Vj; in 1D, replace 17^y with 1;
\Xij\ is the Euclidean distance between x4 and Xj.
Despite outward appearances, the Hl(u) seminorm is equivalent (cf.
Bank and Rose [5]) to the standard Hl(0) seminorm restricted to (S1(T/l/,p)
the space of C piecewise linear polynomial functions on the FVE primal mesh
where H^pu denotes interpolation into S1('Thlp). Then, i71(cu) is the space of
functions that are bounded in the Hl(u) norm:
IMIo,* = (Mo,w + Mi,S) (3.15)
To relate o,w and more directly, we demonstrate a discrete Poincare
inequality.
3.3.3 Discrete Poincare Inequality
Following the continuous case introduced in Section 2.5, let 7 denote
the set of primal mesh points that lie on the domain boundary T. By analogy,
define Hq(u) as the set {u Â£ Hl(u) : u = 0 on 7}. Then the boundedness of
0 translates to a finite number of mesh points in lu and to the following result.
Lemma 3.1 There exists a positive constant C such that
\u\o,u < c \u\iyUJ, V u Â£ Hq(u). (3.16)
Therefore,  ijW is a norm over u) equivalent to  i)C^
Proof: Lemma 2 of Bank and Rose [5] demonstrates that o,S7 A
C in^itkn Since U^pu Â£ S^T11^) C Hq(0), the continuous Poincare in
equality (2.13) yields inj^rtkn < C 11 Therefore, application of the
equivalences (3.12) and (3.14) gives us the desired result (3.16).
As in FE error analysis, the key components of FVE analysis are
boundedness, approximation theory, and ellipticity results: we demonstrate
and discuss each result below. All follow the basic FEstyle framework intro
duced in Sections 2.6.12.6.3 after discrete norm modifications.
64
3.4 Boundedness
To determine an upper bound for the FVE Hform (3.8), we are guided
by the corresponding development for the FE Aform (2.17). In Section 2.6.1,
the upper bound (2.21) for the Aform is given in terms of a product of standard
integerorder Sobolev norms (or sublinear functionals). In the error analysis
of Section 2.6.4, these norms are applied to a continuous trial function and a
discrete test function. Since this upper bound is used in conjunction with the
ellipticity result (2.23), the only theoretical constraint on the bound is that the
norm applied to the test function must match (or can be made to match via
an auxiliary result like the Poincare inequality) the norm used in the ellipticity
resulti.e., I'li^ for FE and for FVE.
In FE analysis, the 77q(0) norm constraint on the test function lead
to the choice of the 7/^(0) norm for the trial functions. However, in FVE
analysis, the Hq(u) norm constraint on the test function leads to the choice
of a problemdependent, bounded functional for the trial functions that cor
responds to the FVE treatment of diffusion in (1.4). In Section 4.2, we will
develop functionals that correspond to FVE treatments of reaction and ad
vection in (1.4). Practically speaking, we want this diffusion functional (and,
later on, reaction and advection functionals) to admit optimalorder or even
superconvergent approximation theorems, since the continuous trial function
in our analysis is the interpolation error for the exact solution.
In terms of the admissible trial space, = W+ Wh, and test space
Xh for FVE, we demonstrate an boundedness result for the Hform (3.8) (cf.
(2.21) of Section 2.6.1).
65
Lemma 3.2 (Diffusion bound)
\B(u,w)\ < T>(u)\w\i^, V u Â£
(3.17)
Z>() is a bounded functional:
() = ( E 4
(3,18)
pj}ew I'd I
dij() is a bounded linear functional arising from the FVE treatment of diffu
sion:
dij(u) = f DVu iiijdS, (3.19)
h,3
where n8J is the unit normal vector pointing outward from Vi into Vr
Proof: By the Sobolev trace theorem, u Â£ Hs+1(0) for S > 1/2
ensures that (3.19), hence (3.18), is bounded. From the definition (3.8) of the
5form and the fact that
Id V;
DVu nidS = dij(u),
jewt
we see that
(3.20)
\B(u,w)\ = J2 d^(u)\>
iei jew,
=  dtJ(u)(wtwJ)j,
{hj Jew
< ( E 4()rr)/2( E
{i,j }ew H'dl {i,j}ew l^dl
where we have used the fact that dij(u) = d^fu) (i.e., n8J = np) in the
second line of (3.21).
66
With an upper bound for the hbform in terms of the nonstandard,
problemdependent, bounded linear functional (3.19), we need an approxima
tion theorem: the theoretical basis for this result is introduced and discussed
below.
67
3.5 BrambleHilbert Lemma
In this section, we state and outline applications of the linear Bramble
Hilbert lemma. This lemma was originally developed by Bramble and Hilbert
[7] for use in approximation theory. Later, the lemma was modified for use
in FE theory by Ciarlet [11] and applied in linear and bilinear forms. In
the standard FE analysis presented in Section 2.6, this lemma only plays a
supporting rolea bilinear BrambleHilbert lemma [11] is employed to analyze
the effects of numerical quadrature. However, in modern FV theory the lemma
plays a central role. In fact, the particular form of the BrambleHilbert lemma
stated here was used in the pioneering FV analysis of Samarskii, Lazarov,
and Makarov [46]. Similarly, the BrambleHilbert lemma plays a central role
in our FVE error analysis: the lemma is used to demonstrate optimalorder
and superconvergence error for bounded linear functionals that arise in the
implementation and analysis of the finite volume element method.
First, we introduce some additional infrastructure for bounded linear
functionals. If / is a bounded linear functional for elements of Hs(0) with
S > 0, / is in the dual space of Hs(0) denoted by H~s(0). The norm of / in
H~s(0) is
II/IIS0 3 sup Mill, (3.22)
uens(n) ps,n
where (, ) denotes the duality pairing between / and u in terms of the standard
L2{0) inner product.
Now we introduce the linear BrambleHilbert lemma.
Lemma 3.3 (BrambleHilbert) Let 0 be a domain in Ud, k > 0
an integer, 6 Â£ (0,1], and / a bounded linear functional on the space Hk+e(0)
68
that vanishes for all polynomials of degree < k on ft:
m = o, v
(3.23)
Then there exists a positive constant C = (7(0) such that
\f(u)\ < C\\f\\(k+e),a \u\k+e,n, VuÂ£ Hk+e(ft). (3.24)
3.5.1 Applications
Now that we have stated the linear BrambleHilbert lemma, we demon
strate its application for some functionals arising from standard hversion finite
element approximations. For simplicity, we work in two dimensions and assume
that the regularity of u is in the interval (l,p + 1],
polant of u on T. Then, rj = IIiU u is the linear interpolation error on T. By
the BrambleHilbert lemma, we have the following error estimate:
where hj is the diameter of T.
Proof: We map the computational element T to a reference element
T by an affine map, x = Mx + c, and we transform our functional:
Let T G Th be an element (i.e., a triangular region in $ft2) in a FE
triangulation of 0, u a function defined on T, and Hiu denote the linear inter
Example 3.1 (Linear interpolation)
S'G (0,1], (3.25)
where
69
and p(x) = 77(x) by the affine map. Since linear interpolation is exact for linear
polynomials, the BrambleHilbert lemma says that
/(i) < C(T) ll/ll_(s+1)if is+1,f, s e (0.1].
Next, we note that 11/ _(5+1) ^ < C\T\1!2 and following Ciarlet [11] we trans
form back to T as in,
M5+1
et M)!/2
where M < C hj/p~ and p~ is the diameter of the circle inscribed in T, and
det M = T/T. Since T is a reference or fixed element, T is chosen so that
prp = 0(1). After transforming back to T and simplifying terms, we obtain the
stated result.
\S+1,T,
3.26
In more generality, we obtain hversion results that are equivalent to
optimalorder L2 results for approximation by polynomials of degree p > 1:
Example 3.2 (pthorder interpolation)
[(ILpuu)dx < Chp1 T1/2u5+it, Ae(0,p]. (3.27)
Jt
Proof: Since pthorder interpolation is exact, IIpu u = 0, when u is
polynomial of degree k < p on T, the BrambleHilbert lemma holds when the
regularity parameter S for u satisfies 1 < A+l < p + 1. We follow the affine
mapping procedure discussed in the previous example and apply Bramble
Hilbert on the reference element. After the inverse mapping procedure to the
original element, we gain S + 1 powers of hj in the final result (3.27).
70
Finally, the estimates for approximation of firstorder derivatives of u
by polynomials of degree p > 1 are reduced by a power of hj (cf. Ciarlet [11],
Theorems 15.12).
Example 3.3 (Firstorder derivatives)
/ dx(Upuu)dx < C hj\T\1/2 \u\s+i}T, S G (0, jo];
JT
(3.28)
Proof: We follow the proof of the previous example with only one
modification. When applying the affine mapping argument, note that we have
the scaling relationship x = C pf1 x. Then after transforming the firstorder
derivative by the chain rule, we find that
dxu = C pTl dgu.
(3.29)
After applying the shape regularity inequality (2.1), the estimate of the previ
ous example is reduced by a power hj.
The hversion estimate (3.28) is equivalent to optimalorder H1 results
for approximation by polynomials of degree p > 1. Since FE uses optimalorder
H1 results to obtain error estimates as in (2.22), use of the BrambleHilbert
lemma does not improve the results of standard FE analysis. However, FVE
uses different functionals than those stated here (i.e., functionals based on
volume elements and not on finite elements) and the BrambleHilbert lemma
can exploit special volumes symmetries to add a power of h to the (p = 1)
approximation result (3.28) for derivatives.
71
3.6 Approximation Theory
As in FE analysis, error estimates for approximation of elements from
fractionalorder Sobolev spaces by elements from piecewise polynomial FE trial
spaces are are at the heart of FVE analysis. By applying the BrambleHilbert
lemma to the bounded linear functional (3.19), we demonstrate approximation
theorems for the diffusion functional Our baseline criterion is that our
results be equivalent to i71(0) optimalorder estimates for polynomial interpo
lation (as stated in (2.22) of Section 2.6.2).
To justify the upper end of the regularity scale (S + 1 >3) associ
ated with highorder polynomial trial spaces (p > 2), we assume that curved
(isoparametricequivalent) elements are being used to approximate the curved
boundaries of smooth domains. However, we ensure the invertibility of the
isoparametric mappings to define these curved elements by requiring a bound
ary discretization sufficiently fine to bound or approximate the isoparametric
map with an affine map (see Ciarlet [11] or Johnson [28] for details). Thus,
it is sufficient to work with the affine mappings that generate affineequivalent
families of polygonal elements. Hence, all of our approximation theory results
in this dissertation are demonstrated for polygonal finite elements under affine
mappings without loss of generality.
3.6.1 Optimalorder Approximation
First, we demonstrate an (i71(0)equivalent) optimalorder result for
pthorder interpolation when the triangulation Th and volumization Vh/p are
regular (defined in Sections 2.1.1 and 2.2.3)cf. Theorem 17.1 of Ciarlet [11],
72
According to the definitions of Section 2.3, the term (h/p)s in (3.30) is under
stood as C(p) hs for hversion FVE and as C(h)p~s for pversion FVE.
Theorem 3.1 Let r]p(u) = H^u u be the interpolation error for
u G W+ into Sp(Th). Then we have the following estimate:
V(r}p(u)) < C (h/p)s\\u\\s+1)m S G (1/2,p]. (3.30)
Proof: As in Cai [10], define a macroelement Eij as the region
formed by joining the endpoints of 77 and Ay, (see Figure 2.6). As noted in
Section 2.2.3, Eij is the union of two subelements (with Ay as a common
edge) in an auxiliary triangulation Ehlp which is assumed to be shape regular.
Thus, Eij = Tfj U Tfj, where each subelement is contained in a single element
of rh/p.
Recalling the nesting of the triangulations Eh/p C Eh^p in (2.6), let
Tp be the element of Eh^p that contains Tfj for k G {1, 2}. In our analysis, poly
nomial functions and their interpolants are defined on Tand then restricted
toTfy
Then we work on each subelement of Eij, as in Example 3 of Section
3.5. First, we split dij and into d\ and d2 which are formed by restricting dij to
each of the subelements: i.e., dj~ is dij restricted to 7^ = 77 ClT^k G {1,2};
dij(nP) = J2 dk(nP), (331)
ke{ 1,2}
dk(r]p) = f DVr]p(u) ntJ dS. (3.32)
Ay*
Second, we affine map from each subelement Tfj to a reference or work
element Tk and transform the bounded linear functionals di and d2: recall
73
that we lose a power of p when we transform firstorder derivatives of it;
\dk(u)\ < Cp,1 j?f}4(h), k G {1,2}, (3.33)
\lij I
where pj~ is the diameter of the circle inscribed in T^. Third, we bound dk(u)
and find that the BrambleHilbert lemma holds when the regularity parameter
S for u satisfies 3/2 < S' + 1 < p + 1: the lower and upper bounds are due
to the Sobolev trace theorem and to the order of the interpolation (i.e., pth
order interpolation is exact for polynomials of degree < p), respectively;
\dk(fjp)\ < C\u\s+1Tk, S G (1/2,p]. (3.34)
Fourth, we transform back to each subelement and find estimates for d\ and
d2 for S G (1/2,p]:
\dk(riP)\
S+l,Tk >
1 ij
(3.35)
where we have used the facts that hk = 0(h/p) is the diameter of TD T^,
\Tj\ > C 7^X4j for some C\ and the shape regularity of 'Thlp (i.e. hk/Pk <
Note that (3.35) is an i71(0)equivalent optimalorder local estimate on the
auxiliary triangulation Thdp. Since Thdp C d'hdp C Th (cf. (2.6) in Section
2.2.3)the triangulations are nested or refinements of one another, this esti
mate on Th/p holds on the FVE triangulation Thlp and the FE triangulation
Th, as demonstrated below.
Combining the results for d\ and d2 from (3.35), we find a local error
estimate for the bounded linear functional (3.19) for S G (1/2,p]:
\dtJ(riP(u))\ < IMVp) l>
ke{ 1,2}
74
1/2
< cWpf( y Mi)1/2 ( y III+,tj)
fce{i,2} MM fce{i,2}
< C(W5fM)/2( E II2S+,t)1/2. (:6)
VMM / fce{i,2}
where we have used the fact that 7b and 7b partition 77.
To obtain the final result (3.30), substitute the local estimate (3.36)
into the diffusion functional (3.18) and simplify: focusing on the Hs+1 semi
norm of it, we have
y~i y~i m5+1,xf ^2 \u\s+i,tp ^2 iuii+i,T5
{i,j}ew ke{i,2} '3 TpeTh/p TeTh
= \u\s+i,Th IMIs+l,m (3.37)
where weve used the nested property of the triangulations (2.6): each element
of Th/p is contained in a single element of 'Thlp which in turn is contained in a
single element of 'Th.
Following Remark 17.2 of Ciarlet [11], we offer an important corollary
to Theorem 3.1.
Corollary 3.1 If u Â£ C(0) fl Hs+1(Th) fl i71(0) for S Â£ (1/2,p],
i.e., u is only piecewise S + 1regular on the triangulation T^, we have the
following estimate:
V(r]p(u)) < C (h/p)s\u\s+ltTh, S Â£ (1/2,p], (3.38)
where
\u\s+i,Th = ( '22 Ms+i,t) (3.39)
TeTh
75
Proof: Follow the proof of Theorem 3.1 up to the penultimate line of
(3.37). We see that u Â£ C'(0) fl Hs+l(Th) fl Hl(Tl) is the necessary for the
proof to hold: (70(O) is the necessary so that Lagrange interpolation is well
defined; Hs+1(Th) fl H1{0) is necessary to ensure that 'D() and the piecewise
seminorm  \s+iTh are bounded.
The result (3.38) of the corollary avoids evaluating highorder deriva
tives of u across the edge Xij which may coincide with dT for an element T in
Th. That is, if it is only piecewise S + 1regular on 'Th, dau may have jump
discontinuities on dT for ct = 1 that generate Dirac delta functions on dT for
M > 2.
We work through an important illustrative application of Corollary
3.1 that will be used to derive a discrete ellipticity result in Section 3.7.
Lemma 3.4 Let u Â£ Sp(fTh) and p > 2. Then the interpolation error
into Sl(Th/p) satisfies the following estimate:
T>(u IIi^pu) < C (h/p)\u\2^h. (3.40)
Proof: For p > 2, u Â£ Sp(Th) C C(Q) D H2(Th) D Hx(0) and the
proof of Corollary 3.1 implies the following estimate for linear interpolation on
rh:
T>(u IlJ'it) < C h\u\2}ph.
Since 'ph!p is a pfold rehnement of Th, the same proof holds for linear inter
polation on Th/p with h replaced by h/p\ i.e.,
V(u IIi^Pu) < C (h/p)\u\2tYh/p,
76
where
\u\2^rh/p = ( y)
TpeTh/p
I 2 V/2
\U\2,TP)
The result (3.40) follows after noting that \u\\T = f2rpeT M^Xp 1=1
3.6.2 Superconvergent Approximation
For linear interpolation in two spatial dimensions, we have a choice
(due to Cai et al. [9, 10] and Siili [50]) between an i71(0)equivalent optimal
order result and a superconvergent result that depends on whether the volu
mization Vh/p is symmetric to the triangulation 'ph!p (see Section 2.2.2). That
is, from the perspective of the FVE diffusion functional, the approximating
power of a linear trial space on a symmetric volumization matches the power
of a quadratic trial space on a regular volumization. At the core of the super
convergence theory for the diffusive flux functional is the 1d fact that centered
differencing is exact for quadratics; we invoke certain volume symmetries, dis
cussed below, to extend this 1d result to multidimensional settings.
Theorem 3.2 (Cai and Siili) Let r]i(u) = II\u u be the linear
interpolation error for u Â£ W+ into Sq('T}i). Then we have the following
estimate:
D(?7i (it)) < C /t5u5__ijii, Se(l/2,P], (3.41)
where P = 1 (optimalorder) or 2 (superconvergent) in the absence or presence
of volume symmetry, respectively.
Proof: A detailed proof is given in Cai [10] that follows the steps
outlined in the proof of Theorem 3.1 for p = 1; instead, we choose to focus on
77
step three where we determine the extent (i.e., the range of regularities for u) to
which the BrambleHilbert lemma can be used to bound the diffusion functional
dij of (3.19). We work on the whole of the macroelement Eij by piecing
 k m ^
together partial results from each T7 in Eij to obtain the superconvergence
result (3.41) for symmetric volumes. That is, if u is quadratic on Tp U Tp
(recall that Tp is the element of Eh^p that contains the element Tj of Ehlp)
and then restricted to Eij, the sum of the linear interpolation errors from T7
2
and Tij (recall that interpolation is defined on Tp and the restricted to Tj in
our analysis) for the normal derivative of u will vanish if the volume is X and
7 symmetric (see Section 2.2.2 for the definition of these symmetries).
The essential ideas behind this result can be seen in the following
argument that holds when the diffusion coefficient is constant on 77. Dropping
the ^notation for convenience, we work on a reference macroelement Eij arising
from a symmetric volumization (see Figure 3.1): specifically, let Xij = [1,1] X
0 and 77 = 0 X [1,1] so that Xij and 77 are perpendicular bisectors of one
another; in this case, the normal derivative of u on 77 is simply dxu.
If u is quadratic, dxu is in the span of {1 ,x,y} on Eij and dxHjU is
in the span of {1} on each T7. Now we work on the common edge Xij of the
T7 subelements, Using the Wsymmetry of the volume, dxTliU is the standard
onedimensional central difference approximation to dxu at (0,0).
Since central differencing is exact for quadratics by a reflection sym
metry, we have the following result:
dxi]i(u)\(o,o) = 0. (3.42)
78
Figure 3.1: Formation of the macroelement Eij = Tl U Tf}.
79
In other words, evaluation of the derivative at the midpoint of an
interval suppresses the midpoint DOF (the midpoint quadratic basis function
has zero slope at the midpoint). Combining (3.42) with the continuity of the
interpolation on Ay,, we find that
^i(^)7ij G span{y}. (3.43)
Therefore, if the volume is also 7symmetric, the functional dij(r] 1) will vanish.
That is, the integrand of (3.19) will be an odd function on a symmetric 77
interval. Alternately, the midpoint or centroid quadrature rule applied to (3.19)
(which gives 0 as a result) will be exact in this case when the integrand is linear.
Therefore, the BrambleHilbert lemma holds when the regularity parameter S
is in (1/2,2] and the result (3.41) holds.
Notice that if bilinear interpolation on a nonuniform tensorproduct
mesh is employed (i.e., we use a bilinear trial space to represent the FVE
solution), as in Siili [50], only Xsymmetry is required for superconvergence.
That is, if u is quadratic, dxTlis in the span of {l,y}. Using the X
symmetry of the volume, we see that the onedimensional central differencing
or reflection symmetry result of (3.42) extends to any y Â£ 77 under bilinear
interpolation:
dxrji,i(u)\(o,y) = 0; (3.44)
i.e., for any fixed value of y the previous result (3.42) holds. Therefore, the
integrand of (3.19) vanishes on any 77 interval. The key to success for bilinears
is that they contain the quadratic xy term that is sufficient to approximate the
normal derivative exactly on 77 when Xsymmetry is present.
80
To demonstrate superconvergence for a variable diffusion coefficient,
we use a piecewise constant approximation D = HqD on 77. Then with the
modified notation di3(u) = di3(u; D)} we consider the splitting induced by D:
dlJ(r]1(u); D) = dtJ(rn(u)] D) + dtJ(rn(u)] r]0(D)). (3.45)
The first term on the righthand side of (3.45) is estimated by the linear
BrambleHilbert lemma as indicated above. However, the second term on the
righthand side of (3.45) requires the bilinear BrambleHilbert lemma of Ciarlet
[11], The bilinear lemma allows us to obtain an 0(h2) estimate by multiplying
the 0(h) estimates (determined by the linear lemma) for linear interpolation
of dxu and constant interpolation of D.
3.6.2.1 Extension to Three Dimensions
The proof of Theorem 3.2 is based on volume symmetries in two spa
tial dimensions (d = 2); however, we will see that the notion of volume sym
metry is dimensionally and geometry dependent. If d = 1, the modifications
in the proof are trivial: 77 is a point in 1d and 17^71 is replaced by 1 as stated
in Section 3.3, so the only volume symmetry is Xsymmetry: i.e., the volume
boundary 77 is at the midpoint of each element edge X7; in this simple case,
the diffusive flux is represented by a centered difference quotient.
In three spatial dimensions, the modifications in the proof are non
trivial: below we outline necessary steps to extend the result (3.41) to this
important case.
First, we upgrade our definition of volume symmetry: for Xsymmetry,
81
we require the planar surface 7 be an orthogonal bisector for X7; for centroid
7symmetry, we require X7 77 and for the intersection point P7 to be the
centroid or center of mass for 77; furthermore, we require that the geometry of
77 is such that the centroid quadrature rule is exact for linear functions: e.g.,
77 is rectangular, triangular, hexagonal, etc. (see Stroud [49] for details).
Second, we consider the use of a C linear trial space on tetrahedra. If
u is quadratic, dxu is in the span of {1, x, y, z} on U7 and dxTl \u is in the span
of {1} on each T7. Using the Xsymmetry of the volume, the onedimensional
central differencing or reflection symmetry (3.42) result becomes:
<9^1 (^)  (0,0,0) = 0. (3.46)
Combining (3.46) with the continuity of the interpolation on element faces that
have as an edge, we find that
dxrji{u)\ltJ G span{y,z}. (3.47)
Therefore, if the volume is also 7symmetric, the functional dij(r] 1) will vanish.
That is, the centroid quadrature rule (which gives 0 as result) will be exact.
Therefore, the BrambleHilbert lemma holds when the regularity parameter S
is in (1/2, 2] and the result (3.41) holds in three dimensions under assumptions
of X and 7 volume symmetries, as it did in two dimensions.
However, for C bilinears (i.e., Hi,iu Â£ span{l,x,y} {1, z}) on
triangular prisms formed in a tensorproduct mesh, volume symmetry requires
Xsymmetry and a 2d 7symmetry: i.e., the x and y crosssections of 77
must possess the centroid symmetry. That is, linearity in x and y requires
82
the 2d 77symmetry in x and y, but bilinearity in z eliminates the 2d 77
symmetry in z. The analog of (3.44) is
dxV i,i(^) (o,o^) = 0; (3.48)
the bilinearity in z eliminates 7 symmetries in the ^direction as in the 2d
bilinear case. However, the residue from the linearity in y remains:
dxr]i,i{u)\ltj G span{y}; (3.49)
in order to eliminate this term, we require that the y crosssection of 77 possess
the centroid symmetry as in the 2d linear case. The same remarks hold for
the analysis of the analysis of dzr]71 follows the 2d bilinear case exactly.
Finally, for C trilinears on rectangular bricks formed by a tensor
product mesh, volume symmetry requires Xsymmetry and a 1d 77symmetry:
i.e., either one of the 1d crosssections of 77 must possess the centroid sym
metry. That is, the quadratic components of a trilinear trial space obviate
77symmetry requirements as in 2d; however, its cubic component invokes a
reduced, 1d 77symmetry. A generalization of (3.44),
(o,j/,o),(o,o,^) = 0, (3.50)
eliminates the quadratic components of the trilinear interpolant. However, the
residue from the cubic component remains:
dxr]i,i,i{u)\lu G span{yz}. (3.51)
In order to eliminate this term, we require that a either a y or z crosssection
of 77 possess the centroid symmetry. So trilinear superconvergence requires a
reduced 7symmetry in addition to Xsymmetry.
83
Unlike standard FE superconvergence results that are usually local
and onedimensional, the FVE superconvergence result (3.41) is always global
and multidimensional. Notice (cf. (3.30)) that the superconvergence is consis
tent with the use of a quadratic, rather than a linear, trial space.
3.6.2.2 Degenerate Volume Paradigm
To help understand the robustness of the superconvergence result
(3.41) in two and three dimensions, we note the the linear FVE on a symmet
ric volumization is identical to a quadratic FVE (using a onepoint centroid
quadrature rule for the diffusion functional (3.19) at the intersection of Xij and
7ij) on a degenerate volumization posed on a uniform tensorproduct mesh.
That is, if we allow the volumes about the midpoint DOF in to degenerate
to sets of measure zero (see Figure 3.2), only vertex DOF in are expressed
in the resulting matrix systemall of the midpoint DOF are systematically
suppressed as can be seen from the onedimensional reflection symmetry re
sults (3.42) and (3.46). Also, the error in the onepoint centroid quadrature is
0(h2) and will not impair the superconvergence, as can be seen in the analysis
of the effects of numerical quadrature in [10].
In this process, the volumization V^2, which is consistent with the
use a quadratic trial space, degenerates to the symmetric volumization Vh,
which is consistent with the use of a linear trial space. As can be easily verified
by direct calculation using the (3.43) and (3.47), the two methodologies yield
identical matrix systems. Furthermore, if a bilinear trial space is used, the
degenerate volume paradigm holds on a nonuniform tensorproduct mesh
(without the use of a centroid quadrature rule) as can be seen from (3.44).
84
Figure 3.2: An example of V^2 degenerating to Vh.
85
Therefore, the symmetric volume method for linear FVE can be thought
of as a degenerate volume method for quadratic FVE. From this perspective,
the robust superconvergence behavior for linear FVE can be understood as a
special case of the ordinary optimalorder behavior for quadratic FVE.
Finally, it is important to emphasize the local character of the er
ror estimates of this chapter, since it may appear that superconvergence (for
linear trial spaces) is strictly relegated to uniform computational meshes and
that only linear convergence is to be expected on nonuniform meshes. Recall
that superconvergence results hold on the macroelement Eij formed by join
ing individual element edges Xij with individual volume interfaces 77. Even
when global symmetry conditions are violated, local (secondorder) supercon
vergence holds when local symmetry conditions for Xij and 77 are satisfied.
Qualitatively, global superconvergence can be seen as the cumulative effect of
local superconvergence behavior. Quantitatively, local superconvergence leads
to superlinear global convergence rates: the precise rate depends upon the
prevalence of local symmetry or local uniformity on the global computational
mesh.
Therefore, on general nonuniform meshes, partial superconver
gence can exist globally even when full secondorder superconvergence is lost
(see [9] for details). Even on highly nonuniform meshes, there are usually
significant regions of uniformity and symmetry. In these local regions, the
solution converges quadratically and remediates the global solution: super
linear, rather than linear, convergence is observed globally.
86
3.7 Ellipticity
In Section 2.6.4, use of the ellipticity result (2.23) in the FE error
analysis required that the trial space representation error be used as a test
function. In FE, test and trial spaces coincide so this was straightforward.
In FVE, test and trial spaces differ, and we must describe the test function
representation of a trial function.
Let u Â£ Wh with nodal representation u = J2iei uiwhere {8}8e/
is the DOF basis for Wh. Then u also has a corresponding lumped or test space
representation, Tt, in Xh:
u = J2u^Xt (3.52)
iei
In other words, if it Â£ Sp(Th), u Â£ S(Vh/p) is a volumewise constant inter
polant of it: u = Tl^pu. Our ellipticity result for FVE will use the lumped
representation defined by (3.52) to characterize test functions.
As in FE, a discrete ellipticity result is perhaps the most important
technical component of FVE analysis. Once ellipticity is established, FVE error
analysis for diffusion equation reduces to the approximation theory developed
in the previous section.
Using the optimalorder approximation result (3.40) of Lemma 3.4
from Section 3.6, we derive (for p > 2) an asymptotic ellipticity result for
the diffusion bilinear form in terms of our discrete FVE trial space Sp(l~h) and
our discrete H1 seminorm. That is, for (h/p) sufficiently small, the ellipticity
constant (3 is independent of the computational mesh. This type of argument is
unavoidable if we wish to maintain the FE formalism in our FVE error analysis
87
for higherorder (p > 2) polynomial trial spaces, since our discrete H1 semi
norm is defined as the restriction of the continuous H1 seminorm to piecewise
linear polynomials on the FVE triangulation 'Thlp.
Lemma 3.5 (Discrete Ellipticity) There exists a positive constant
/3 = /3(Dm; AÂ£>), for h/p < when p > 2, such that
F?(u, u) > /3 V u Â£ W^, (3.53)
where is the positive lower bound on the diffusion coefficient and AÂ£> =
lDmC11 wi,w/^2,Tfe is an asymptotic ellipticity parameter.
Proof: For u Â£ ^(T^), we follow Cai[10] for the following sub
result. By the definition of the Fbform (3.8), we have
B(u, u) = dij(U)i
iei jew,
= (uiujjdijtu),
{hjjew
since dij(u) = d^fu) (i.e., n4j = ny).
For any u Â£ (T^), we recognize that
n,; u,
I*
Â£ I J'Hj
Then, it follows that
B(u,u) = v ^ I DdS>
{i,j }ew I** I
p Dm 'y' (ui Uj)
2 Ihol
{i,j}ew
2
\Xii I
Dm \u\1}UJ.
(3.54)
Thus, the lemma is proved for the p = 1 case with (3 = Dm.
For u Â£ Sp(Th) with p > 2, define u\ = IIh^pu. Now, analyze the
following splitting induced by u\\
B(u}u) = B(ui, Hi) + B(u Ui, u),
> Dm \u\\^ C (h/p) \u\2tTh \u\1:LU,
> \Dm \u\l, (3.55)
where h/p < AÂ£> = \DmC~1 \u\i^/\u\2^h.
In the first three lines of (3.55), we use the fact that u = Tq, the
subresult (3.54), the fact that = iti)C, the bound (3.17), and the
approximation result (3.40), respectively. Thus, the lemma is proved for the
p > 2 case with (3 = Dm.
In addition, we see that B(u,u) = uTBu = Qb(u) V u Â£ Wh, where
u is the vector of nodal values for u Â£ B is the FVE matrix generated
by the Hform, and Qb(') is the associated quadratic form. Therefore, (3.53)
simply states that B is positive definite:
Qb(u) >0, V u, (3.56)
with equality holding in (3.56) only when u (i.e., u Â£ Wh) is identically zero.
Hence, (3.53) guarantees the existence and uniqueness of the numerical solution
uh Â£ Wh.
The FVE ellipticity result (3.53) clearly has the same form as the
corresponding FE result (2.23) (when restricted to Wh). In fact, the two
89
