Spectral Exterior Calculus and its Imple

Spectral Exterior Calculus and its Implementation - CaltechTHESIS
CaltechTHESIS
A Caltech Library Service
About
Browse
Deposit an Item
Instructions for Students
Spectral Exterior Calculus and its Implementation
Citation
Rufat, Dzhelil Sabahatin
(2017)
Spectral Exterior Calculus and its Implementation.
Dissertation (Ph.D.), California Institute of Technology.
doi:10.7907/Z9VX0DKV.
Abstract
Preserving geometric, topological and algebraic structures at play in partial differential equations has proven to be a fruitful guiding principle for computational methods in a variety of scientific fields. However, structure-preserving numerical methods have traditionally used spaces of piecewise polynomial basis functions with local support to interpolate differential forms. When solutions are known to be smooth, a spectral treatment is often preferred instead as it brings exponential convergence. While recent works have established spectral variants of discrete exterior calculus, no existing approach offers the full breadth of exterior calculus operators and a clear distinction between vectors and covectors. We present such a unified approach to spectral exterior calculus (SPEX) and provide detail on its implementation. Notably, our approach leverages Poincare duality through the use of a primal grid and its dual (with a natural handling of boundaries to facilitate the treatment of boundary conditions), and uses a twin representation of differential forms as both integrated and pointwise values. Through its reliance on the fast Fourier transform, the resulting framework enables computations in arbitrary dimensions that are both efficient and have excellent convergence properties.
Item Type:
Thesis (Dissertation (Ph.D.))
Subject Keywords:
spectral exterior calculus discrete differential geometry numerical methods
Degree Grantor:
California Institute of Technology
Division:
Engineering and Applied Science
Major Option:
Applied Physics
Thesis Availability:
Public (worldwide access)
Research Advisor(s):
Desbrun, Mathieu
Thesis Committee:
Desbrun, Mathieu (chair)
Colonius, Tim
Porter, Frank C.
Bellan, Paul Murray
Defense Date:
5 April 2017
Record Number:
CaltechTHESIS:05302017-094600781
Persistent URL:
DOI:
10.7907/Z9VX0DKV
Related URLs:
URL
URL Type
Description
DOI
Article adapted for Ch. 3
ORCID:
Author
ORCID
Rufat, Dzhelil Sabahatin
0000-0001-8766-2338
Default Usage Policy:
No commercial reproduction, distribution, display or performance rights in this work are provided.
ID Code:
10219
Collection:
CaltechTHESIS
Deposited By:
Dzhelil Rufat
Deposited On:
02 Jun 2017 20:10
Last Modified:
28 May 2025 22:38
Thesis Files
Preview
PDF
- Final Version
See Usage Policy.
39MB
Repository Staff Only:
item control page
CaltechTHESIS is powered by
EPrints 3.3
which is developed by the
School of Electronics and Computer Science
at the University of Southampton.
More information and software credits
Spectral Exterior Calculus
and Its Implementation

Thesis by

Dzhelil S. Rufat

In Partial Fulfillment of the Requirements for the
Degree of
Doctor of Philosophy

CALIFORNIA INSTITUTE OF TECHNOLOGY
Pasadena, California

2017
Defended April 5, 2017

Dzhelil S. Rufat
ORCID: 0000-0001-8766-2338

iii

To my mother Bakie.

ACKNOWLEDGEMENTS
First of all, I would like to thank my advisor Mathieu Desbrun for his mentorship
and guidance, without which this piece of work would not have been possible.
I would also like to acknolwedge the contribution of the late Jerrold Marsden
who first introduced me to the wonders of differential geometry and symmetry in
mechanics, first as a teacher and then briefly as an advisor.
I am grateful to Tim Colonius for his guidance in numerical fluid mechanics.
I am thankful to Yiying Tong, Joris Vankerschaver, Patrick Mullen, and Gemma
Mason for their suggestions and collaboration.
Additionally, the support and encouragement of my friends Daniel Sandoval, Vaclav
Cvicek, Nikola Kamburov, Silviu Pufu, Ana Caraiani, and Andreas Hoenselaar have
been invaluable to me during the years.
I also want to extend my deepest appreciation to staff Kevin Austin, Linda Schutz,
Alice Sogomonian, Divina Bautista, Sheila Shull, and Maria Lopez.
Last but not least, I want to say thank you to my mother Bakie, my father Sabahatin,
and my sister Refie for their love and unwavering dedication to me.

ABSTRACT
Preserving geometric, topological and algebraic structures at play in partial differential equations has proven to be a fruitful guiding principle for computational methods
in a variety of scientific fields. However, structure-preserving numerical methods
have traditionally used spaces of piecewise polynomial basis functions with local
support to interpolate differential forms. When solutions are known to be smooth,
a spectral treatment is often preferred instead as it brings exponential convergence.
While recent works have established spectral variants of discrete exterior calculus,
no existing approach offers the full breadth of exterior calculus operators and a clear
distinction between vectors and covectors. We present such a unified approach to
spectral exterior calculus (Spex) and provide detail on its implementation. Notably,
our approach leverages Poincaré duality through the use of a primal grid and its dual
(with a natural handling of boundaries to facilitate the treatment of boundary conditions), and uses a twin representation of differential forms as both integrated and
pointwise values. Through its reliance on the fast Fourier transform, the resulting
framework enables computations in arbitrary dimensions that are both efficient and
have excellent convergence properties.

PUBLISHED CONTENT AND CONTRIBUTIONS

[1] Dzhelil Rufat et al. “The chain collocation method: A spectrally accurate
calculus of forms”. In: Journal of Computational Physics 257 (Jan. 2014),
pp. 1352–1372. doi: 10.1016/j.jcp.2013.08.011.
Rufat participated in the conception of the project, derived the spectrally
accurate basis functions, formulated the corresponding discrete operators,
implemented them in code, and validated the results via extensive unit testing.

vii

CONTENTS

Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Published Content and Contributions . . . . . . . . . . . . . . . . . . . . . . . . vi
Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
Chapter I: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Organization of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 2
Chapter II: Differential Forms, Vector Fields, and Associated Operators . . . . 3
2.1 Vectors and Forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Exterior Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Wedge Product . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Hodge Star . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Contraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.6 Lie Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.7 Laplace-Beltrami . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.8 Flat and Sharp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
Chapter III: From Point Collocation to Chain Collocation . . . . . . . . . . . . 9
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Discrete Differential Forms and Operators . . . . . . . . . . . . . . . . 13
3.3 Basic Spectral Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.4 Spectrally Accurate Discrete Operators . . . . . . . . . . . . . . . . . . 24
3.5 Implementation in Fourier Space . . . . . . . . . . . . . . . . . . . . . 28
3.6 Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . 38
Chapter IV: Spectral Exterior Calculus . . . . . . . . . . . . . . . . . . . . . . . 41
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2 Context and Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Spatial Discretization and Associated Basis Functions . . . . . . . . . 53
4.4 Discretization of Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.5 Operators on Discrete Forms and Vector Fields . . . . . . . . . . . . . 69
4.A Transformations under Coordinate Change . . . . . . . . . . . . . . . . 74
4.B Chebyshev Identities and Limits . . . . . . . . . . . . . . . . . . . . . . 76
Chapter V: Implementation, Validation, and Results . . . . . . . . . . . . . . . 78
5.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.2 Convergence and Validation . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.A Auxilary Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.B Boundary Conditions for Laplace’s Equation . . . . . . . . . . . . . . 106

viii
Chapter VI: Programming Contributions . . . . . . . . . . . . . . . . . . . . . . 108
6.1 Spexy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
6.2 LicPy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.3 PyBindCpp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.4 Third Party Modules . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
Chapter VII: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
7.1 Review of Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 124
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Appendix A: Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
A.1 Invariance of the Discrete Wedge Product . . . . . . . . . . . . . . . . 127
A.2 Involutivity of the Discrete Hodge Star Operator . . . . . . . . . . . . 128
A.3 Types of Hodge Star Matrices . . . . . . . . . . . . . . . . . . . . . . . 130
A.4 In Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
A.5 Variational Derivation of the Equation of Motion of an Ideal Fluid . . 137

LIST OF FIGURES

Number
Page
2.1 A manifold with charts. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Tangent space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.1 Logically rectangular grids . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.2 Illustration of a regular 2D grid . . . . . . . . . . . . . . . . . . . . . . 14
3.3 Examples of periodic interpolators and histopolators . . . . . . . . . . 19
3.4 The map ϕ mapping the canonical interval to the bottom unit hemicircle 21
3.5 Chebyshev primal basis functions for a grid with N = 7 . . . . . . . . . 23
3.6 Chebyshev dual basis functions for a grid with N = 7 . . . . . . . . . . 24
3.7 Convergence graphs for a 1D Poisson equation . . . . . . . . . . . . . 37
3.8 Convergence graphs for a 2D Poisson equation . . . . . . . . . . . . . 39
3.9 Laplace’s equation for 1-form . . . . . . . . . . . . . . . . . . . . . . . 40
4.1 Computing an orthonormal frame for a surface embedded in 3D . . . 52
4.2 Computing an orthonormal frame for a surface embedded in 2D . . . 53
4.3 Schematic representation of the computational grid in 1D . . . . . . . 55
4.4 Cell complex with clipped edges . . . . . . . . . . . . . . . . . . . . . . 55
4.5 Computational grids of domain D . . . . . . . . . . . . . . . . . . . . . 60
4.6 Enumeration of the cells of a 4 × 4 grid . . . . . . . . . . . . . . . . . . 62
4.7 Enumeration of the cells of a 3 × 3 grid with clipped edges . . . . . . 63
4.8 Two discrete representations of forms. . . . . . . . . . . . . . . . . . . 64
4.9 Cell complex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.10 Schematic description of the spaces of forms and their various operators. 69
4.11 2 × 2 and 3 × 3 grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.1 Convergence for the Laplace operator in a square domain . . . . . . . 83
5.2 Convergence of the Laplace operator in a non-square domain. . . . . 84
5.3 Laplace’s equation for a 0-form f . . . . . . . . . . . . . . . . . . . . . 86
5.4 Laplace’s equation for a 1-form f . . . . . . . . . . . . . . . . . . . . . 87
5.5 Solution to Laplace’s equation for 1-forms on curved domains . . . . 88
5.6 Laplace’s equation for a uniform field . . . . . . . . . . . . . . . . . . . 89
5.7 Laplace’s equation for a rotational field . . . . . . . . . . . . . . . . . . 90
5.8 Velocity reconstruction from vorticity . . . . . . . . . . . . . . . . . . . 92
5.9 Diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

5.10
5.11
5.12
5.13
5.14
5.15
6.1
6.2

One-dimensional wave propagation . . . . . . . . . . . . . . . . . . . . 95
Wave propagation with no initial velocity . . . . . . . . . . . . . . . . . 96
Wave propagation with initial velocity along the x-axis . . . . . . . . . 97
Wave propagation with initial radial velocity . . . . . . . . . . . . . . . 98
Wave propagation over a hole . . . . . . . . . . . . . . . . . . . . . . . 100
Lie advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Streamline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Comparison of quiver plots with line integral convolution plots . . . . 118

LIST OF TABLES
Number
Page
3.1 Meaning of the basic notations. . . . . . . . . . . . . . . . . . . . . . . 13
4.1 List of basic and a few compound operators of exterior calculus, with
their dependence on the metric being indicated in the last column. . . 43
4.2 Meaning of the basic notations used throughout this document. . . . . 46
4.3 List of transformation for vectors and forms under a basis change in
2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.4 Number of primal and dual elements in 1D grids as a function of N. . 55
4.5 Correspondence between notation for periodic and bounded domains. 56
5.1 List of auxilary operators . . . . . . . . . . . . . . . . . . . . . . . . . . 80

Chapter 1

INTRODUCTION
The laws of physics are coordinate-independent as a consequence of general covariance. Many discretization schemes, however, rely on an arbitrary choice of a
coordinate system to describe and simulate these laws. Consequently, such discretizations often suffer from spurious numerical modes and other numerical shortcomings that are artifacts of the particular computational method and do not exist
in the underlying physical system.
Using the language of exterior calculus and differential geometry one can express
the laws of physics in a coordinate-independent manner. Whether studying fluid
mechanics, electromagnetism, wave propagation, diffusion, or advection, one can
readily cast the governing laws and equations using the objects and operators of exterior calculus and differential geometry. Then, rather than discretizing the coordinate
dependent equations of motion, one can instead discretize directly the intrinsic operators and objects of exterior calculus, which is precisely the approach taken by
discrete exterior calculus (Dec).
In addition to coordinate independence, another important goal of Dec is structure
preservation. Dec strives to preserve in the discrete setting as many as possible
of the properties of the smooth operators from the continuous realm. This allows
for important features of the continuous theory to be true in the discrete world by
construction rather than as an ad hoc feature. However, many of the previous works
on Dec are incomplete in terms of objects or operations that they consider, or they
are low order approximations.
In this thesis, we introduce spectral exterior calculus (Spex), which provides a
complete hierarchy of operators that have a high degree of accuracy and that can
be used to simulate a wide variety of physical problems. Spex is a merger between
well known spectral numerical methods and discrete exterior calculus in order
to create a formulation that has excellent convergence properties and is fast and
efficient in implementation, with operators that are discretized in a coordinate-free
and structure-preserving manner.

1.1

Organization of this Thesis
• Chapter 2 reviews manifold theory, differential forms, vector fields, the operators associated with them, and their properties. All these objects and operators
will be subject to discretization in the chapters that follow.
• Chapter 3 describes our first attempt at defining a spectrally accurate calculus
of forms. We present some discrete spectral operators such as the Hodge star
and the wedge product. The framework presented is limited to flat domains
defined over either periodic or bounded intervals.
• Chapter 4 presents spectral exterior calculus (Spex) with a complete hierarchy
of operators defined on arbitrary curved domains that are logically rectangular.
A twin representation of discrete forms is introduced, both as cochains and
components, along with the ability to seamlessly map between them. Metric
dependent operators are performed in the component representation by first
transforming into an orthonormal frame, in which the operators have a very
simple form.
• Chapter 5 provides detail of the Spex implementation in terms of the fast
Fourier transform and elementary array operations, describes how said implementation is validated, demonstrates its convergence, and presents results
and applications to a variety of problems in physics.

Chapter 2

DIFFERENTIAL FORMS, VECTOR FIELDS, AND
ASSOCIATED OPERATORS
In this chapter, we review basic manifold theory in order to introduce the objects
that live on manifolds and their associated operators, that in later chapters will be
subject to discretization. In particular, we will list the properties of these operators,
all of which will naturally carry over to our spectrally accurate discretization.
Roughly speaking manifolds are topological spaces that locally resemble Euclidean spaces. Given a set M, a chart
of M is a bijective map ϕ ∶ U → ϕ(U) ⊂
Rn that takes a subset U of M and maps
it to (x 1, . . . , x d ). We call the scalar
values x i the coordinates of the point
m ∈ U ⊂ M. A manifold is then a set
M for which every point is covered by at
least one chart, and charts are compatible with each other to ensure smoothness [1].
2.1

Figure 2.1: A manifold with charts.

Vectors and Forms

To each point m ∈ M we can associate a tangent space denoted by Tm M. A tangent
vector v to a manifold M at point m is an element of Tm M. The values (v 1, . . . , v d )
are called the components of v with respect to some local coordinate chart. The
collection of all tangent spaces forms a tangent bundle denoted by TM:
TM = ∪ Tm M.
m∈M

To each tangent space Tm M we can associate a cotangent space Tm∗ M whose
elements are the linear functions α ∶ Tm M → R that are known as covectors or
1-forms. The values (α1, . . . , αd ) are known as the components of the 1-form α.

↵ ^ ? = (↵ · ) µ

• symmetric

µ = ?1

Basic Manifold↵ ·Theory
= ·↵

Basic Manifold
1 BasicTheory
Manifold Theory
• bilinear

• Manifold M

(a↵ + b ) ·
• Manifold M

= a(↵ · ) + b( · )

• Manifold M
· (b: U
+c!
) = '(U
b(↵ · )
) +⇢
c(↵R· d) U ⇢ M 7! (x1
• Coordinate Charts↵'

A vector field X on the manifold M is a map X ∶
• Coordinate Charts ' : U ! '(U ) ⇢ Rd U ⇢ M 7!
1 · 1 = 1, dx · dx = 1, dx · dy = 0, dy · dy =d1, (dx
1 1
1 ^ dy) · (dx
• Charts
Tangent Space
v2
. . ,1M
v d )^ dy)d=(x
mM
Coordinate
' : UT!
'(U
) T⇢mvM
R2 Tm(vM
U, .⇢
M → TM that assigns a vector•X(m)
to each point
• Tangent
Space
Tm M
(v , . . . 7!
,v )
↵^
⇤ ? = (↵ · ) µ
⇤µ = ?1
m of the manifold. The space of all vector fields
• Cotangent Space Tm M ↵⇤ 2 Tm1M ⇤(↵1 ,d. . . , ↵d )
Cotangent
Tm M(v↵ ,2. T. m
• Tangent Space Tm•M
v 2 Space
Tm M
.M
, v )(↵1 , . . . , ↵d )
is denoted by X(M). Similarly, a 1-form field•isTangent
a1 Basic
Manifold
Theory
M = TM
[ =Tm M
• Bundle
Tangent T
Bundle
[ Tm M
m2M
⇤ M
map α ∶ M → T ∗ M that assigns
covector
α(m)
Manifold
• Cotangent Space Tm M ↵ 2⇤ Tm M m2M
(↵
1 , . . . , ↵d )
⇤1
• Cotangent
Bundle
M'(U=T) ⇢⇤ M
[d =
T⇢m
M7!T(xm
Cotangent
Bundle
Coordinate
Charts
, . . . , xd )
to each point of the manifold. The space of all 1m2M m2M
2 Tm
• Tangent Bundle• Tangent
T MSpace
= Tm M
[ vT
M(v1, . . . , vd)
mM
form fields is denoted by Λ1 (M). Scalar fields• are
•Field
Vector
Field
!T
M X(m)
2 Tm M
m2M
Vector
X(m)
2 Tm M
• Cotangent Space Tm
M ↵ 2 Tm
M (↵1 , . . . , ↵d )
0-forms, and are denoted by Λ (M), whereas k- • Tangent
k ⇤Bundle !
k-Form
· ·M
⇥ T!
[: T⇥
mM
mT
k-Form•!T
(m)
T=mm2M
·⇤·M
· ⇥⇥T·m
R} ! R
• Cotangent•Bundle
M: T=M(m)
|T
|m2M
{zM {z
Figure 2.2: Tangent
space.
forms, which are antisymmetric tensors acting on k
• Cotangent Bundle T M = [ Tmk
m2M

vectors, are denoted by Λ k (M), i.e., if ω k ∈ Λ k (M) • Vector
• Space
ofMVector
Fields
X(M)
Field
! T MX(m)
X(m) 2 T2
m MT
• Space
of Vector
Fields
X(M)
• Vector Field
X:M
! XT: M
mM
then
• k-Form
! (m) of
: k-Forms
Tm M ⇥ · · · ⇥
Tm(M)
!R
• Space
{z
• Space of k-Forms |⇤k (M)
(m)
Tm M ! R
ω k (m) •
∶ k-Form
Tm M × ⋅ ⋅!
⋅ ×k T
M:→ R.Tm M ⇥ · · · ⇥
X(M)
´¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¶ • Space| of Vector Fields{z

• Space of k-Forms ⇤k (M)

Having defined the spaces of vector fields X(M) and form fields Λ k (M), we now
• Space of Vector Fields X(M)
turn our attention to the operators that act upon them.
2.2

• Space of k-Forms ⇤k (M)

Exterior Derivative

The continuous exterior derivative,
d∶

Λ (M) → Λ

k+1

(M),

extends the notion of derivative to differential forms, turning a k-form ω k to a
(k +1)-form dω k . This metric-independent operator satisfies:
• linearity
d(c1 ω1 +c2 ω2 ) = c1 dω1 + c2 dω2

• nilpotency
d2 = 0
• locality
d(ω∣U ) = (dω)∣U .
The exterior derivative encompasses the usual linear operators in vector calculus
(gradient, curl, and divergence); its nilpotency reflects the classical identities ∇ ⋅
(∇× v⃗) = 0 and ∇× (∇ f ) = 0. Moreover, a number of theorems (the fundamental
theorem of calculus, as well as Green’s, Kevin’s, and Gauss’s theorems) can be seen

as special cases of the general Stokes’ theorem which states that d and the boundary
operator ∂ are dual operators, in the sense that:
dω =
ω.

∂σ

Fortunately, one can use the notion of chains and cochains from algebraic topology to
provide discrete realization of this derivative operator [73, 49, 8, 21] that captures all
its defining properties. The exterior derivative is invariant under diffeomorphisms
ϕ∗ dα = dϕ∗ α.
Given its purely topological (metric-independent) definition, no special care is necessary to adapt the traditional discrete exterior derivative to our spectral framework
(see Section 3.4).
2.3

Wedge Product

The continuous wedge product is also a metric-independent operator, assembling a
(k +l)-form from a k-form and an l-form; it generalizes the notion of cross product
to differential forms. The wedge operator
∧∶

Λ k (M) × Λl (M) → Λ k+l (M)

has the following properties:
• associativity
α ∧ (β ∧ γ) = (α ∧ β) ∧ γ
• bilinearity
(aα + bβ) ∧ γ = a(α ∧ β) + b(β ∧ γ)
α ∧ (bβ + cγ) = b(α ∧ β) + c(α ∧ γ)
• anti-commutativity
α k ∧ βl = (−1) kl βl ∧ α k .
The wedge product also satisfies a Leibniz rule (also known as a product rule) with
the exterior derivative
d(α k ∧ βl ) = dα k ∧ βl + (−1) k α k ∧ dβl .

(2.1)

Various discretizations of this operator have been proposed over the years, from
primal-dual versions [35], to metric-independent definitions based on cup product [22, 7]. While most properties of the wedge product are retained by these
discrete counterparts, associativity is only satisfied for closed forms, or in the limit
of mesh refinement [40]. Additionally, the wedge product is a purely topological
(metric-independent) operator that is invariant under a diffeomorphism:
ϕ∗ (α ∧ β) = (ϕ∗ α) ∧ (ϕ∗ β).
This property is preserved in the discrete setting (see A.1).
2.4

Hodge Star

The continuous Hodge star operator, denoted by ⋆, is a metric-dependent operator
that maps k-forms on a d-dimensional manifold to (d−k)-forms:
⋆∶

Λ k (M) → Λd−k (M).

The Hodge star is an (anti-)involution in that if it is applied twice, the result is the
identity up to a sign; more precisely, on a d-dimensional manifold:
⋆ k ⋆d−k = (−1) k(d−k) .

(2.2)

Under certain conditions, the discrete equivalent is also an involution (see A.2). In
dimension two, the Hodge star operator in coordinates satisfies:
⋆1 = dx ∧ dy,

⋆dx = −dy,

⋆dy = dx,

⋆(dx ∧ dy) = 1.

Along with the wedge product, the Hodge star defines a natural inner product on
k-forms:
α k ⋅ β k = ⋆(α k ∧ ⋆β k ).
(2.3)
2.5

Contraction

The contraction (sometimes referred to as the interior product) is a metric independent operator
⌟ ∶ X(M) × Λ k (M) → Λ k−1 (M)
that arises from the natural pairing of vectors and forms. It turns k-form ω k into a
(k − 1)-form X ⌟ ω k . Since a k-forms are linear maps that take k vectors as inputs,
the contraction simply places the vector field X into the first slot of the k-form to
reduce its degree by one:
(X ⌟ α k )(v2, . . . , vk ) = α k (X, v2, . . . , vk ).

The contraction is left-distributive over the wedge product
X ⌟ (α ∧ β) = (X ⌟ α) ∧ β + (−1) k α ∧ (X ⌟ β).
2.6

(2.4)

Lie Derivative

The notion of a Lie derivative L X extends the usual concept of the derivative of a
function along a vector field X to k-forms:
L∶

X(M) × Λ k (M) → Λ k (M).

Although a dynamic definition of this operator can be given in terms of extrusions
along the flow of a vector field [48], in this work we will rely on an algebraic
definition given via Cartan’s Magic formula:
L X α = dX ⌟ α + X ⌟ dα.
The Lie derivative can be directly defined as a composition of the metric independent
operators contraction ⌟ and exterior derivative d, without the need for its own
separate discrete definition. An important property, which will also hold in the
discrete setting, is that the Lie derivative commutes with the exterior derivative
L X d = dL X , which can be easily shown to hold true because of the fact that d2 = 0:
dL X α = d(dX ⌟ α + X ⌟ dα)
= d(X ⌟ dα)
= (dX ⌟ +X ⌟ d)dα
= L X dα.
A useful consequence of this is that discrete Lie advection of closed forms will
remain closed by construction. As a consequence of Eq. (2.1) and Eq. (2.4), the Lie
derivative, like the exterior derivative, also satisfies a Leibniz rule over the wedge
product:
L X (α ∧ β) = (L X α) ∧ β + α ∧ (L X β).
2.7

Laplace-Beltrami

The Laplace-Beltrami operator (sometimes referred to only as the Laplacian) is an
elliptic operator that is ubiquitous in physics and mathematics. For example, it
arises in the diffusion and wave equations:
∆∶

Λ k (M) → Λ k (M).

The Laplace-Beltrami operator acts on differential forms to produce forms of the
same degree. Explicitly, it is a composition of Hodge stars and exterior derivatives:
∆ = ⋆d ⋆ d + d ⋆ d⋆
and its discrete version can be expressed as compositions of the discrete counterparts
of the above operators. The Laplace-Beltrami operator is metric dependent because
of the presence of the Hodge star in its definition.
2.8

Flat and Sharp

The flat and sharp are metric dependent operators that convert vectors to 1-forms
and vice versa:
♭ ∶ X(M) → Λ1 (M),
♯∶

Λ1 (M) → X(M).

Known also as the musical operators, they are inverses of each other:
♭ ♯ = ♯ ♭ = id
In indicial notation the flat raises the indices of a vector to convert it to a form,
whereas the sharp lowers the indices of a form to create a vector. The musical
operators along with the contraction define an inner product (or metric) on both
vectors (X, Y ) and forms (α, β) that is given by
X ⋅ Y = X ⌟ Y ♭,

α ⋅ β = α♯ ⌟ β.

To illustrate that the operators considered in this chapter form a complete set, we
will cast some of the familiar operators of vector calculus into the ones that we have
defined so far. Operators such as the gradient, curl, and divergence can be readily
expressed in 3D as
∇ f = (d f )♯

(gradient),

∇ × v = [⋆(dv♭ )]♯

(curl),

∇ ⋅ v = ⋆d(⋆v♭ )

(divergence).

Additionally, the vector cross product in 3D can be written as
u × v = [⋆(u♭ ∧ v♭ )]♯
. Keep in mind that associativity of the wedge product does not imply associativity
of the cross product.

Chapter 3

FROM POINT COLLOCATION TO CHAIN COLLOCATION

3.1

Introduction

Recent years have seen the development of novel discretizations for a wide variety
of systems of partial differential equations. In particular, preserving in the discrete
realm the underlying geometric, topological, and algebraic structures at stake in differential equations has proven to be a fruitful guiding principle for discretization [2,
65, 3, 58]. This geometric approach has led to numerical methods, analyzed in,
e.g., [36, 3], that inherit a variety of properties from the continuous world. However,
geometric discretizations of elasticity, electromagnetism, or fluid mechanics have
mostly been demonstrated using spaces of piecewise polynomial differential forms.
Many problems where solutions are smoothly varying in space call for a spectral
numerical treatment instead, as it produces low-error, exponentially converging
approximations by leveraging fast implementations of transforms such as the fast
Fourier transform. In an effort to provide structure-preserving numerical tools with
spectral accuracy on logically rectangular grids over periodic or bounded domains,
we present a spectral extension of the discrete exterior calculus described in [8,
21, 7, 31]—and point out that the resulting computational tools extend well-known
spectral collocation methods.
Review of Previous Work
Computational methods preserving geometric structures have become increasingly
popular over the past few yeears, gaining acceptance among both engineers and
mathematicians [4]. Computational electromagnetism [8, 65], mimetic (or natural)
discretizations [52, 7], finite-dimensional exterior calculus (including discrete exterior calculus (Dec, [20, 21]), and finite element exterior calculus (FEEC, [2, 3]))
have all proposed discretizations that preserve vector calculus identities in order to
improve numerics. In particular, the relevance of exterior calculus (Cartan’s calculus of differential forms [18]) and algebraic topology [49] to computations came to
light.
Exterior calculus is a concise mathematical formalism to express differential and
integral equations on smooth and curved spaces, while revealing the geometric

10
structures at play and clarifying the nature of the physical quantities involved. At the
heart of exterior calculus is the notion of differential forms, denoting antisymmetric
tensors of arbitrary order. As integration of differential forms is an abstraction of
the measurement process, this calculus of forms provides an intrinsic, coordinatefree approach particularly relevant to the neat description of a multitude of physical
models making heavy use of line, surface and volume integrals [1, 15, 25, 26, 60].
Moreover, physical measurements, such as fluxes, represent local integrations over a
small surface of the measuring instrument. Pointwise evaluation of such quantities
does not have physical meaning; instead, one should manipulate these quantities
only as geometrically-meaningful entities integrated over appropriate submanifolds.
Algebraic topology, specifically the notion of chains and cochains [73, 49] has been
used to provide a natural discretization of differential forms and to emulate exterior
calculus on finite grids: a set of values on vertices, edges, faces, and cells are proper
discrete versions of respectively pointwise functions, line integrals, surface integrals, and volume integrals [8]. This point of view is entirely compatible with the
treatment of volume integrals in finite volume methods, or scalar functions in finite
element methods; however, it also involves the “edge elements” and “facet elements”
(as first introduced in computational electromagnetism) as special Hdiv and Hcurl basis elements [51]. Equipped with such discrete forms of arbitrary degree, Stokes’
theorem connecting differentiation and integration is automatically enforced if one
thinks of differentiation as the dual of the boundary operator—a particularly simple
operator on meshes. With these basic building blocks, important structures and
invariants of the continuous setting directly carry over to the discrete world, culminating in a discrete Hodge theory [21, 3]. As a consequence, such a discrete exterior
calculus has already proven useful in many areas such as electromagnetism [8, 65],
fluid simulation [58], (re)meshing of surfaces [68, 47], and graph theory [31] to
mention a few. So far, only piecewise polynomial basis functions [2, 72] have been
employed in these applications, thus limiting their computational efficiency in terms
of convergence rates.
Spectral Methods
Spectral methods are a class of spatial discretizations of differential equations widely
recognized as crucial in fluid mechanics, electromagnetics, and other applications
where solutions are expected to be smooth. Central to the efficiency of this large
family of numerical methods is the fact that the approximation of a periodic C ∞
function by its trigonometric interpolation over evenly spaced points converges faster

11
than any polynomial order of the step size. This is sometimes referred to as “spectral
accuracy” or “super-convergence.” In practice, spectral accuracy can be achieved
for bounded domains through continuation methods [14] or using Gauss-Lobatto
quadrature on Legendre or Chebyshev grids [17]. A larger number of spectral methods have been designed, varying in the mesh they consider (primal grids only, or
staggered grids [39]), and the locations at which they enforce partial differential
equations (PDE). Be it for Galerkin, Petrov-Galerkin, or collocation-based spectral
schemes, it has however been noticed that besides constructing spectrally accurate
approximations of the relevant fields and their derivatives involved in a PDE, numerically preserving conservation properties helps in obtaining stable and/or physically
adequate results [17]. Yet, numerical schemes are often proven conservative a posteriori, as a formal approach to guarantee conservation properties by design remains
elusive.
Motivations and Contributions
Despite an increasingly large body of work on numerical approaches based on
exterior calculus, developing a spectrally accurate calculus of discrete forms has
received very little attention—with a few recent exceptions [59, 10, 28, 55] that we
will build upon. We present a discrete exterior calculus of differential forms on
periodic or bounded domains, including wedge product, Hodge star, and exterior
derivative, all of which converge spectrally under grid refinement while utilizing
fast Fourier methods to remain computationally efficient. In order to construct
a spectral representation of the operators on differential forms, we expand the
conventional tools of spectral methods to give spectrally accurate approximations
of fields for which integral values over specified domains are known—a process
referred to as histopolation [59, 72]. In this chapter we construct a histopolation
using trigonometric polynomials on periodic domains, and consider the extension
to bounded domains using a Chebyshev grid, thereby allowing the use of the fast
Fourier transform for efficient calculation.
Our work lays out a set of spectral, structure-preserving computational tools with
the following distinguishing features:
• We leverage existing work in algebraic topology to discretize space through
chains (linear combination of mesh elements) and differential forms through
cochains (discrete forms). The resulting discrete de Rham complex, that
by construction satisfies Stokes’ theorem, offers a consistent, “structure-

12
preserving” manipulation of integrals and differentials which respect important conservation laws. This approach, used mainly so far in non-spectral
computations [2], was identified in [10, 28] as a significant departure in the
construction of conservative schemes from traditional spectral methods, since
divergences, gradients, and curls are no longer computed through derivation
but directly evaluated via a metric-independent exterior derivative without
having recourse to approximations—thus exactly enforcing the divergence
theorem, Green’s theorem, etc.
• We extend the basic Dec operators to ensure spectral accuracy. Using a
framework with some similarities to that of Bochev [7] and Dec [21], we
construct operations on discrete differential forms by interpolating their finitedimensional (point-sampled or integrated) representation, performing exterior
calculus operations in continuous space, and point-sampling or integrating the
results to project back to their finite-dimensional representation. We differ
from previous spectral Dec approaches [55, 10, 28] in that we provide a
natural Hodge star leveraging mesh duality, introduce a wedge product, satisfy
consistency conditions between the various operators and basis functions, and
provide their explicit expressions to facilitate implementation.
• Our spectrally accurate operators result in a spectral chain collocation method
to solving partial differential equations: a differential equation is not enforced
at a series of points, but at a series of mesh elements (points, edges, faces,
etc), extending point collocation methods to respect the natural geometric
structure of the continuous equation. We also provide a discussion on the
implementation of this spectral computational framework, and show results
on Poisson problems with various boundary conditions.
• Finally, our approach is versatile: since it builds on the notion of integrated
values (cochains to be precise), one can use the same computational tools on
curved grids with only minor changes.
For simplicity of presentation, we will first treat regular and Chebyshev onedimensional (1D) grids before extending our approach to grids of arbitrary dimension via tensor products. Other Chebyshev-like grids are easy to derive as well,
since integrals of differential forms (the building blocks of our approach) are unchanged by pushforward or pullback. For the reader’s convenience, Table 3.1 lists
the main notations we will use throughout this chapter.

13
Symbol

Meaning

Dimension of the domain

Spatial domain in Rd

Computational grid of Ω

σk

Primal k-dim element of grid K
Dual k-dim element of grid K

Dual grid to K

Duality operator (∗σ k = ̃
σ n−k )

Boundary operator on mesh elements

Λk

Space of differential k-forms ωk over Ω

Space of discrete k-forms ω̄k on K

Discrete exterior derivative operator d

Discrete wedge operator ∧

Discrete Hodge star operator ⋆

Discrete Fourier transform

.̂

Quantities expressed in Fourier space

Table 3.1: Meaning of the basic notations.
3.2

Discrete Differential Forms and Operators

Differential forms were pioneered by Cartan [18] in an effort to provide a unified
approach to defining differentiation and integration over curves, surfaces, volumes,
and higher-dimensional manifolds. Since then, forms have been shown to be crucial
in elucidating the structures and invariants in physics [26]. Since most of our
measurements of the world are of integral nature, forms have been also at the
core of a number of numerical approaches targeting the definition of numerical
counterparts of differential operators that are compatible with the geometric and
topological structures underlying well-posed partial differential equations.
Most computational methods based on exterior calculus use the notion of simplicial
k-cochain as a fundamental object that assigns a number to each simplex of dimension k of the simplicial complex. Cochains are, in a sense, a natural discretization of
differential forms of degree k, which are objects to be integrated over k-dimensional
domains to return a number; in fact, de Rham’s theorem states that the (discretization) map from the de Rham complex to the simplicial cochain complex induces
an isomorphism on cohomology. Conversely, one can go from cochains back to
differential forms using Whitney forms [73, 8] or higher-order piecewise polynomial finite element spaces [2, 72]. These mimetic properties of finite-dimensional

14
approximations to differential forms have led to the term discrete differential forms
to describe cochains [21].
In this section, we review the principles and approaches used in discrete versions
of exterior calculus in preparation for extending them to spectral computations. By
construction, the calculus of discrete differential forms automatically preserves a
number of important geometric structures, including integration by parts (with a
proper treatment of boundaries), Stokes’ theorem, the de Rham complex, Poincaré
duality, Poincaré’s lemma, and Hodge theory. Therefore, it provides a suitable foundation for a coordinate-free discretization of geometric field theories and associated
computations, with applications including electromagnetism [65], incompressible
Euler and complex fluids [58], and meshing algorithms [68, 19]. The particular
“flavor” of discrete differential forms and operators we will be leveraging is known
as discrete exterior calculus, or Dec for short; see [35, 43]. (For related work in
this direction, see also [33] and [2]; Dec most significantly differs from these other
discrete calculus approaches in the way the Hodge duality is defined through the use
of a dual mesh.) We refer the reader to [21] for a brief introduction to Dec.

(a) Periodic

(b) Chebyshev

Figure 3.1: Logically rectangular grids such as (a) regular grids over periodic domains,
(b) Chebyshev grids over bounded domains, or (c) an arbitrarily mapped rectangular grid.
Primal elements are displayed in black, while their duals are in red.

(a) Primal/dual grid

φ̄1y
mn

φ̄2mn

φ̄0mn

φ̄1x
mn
(b) Primal grid

φ̃1x
mn

φ̃0mn

φ̃2mn

φ̃1y
mn

Figure 3.2: Illustration of a regular 2D grid (in black) along with its dual (in red), along
with our mesh orientation and index convention.

15
Mesh and Dual Mesh as Spatial Discretization
In our approach, the continuous domain (usually a smooth d-dimensional manifold
Ω) is approximated by a mesh—more precisely, by a cell complex that is manifold,
admits a metric, and is orientable. We will generally denote the complex by K, and a
cell in the complex by σ. Superscripts will be used at times to denote the dimension
of a cell; i.e., σ 0 represents a vertex of K, σ 1 represents an edge, etc. Moreover, we
restrict our study to the case of logically rectangular grids over periodic and bounded
domains that include regular grids and Chebyshev grids in Rn (see Figure 3.1); that
is, each k-cell for k ≥ 1 has 2k non-degenerate faces, and the mapping function
between the rectangular grid of a rectangular domain and the actual spatial grid is
assumed to be bijective and C ∞ . Consequently, the mesh topology is easily encoded
using indices as indicated in Figure 3.2; we will also assume prescribed orientations
as the figure indicates as well.
From a mesh K, one can construct a dual mesh—where duality is to be understood
in the sense of graph duality. The dual mesh is defined via the location of its dual
vertices, whereas its connectivity is directly induced by the primal connectivity.
Once a primal mesh and its dual are defined, one can define the duality operator ∗
that maps each k-cell σ k of the primal mesh K to a dual (n−k)-cell ∗σ k in its dual cell
complex ∗K. Note that a refined definition of the dual mesh for bounded domains
(where dual cells at the boundary are restricted (“clipped”) to K to allow proper
enforcement of boundary conditions) is assumed as well (see, for instance, [22, 65]).
The canonical placement of primal and dual nodes for Chebyshev grids will be given
in Section 3.3; but arbitrarily distorted grids will be straightforward to handle (see
Section 3.6), owing to the fact that the exterior derivative and the pullback operator
commute.
Discretization of Differential Forms
The fundamental objects of Dec are discrete differential forms. A discrete k-form
ω̄ k , the discrete counterpart of a continuous k-form ω k , assigns a real number ω̄ik
to each oriented k-dimensional cell σik in the mesh K—and zero for all other mesh
elements. (The superscripts k are not actually required by the notation, but they are
often useful as reminders of what dimension of form or cell we are dealing with.)
The value of the form on σik is sometimes denoted by ⟨ω k , σik ⟩, and represents the
value of the continuous form ω k integrated over the mesh element σik , i.e.,
ωk .
ω̄i ≡ ⟨ω , σ ⟩ =
σk

16
For implementation purposes, a discrete k-form can therefore be thought of as an
array of values, one per k-dimensional cell in the mesh K. For example, discrete
0-forms sample values of an associated continuous 0-form at vertices, discrete 1forms sample values of a continuous 1-form on all edges, etc. Note that we can
even integrate these discrete forms over discrete paths by linearity: simply add the
form’s values on each cell in the path, taking care to flip the sign if the path is
oriented opposite the cell. Formally, these “paths” of k-dimensional elements are
called chains, and discrete differential forms are cochains, where ⟨⋅, ⋅⟩ is the pairing
between cochains and chains [49].
Discrete differential forms can be defined either on the mesh K or on its dual ∗K. We
will refer to these as primal forms and dual forms respectively, and denote the space
̃ Note that
of primal discrete forms by Λ̄ and the space of dual discrete forms by Λ.
there is a natural correspondence between primal k-forms ω̄ k and dual (n−k)-forms
̃n−k , corresponding to the mesh duality discussed in Section 3.2, for which each
primal k-cell σ k has an associated dual (n−k)-cell ∗σ k (Figure 3.2). This important
property will be used below to define the discrete Hodge star operator.
Reduction and Reconstruction Maps
The reduction and reconstruction maps provide a way to go back and forth between
continuous forms and their discrete realizations.
Reduction. The reduction map (also called the de Rham map) is a linear operator
P that projects a continuous form to its discrete realization on the grid through
integration over mesh elements:
P∶

Λ k → Λ̄ k
ω → ω̄

with

ω̄ik = (Pω k )i ≡

σik

ωk .

̃ the analogous operator mapping continuous forms to their dual
We will denote by P
discrete counterparts in a similar fashion:
̃∶

̃k
Λk → Λ
̃,
ω →ω

with

̃ k )i ≡
̃ik = (Pω

σik

ωk .

Note that this definition of reduction extends the notion of point sampling: while
the reduction of a 0-form is found by simply point-sampling its value at each vertex
of the grid, the reduction of a general k-form is its evaluation (i.e., integral) on all

17
the k-dimensional elements (vertices for k =0, edges for k =1, faces for k =2, etc) of
the grid.
Reconstruction. Conversely, the reconstruction map (R) is a map which reconstructs a continuous k-form from its discrete realization by interpolation for k = 0,
and by histopolation otherwise:
R∶

Λ̄ k → Λ k
ω̄ k → ω k .

̃ the analogous operator mapping dual discrete forms to continWe will denote by R
uous forms in a similar fashion. Note that we will sometimes omit the dual sign ̃.
for clarity, as which reconstruction operator is meant is unambiguously implied by
the (primal or dual) nature of the discrete form it is applied on.
If one has a set of basis functions {φ0k (x), φ1k (x), . . . } for k-forms that satisfy the
property
∀i, j,

σik

φ kj = δi j

(3.1)

then R can trivially be defined as
ω k = Rω̄ k ≡ ∑ ω̄ik φik .

One can readily verify that this reduction map is a left inverse of the reconstruction
map:
PR = Id.
However, the converse is not true; the reconstruction map is only approximately the
right inverse of the reduction map, with equality in the limit when the mesh element
size h approaches 0:
∥ω − RPω∥ ÐÐ→ 0,
h→0

with a rate of convergence determined by the chosen norm on forms and the degree
of the basis functions. While Whitney first introduced a one-sided inverse of the
de Rham map using what amounts to piecewise linear basis functions [73], we will
instead use global basis functions satisfying Eq. (3.1) to provide spectrally accurate
reconstructions.

18
3.3

Basic Spectral Tools

Before delving into the design of spectrally accurate discrete operators, we must
define a series of basic tools and conventions which will be particularly useful for
our task. We start by defining proper periodic basis functions on regular grids,
before describing the case of Chebyshev grids.
1D Periodic Interpolator and Histopolator Functions
To build our spectral wedge and Hodge star operators to work on discrete forms
of arbitrary degree, we will need not only spectral interpolating basis functions,
but also spectral histopolating basis functions, i.e., basis functions which integrate
to one over assigned intervals [10]. To this end, we consider a one-dimensional
periodic domain of width 2π with N regularly-spaced nodes, and define over this
canonical domain two scalar functions αN and βN —that we will respectively call
interpolator and histopolator—as follows:
Nx
1⎪
⎪cot sin 2
αN (x) = ⎨ 2
Nx
N⎪
⎩csc 2 sin 2
and

if N even,

(3.2)

if N odd,

N/2
Nx
n cos nx
cos
N ∑ sin nπ
⎪ 2π 4
n=1
βN (x) = ⎨
(N−1)/2
cos nx
nπ
∑ nsin
2π + N
n=1

if N even,
(3.3)
if N odd.

For convenience, let us denote any translation of the above functions using the
notation below:
αN,n (x) = αN (x − hn),

and

βN,n (x) = βN (x − hn),

where h = 2π/N is the mesh’s edge element size and the nodes are enumerated by
xn = nh. For any given number N of nodes, these two functions satisfy the following
important properties mentioned in Eq. (3.1) (where δmn refers to the Kronecker
delta):
ˆ xm + h
βN,n (x) dx = δmn .
αN,n (xm ) = δmn, and
xm − h2

In other words, αN provides a smooth interpolation of a discrete function with 1
at node x0 and 0 at every other node, while βN integrates to 1 over the interval
[x0 − 2h , x0 + 2h ], and to 0 over all other intervals (see Figure 3.3). Notice that these

19
functions are the only Fourier series with N sinusoidal components satisfying the
point-wise (resp., interval-wise) constraints. Note also that αN and βN are related
to each other:
ˆ x+ h
βN′ (ξ)dξ = αN (x), or equivalently, βN (x + ) − βN (x − ) = αN′ (x).
x− h

Next, we show that these two functions provide building blocks to construct basis
functions for arbitrary forms on regular, periodic grids—from which we will derive
basis functions for bounded domains via pushforward.

(a) α6 (x)

(b) α7 (x)

(d) hβ7 (x)

Figure 3.3: Examples of periodic interpolators for N = 6 (a) and 7 (b), and corresponding
periodic histopolators (c) and (d), scaled by h = 2π/N for clarity. While the interpolator α N
satisfies α N (nh mod 2π) = δ0n ∀n, the histopolator β N integrates to 1 over the dual cell
straddling x = 0, and to 0 over other dual cells in the range [0, 2π]. Note that the alternating
red and green colors are used to mark out dual cells, and to illustrate that the integral of β N
over each of these dual cells sums to zero or one.

Spectral Basis Functions on Regular Grids
Basis functions for periodic grids in arbitrary dimensions can be easily built through
tensor products of (translated) interpolators αN and histopolators βN , where the
number of βN used in the tensor product is equal to the degree of the form. To make
this point clear, we will now introduce our notation for 1-, 2-, and 3-dimensional
basis functions, with the extension to higher dimensions being straightforward. We
p[comp]
will denote the basis functions as φ[ind] with superscripts used to indicate the
degree p, followed by the component of the form when appropriate (e.g., x, x y,
etc), and subscripts used for grid indices; for instance, φ1y
mnk (.) is the function of R
representing the dy-component of the 1-form basis, located on the edge parallel to
the y axis indexed by (m, n, k). See Figure 3.2 for an illustration of the 2D notations.
Primal and Dual Nodes. On a 1D regular grid the N primal nodes are equally
spaced, and the dual nodes are placed in the middle between two adjacent primal
nodes. Consequently, the primal nodes are: xn = 2πn/N, and the dual nodes are at
xn = 2π(n + 12 )/N.

20
Basis Functions for Forms in 1D. Because of their interpolation (resp., histopolation) property, translated versions of the functions αN and βN can directly be used
as basis functions for 0- and 1-forms respectively as follows:
φ0N,n (x) = αN,n (x),

and

φ1N,n (x) = βN,n+ 1 (x) dx,

where n ∈ {0, 1, . . . , N −1}. Indeed, a discrete 0-form ω̄0 (resp., a discrete 1-form ω̄1 )
can be reconstructed as a smooth form through ω0 = ∑i ω̄i0 φi0 (resp., ω1 = ∑i ω̄i1 φi1 );
the reconstructed form then satisfies Pω0 = ω̄0 (resp., Pω1 = ω̄1 ). Similarly, dual
basis functions are easily designed as well through
φ̃0N,n (x) = αN,n+ 1 (x)

and

φ̃1N,n (x) = βN,n (x) dx.

Basis Functions for Forms in 2D. In two dimensions, tensor products of the onedimensional bases provide basis functions for 0-, 1-, and 2-forms on a regular M×N
grid. These functions are expressed as follows:
φ0MN,mn (x, y) = φ0M,m (x) ∧ φ0N,n (y),
φ1x
MN,mn (x, y) = φ M,m (x) ∧ φ N,n (y),
φ1y
MN,mn (x, y) = φ M,m (x) ∧ φ N,n (y),

φ2MN,mn (x, y) = φ1M,m (x) ∧ φ1N,n (y).
One can easily check that these functions are 1 on their associated degree of freedom
and 0 on all others (vertices for 0-forms, edges for 1-forms, and faces for 2-forms),
thus offering a proper set of bases for smooth reconstructions of discrete forms.
Formulas for the dual basis functions are strictly analogous, where φ̃ is used in lieu
of φ.
Basis Functions for Forms in 3D. The same construction can be used in three
or higher dimensions. For completeness, we describe the three-dimensional basis
functions for primal forms:
φ0MNK,mnk (x, y, z) = φ0M,m (x) ∧ φ0N,n (y) ∧ φ0K,k (z),
φ1x
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),
φ1y
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),
φ1z
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),
φ2x
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),
φ2xz
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),
φ2yz
MNK,mnk (x, y, z) = φ M,m (x) ∧ φ N,n (y) ∧ φK,k (z),

φ3MNK,mnk (x, y, z) = φ1M,m (x) ∧ φ1N,n (y) ∧ φ1K,k (z).

21
Note that since the wedge products above are the continuous (as opposed to discrete) ones, associativity holds, and no ambiguity is introduced by the omission of
parentheses. Here again, it is easy to check that these functions provide smooth
reconstructions from values of vertices, edges, faces, or volumes on a regular grid
indexed by m, n, and k, in a periodic domain. Formulas for the dual basis functions
are strictly analogous, where φ̃ is used in lieu of φ.
Chebyshev Grids over Bounded Domains
For bounded domains, a popular choice of spatial discretization in spectral methods
is the use of Chebyshev computational grids [17]. It is well known that Chebyshev
polynomials can be derived as the pullback of Fourier basis functions from a circle
onto its diameter; see Figure 3.4. We can, in fact, perform the same pullback of
our regular-grid basis functions of forms to obtain new spectral bases of forms
applicable for Chebyshev grids over non-periodic domains.
x = −1

x=1

θ = 0, 2π

θ=π

Figure 3.4: The map ϕ mapping the canonical interval [−1, 1] to the bottom unit hemicircle

(0 ≤ θ ≤ π).

Denote by ϕ the map from the interval between −1 and 1 on the real line (with
Cartesian coordinate x) to the unit semicircle (with polar angle θ, see Figure 3.4):
ϕ∶

[−1, 1] ⊂ R → [0, π] ⊂ S 1
x ↦ θ = acos(−x).

Primal and Dual Nodes. Let the 1D grid consist of N points (including the
boundaries) in the interval [−1, 1]. The corresponding unit circle will then have
2N − 2 points (see Figure 3.4). The primal nodes {xn } (n = 0..N−1) and dual nodes
{̃
xn } (n = 0..N −2) of the Chebyshev grid are thus
nπ
xn = − cos
N −1

and

(n + 21 )π
xn = − cos
N −1

22
Basis functions. To design our basis functions, we first define functions on the
semicircle (θ ∈ [0, π]) by mirroring/antimirroring the regular basis functions φ on
the whole circle (see [69] for the usual case of primal 0-forms), in order to satisfy the
interpolation/histopolation properties (Eq. (3.1)) on the semicircle. The resulting
functions κ are expressed as:

⎪φ
κN,n
= ⎨ 2N−2,n
⎩φ2N−2,n + φ2N−2,2N−2−n,

n = 0 or n = N − 1 (endpoints)
n ∈ {1, . . . , N − 2} (midpoints)

for primal 0-forms,
κN,n
= φ12N−2,n − φ12N−2,2N−3−n,

n ∈ {0, . . . , N − 2}

for primal 1-forms,
κ̃N,n
= φ̃02N−2,n + φ̃02N−2,2N−3−n

n ∈ {0, . . . , N − 2}

for dual 0-forms, and
δN (θ)
κ̃N,n
(θ) = ⎨[φ̃02N−2,n − φ̃02N−2,2N−2−n − ρN,n ] (θ)
⎩δN (π − θ)

n=0
n ∈ {1, . . . , N − 2}
n= N −1

for dual 1-forms, where
cos ((N − 1)θ) ) sin θ dθ,
ρN,n (θ) = 2 (γ2N−2,n δN (θ) + γ2N−2,N−n−1 δN (π − θ)) dθ,

δN (θ) = ((N − 1)2 α2N−2,0 (θ) +
and

ˆ π
with: γN,n =

βN,n (θ) dθ

N/2
1−(−1)k
sin(2knπ/N)−sin((2k−1)nπ/N)
+ N1 ∑
2N
sin nπ
n=1
= ⎨
(N−1)/2
sin(2knπ/N)−sin((2k−1)nπ/N)
+ N1 ∑
2N
sin nπ
n=1

if N even,
if N odd.

The function δ is used to deal with the special case of the two boundary dual
(half)edges, and the function ρ adds contributions to intermediate basis functions

23
so that they integrate to zero at both boundary dual edges. Finally, we pull back
the functions κ by ϕ to obtain the form basis functions ψ on the Chebyshev 1D grid
(ψ = ϕ∗ κ):
ψN,n
(x) = κN,n
(acos(−x)),

dx
ψN,n
(x) = κN,n
(acos(−x)) √
1 − x2
ψ̃N,n
(x) = κ̃N,n
(acos(−x)),
dx
ψ̃N,n
(x) = κ̃N,n
(acos(−x)) √
1 − x2
One can easily check that the functions above satisfy the property in Eq. (3.1)
required for basis functions, as the functions κ were designed to satisfy these properties, and the pullback ϕ∗ commutes with integration.

(a) ψ00 (x)

(b) ψ10 (x)

(d) ψ30 (x)

(e) (x1 − x0 )ψ01 (x) (f) (x2 − x1 )ψ11 (x) (g) (x3 − x2 )ψ21 (x)

Figure 3.5: Chebyshev primal basis functions for a grid with N = 7. We normalize
the one-form basis functions by xn − xn−1 to have approximately the same scale in our
visualizations.

Note that the basis functions for primal zero-forms turn out (unsurprisingly) to be
the Lagrange polynomials of order N, i.e.,
N−1 x − x

ψN,n
(x) = ∏

m=0 xn − xm

m≠n

where {xn }n=0,...,N−1 represent the coordinates of the primal points. All other basis
functions of forms are also polynomials for any choice of N. Examples for the primal
and dual basis functions for 0- and 1-forms are provided in Figures 3.5 and 3.6 for
N = 7.

(a) ψ00 (x)

(b) ψ10 (x)

(d) ψ30 (x)

(e) (x1 − x0 )ψ01 (x) (f) (x2 − x1 )ψ11 (x) (g) (x3 − x2 )ψ21 (x)

Figure 3.6: Chebyshev dual basis functions for a grid with N = 7
Finally, basis functions in higher dimensions for arbitrary forms are derived using
tensor products of these two primal and two dual 1D basis functions, just as we
explained in Section 3.3.
3.4

Spectrally Accurate Discrete Operators

Equipped with these basis functions on regular and Chebyshev grids, we can now
derive the discrete, spectrally accurate version of the exterior derivative d, the wedge
product ∧, and the Hodge star ⋆. Emphasis is placed on the mathematical definitions;
fast evaluations of these operators will be discussed subsequently in Section 3.5.
Discrete Exterior Derivative D
The most common discrete realization of the exterior derivative based on algebraic
topology [49] using chains and cochains is a linear operator D:
D∶

Λ̄ k → Λ̄ k+1,

which is the (signed) incidence matrix between (k + 1) elements and k elements of
the grid, with a sign determined by the relative orientation of the elements. In other
words,
D = ∂ t,
where ∂ refers to the boundary operator acting on chains (see [49, 21]). Note that
the operator D is thus implemented via a sparse matrix whose non-zero elements
have only +1 and −1 values to indicate incidence between mesh cells. It also satisfies
D2 = 0 like its continuous equivalent (since the boundary of a boundary is always

25
the empty set), and that a discrete Stokes’ theorem is also enforced on chains since
the very definition of D is tantamount to enforcing
dω =

∂σ

for every mesh element σ. This discrete realization of the exterior derivative is
exact, in the sense that the operator D commutes with the reduction operator:
dP =PD.

(3.4)

This implies that the following diagram fully commutes:
̃k o

/ Λk o

̃ k+1 o
Discrete
Dual Forms

/ Λ̄ k

Λ k+1
Differential
Forms

/ Λ̄ k+1

Discrete
Primal Forms

Therefore, as mentioned in Section 2.2, the classic discrete exterior derivative used
in mimetic methods, Dec, or finite-dimensional exterior calculus needs no special
treatment to obtain spectral accuracy. However, its use for spectral computations
represents a clear departure from the conventional spectral methods [10], since all
operators of classical vector calculus (divergence, gradient, and curl) can be achieved
exactly via D [28]. In our framework, we also use the exterior derivative on dual
forms:
̃k → Λ
̃ k+1 .
̃∶ Λ
Its implementation and properties are no different from its primal version, since the
adjacency of the dual mesh is directly derived from the adjacency on the primal:
̃t .
̃ =∂
We now turn to the definition of a discrete wedge product and discrete Hodge star
operator.
Discrete Wedge Product W
Various discrete definitions of the wedge product have been proposed, sometimes
mixing primal and dual elements [35]. We instead follow Bochev and Hyman’s
treatment [7]: by utilizing the de Rham and reconstruction maps, one can define a

26
discrete wedge product W(ᾱ, β̄) in a general way that applies to any pair of discrete
(primal or dual) forms ᾱ and β̄ through:
W(ᾱ, β̄) = P(Rᾱ ∧ R β̄) .

(3.5)

That is, we first reconstruct the two discrete forms with our trigonometric or Chebyshev basis functions, we then apply the wedge product to the two continuous forms
that we have created, and finally, integrate (or, in the case of 0-forms, point-sample)
the result to get the final discrete form W(ᾱ, β̄). The resulting operator satisfies a
discrete Leibniz rule, equivalent to Eq. (2.1):
DW(ᾱ, β̄) = W(D(ᾱ), β̄) + (−1) k W(ᾱ, D( β̄)),
where ᾱ is a discrete k-form, and β̄ is an arbitrary discrete form. Indeed, because
the exterior derivative commutes with the reduction and reconstruction maps, this
discrete Leibniz rule can be derived directly from the continuous Leibniz rule.
Our spectrally accurate definition of the discrete wedge product is bilinear and
anticommutative, but it is not associative. It is, however, associative in the limit as the
mesh size approaches zero. That is, the associator W(ᾱ, W( β̄, γ̄)) − W(W(ᾱ, β̄), γ̄)
(i.e., the measure of nonassociativity) will approach zero exponentially fast as the
step size is reduced. The difficulty of creating an associative discrete wedge product
has been previously noted by, e.g., Kotiuga [40], who refers to it as the “commutative
cochain problem” and discusses some deeper topological reasons behind it. The
spectral accuracy of our wedge operator will mitigate this lack of associativity
exponentially fast under mesh refinement.
Note that the discrete wedge product of a (primal or dual) p-form and a (primal or
dual) q-form can also be seen as a three tensor
p,q
Wi j k =
φ pj ∧ φqk,
σip+k

since we can write the value of the wedge product on element σip+q as :
[W(ᾱ, β̄)]i = ∑ Wi j k ᾱ j β̄k .
j,k

Discrete Hodge Star H
The different methods to discretize the Hodge star operator fall mostly in two
categories. A majority of approaches start from a discretization of the inner product
of forms (Eq. (2.3)), then derive a discrete Hodge star from it. This leads to

27
what is called the “derived discrete ⋆” in [7]. A few authors, instead, use a direct
(called “natural” in [7]) definition of the Hodge star, which exploits the existence of
primal/dual structures on meshes; i.e., the discrete Hodge star H now maps discrete
primal forms to discrete dual forms, and vice-versa through the inverse of the Hodge
star. Discretizations of the Hodge star typically rely on local reconstructions and
integration, the simplest one being what is often referred to as the diagonal Hodge
star [8], corresponding to a piecewise constant reconstruction. Its implementation
involves each primal value being multiplied by the ratio of the measure of the primal
and dual mesh elements to obtain the dual value. It is hence represented by a
diagonal matrix. Higher-order Galerkin Hodge stars have also been considered,
resulting in better accuracy at the cost of denser matrices. As we will see, by
choosing proper basis functions, we will be able to construct a spectrally accurate
natural discrete Hodge star that is computationally efficient, as the star will become
a Toeplitz matrix in the periodic case.
As explained earlier, our discrete Hodge star H exploits the notion of mesh duality
in that the discrete Hodge star of a primal k-form is a dual (d − k)-form, and viceversa. The discrete Hodge star for a discrete form is realized conceptually by first
reconstructing the continuous form with our primal (resp., dual) spectral bases,
applying the continuous Hodge star to this form, and then projecting this form back
to the dual (resp., primal) grid. In our notation, this can be written as
̃ ⋆ k R,
Hk = P

(3.6)

and is easily understood from the following diagram:
Λk o

Λ̄ k

⋆k

Λn−k

/Λ
̃ n−k

The operator H can be built as a matrix with easily precomputed entries. For
instance, for k-forms, this matrix contains the terms
Hi j =
⋆φ kj .
(3.7)
σin−k

Similarly, the dual Hodge star operator is given by:
̃ =
⋆φ̃kj .
ij
σin−k

(3.8)

28
Notice that if we require our discrete Hodge stars to satisfy the discrete equivalent
of Eq. (2.2):
̃ n−k = H
̃ n−k H k = (−1) k(n−k) Id,
Hk H
(3.9)
this imposes a constraint on dual basis functions. Indeed, for Eq. (3.9) to hold true,
they must be a linear combination of the primal basis functions:
φ̃ik = ∑[(Hn−k )−1 ] ji ⋆ φn−k
j .

(3.10)

Both the regular (φ) and Chebyshev (ψ) basis functions that we defined above satisfy
this constraint.
Finally, our Hodge matrices are dense; however, each matrix is circulant for a regular
grid on periodic domains (due to the invariance of the grid under translation),
and centrosymmetric for a Chebyshev grid (due to mirror symmetry around its
center). While circulant matrices are special cases of Toeplitz matrices, for which
multiplications by vectors can be done O(N log N), we will show that all discrete
Hodge stars for arbitrary dimensions can be efficiently implemented in Fourier space
in O(N log N) time complexity where N is the total number of spatial points.
Other Operators
There are important additional operators, including contraction by vector field, and
Lie derivative (extending directional derivatives to differential forms). While these
operators play a crucial role in certain applications such as computational fluids and
other advection problems, their extension to our spectral framework brings additional
difficulty due to their dependence on a vector field. They can, however, be handled
in our framework, but only through the use of coordinates, whereas in the continuous
world the application of a Lie derivative on a differential form is a coordinate-free
metric-independent operation. A first step towards a geometric discretization of
contractions and Lie derivatives was proposed in [48]; but a spectrally accurate
discretization preserving key properties of these other operators is left for future
work.
3.5

Implementation in Fourier Space

We now describe how to efficiently implement the operators described in the last
section. We will see that the Hodge star and wedge operators can have very fast
implementations as they can be expressed as the product of fast Fourier transforms
and sparse matrices. On periodic, regular grids of size N d , the basic Dec operators

29
can be expressed in terms of the discrete equivalents of two elementary operators,
namely a convolution with a box function over [a, b] (denoted Ia,b ) and a translation
by a (denoted Ta ). Both of these operations are expressed as diagonal matrices
in Fourier space. The wedge product and Hodge star for Chebyshev grids require
additional sparse operators that we describe next.
Sparse Matrices for Periodic 1D Domains
We define our convolution and translation operators to act on an array of function values f¯i at evenly spaced points xi on a periodic grid by first constructing a
trigonometric interpolation R f¯ as described in Section 3.2, then by transforming it
as follows:
ˆ xi +b
a,b ¯
(R f¯) (ξ) dξ,
I ∶ fi ↦
(integration)
(3.11)
xi +a

Ta ∶ f¯i ↦ (R f¯) (xi + a).

(translation)

(3.12)

These two operators are implemented through discrete Fourier transforms as follows:
Ia,b = F −1̂Ia,b F,

̂ a,b F .
Ta,b = F −1 T

The operator F represents the shifted discrete Fourier transform and F −1 represents
its inverse: rather than placing the zero-frequency term (the mean of the signal)
first, it is shifted to the middle with the negative Nyquist frequency becoming the
first element, and the positive Nyquist frequency becoming the last.
In Fourier space (denoted using a hat ̂. ), each operator is simply encoded by a
diagonal matrix with the following coefficients:
ikb − eika

̂Ia,b = e
kk

̂ a = eika,
, and T
kk
ik
where the special case k = 0 is treated as:
eikb − eika
= b − a.
k→0
ik
lim

Note that since these two matrices are diagonal, the order in which they are applied
does not matter.
Sparse Matrices for 1D Chebyshev Grids
For conciseness in later expressions, we also define a number of sparse matrices
representing simple linear operators that will be particularly useful for the Hodge
stars and wedge products on Chebyshev grids.

30
Mirroring The mirroring matrices will allow us to map quantities from the Chebyshev grid to the unit circle and back (see Fig 3.4). We will need four versions of
mirroring: a mirroring M+0 and anti-mirroring M−0 operator for 0-forms mapping
N nodal values of the Chebyshev grid to 2N − 2 values on the circle; and similarly, a mirroring M+1 and anti-mirroring M−1 operator mapping N nodal values of
the Chebyshev grid to 2N values on the circle will be useful for 1-forms. Their
expressions are:
⎡ x ⎤
⎢ 0 ⎥
⎡ x ⎤ ⎢⎢ x1 ⎥⎥
⎢ 0 ⎥ ⎢
⎢ x ⎥ ⎢⎢ ⋮ ⎥⎥
⎢ 1 ⎥ ⎢
⎥ ⎢ xN−2 ⎥⎥
⎥,
M±0 ∶ ⎢⎢ ⋮ ⎥⎥ ↦ ⎢
⎥ ⎢⎢ xN−1 ⎥⎥
⎢ xN−2 ⎥ ⎢
⎥ ⎢
⎥ ⎢±xN−2 ⎥⎥
⎢ xN−1 ⎥ ⎢
⎦ ⎢ ⋮ ⎥⎥
⎢ ±x1 ⎥

⎢ x0 ⎥
⎢ x ⎥
⎢ 1 ⎥
⎡ x ⎤ ⎢⎢ ⋮ ⎥⎥
⎢ 0 ⎥ ⎢
⎢ x ⎥ ⎢⎢ xN−2 ⎥⎥
⎢ 1 ⎥ ⎢
⎥ ⎢ xN−1 ⎥⎥
⎥.
M±1 ∶ ⎢⎢ ⋮ ⎥⎥ ↦ ⎢
⎥ ⎢⎢±xN−1 ⎥⎥
⎢ xN−2 ⎥ ⎢
⎥ ⎢
⎥ ⎢±xN−2 ⎥⎥
⎢ xN−1 ⎥ ⎢
⎦ ⎢ ⋮ ⎥⎥
⎢ ±x1 ⎥
⎢ ±x0 ⎥

The left inverses of the above operators, satisfying M†0 M±0 = Id and M†1 M±1 = Id, are
then given by

M†0 ∶

⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥ ⎡ x ⎤
⎢ 1 ⎥ ⎢ 0 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢ x ⎥
⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ xN−1 ⎥ ↦ ⎢ ⋮ ⎥,
⎥ ⎢
⎥ ⎢
⎢ xN ⎥ ⎢ xN−2 ⎥
⎥ ⎢
⎥ ⎢
⎢ ⋮ ⎥ ⎢ xN−1 ⎥
⎥ ⎣
⎢ x2N−1 ⎥

M†1 ∶

⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥ ⎡ x ⎤
⎢ 1 ⎥ ⎢ 0 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢ x ⎥
⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ xN−1 ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎢ xN ⎥ ⎢ xN−2 ⎥
⎥ ⎢
⎥ ⎢
⎢ ⋮ ⎥ ⎢ xN−1 ⎥
⎥ ⎣
⎢ x2N ⎥

Extensions Extension matrices will allow us to resample forms at higher or lower
resolutions, an operation needed in the wedge product as we will mention in Seĉ n that will act in Fourier
tion 3.5. We define the upsampling extension matrix E

31
space to add 2n samples as follows:

∀n > 0,

̂n ∶

⎡ ⎫
⎢ 0 ⎪
⎢ ⎪
⎢ ⋮ ⎪
⎢ ⎬ n⎥
⎢ ⎪
⎢ 0 ⎪
⎢ ⎪
⎢ ⎭ ⎥
⎡ x ⎤
⎢ x ⎥
⎢ 0 ⎥
⎢ 0 ⎥
⎢ x ⎥
⎢ x ⎥
⎢ 1 ⎥
⎢ 1 ⎥
⎢ ⋮ ⎥ ↦ N + 2n ⎢ ⋮ ⎥.
N ⎢
⎢ xN−2 ⎥
⎢ xN−2 ⎥
⎢ xN−1 ⎥
⎢ xN−1 ⎥
⎢ ⎫ ⎥
⎢ 0 ⎪
⎢ ⎪
⎢ ⎪
⎪ ⎥
⎢ ⋮ ⎬ n⎥
⎢ ⎪
⎢ ⎪
⎢ 0 ⎪
⎣ ⎪
⎭ ⎦

̂ n† such that
The left inverse of this extension matrix is a downsampling matrix E
̂ n† E
̂ n =Id, and is expressed as

̂ n† ∶

⎡ x ⎫
x0 + xN
−n ⎪
⎬ n⎥⎥
⋮ ⎬ n ⎥⎥
⎢ x +x
⎢ x ⎪
n−1
N+n−1
−1 ⎪
⎭ ⎥⎥
⎭ ⎥
xn
x0
⎥ N − 2n ⎢
⎥.
⎥↦ N ⎢
N−n−1
N−1
⎢ xN ⎫
N−n
−n ⎪
⎬ n ⎥⎥
⎬ n⎥
⎢ xN−1 + x−1 ⎪
⎢ xN+n−1 ⎪
⎭ ⎦

With these extensions matrices, we can upsample a discrete form ω̄ from a N-point grid to a
(N+2n)-point grid by first performing a Fourier
̂n , and applying
transform FN , then applying E
−1
an inverse Fourier transform FN+2n
(see inset
for an illustration of E and its left inverse):
−1 ̂ n
En ω̄ = FN+2n
E FN ω̄.

Downsampling is done similarly through En† , with now FN−1 and FN+2n .

32
Scaling The metric dependent matrix ΩN is a diagonal matrix which has, as its
diagonal elements, the weight by which each nodal value needs to be scaled in order
for the integration over its corresponding dual region on the unit circle to faithfully
represent integration over its corresponding dual region on the original Chebyshev
grid. Unlike previous matrices in this section, ΩN is metric-dependent in the sense
that it does not only depend on the grid topology, but on the volume change between
the grid on the circle and the Chebyshev grid:
2π
Ω k k = sin ( (k + 2 ))) = 1 − ̃
x k2 .
Other Matrices Finally, we will need simple, additional matrices denoted A, B,
B† , and C; they respectively lump every consecutive pair of non-boundary dual edge
values, evaluate the nodal values of the 1-form basis functions of the left and right
boundary dual edges (i.e., δN (θ) and δN (π − θ) in Section 3.3), and remove the end
node values. Their expressions are thus:
⎛1
⎜ 1 1
1 1
A = ⎜⎜
⋱ ⋱
1 1⎟
1⎠

´¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¶
N×(2N−2)

⎛(N − 1)2 + 1/2
−1/2
1/2
B = ⎜⎜
−1/2
1/2
⎝ (−1)N+1 /2

´¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¶
N×N

⎛0
⎜ 1
⎟,
C = ⎜⎜
1 ⎟
0⎠

´¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¶
N×N

† ⎜
B = ⎜⎜

(−1)N+1 /2 ⎞
1/2
⎟.
−1/2
1/2
−1/2
(N − 1) + 1/2⎠

Wedge Product W
Let ⊙ denote the element-wise multiplication (typically used in Fourier space) such
that, for two vectors a and b, a ⊙ b is the vector with components
(a ⊙ b)i = ai bi .
Remember that we are (conceptually) working with bandlimited functions f (x)
whose frequencies are below the Nyquist rate for the spatial grid resolution; i.e., with

33
functions that can be reconstructed exactly from their sample values {ai = a(xi )} on
the grid. However, when the element-wise product ⊙ is applied to such functions, the
result is no longer necessarily below the Nyquist rate, as the highest frequency can
potentially be doubled by the multiplication. In order to compute wedge products
exactly within the Nyquist range, we will have recourse to a temporary upsampling
of the operands to a higher grid resolution, perform the multiplication at the higher
resolution, and downsample back to the original grid resolution.
1D Wedge Based on our definition in Eq. (3.5), a wedge product of two 0-forms
in 1D is expressed as
ω̄0 = W(ᾱ0, β̄0 ) = ᾱ0 ⊙ β̄0 .
Note that there is no reason to use the up/downsampling operators here: while we
obviously cannot capture the wedge product content past the Nyquist rate, all the
frequencies below are exactly evaluated, thus leading to spectral accuracy. However,
the wedge product of a 0-form and a 1-form differs, in the sense that to avoid losing
accuracy, we must use the following definition:
ω̄1 = W(ᾱ0, β̄1 ) = EN† J ((EN ᾱ0 ) ⊙ (EN J−1 β̄1 )),
where J maps from 0-forms to 1-forms (on the periodic domain) through
J = I− 2 , 2 T 2 .
h h

Indeed, J being an integration, the wedge without EN would lose accuracy at the
Nyquist rate. Should dual forms be used, they could easily be converted to primal
forms via the translation operator T before the wedge product is applied. Therefore,
wedge products in 1D can always be implemented with efficient evaluation through
discrete fast Fourier transforms, both for periodic and Chebyshev grids, since their
expressions hold in each case.
2D and Beyond The various wedge products between forms can be expressed in
terms of the following functions:
w00 (a, b) = a ⊙ b,
w01 (a, ( bbyx )) = ( a⊙b
a⊙by ),

w11 (( aayx ), ( bbyx )) = a x ⊙ b y − a y ⊙ b x,
and similarly for higher dimensions.

34
Now, let:
Jx ≡ J ⊗ Id,

Jy ≡ Id ⊗ J,

E ≡ EN ⊗ EN ,

Jxy ≡ J ⊗ J,

E† ≡ EN† ⊗ EN† .

Then we can define the wedge product in 2D for all possible pairs of forms as:
W(ᾱ0, β̄0 ) = w00 (ᾱ0, β̄0 ),
β̄1

W (ᾱ0, ( β̄x1 )) =

⎛E† Jx
⎞ 01
EJ −1 β̄1
w (Eᾱ0, ( EJx −1 β̄x1 )),
E Jy ⎠

ᾱ1

β̄1

EJ −1 ᾱ1

EJ −1 β̄1

W (( ᾱx1 ), ( β̄x1 )) = E† Jxy w11 (( EJx −1 ᾱx1 ), ( EJx −1 β̄x1 )),

W(ᾱ0, β̄2 ) = E† Jxy w00 (Eᾱ0, EJxy −1 β̄2 ) .
Here again, the use of up/downsampling is purely to ensure no loss of accuracy
for frequencies at or below the Nyquist rate. Extension of the wedge product to
dual forms or in higher dimensions is straightforward as the same principles apply.
While the above derivation was for periodic domains, the wedge product is preserved
under pullback; so the equations hold for any topologically equivalent space. For
Chebyshev grids, J is no longer square, and there are only N − 1 edges—but the
formulas are basically unchanged.
Hodge Star H
We can also express our discrete Hodge stars in a way that lends itself to fast
evaluations.
Periodic Grid in 1D The Hodge star for a periodic grid is very simple, as it is
simply expressed as
̃ 0 = I− h2 , h2
H0 = H

̃ 1 = I− h2 , h2
H1 = H

−1

(3.13)

Chebyshev Grid in 1D To define the Hodge-Star operator on a Chebyshev grid,
we make use of the matrices we defined earlier in Section 3.5. In particular, we get
̃ 0 = M1 † I− h2 , h2 ΩM+,
h h −1

H1 = M1 † Ω−1 I− 2 , 2

M−1 .

The expressions for the other two stars are slightly more complicated because of the
special treatment needed for the dual (half-)edges at the two ends of the domain:

H0 = AM†0 I0,+ 2 ΩEN−1 M+0
̃ 1 = M† (T− h2 Ω−1 T h2 − BI0, h2 − B† I− h2 ,0 ) I− h2 , h2

−1

M−0 C + B + B† .

35
2D and beyond The extension of the fast 1D implementation to arbitrary dimension is rather direct: we construct Hodge stars on higher-dimensional grids by taking
the Cartesian products of stars on one-dimensional grids. Using 2D for illustration
purposes, we first introduce the following operators
Hx ≡ H ⊗ Id,

Hy ≡ Id ⊗ H.

where Hx acts along the first index (x-direction), and Hy acts along the second index
(y-direction). Then from these two 1D Hodge stars, the 2D Hodge star operators
can be expressed as
H0 = H0x H0y,
H1 =

⎛ 0
−H0x H1y, ⎞
⎝ H1x H0y

H2 = H1x H1y .
Equivalent equations hold for the dual Hodge stars. Note that we have H1 H1 = −Id
as desired. For completeness, we also provide the 2D exterior derivative in terms of
the Cartesian products of the 1D Dx and Dy operators:
x ᾱ
),
D0 ᾱ0 = ( D
Dy ᾱ0

ᾱ1

D1 ( ᾱx1 ) = −Dy ᾱx1 + Dx ᾱy1

3.6

Numerical Tests

We now discuss how our spectral operators can be used for practical computations.
In particular, we present numerical tests involving Poisson equations on 0- and
1-forms in 1D and 2D.
Solving Differential Equations
Our spectral framework can be used in various ways, either as is, or with minor
modifications.
Chain Collocation Our computational framework lends itself naturally to a spectral method that we call the Chain Collocation Method. Just as the point collocation
method enforces a partial differential equation on all grid nodes, we now enforce
partial differential equations on all appropriate chains. If the unknown is a scalar
function, then the PDE is enforced at every node as usual; but for a k-form, the PDE

36
will now be enforced on each k-chain. Since any differential equation can be rewritten with exterior forms, our operators can be used as is on discrete approximations of
these forms, and the conversion of the continuous operators into discrete operators
combining D, H, and W will enforce the differential equation on each appropriate
mesh element. For instance, a 2D Poisson equation ∆v = f for a vector field v can be
rewritten with, for instance, a 1-form ω as ∆ω = γ where ω is now the circulation
of the continuous field v on all edges, while γ is the circulation of f—leading to
a coordinate-free setup, drastically different from the usual coordinate-based treatment. Modulo boundary conditions, a direct discretization of the Laplace operator
∆ = d ⋆ d ⋆ + ⋆ d ⋆ d will now enforce that a discrete form ω̄ satisfies the PDE in a
discrete sense, i.e., on each edge. (Various 1D and 2D examples of this approach are
shown at the end of this section.) Because of the algebraic topology properties of
our framework, conservative differential equations will lead to conservative discrete
equations, and structure-preservation will be enforced.
Other Numerical Approaches Our framework can be easily altered to accommodate other numerical solves. First, different grids can be used if necessary, as we
spelled out the constraints that one needs to enforce to get proper discretizations.
Second, one can modify our Hodge stars by deriving them from the inner product
instead: this is the approach spelled out in [55, 28]. While the Hodge star operators
are no longer converting primal to dual forms (and vice-versa), the advantage of this
approach is that the resulting operators inherit the symmetry of the inner product.
Logically rectangular grids Finally, modifying our approach for arbitrary logically rectangular computational grids is simple. First, both D and W are unchanged
since the integral of forms are unchanged by the pullback ϕ∗ and our numerical
treatment of these operators is metric independent. The Hodge star operator needs,
however, to be modified. In fact, only the matrices ΩN require modification: one
simply has to rescale the diagonal elements of ΩN by how much the orthogonal
Chebyshev grid is deformed into the computational grid. In this sense, our spectral
framework is reminiscent of the regular Dec framework, as the diagonal Hodge
stars are also subject to be changed based on the volume forms of the computational
domains [21].

37
Primal
Dual

1e0

f−∆−1 q ∞

1e-4

f−∆−1 q ∞

1e-4

Primal, Dirichlet bc
Dual, Neumann bc

1e0

1e-8

1e-12

100

(a) Periodic Domain

100

(b) Bounded Domain

Figure 3.7: Convergence graphs for a 1D Poisson equation: (a) we solve ∆ f =

esin x (cos2 x − sin x) on a periodic domain for either a primal, or dual 0-form f ; (b) we
solve ∆ f = e x on a Chebyshev grid, for either a primal 0-form with Dirichlet boundary conditions f (−1) = e−1 and f (1) = e, or for a dual 0-form with Neumann boundary conditions
f ′ (−1) = e−1 and f ′ (+1) = e. All of our results exhibit spectral convergence (measured
through the L ∞ error ∥ f − ∆−1 q∥∞ ), with the conventional plateau when we reach the limit
of accuracy of the representation of floating point numbers.

Convergence Tests
In order to provide numerical evidence of spectral behavior, we tested our operators
on a variety of examples. We show below a series of Poisson equations with Dirichlet
and Neumann conditions in 1D and 2D.
1D Poisson The traditional Poisson equation for a scalar function f reads:
∇2 f = q,
which can be rewritten in exterior calculus as an equation on a 0-form f :
⋆d ⋆ d f = q.
Using a periodic domain (where a regular grid is placed at [0, 2π]) or a bounded
domain (chosen to be [−1, +1], discretized by a Chebyshev grid), we can directly
discretize this last equation. The only choice we are left with is to decide whether
to think of f as a primal ( f¯) or a dual ( f˜) 0-form. Thus, the 1D Poisson equation is
implemented as either
̃ 1D
̃ 1 H1 D0 f¯0 = q̄0
or
̃ 1D
̃ 0 f˜0 = q̃0 .
H 1 D1 H

38
Note that the boundary conditions (Dirichlet or Neumann) are simply incorporated
into the D operators: prescribed primal values or dual values at the boundary will be
used as is instead of being considered unknowns. This is substantially simpler than
the traditional use of, e.g., generalized Lagrange polynomials for the treatment of
boundary conditions. The convergence plots are shown in Fig 3.7, proving spectral
convergence.
2D Poisson We also tested our approach on 2D problems. This time, we computed
Poisson equations not only for scalar functions, but also for vector fields—or in
exterior calculus jargon, we can compute the Laplacian ∆ = ⋆d ⋆ d + d ⋆ d⋆ for 0-,
1-, and 2-forms. Given that these forms can also be performed on the primal or the
dual, we now have six discrete versions of the Laplace operator:
̃ 1D
̃ 1 H1 D0
̃ 1D
̃0
H1 D1 H
̃ 1D
̃ 0 H2 D1 + D0 H
̃ 2D
̃ 1 H1
⎪H
P(⋆d ⋆ d + d ⋆ d⋆)R = ⎨
̃ 2D
̃1 + D
̃ 0 H2 D1 H
̃1
H1 D0 H
̃ 1D
̃ 0 H2
D1 H
̃1 1 0 ̃ 2
⎩D H D H

for primal 0-forms
for dual 0-forms
for primal 1-forms
for dual 1-forms
for primal 2-forms
for dual 2-forms.

The convergence plots, shown in Fig 3.8, demonstrate spectral convergence again.
Finally, to distinguish our approach from other conventional discretization techniques and reinforce the point that our spectral method is not solely limited to scalar
functions, we depict a solution of the Poisson equation for primal 1-forms (i.e., the
equivalent of a vectorial Poisson equation) on Fig 3.9.
3.7

Conclusions and Future Work

In this chapter, we laid out the foundations for a spectral instance of discrete exterior
calculus. Using chains and cochains as discretization for differential forms, we
provide a discrete exterior derivative, a discrete Hodge star, and a discrete wedge
product that are spectrally accurate. One can use the resulting calculus as a chain
collocation method to solve differential equations.
While the three exterior operators we treated cover already a large range of differential operators (allowing a number of computational tasks including Hodge
decomposition), contractions (exterior product) and Lie derivatives still need to be
derived in a similar geometric, coordinate-free manner—currently, they can only

39
Primal
Dual

1e0

1e-4

∆f−q ∞

1e-4

Primal
Dual

1e0

1e-8

1e-12

1e-8

1e-12

100

(a) Periodic 0-forms

Primal
Dual

1e0

1e-4

∆f−q ∞

1e-4

∆f−q ∞

100

(b) Periodic 1-forms
Primal
Dual

1e0

1e-8

1e-12

1e-8

1e-12

100

(d) Chebyshev 1-forms

Figure 3.8: Convergence graphs for a 2D Poisson equation; (a) we solve ∆ f =

esin x (cos2 x − sin x) + esin y (cos2 y − sin y) on a periodic domain for either a primal, or
dual 0-form f ; (b) now for ∆ f = esin x (cos2 x − sin x)dx + esin y (cos2 y − sin y)dy; (c) we
solve ∆ f = e x + ey on a Chebyshev grid, for either a primal 0-form with Dirichlet boundary
conditions f (x, y) = e x + ey , or for a dual 0-form with Neumann boundary conditions
∇ f (x) ⋅ n = (e x ey )t n; (d) now for ∆ f = e x dx + ey dy. All of our results exhibit spectral
convergence (measured through the L ∞ error ∥∆ f − q∥∞ ), with the conventional plateau in
accuracy for fine meshes.

be handled using coordinates. Note that defining the contraction operator would
be enough, since one can then define the Lie derivative algebraically with Cartan’s
“magic formula" [22]. But a direct dynamic (“flowed-out") definition as presented
in [48] may instead prove beneficial.

Figure 3.9: Plot for the solution of (⋆d⋆d+d⋆d⋆) f = 0 with the boundary conditions of
⋆ f ∣B = 0 (vector field is tangent to the boundary) and ⋆d⋆ f ∣B = 1 (the curl at the boundary
is equal to 1). The domain is [-1, 1]2 , discretized by a 10×10 Chebyshev grid.

41
Chapter 4

SPECTRAL EXTERIOR CALCULUS

4.1

Introduction

Continuum mechanics formulates the governing equations of deformable bodies and
fluids using operators from integral and differential calculus on Riemannian manifolds. While the mathematics of Riemannian geometry has been long established,
computational methods require the discretization of domains, fields, and operators
since numerical simulation necessitates finite-dimensional approximations—and
many questions remain open about what are the best approaches to discretize the
aforementioned continuum theories. Discretization of the domain geometry is often
done through the design of a simplicial or cell complex, while fields and operators are approximated using finite element basis functions and weak formulations.
The locality of the basis functions used in computations directly influences both
the computational complexity (i.e., simulation times and memory usage) and the
accuracy (i.e., order of approximation) of the resulting computations.
Spectral methods have gained popularity in the computation of continuum mechanics
problems whose solutions are sufficiently smooth: their use leads to exponential rates
of convergence (referred to as spectral convergence [17, 69]). Among the various
spectral methods, spectral collocation that uses global smooth trial functions and
Dirac test functions at collocation points is particularly attractive due to its simplicity
and its computational efficiency.
Recent years have also seen the development of novel discretizations aimed at
preserving in the discrete realm the underlying geometric, topological, and algebraic
structures at stake in differential equations. This has proven to be a fruitful guiding
principle for discretization [7, 2, 21, 3]. This geometric approach has led to new
structure-preserving numerical methods, analyzed in, e.g., [3, 36], that inherit a
variety of properties from the continuous world.
Yet, combining structure preservation with spectral convergence is far from a simple
matter. Only a handful of works have been proposed [28, 56, 57, 62, 29], in which the
foundations of exterior calculus [1] are adapted to spectral methods to offer spectrally
accurate treatments of various operators. To this day, there is no unifying framework

42
covering all the exterior calculus operators needed for computational modeling that
offers an efficient implementation of generic computations with structure-preserving
and spectral properties.
In this paper, we present Spex, an exterior calculus based spectral approach to computations on simply-connected domains of arbitrary dimension, equipped with an
arbitrary Riemannian metric. We leverage two complementary discrete representations of differential forms to offer an efficient evaluation of both metric-dependent
and metric-independent operators: namely, we represent k-forms via their components at every mesh node and via their integral over mesh k-dimensional elements.
The resulting spectral calculus extends previous approaches [57, 62] and provides
efficient evaluations of all the exterior operators (d, ⋆, ∧, ⌟, L) and the musical isomorphisms ♭ and ♯ (to convert vector fields to one-forms and vice-versa), while
preserving many of the topological properties and features of the continuum theory.
4.2

Context and Motivations

Computational methods that inherently preserve geometric structures have become
increasingly popular over the past few years. They have gained acceptance among
engineers and mathematicians [4] by deriving numerical evaluations based on wellestablished geometric foundations of calculus on smooth manifolds as we briefly
review next.
Exterior Calculus Based Computations
Exterior calculus [1] is a concise mathematical formalism that expresses differential
and integral equations on smooth and curved spaces while revealing important
geometric and topological structures behind integration and derivation. At the
heart of exterior calculus is the notion of differential forms (pioneered by Cartan),
denoting antisymmetric tensors of arbitrary rank. Formally, a k-form on a manifold
M is a multi-linear map that takes, at every point, k tangent vectors as input and
returns a real number; in this paper, we will denote the space of k forms as Λ k
and the space of vector fields as X. Note that consequently, Λ0 is simply the space
of scalar functions on M, and that Λ1 and X are dual in the algebraic sense, i.e.,
one-forms can be understood as covectors since a one-form applied to a vector
field returns a scalar function. As integration of differential forms is an abstraction
of the process of measurement, this form-based calculus provides an intrinsic,
coordinate-free approach to describing a multitude of physical models that make
heavy use of line, surface and volume integrals [15, 25, 26]. Moreover, physical

43
measurements (such as fluxes or currents) typically represent local integrations over
a small line, surface or volume element of the measuring instrument: pointwise
evaluation of such quantities does not have physical meaning. Thus, one should
manipulate these quantities as geometrically-meaningful entities integrated over
appropriate submanifolds. With differential forms as its building blocks, exterior
calculus consists of seven operators as listed in Table 4.1, and compositions thereof
such as the inner product of forms and the Laplace-Beltrami operator. (Note that
the contraction of a one-form α by a vector field X, X ⌟ α, is sometimes written as
⟨α, X⟩, or as i X α, or as α(X) depending on the authors.)
Name

Symbol

Derivative

Λk → Λk+1

Wedge product

Λk × Λl → Λk+l

Contraction

Domain → Range Metric

X×Λ

→ Λ

k−1

→ Λ

n−k

yes

Hodge star

Flat

X → Λ1

yes

Sharp

Inner product

_ ⋅ _ ≡ ⋆(_ ∧ ⋆_)

Λ k × Λ k → Λ0

yes

→ Λ

yes

X×Λ

→ Λ

Laplace-Beltrami ∆ ≡ d ⋆ d ⋆ + ⋆ d ⋆ d
Lie derivative

L = _ ⌟ d_ + d_ ⌟ _

→ X

Table 4.1: List of basic and a few compound operators of exterior calculus, with their
dependence on the metric being indicated in the last column.

(Note that the contraction of a one-form α by a vector field X, X ⌟ α, is sometimes
written as ⟨α, X⟩, or as i X α, or as α(X) depending on the authors.) This exterior
calculus notation is very convenient and universal, as it allows us to write a number
of key theorems and physical models on arbitrary (flat or curved) domains in a very
compact manner that is independent of any particular choice of coordinate system.
The Laplace operator ∆ for instance encapsulates both the scalar Laplacian (used in
diffusion) when applied to zero-forms, and the vector Laplacian (used in the wave
equation for the electric field in E&M) when applied to one-forms. Stokes’ theorem
(which supersedes Green’s, Kevin’s, and Gauss’s theorems) states that the integral
of the derivative d of an arbitrary form α over a domain D is equal to the integral
of the said form over the boundary ∂D of D, i.e., D dα = ∂D α; that is, d and the
boundary operator ∂ are adjoint with respect to integration. In fluid dynamics, the
incompressible Navier-Stokes equations on arbitrary domains in the case of a fluid

44
density being constant in time and uniform in space are simply
⎪vt + Lv♯ v = −dp + Re ∆v,
⎩d ⋆ v = 0,
where v is the velocity one-form, p is the pressure zero-form, and Re is the dimensionless Reynolds number measuring the importance of viscous forces relative to
inertial forces. Maxwell’s equations of electromagnetism can also be formulated
very compactly using the operators of exterior calculus:
dF = 0,
d ⋆ F = J,
where F is the Faraday 2-form (incorporating the electric and magnetic fields)
and J is the electric current 3-form (incorporating the current and electric charge
densities) [64]. Typical boundary conditions (Dirichlet, Neumann) are also trivially
expressed with exterior operators and inclusions maps.
Finite Dimensional Exterior Calculus
Discrete exterior calculus endeavors to extend exterior calculus to discrete spaces,
be they graphs or manifold meshes. Behind the various efforts in this direction [8,
22, 2, 31, 60, 30] is the desire to preserve as many properties as possible of the
continuous forms and exterior calculus operators in the finite-dimensional case for
computational purposes. For instance, enforcing Stokes’ theorem at the discrete
level guarantees various conservation properties, while a proper notion of a discrete
Hodge decomposition prevents spurious modes in solutions. Most of these efforts
use the notion of chains and cochains from algebraic topology to discretize manifolds
and differential forms, respectively. Integration is thus seen as the pairing between
chains (linear combination of mesh elements) and cochains, and Stokes’ theorem
is enforced by construction. This approach, creating an entire discrete de Rham
complex, extends traditional finite element methods with now (low- or high-order)
Whitney basis functions associated to each element of a (simplicial or cell) complex
that represents the computational domain: k-forms can thus be histopolated (i.e.,
reconstructed from their integrated values on chains) over the whole domain.
Spectral Chain Based Methods
While finite element methods use locally supported basis functions, spectral methods favor globally supported basis functions instead, combined with the use of the

45
fast Fourier transform for efficiency of computation. They are often a preferred
approach for problems where the solution is smooth because their global nature
offers exponential convergence under refinement. However, structure preservation
was initially not a topic of interest for spectral methods. In fact, the point collocation method, one of the most common spectral approaches, computes derivatives
nodewise in Fourier space, thus failing to satisfy Stokes’ theorem.
A decade ago, Mai-Duy et al. [44] proposed to modify the spectral collocation
method by using integrated Chebyshev polynomials to numerically solve biharmonic
differential equations with easier enforcement of boundary conditions. Shortly after,
Gerritsma [28] described a method for constructing higher order edge basis functions, a staple of discrete counterparts to exterior calculus, to interpolate 1-forms.
From there on, a few approaches were recently formulated that marry structure
preservation properties and spectral accuracy [56, 57, 62, 29]. The ingredients at
the core of these methods are quite similar: all rely on the use of cochains as a finitedimensional representation of forms and they introduce spectral approximations of
Hodge stars. Most related to our contribution, the work of [62] presents a spectrally
accurate calculus of discrete differential forms described as a “chain collocation
method” as it extends the traditional point collocation method to collocations over
any type of mesh elements.
Limitations of Previous Work
While discrete, structure-preserving versions of exterior calculus offering spectral
accuracy are readily available, the current state of the art is limited in various ways.
• While exterior derivative, Hodge star, and wedge product have been formulated in these frameworks, there are currently no associated contraction and
Lie derivative that fit the same computational framework: previous approaches
rely on high-order accurate extrusion of chains [47] or invoke pointwise notions of forms [29]) instead of cochains—both incompatible with existing
spectral cochain-based methods.
• In existing works, vector fields and 1-forms are usually not differentiated: vector fields are simply treated through their corresponding one-forms. Yet these
two entities are clearly different from a geometric perspective (Lie derivatives
of forms vs. vector fields, for instance, are clearly distinct operations), as they
do not even transform similarly under a change of basis.

46
Symbol

Meaning

Dimension of the manifold

Reference domain of dimension d

Manifold of dimension d

Diffeomorphism ϕ ∶ D → M

Σ̄ k
Σk

Primal k-dim cells of M

Σ̊

Component 0-dim cells (vetices) of M

σ̄ik
σik
φ̄ik
φ̃ik

Primal k-dim cell of grid Σ̄ k with index i
Dual k-dim cell of grid ̃
Σ k with index i

Space of differential k-forms ωk on M

Λ̄k
̃k

Space of primal discrete k-forms ω̄k on C̄M
̃k on CM
Space of dual discrete k-forms ω

Λ̊k

̃k on C̄ 0 M ∪ C̃0 M
Space of component k-forms ω

Space of vector fields on M (dual of Λ1 )

X̊

Space of discrete component vector fields on M

Poincaré duality operator (∗σ̄ik = ̃
σid−k )

Boundary operator on grid elements

Dual k-dim cells of M

Primal basis function for k-forms associated with σ̄ik
Dual basis function for k-forms associated with ̃
σik

Table 4.2: Meaning of the basic notations used throughout this document.
• The notion of discrete metric is also rarely discussed in existing works. In
particular, the derivation of Hodge stars from arbitrary metrics is not available,
and neither are the musical isomorphisms ♭ and ♯ which lower and raise indices
for one-forms and vector fields.
Overcoming these limitations will require enriching the cochain-based representation that is the basis of all spectral calculus approaches thus far to offer a fast, and
spectrally accurate formulation of all the exterior calculus operators.
Spex Overview
We now formulate the basic design principles that we will use to derive our spectral
exterior calculus.
Pointwise vs. Integral Representations The operators derived in the context of
existing spectral calculus frameworks have been mostly unary. However, binary

47
operators such as the wedge product ∧ and the contraction ⌟, which have two inputs,
require combining components through multiplication. Therefore any discrete form
that is not associated with a vertex (i.e., with a degree higher than zero) must be
converted to a component representation at vertices that are collocated to allow for
direct pointwise multiplication. Additionally, since multiplication can double the
spectral band, the multiplication has to be performed on a refined grid with twice
as many vertices along each dimension in order to properly resolve the result. This
necessitates the introduction of a component-based representation of differential
forms, which we will call Component Forms, for which the space of all component
forms will be denoted Λ̊. Additionally, a space of Component Vector Fields (denoted
X̊) will also be similarly defined, so that the operators ♭ and ♯ providing the mapping
to and from Component One Forms Λ̊1 can be easily implemented. Importantly,
we will show that being able to easily convert between these twin (chain-based and
component-based) representations will be crucial to the efficient implementation of
our operators. The various maps between the spaces of continuous and discrete
forms are shown in Figure 4.10.
Chain Based Exterior Derivative Discrete and finite element exterior calculus
strive to mimetically reproduce exterior calculus on finite degrees of freedom. With
this goal in mind, the use of chains and cochains was shown to be the most natural
approach for constructing a discrete counterpart of the continuous d operator: its
topological nature and adjointness to the boundary operators make it best described
in the discrete sense purely based the incidence of elements of the computational
grid, as it enforces a simple discrete version of Stokes’ theorem for all chains of the
grid [35]. Spex will proceed similarly, with the discrete exterior derivative acting
naturally on Λ̄ and Λ.
Hodge Operator and Poincaré Duality If M is a d-dimensional manifold, then
there is a natural isomorphism between k-forms Λ k and (d −k)-forms Λd−k , called
Poincaré duality. The duality map between the two spaces is through the Hodge
star operators ⋆ k and ⋆d−k , acting on k- and (d−k)-forms respectively:
⋆k ∶

Λ k → Λd−k ,

⋆d−k ∶

Λd−k → Λ k .

In the discrete settings, this duality calls for the construction of a dual grid: the
discrete star operators will then provide an isomorphism between primal (Λ̄) and
̃ since the dimension of the space of primal (resp., dual) k-forms
dual cochains (Λ)
matches the dimension of the space of dual (resp., primal) (d−k)-form [35].

48
High-order Operators Previously, the Hodge star has been evaluated using zeroth
or first order basis functions, which while sparse and efficient, have poor convergence properties that do not permit its use in expressions involving combinations
of multiple operators such as the composite Poisson operator ⋆d ⋆ d + d ⋆ d⋆. For
lower order ⋆, the result of ⋆d ⋆ f is not meaningful since the error of ⋆ f ∼ O(h),
where h is the resolution of the discrete mesh, would become d ⋆ f ∼ O(1) after
differentiation. Solving this issue can be achieved through a weak formulation;
but with higher order methods, where ⋆ f ∼ O(h k ), taking the derivative returns
directly a meaningful result d ⋆ f ∼ O(h k−1 ) for sufficiently high grid resolution,
and therefore the Poisson operator, or even higher-order operators, can be computed
directly.
Additionally, since multiplication doubles the spectral band, in order to properly
resolve the result, the multiplication has to be performed on a refined grid with
twice as many points along each dimension.
Orthonormal Frames The metric will be stored in this representation per node
and will be used to locally construct an orthonormal frame, in which the ⋆, ♭ and ♯
operators will have will be easy to compute.
In Riemannian geometry and manifold theory, orthonormal frames are useful for
studying the structure of a differentiable manifold equipped with a metric g. If M
is a manifold and g is its metric, then an orthonormal frame at a point m of M
is an ordered basis of elements of the tangent bundle Tm M consisting of vectors
which are orthonormal with respect to the bilinear form gm . The metric provides
a mapping to a basis set of the dual tangent space Tm∗ M given by the flat operator.
To make that clear, consider a vector v ∈ Tm M; then its corresponding one-form is
v♭ (⋅) = gm (v, ⋅).
The components of forms or vectors are always expressed relative to a basis frame.
A vector frame B is an ordered set of vectors {ei } for each point of the manifold,
whereas a co-vector frame B ∗ is an ordered set of co-vectors {ei }. Formally ei is an
assignment of a vector to each point, i.e. ei ∶ M → TM, whereas ei is the same on
the cotangent bundle ei ∶ M → T ∗ M.
The reference coordinate basis for a rectilinear d-dimensional domain is
B={

, . . . d },
∂x ∂x
∂x

B ∗ = {dx 1, dx 2, . . . dx d },

49
where the coordinates {x 1, x 2, . . . , x d } span the manifold M. This basis set is
holonomic because for any differentiable manifold the set of basis vectors {ei } are
integrable and they are such that there exists a coordinate system {x i } for which
ei =

∂ xi

ei = dx i,

(4.1)

The tangent and cotangent basis sets satisfy the contraction rule
ei ⌟ e j = δij .
The contraction pairs forms with vectors. If we want to pair vectors with vectors or
forms with forms, we need the metric. The coefficients of the metric are expressible
in terms of the basis elements as follows:
gi j = ei ⋅ e j ,

gi j = ei ⋅ e j ,

and the metric is related to its inverse by
gik g k j = δij ,

gik gk j = δij .

Unless the metric is flat everywhere, this basis set will not be orthonormal in the
general case, i.e. ei will not be perpendicular to e j .
The metric establishes a duality between ei and ei by lowering and raising their
indices:
ei = gi j e j ,
ei = gi j e j .
The reason orthonormal frames are useful is because it is easy to express the Hodge
star and the inner product relative to an orthonormal frame of reference. In such a
frame, the Hodge star operator acts on a basis form of degree r as follows:
⋆ (ei1 ∧ ei2 ∧ ⋅ ⋅ ⋅ ∧ eir ) = i1i2 ...ir jr+1 ... jm e jr+1 ∧ e jr+2 ∧ ⋅ ⋅ ⋅ ∧ e jd ,

(4.2)

and the inner product of basis forms of degree r is given by
(ei1 ∧ ei2 ∧ ⋅ ⋅ ⋅ ∧ eir ) ⋅ (e j1 ∧ e j2 ∧ ⋅ ⋅ ⋅ ∧ e jr ) = δi1 j1 δi2 j2 . . . δir jr

(4.3)

assuming the basis indices are ordered i1 < i2 < ⋯ < ir and j1 < j2 < ⋯ < jr . For
example, in two dimensions, for a flat metric, the Hodge star and the inner product
are completely defined by the equalities
⋆ 1 = e1 ∧ e2,

⋆e1 = e2,

⋆e2 = −e1,

⋆(e1 ∧ e2 ) = 1,

(4.4)

50
1 ⋅ 1 = 1,

e1 ⋅ e1 = 1,

e1 ⋅ e2 = 0,

e2 ⋅ e2 = 1,

(e1 ∧ e2 ) ⋅ (e1 ∧ e2 ) = 1. (4.5)

In order to get to an orthonormal frame, vectors and forms must undergo a basis
transformation. To avoid ambiguity, let us denote the orthonormal frame with a hat:
B̂ = {ê1, ê2, . . . , êd },

B̂ ∗ = {ê1, ê2, . . . , êd },

where êi ⊥ ê j for any pair of basis elements with respect to the relevant metric. The
orthonormal frame is defined on each tangent space of the manifold and it may be
non-holonomic as there does not necessarily exist an integrable global coordinate
system that is linked to the basis vectors by Eq. (4.1).
The transformation between the two basis sets êi and ei will be given by
êi = J ij e j

êi = e j (J −1 )ij .

Notice that vectors and forms transform in such a way that their contraction remains
invariant:
êi ⌟ ê j = ei ⌟ e j = δi j .
Namely, if vectors transform as v ↦ Jv, then one-forms must transform as α ↦ αJ−1 ,
where the contraction corresponds to the matrix product of a row with a column
vector v ⌟ α = αv. Section 4.A provides a detailed overview of how different objects
transform under a coordinate change. The results relevant to the two-dimensional
case are summarized in Table 4.3. In matrix notation, vectors are treated as a column
matrix with one single column and differential forms are treated as a row matrix
with a single row.
Name

Element Transformation

vector

X ∈X

X ↦ JX

f ↦ f

α ↦ αJ−1

0-form

f ∈Λ

1-form

α∈Λ

2-form

ω ∈ Λ2

ω ↦ ω det(J)−1

Table 4.3: List of transformation for vectors and forms under a basis change in 2D.

51
Induced Metric The reference computational grid will be denoted by D, which
for the periodic case will be [0, 2π)2 and for the bounded cases will be [−1, 1]2 .
The reference domain D will be mapped to the physical domain M by the map
ϕ ∶ D → M.

Typically, the manifold M will be embedded in either two- or three-dimensional
Eucledian space, i.e. M ⊂ R{2,3} , and it will have an associated metric gM that will
be the standard Eucledian metric for the embedding space. The map ϕ induces a
metric gD on the reference domain D given by the pullback of the metric from the
physical manifold M:
gD = ϕ∗ (gM ).
The components are given by its application to the basis elements
gi j = gD (ei, e j )
Given this induced metric, we want to compute a transformation J ij such that êi =
e j (J −1 )ij is an orthonormal frame with respect to gD :
ĝi j = gD (êi, ê j ) = δi j .
Because the expression above is invariant under rotations,
gD (êi, ê j ) = gD (Rêi, Rê j ),
there are many choices for the gauge. The matrix square root is one such choice. Let
us denote the elements of the transformation J ij using the matrix J, and the elements
of the tensor gi j using the matrix g. For example, in two dimensions ,
⎡ 1 1 ⎤
⎢ J J ⎥
J = ⎢⎢ 12 22 ⎥⎥ ,
⎢ J1 J2 ⎥

⎢ g11 g12 ⎥
⎥.
g=⎢
⎢ g21 g22 ⎥

(a) {e1, e2 }

(b) {ϕ∗ e1, ϕ∗ e2 }

(d) {ϕ∗ ê1, ϕ∗ ê2 }

Figure 4.1: Computing an orthonormal frame for a surface embedded in 3D Euclidean
space. The black arrows represent e1 , and the red arrows represent e2 . (a) and (b) illustrate
the coordinate basis B in D and M respectively, and (c) and (d) illustrate the orthonormal
basis B.

For J to map to an orthonormal frame, it must be related to g by
g = JT J.
To see that this does indeed result in an orthonormal frame:
ĝi j = gD (êi, ê j )
= gD (em, en )(J −1 )im (J −1 )nj
= gmn (J −1 )im (J −1 )nj
= (J−T g J−1 )i j
= δi j

53
as required. One possibility for a J that satisfies the above requirement is given by
the matrix square root, which results in a symmetric J since g is symmetric:
J=

(a) {e1, e2 }

(b) {ê1, ê2 }

(d) {ê′′1 , ê′′2 }

Figure 4.2: Computing an orthonormal frame for a surface embedded in 2D Euclidean
space. (a) illustrates the standard reference basis and its push forward, (b) - an orthonormal
basis by taking the square root of the metric, (c) - an orthonormal basis using the GranSchmidt process starting with e1 , and (d) - starting with e2 .

Figure 4.1 illustrates the results of constructing of an orthonormal frame using the
above square root method for a manifold embedded in R3 . Note that the orthonormal
frame is by no means unique and there are infinitely many equally valid choices for
constructing J from g. For example, given any set of linearly independent vectors
{e1, e2, . . . , ed }, instead of taking the square root as above, one could use the GramSchmidt process to compute an orthonormal set {ê1, ê2, . . . , êd }. Figure 4.2 compares
the matrix square root and the Gram-Schmidt orthogonalization.
4.3

Spatial Discretization and Associated Basis Functions

We begin our exposition by discussing how the computational domain D can be
discretized into a cell complex, and how the basis functions (each associated with
a cell of the computational grid) are expressed accordingly. We will consider two
types of domains - periodic and bounded.

54
One-dimensional Cell Complexes
For a one-dimensional domain, its discretization into a computational cell complex
consists of only vertices (also called nodes) and edges. We will use N to refer to
the grid size, defined by the number of primal edges. As noted earlier, our discrete
calculus setup exploits Poincaré duality by using both a primal and dual grid between
which the Hodge stars can act as a duality operator that converts 0-forms to 1-forms
and vice versa. Consequently, a one-dimensional grid has four different cells: primal
̃ and dual edges (E).
vertices (V̄), primal edges (Ē), dual vertices (V),
Periodic Domain When the 1D domain is periodic, we consider D = [0, 2π)
without loss of generality. The positions of primal and dual vertices are set to
V̄N = {θ̄ n =

2π
n,

n = 0, . . . , N − 1},

̃N = {θ̃n = 2π (n + 1 ),

n = 0, . . . , N − 1},

offering a regular sampling of the domain. Primal (resp., dual) edges are between
consecutive primal (resp., dual) vertices. The resulting cell complex will be referred
to as the periodic cell complex. Its cells will be given by
σ̄n0 = [θ̄ n ],

σn0 = [θ̃n ],

σ̄n1 = [θ̄ n, θ̄ n+1 ],

σn1 = [θ̃n−1, θ̃n ].

Bounded Domain Without loss of generality, we consider D = [−1, 1] in the
bounded case. Positions of the primal and dual vertices are set to
n+1
π),

n = 0, . . . , N − 2},

n + /2
̃N = {̃
π),
xn = − cos (

n = 0, . . . , N − 1}.

V̄N = { x̄n = − cos (

Just like in the periodic case, primal (resp., dual) edges are between consecutive
primal (resp., dual) vertices; but in this bounded domain case, two boundary vertices
are also added at each end of the domain D in order to ensure the closure of the
boundary operator ∂. The resulting cell complex will be called the bounded cell
complex. The cells, as in the periodic case, are given by
σ̄n0 = [ x̄n ],

σn0 = [̃
xn ],

σ̄n1 = [ x̄n, x̄n+1 ],

σn1 = [̃
xn−1, ̃
xn ].

Figure 4.3 illustrates the relative topological arrangement of its elements.

V̄ 0
Ē 0
V0

V̄ 1
Ē 1
V1

V̄ 2
Ē 2

Ē 3
V3

Figure 4.3: Schematic representation of the computational grid in the one-dimensional

̃ and Ẽ are used to store the discrete forms in Λ̄0 , Λ̄1 ,
bounded case. The simplices V̄, Ē, V,
Λ , and Λ respectively. Vertices are drawn equally spaced for simplicity.

V̄ 0

V̄ 1
E0

Ē 0
V0

V̄ 2
Ē 1
V1

V̄ 3
Ē 2

V̄ 4
Ē 3
V3

Ef4

Figure 4.4: Cell complex with clipped edges at the boundary (old setup).
Comparison Note that the cell complex shown in Figure 4.3 is different from the
one used in the previous chapter [62]. In the new grid there are no longer any
clipped edges. All the edges that form the cell complex are whole, whereas before
they were clipped at the boundary (see Figure 4.4). While this difference is rather
insignificant from a geometric point of view, the use of half-edges complicates the
analytical expressions for histopolation as well as the implementation of exterior
calculus operators. We thus default to using the simpler arrangement of edges in
Figure 4.3: the new arrangement for the dual edges no longer fills the entire domain,
and only the primal edges fill up the domain completely. Boundary vertices are
considered separately, as they will allow us to prescribe boundary conditions.
1D Grid
Periodic
Bounded
Bounded (old)

∣V̄N ∣
N −1
N +1

∣ĒN ∣

̃N ∣
∣V

∣ẼN ∣
N −1
N +1

∣V̊N ∣
2N
2N − 1
2N + 1

Table 4.4: Number of primal and dual elements in 1D grids as a function of N.
Component vertices For any cell complex, the component vertices (V̊) are then
defined as the union of the primal and dual vertices:
̃.
V̊ ≡ V̄ ∪ V

56
The total number of vertices and edges for the one-dimensional grid for a given
discretization level N is shown in Table 4.4. Notice that our construction implies
V̊N = V̄2N ,
that is, the component vertices are identified to the primal vertices of the twice-finer
grid V̄2N , a property that we will leverage to seamlessly upsample or downsample a
field later on.
One-dimensional Basis Functions
Having defined the cell complexes, we now proceed to describe a general approach
for deriving the basis functions associated with them. We need to define the
following set of basis functions for the cell complexes described in the preceding
section
φ̄0n

φ̃0n

φ̄1n

φ̃1n

The topic of constructing cardinal basis functions associated with the vertices of
a computational grid has been well studied in the literature [11, 42] and can be
used to interpolate vertex-based functions (0-forms). Additionally, we can construct
histopolation basis functions that can be used to do the same for edge-based functions
(1-forms).

position variable
vertex positions
grid function
vertex functions
edge functions

periodic
θn
ΦN (θ)
V[ΦN , n](θ)
E[ΦN , n](θ)

bounded
xn
PN (x)
V[PN , n](x)
E[PN , n](x)

Table 4.5: Correspondence between notation for periodic and bounded domains.
Table 4.5 illustrates the correspondence between the notation for periodic and
bounded domains. The grid function is a function whose zeroes correspond to
the vertices of the grid, hence its name. Its shape defines the grid. For the periodic
case, the grid function will be denoted by ΦN (θ) and in the bounded case it will
be PN (x). From the grid functions we can generate the cardinal vertex and edge
functions, which we will denote using V and E respectively.

57
Vertex Interpolation
Bounded For bounded domains we use polynomial interpolation. The Lagrange polynomials, which correspond to the vertex basis functions, are given by
x − xm
m≠n xn − xm

pn (x) = ∏

The polynomial pn (x) assumes the value one only at the vertex xn , and is zero
at every other vertex xm (m ≠ n). If xn are the roots of the polynomial PN , i.e.
PN (xn ) = 0 for ∀n = 0, . . . , N − 1, then the Lagrange polynomials can be written in
the equivalent form
1 PN (x)
pn (x) =
N x − xn
where N is a normalization factor chosen to enforce the condition
pn (xn ) = 1.
By L’Hospital’s rule, we can compute the limit as
1 ′
P (xn ),
x→xn
N N
where the prime denotes differentiation. Thus the appropriate value for the normalization constant is N = PN′ (xn ) and the formula for the cardinal functions becomes
lim pn (x) =

V[PN , n](x) ≡ pn (x) =

1 PN (x)
PN (xn ) x − xn

We will say that these are the vertex functions associated with the roots of PN , and
we will write them as V[PN , n](x) in order to underline the fact that the vertex
functions can be computed solely from PN and n. Note that V[PN , n] as a function
of PN is homogeneous in PN , i.e. for any nonzero constant a it holds true that
V[aPN , n] = V[PN , n].
Periodic A similar formula holds true for trigonometric interpolation over
periodic domains.
sin 12 (θ − θ m )
φn (θ) = cN (θ) ∏
m≠n sin 2 (θ n − θ m )
where cN is a correction factor for when N is odd/even.1 If θ n are the roots of the
trigonometric function ΦN , i.e. ΦN (θ n ) = 0, then the periodic cardinal functions

For equidistant points c N is given by c N (θ)

1+(−1) N
cos θ2

= cos2 θ4 − (−1) N sin2 θ4 .

⎪1
= ⎨
cos θ2

if N odd
if N even

1−(−1) N

58
can be written as
V[ΦN , n](θ) ≡ φn (θ) =

ΦN (θ)
ΦN (θ n ) 2 sin 12 (θ − θ n )

which emphasizes the fact that the above cardinal functions can be solely derived
from ΦN and n. It can be easily checked by applying L’Hospital’s rule that the above
expression does indeed constitute a cardinal function by verifying that φn (θ m ) = δnm .
Edge Histopolation
Bounded In order to compute the edge functions, which are needed for histopolation, we follow the approach taken by Gerritsma [28]. Namely,

E[PN , n](x) = − ∑ V[PN , m]′ (x).
m=0

The above expression corresponds to taking the derivative of a “step function”. Only
the edge across which the step shifts from 0 to 1 will have a non-zero value for the
integral, whereas for the rest of the edges that value will be equal to zero.
The derivative is explicitly given by
V[PN , n]′ (x) =

(x − xn )PN′ (x) − PN (x)
(x − xn )2 PN′ (xn )

Sometimes we will need to evaluate the above expression at xn , and its limit there is
V[PN , n]′ (xn ) =

PN′′ (xn )
2PN′ (xn )

Explicit expressions for the limits at the boundaries for Chebyshev polynomials are
given in 4.B.
Periodic To go from vertex functions to edge functions, we first note that the
convolution of the edge function with a box function that spans the corresponding
edge must be equal to the vertex function [62]. Given a function f (θ) we define the
following functional transformation on it:
ˆ θ+ h
2π
T [ f ](θ) =
f (θ ′ ) dθ ′ with h =
θ− h

The cardinal edge functions then can be expressed in terms of the inverse of the
above transformation:
E[ΦN , n](θ) = T −1 [V[ΦN , n]](θ).

59
Basis Functions Explicitly
Periodic The primal and dual vertices are the roots of the trigonometric functions below:
N odd

Φ̄N (θ) = sin ( Nθ)

V̄ 0

V̄ 1
V0

̃ N (θ) = cos ( 1 Nθ)

if N odd.

V̄ 2

V̄ 4

V̄ 3
V2

V̄ 0

2π

N even

Φ̄N (θ) = sin ( Nθ) cos

V̄ 0

̃ N (θ) = cos ( 1 Nθ) cos θ

V̄ 1
V0

V̄ 2
V1

if N even.

V̄ 3
V2

V̄ 0
V3

2π

Then the basis functions in the periodic case become
φ̄0n (θ) = V[Φ̄N , n](θ),

̃ N , n](θ),
φ̃0n (θ) = V[Φ

̃ N , n](θ),
φ̄1n (θ) = E[Φ

φ̃1n (θ) = E[Φ̄N , n](θ).

Bounded One distinguishes between Chebyshev polynomials of the first kind
which are denoted by Tn and Chebyshev polynomials of the second kind which
are denoted Un . The primal vertices for the Chebyshev grid are the roots of the
polynomials of second kind, whereas the dual vertices are the roots of the first kind.
Hence, the primal and dual grid functions can be written as

V̄ 0
V0

−1

V̄ 1

V̄ 2

V̄ 3

For the primal edges, we need to clamp the polynomial at the endpoints by multiplying by 1 − x 2 :
clamp[PN ](x) ≡ (1 − x 2 )PN (x).
Then basis functions are
φ̄0n (x) = V[UN−1, n](x),

φ̃0n (x) = V[TN , n](x),

φ̄1n (x) = E[clamp[UN+1 ], n](x),

φ̃1n (x) = E[TN , n](x).

Higher-dimension Grids and Associated Basis Functions
The case of grids in higher dimensions are simply handled via Cartesian product
of the 1D grids we just described: Figure 4.5) shows the case of 2D grids for the
periodic domain [0, 2π)2 and the bounded domain [−1, 1]2 . Note that the component
vertices are also obtained by Cartesian product of the 1D component vertices: they
are no longer just the union of primal and dual vertices, but also contain intersections
of primal and dual elements (see example in Figure 4.9). Finally, the associated
basis functions to each of the cells of these grids in nD are simply Cartesian products
of the one-dimensional ones.

(a) Periodic

(b) Bounded

Figure 4.5: Computational grids of domain D for the periodic and non-periodic cases.
Hollow nodes denote boundary vertices, whereas dashed edges denote boundary edges.

Note that while the use of vertices (V) and edges (E) were particularly convenient
for the 1D case, we will now refer to arbitrary d-dimensional cells as σ̄ for primal

61
ones and σ̃ for dual ones, with an upper index indicating its dimension, and a lower
index indicating its label; i.e., in 1D, we will now use
V̄k ≡ σ̄k0,

Ṽk ≡ σ̃k0,

Ēk ≡ σ̄k1,

Ẽk ≡ σ̃k1,

and higher-dimensional cells will be similarly denoted with higher upper indices.
We will denote by φ̄ip (resp., φ̃ip ) the p-form primal (resp., dual) basis function of
the i-th p-dimensional element σ̄ip (resp., σ̃ip ); by construction, they satisfy:
φ̃ip = δi j .
φ̄i = δi j and
σ̄jp

σ̃jp

We also use a systematic enumeration of cells in our code, an example of which is
presented in Figure 4.6. In the remainder of this paper, we will use Σ0 ≡ {σk0 } k to
refer to the set of all nodes, Σ1 for all the edges, etc.
4.4

Discretization of Fields

Equipped with the periodic and bounded computational grids, we can now discuss
how the main fields (vector fields and k-forms) on a smooth manifold M are stored
on this grid. In particular, we will use the property that integration of differential
forms is preserved by the pullback of a smooth mapping function between two
manifolds to motivate a cochain representation of arbitrary forms. A pointwise
coordinate representation of forms, vector fields, and metric will also be leveraged
to offer flexibility in the type of computations that one can achieve with Spex.
We introduce spaces of discrete forms and we will be denote them by adding a
superscript on top of the continuous symbol. Primal and dual cochain forms will be
̃ respectively.
denoted by Λ̄ and Λ,
The space of pointwise-component forms will be denoted by Λ̊.
Discrete Forms As Cochains
As discussed in Section 4.2, cochains are a very powerful discretization of differential
forms. For an arbitrary p-form ω p on the original manifold M, its integrations over
the submanifolds formed by the grid elements of order p of the domain D via
the map φ provide a simple finite-dimensional approximation of the original form:
these integrated form values of (primal or dual) cells of D can be histopolated using
the (primal or dual) basis functions we described in the section above, in order to
reconstruct a smooth approximation of the pull-back of ω on D. For illustration
purposes, we show the cells associated with our cochain representation of forms of

V̄ 2

V̄ 5

V̄ 8

V̄ 1

V̄ 4

V̄ 7

V̄ 0

V̄ 3

V̄ 6

Ē 15

Ē 19

Ē 23

Ē 2

Ē 5
Ē 14

Ē 1

Ē 8
Ē 18

Ē 4
Ē 13

Ē 0

Ē 22
Ē 7

Ē 17
Ē 3

Ē 12

Ē 11
Ē 10
Ē 21

Ē 6
Ē 16

Ē 9
Ē 20

F̄ 3

F̄ 7

F̄ 11

F̄ 15

F̄ 2

F̄ 6

F̄ 10

F̄ 14

F̄ 1

F̄ 5

F̄ 9

F̄ 13

F̄ 0

F̄ 4

F̄ 8

F̄ 12

Ef2
Ef1
Ef0

Ef15

Ef19

Ef23

Ef14
Ef13
Ef12

Ef5
Ef4
Ef3

Ef18
Ef17
Ef16

Ef8
Ef7
Ef6

Ef22
Ef21
Ef20

Ef11
Ef10
Ef9

V 11

V 15

V 10

V 14

V 13

V 12

Figure 4.6: Enumeration of the cells of a 4 × 4 grid. Boundary vertices and edges are
excluded from the figure.

V̄ 3

V̄ 7

V̄ 11

V̄ 15
F3

F 11

F 15

F 10

F 14

F 13

F 12

V̄ 2

V̄ 6

V̄ 10

V̄ 14

V̄ 1

V̄ 5

V̄ 9

V̄ 13

V̄ 0

V̄ 4

V̄ 8

V̄ 12

Ē 23

E 14

Ē 3
Ē 14

Ē 7
Ē 17

Ē 2
Ē 13

Ē 6

Ē 1
Ē 12

Ē 10
Ē 19

Ē 5
Ē 15

Ē 11
Ē 20

Ē 16

Ef13

Ē 9

Ē 4

Ē 8

F̄ 2

F̄ 5

F̄ 8

F̄ 1

F̄ 4

F̄ 7

F̄ 0

F̄ 3

F̄ 6

Ef12

Ef7

Ef16

E 20

Ef6

Ef1
Ē 21

Ē 0

E 17

Ef2
Ē 22

Ē 18

Ef3

Ef19

Ef22
Ef9

Ef18

Ef23

Ef10

Ef5
Ef15

Ef11

Ef21

Figure 4.7: Enumeration of the cells of a 3 × 3 grid. This corresponds to the old setup
with clipped edges at the boundary (e.g., the half-edges Ẽ12 ) visible in the dual grid.

(a) Cochain
Representation.

(b) Pointwise Component
Representation.

Figure 4.8: Two discrete representations of forms.
various orders for a 2D grid on the first two rows of Figure 4.9. We now make these
concepts of discretization and histopolation more formal by introducing cochain
spaces, along with a reduction operator and a reconstruction operator.
Spaces of Primal and Dual Cochains In exterior calculus, differential forms of
order p form a vector space Λ p . As exploited in various efforts to construct a
discrete equivalent of exterior calculus [35, 7, 3], the algebraic topology notion
of cochain is a natural discretization of differential forms of degree p: just as
differential forms are objects to be integrated over p-dimensional domains to return
a scalar, cochains are fundamental objects that assign a number to each cell of
dimension p of the grid. Formally, the space of cochains is the dual linear space to
the space of chains, i.e., linear mappings from chains to real numbers. For instance,
a primal p-cochain is a linear function defined on the primal p-chain space, which
consists of linear combinations of p-dimensional primal mesh elements. As a linear
function, a primal p-cochain can be determined by the value that it assigns to each
p-dimensional primal mesh element.
In this work, we also use cochains as one discrete representation of differential
forms. We will denote the space of primal p-cochains as Λ̄ p , and the space of dual
̃ p . Note that we may sometimes omit the superscript (p) for brevity,
p-cochains as Λ
when the order of the cochain is obvious from context.

(a) B̄0

(b) Σ̄0

(d) B̄1

(e) Σ̄1

(f) ̃
Σ1

(g) Σ̊

(h) Σ̄2

(i) ̃
Σ0

Figure 4.9: Boundary elements are shown in the first column. Boundary vertices (a) are
marked as hollow points, and boundary edges (d) as dashed line segments. The component
vertices (g) Σ̊ are the points at which component forms are collocated regardless of their
degree. The primal cell complex Σ̄ k of degree k is shown in the middle column. The dual
cell complex ̃
Σ k is shown in the last column. There is a correspondence between Σ̄ k and
Σ n−k given by the Poincaré duality operator ∗.

Reduction Operators We denote the linear map from a continuous form to its
corresponding, finite-dimensional primal (resp., dual) cochain representation as the
̃ This reduction is achieved through integration over
reduction operator P (resp., P).
each of the grid elements σ k (or equivalently, on their mapped submanifolds φ(σ k )
onto M) of the same order k as the differential form; i.e.,
P̄ ∶

Λ k → Λ̄ k
ω → ω̄ ,

with

ω̄ik =

σik

ω k,

66
and its dual version
̃∶

̃k
Λk → Λ
̃,
ω →ω

with

̃ik =

σik

ωk .

This reduction operator commutes with the exterior derivative (defined as the adjoint
of the boundary operator on primal or dual chains) and creates an isomorphism
between de Rham cohomology and singular cohomology, two particularly useful
properties for computations [7, 21, 13]. Note that an obvious alternative is through
L 2 projection onto the basis functions of forms.
Reconstruction Operators In order to go from primal or dual cochains to continuous forms, we simply apply interpolation/histopolation based on the basis functions
defined in Section 4.3: with the value of each cell as the coefficient for the basis function associated to the cell, the resulting linear combination reconstructs a continuous
form, i.e.,
R̄ ∶

Λ̄ k → Λ k
ω̄ k → ω k ,

with ω k ≡ ∑ ω̄ik φ̄ik

as well as
̃∶

Λ̃ k → Λ k
ω̃ k → ω k ,

with ω k ≡ ∑ ω̃ik φ̃ik .

For simplicity, we will omit the bar or the tilde when it can be trivially determined
by the operand that the operator acts upon.
Mimetic Projection Maps The primal (resp., dual) mimetic projection maps πh̄
(resp., π̃h ) combine reduction and reconstruction [54] through
πh̄ ∶= R̄P̄,

̃P.
π̃h ∶= R

It can be readily verified that the above maps are projections since they are idempotent:
πh πh = πh .

67
Any differential form originally inside the subspace spanned by the form basis
functions can be discretized and reconstructed without any loss. Moreover, as
the resolution of the meshes increases, the form subspaces become dense, and the
projection onto these finite-dimensional subspaces converges to the identity in the
continuous limit, i.e.,
lim ∥ω − πh ω∥ = 0.
h→0

Discrete Forms as Pointwise Components
While cochains are a convenient representation of integrals of forms (which will
allow us, for instance, to enforce Stokes’ theorem on the computational grid), we
will also make use of a pointwise, component-based representation of forms. Our
“pointwise component” representation of forms is defined as a set of discrete zero
forms stored on the component (i.e., refined) grid Σ̊ (see Figure 4.9). Each zero
form represents one of the coordinates of the form on a chosen basis; for instance,
a 1-form on a 2-manifold parmaterized through two coordinates u and v can be
represented as f (u, v)du + g(u, v)dv, so its point component representation will be
given as the values of f and g sampled at the vertices of the component grid. More
formally, the space of discrete “pointwise component” k-forms will be denoted by
d!
Λ̊ k , with a number of components in d dimensions given by (dk ) = k!(d−k)!
Λ̊Nk ≅ Λ̄02N × ⋯ × Λ̄02N .
´¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹¸¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¹ ¶
(dk)

For example, on a two-dimensional manifold, a 0-form has a single component, a
1-form has two components (e.g., du and dv), and a 2-form has a single component
(e.g., du ∧ dv).
Conversion Between Discrete Forms
To enforce consistency between our twin representations of forms in our finitedimensional setting, we need to be able to easily convert from one representation to
another. For this purpose, we define two operators: a (primal or dual) upsampling
operator, and a (primal or dual) downsampling operator.
From Cochains to Pointwise Components The upsampling operator P̄↑ (resp.,
̃↑ ) converts primal (resp., dual) k-cochains to (d ) components on the component
grid Σ̊:
̃ k → Λ̊ k .
̃↑ ∶ Λ
P̄↑ ∶ Λ̄ k → Λ̊ k

68
Notice that their implentations are simple: one only has to apply the reconstruction
operator RN , followed by a pointwise sampling of the resulting k-form over the
vertices of the component grid.
From Pointwise Components to Cochains Conversely, the downsampling oper̃↓ ) converts a k-form given via its components on the grid Σ̊ into
ator P̄↓ (resp., P
primal (resp., dual) k-cochains on Σ̄ (resp., ̃
Σ):
P̄↓ ∶

̃↓ ∶

Λ̊ k → Λ̄ k

̃k
Λ̊ k → Λ

Note that their implementation is rather simple as well: from the reconstruction of
each component using R2N , we simply have to apply the projection operator Pn to
integrate the resulting continuous form onto each of the k-dimensional elements of
Σ̄N or ̃
ΣN .
Mimetic Conversion Finally, we point out that our downsampling is the left
inverse of upsampling, i.e., one has
P̄↓ P̄↑ = id

and

̃↓ P
̃↑ = id.

Thus, as long as we apply first downsampling then upsampling, any differential form
originally inside the subspace spanned by the form basis functions of a grid stays
unaltered.
Choice of Pointwise Basis for Components Finally, note that a local form basis
needs to be defined for the components to be meaningful. To denote in which basis
the components are expressed, we will add a subscript to the pointwise component
k will denote k-forms with components given in a basis B
spaces. For example, Λ̊B
defined at each vertex of Σ̊. Furthermore, a given basis B can be transformed into
another basis B̂ by a non-singular basis change JBBˆ , thus offering a map between
Λ̊Bˆ and Λ̊B . This will be particularly convenient later on as we will see that
choosing an orthonormal basis with respect to the local metric simplifies a number
of metric-dependent operators like the Hodge stars.
Figure 4.10 summarizes all the spaces of forms, both continuous and discrete, and
the maps between them. Dashed lines signify maps that are surjections (meaning information may be lost), whereas solid lines indicate injections (where no information
is lost).

69
4 Λ̄ j

P¯

P̄↓

P̄↑

P↑

*̃t

JBB̂

4 Λ̊B k

Λ̊Bˆ

JB̂B
P↓

Figure 4.10: Schematic description of the spaces of forms and their various operators.
Discrete Metric
The metric from the original manifold M can be pulled back onto the computational
grid D, and sampled at the component vertices Σ̊. For a given basis B, the resulting
discrete metric can be represented as three components per node, since the metric
is always symmetric. Note that a continuous metric field over D can be easily
reconstructed through R2N .
Discrete Vector Fields
Finally, discrete vector fields are also represented as d components on Σ̊ for a given
vector basis per node. The resulting space of discrete vector fields will be denoted
by X̊.
4.5

Operators on Discrete Forms and Vector Fields

On Cochains
Since the exterior derivative is a fully topological operator, its most natural implementation is done directly using the cochain representation so as to enforce, by
definition, Stokes’ theorem,

dω =

∂σ

on each mesh element σ. Consequently, the discrete exterior derivative is the signed
incidence matrix between (k + 1) elements and k elements of the (primal or dual)
grid, whose non zero elements are only +1 and −1 values to indicate incidence
between mesh cells with the sign determined by the relative orientation of the
elements:
̃T ,
̃ =∂
D̄ = ∂¯ ,
̃) refers to the primal (resp., dual) boundary operator acting
where ∂¯ (resp., ∂
on primal (resp., dual) chains (see [49, 21]). Note that these operators are thus
̃ 2 = 0 like
implemented via sparse matrix multiplication. They also satisfy D̄2 = D

70
their continuous equivalent (since the boundary of a boundary is always the empty
set). These discrete realizations of the exterior derivative are exact, in the sense that
the operators commute with their reduction operators:
D̄P̄ = P̄d

̃ = Pd.
̃P

and

(4.6)

Therefore, these operators need no special treatment to ensure spectral accuracy.
Note that their use for spectral computations represents a clear departure from the
conventional spectral methods [17, 69], which evaluate derivatives pointwise and in
Fourier space. Based on the enumerations in Figure 4.11, the derivative matrices

Ē 3
V̄ 0

Ē 0

V̄ 3

V̄ 0

V̄ 2

Ē 8
Ē 1

Ē 11
Ē 3

Ē 7
Ē 0

Ē 2
Ē 9

will be

D̄ =

⎡1⎤
⎢ ⎥
⎢ ⎥
⎢−1⎥
⎢ ⎥
⎢ ⎥
⎢1⎥
⎢ ⎥
⎢ ⎥
⎢−1⎥
⎣ ⎦

F̄ 2

D̄ =

F̄ 2

F̄ 5

F̄ 8

F̄ 1

F̄ 4

F̄ 7

F̄ 0

F̄ 3

F̄ 6

Ē 4

Figure 4.11: 2 × 2 and 3 × 3 grids

F̄ 0

Ē 5
Ē 10

Ē 6

2 × 2 Grid

F̄ 3

Ē 1
Ē 2

V̄ 1

F̄ 1

⎡−1 0 1 0 ⎤
⎢1 0 0 1⎥
⎢ 0 −1 −1 0 ⎥
⎢ 0 1 0 −1⎥

71
3 × 3 Grid

D̄0 =

⎢1 0 0 0⎥
⎢0 1 0 0⎥
⎢−1 0 1 0 ⎥
⎢ 0 −1 0 1 ⎥
⎢ 0 0 −1 0 ⎥
⎢ 0 0 0 −1⎥
⎢1 0 0 0⎥
⎢−1 1 0 0 ⎥
⎢ 0 −1 0 0 ⎥
⎢0 0 1 0⎥
⎢ 0 0 −1 1 ⎥
⎢ 0 0 0 −1⎥

D̄1 =

⎡ −1 0 0 0 0 0 1 0 0 0 0 0 ⎤
⎢ 1 −1 0 0 0 0 0 1 0 0 0 0 ⎥
⎢ 0 1 0 0 0 0 0 0 1 0 0 0 ⎥
⎢ 0 0 −1 0 0 0 −1 0 0 1 0 0 ⎥
⎢ 0 0 1 −1 0 0 0 −1 0 0 1 0 ⎥
⎢ 0 0 0 1 0 0 0 0 −1 0 0 1 ⎥
⎢ 0 0 0 0 −1 0 0 0 0 −1 0 0 ⎥
⎢ 0 0 0 0 1 −1 0 0 0 0 −1 0 ⎥
⎢ 0 0 0 0 0 1 0 0 0 0 0 −1 ⎥

and the corresponding matrices on the dual grid will be
̃ 0 = D̄1 T

̃ 1 = −D̄0 T .

Boundary Conditions
With our choice of primal and dual Chebyshev grids as defined in Section 4.3, the
Poincaré duality operator ∗ is closed under the primal and dual cells
∗∶

Σ̄ k → ̃
Σ d−k ,
Σ k → Σ̄ d−k .

However, the boundary operator ∂ is not closed under the Σs. We need to introduce
boundary cells to complete its range. Note that for a boundary cell σ ∈ B k we have
no notion of a dual ∗σ.
∂ ∶ Σ k+1 → Σ k ⊕ B k ,
Because D is the adjoint of ∂, this leads to a boundary term in its definition:
D∶

Λ k ⊕ B k → Λ k+1 .

Thus, if λ ∈ Λ k and β ∈ B k , then because D is a linear operator, it can be decomposed
D(λ + β) = D(λ) + D(β) into its restrictions to the internal Λ̄ k and boundary B k
cochain subspaces. In our present setup for the bounded domain the boundary cells
are associated with the primal complex, and therefore only the primal D will have a
boundary term.
Let us denote the restriction of D to Λ̄ k with the symbol D̄, and to B k with the
symbol D:
D̄ ∶ Λ̄ k → Λ̄ k+1,
D ∶ B k → Λ̄ k+1 .

72
Given a k-cochain λ in Λ̄ k , and given a boundary condition β in B k , we will say
that D̄ β λ is the discrete derivative of λ with boundary condition β:
D̄ β λ = D̄λ + Dβ.
̃ will have no associated boundary conditions because there are no
The dual D
boundary vertices or edges associated with the dual grid.
On Components
All the operators that involve multiplication need to be performed in the pointwise
representation. The pointwise representation is useful for implementing the discrete
Hodge star, inner product, sharp, and flat (⋆ ⋅ ♯ ♭), which depend on the metric, and
the contraction and wedge product (∧⌟) which do not depend on the metric. The
discrete metric is also stored in this representation as a matrix per point. There is
no distinction between primal and dual forms in this representation.
Let us denote the components of the discrete component forms (Λ̊) using square
brackets [_]. To distinguish vectors from forms, we denote the vector components
(X̊) using n_o. The components will simply be numbers expressing the coordinates
of the discrete vector or forms relative to a predefined frame of reference. When the
frame is orthonormal, the metric dependent operators are very easy to express. For
the 2D case, which we implement in this work, the operators are given by:
Hodge Star
H̊0 ([a]) = [a]
⎡ ⎤
⎢ ⎥
1 ⎛⎢a x ⎥⎞ ⎢−a y ⎥
H̊ ⎢ ⎥ = ⎢
⎝⎢a y ⎥⎠ ⎢ a x ⎥⎥
⎣ ⎦
H̊ ([a x y ]) = [a xy ]

(4.7)
(4.8)
(4.9)

Inner Product
I̊00 ([a], [b]) = [a ⊙ b]
⎡ ⎤ ⎡ ⎤
⎛⎢a x ⎥ ⎢b x ⎥⎞
I̊11 ⎢⎢ ⎥⎥ , ⎢⎢ ⎥⎥ = [a x ⊙ b x + a y ⊙ b y ]
⎝⎢a y ⎥ ⎢b y ⎥⎠
⎣ ⎦ ⎣ ⎦
22
I̊ ([a xy ], [b xy ]) = [a x y ⊙ b xy ]

(4.10)
(4.11)
(4.12)

The general expressions for dimensions higher than two can be found in Eq. (4.2)
and Eq. (4.3).

73
Flat & Sharp The expressions for the discrete flat and sharp acting on component
forms in the orthonormal frame will be given by

⎡ ⎤
⎡ ⎤
⎛⎢⎢a x ⎥⎥⎞
⎛ a x ⎞ ⎢⎢a x ⎥⎥
ax
(4.13)
= ⎢ ⎥ ♯˚ ⎢ ⎥ =
♭˚
⎝⎢a y ⎥⎠
⎝ a y ⎠ ⎢a y ⎥
⎣ ⎦
⎣ ⎦
For the wedge product and contraction one does not need an orthornormal frame,
because they are metric independent and invariant under a basis change. Therefore,
these operators can be evaluated using the expressions below for any basis frame
that is composed of a linearly independent set of vectors.
Wedge Product
W̊00 ([a], [b]) = [a ⊙ b]
⎡ ⎤
⎢b x ⎥⎞ ⎢a ⊙ b x ⎥
01 ⎛
W̊
[a], ⎢ ⎥ = ⎢
⎢b y ⎥⎠ ⎢a ⊙ b y ⎥
⎣ ⎦
02
W̊ ([a], [b xy ]) = [a ⊙ b x y ]
⎡ ⎤ ⎡ ⎤
⎢ ⎥ ⎢ ⎥
11 ⎛⎢a x ⎥ ⎢b x ⎥⎞
W̊
= [a x ⊙ b y − a y ⊙ b x ]
⎝⎢⎢a y ⎥⎥ ⎢⎢b y ⎥⎥⎠
⎣ ⎦ ⎣ ⎦
Contraction

⎡ ⎤
a x ⎢⎢b x ⎥⎥⎞
C̊1
= [a x ⊙ b x + a y ⊙ b y ]
⎝ a y ⎢⎢b y ⎥⎥⎠
⎣ ⎦
⎞ ⎢⎢−a y ⊙ b xy ⎥⎥
2 ⎛ ax
C̊
, [b xy ] = ⎢
⎝ ay
⎠ ⎢ a x ⊙ b x y ⎥⎥

(4.14)
(4.15)
(4.16)
(4.17)

(4.18)

(4.19)

On Integrated Forms
The corresponding operators on Dec forms can be easily constructed by mapping
from and to the component space via upsampling and downsampling. For example,
for the wedge product,
W̄(a, b) = P̄↓ W̊(P̄↑ a, P̄↑ b).
For the Hodge star one needs to transform to an orthonormal frame:
Λ̄ k
H̄k

P̄↑k

/ Λ̊ k

BB̂

/ Λ̊ k

Bˆ

H̊k

n−k o
̃ n−k o
Λ̊B
Λ̊n−k
n−k
̃n−k
P↓

B̂B

(4.20)

H̄(a) = P̄↓ JBB
ˆ H̊ JBBˆ P̄↑ a.

(4.21)

Similarly, for the inner product,
Ī(a, b) = P̄↓ JBB
ˆ I̊(JBBˆ P̄↑ a, JBBˆ P̄↑ b).

(4.22)

Analogous relationships hold for the dual operators with the bar replaced by a tilde.
4.A

Transformations under Coordinate Change

Let us consider how different objects transform under a coordinate change. Forms
and vectors transform as follows
êi = J ij e j ,

êi = e j (J −1 )ij ,

where

⎡ J1 J1 ⋯ ⎤
⎢ 1 2
⎢ 2 2
J = ⎢ J1 J2 ⋯ ⎥⎥ .
⎢ ⋮
⋮ ⋱ ⎥⎦
Vectors will be represented by column arrays; forms will be represented by row
arrays.
⎡ X1 ⎤
X = ⎢⎢ X 2 ⎥⎥ ,
⎢ ⋮ ⎥

α = [ α1 α2 ⋯ ] .

Vectors X
Vector components transform as
X̂ i = J ij X j
so that
X = X̂ i êi = X i ei
remains invariant under a transformation.
One-Forms (Covectors) Λ1
One form components transform as
α̂i = α j (J −1 )ij

75
so that
α = α̂i êi = αi ei
remains invariant under such transformation. Notice also the contraction between a
vector and a form is independent of the choice of basis:
X ⌟ α = ( X̂ i êi ) ⌟ (α̂ j ê j ) = X̂ i α̂i = X i αi .
Zero-Forms Λ0
fˆ = f
Two-Forms Λ2
Given a two-form
ω = ωi j ei ∧ e j ,
its components transform as
ω̂i j = Ji m J j n ωmn .
Because of anti-symmetry ωi j = −ω ji , the only component in two dimensions will
be ω12 , and the transformation simply becomes
ω̂12 = det(J)−1 ω12
since
ω̂12 = ω12 (J −1 )11 (J −1 )22 + ω21 (J −1 )21 (J −1 )12
= ω12 [(J −1 )11 (J −1 )22 − (J −1 )21 (J −1 )12 ]
= ω12 det (J)−1 .
Again under a coordinate change the contraction will remain invariant:
ω̂( X̂, Ŷ ) = ω(X, Y ),
where
ω(X, Y ) = Y ⌟ X ⌟ ω = X iY j ωi j .
Metric X × X → R
g = gi j ei ⊗ e j ,
ĝi j = gmn (J −1 )im (J −1 )nj,
or using matrix notation
ĝ = J−T gJ−1 .

76
Contraction ⌟ ∶ Λ1 × X → R
X ⌟ α = αX = αi X i
Flat ♭ ∶ X → Λ1
αi = gi j X j
α = X ♭ = XT g
Sharp ♯ ∶ Λ1 → X
X i = gi j α j
X = α♯ = g−1 αT
Inner Products
X×X → R
X ⋅ Y = X i gi j Y j
X ⋅ Y = X T gY
X ⋅ Y = X♭ ⋅ Y ♭ = X ⌟ Y ♭
Λ1 × Λ1 → R
α ⋅ β = αi gi j β j
α ⋅ β = αg−1 βT
α ⋅ β = α♯ ⋅ β♯ = α♯ ⌟ β
4.B

Chebyshev Identities and Limits

The higher order derivatives of Chebyshev polynomials can be expressed in terms
of their lower order counterparts. The basis functions used in this work require the

77
derivatives up to order two.
T0 (x) = 1

U0 (x) = 1

T1 (x) = x

U1 (x) = 2x

TN (x) = 2xTN (x) − TN−1 (x)

UN (x) = 2xUN (x) − UN−1 (x)

TN′ (x) = NUN−1 (x)
TN′′ (x) = NUN−1
(x)

xUN (x) − (N + 1)TN+1 (x)
1 − x2
3xUN (x) − N(N + 2)UN (x)
UN′′ (x) =
1 − x2
UN′ (x) =

For some of the expressions, it is necessary to evaluate TN and UN at the endpoints
x = ±1. While formulas exist for TN in the literature, no similar expressions were
found for UN .
TN (±1) = (±1)N
TN′ (±1) = N 2 (±1)N
TN′′ (±1) = N 2 (N − 1)(N + 1)(±1)N

The expressions below were derived empirically, and they have been checked to
hold true for up to a very large N using a modern computer.
UN (±1) = (N + 1)(±1)N
UN′ (±1) = N(N + 1)(N + 2)(±1)N
UN′′ (±1) = (N − 1)N(N + 1)(N + 2)(N + 3)(±1)N
15

78
Chapter 5

IMPLEMENTATION, VALIDATION, AND RESULTS

5.1

Implementation

The Spex operators are constructed as compositions of the one-dimensional elementary operators that we consider here. For the periodic and bounded domains,
we need to define the following Spex operators for k = {0, 1}:
D̄ k

H̄ k

P̄↑k

P̄↓k

̃k

̃k .

Because of the two possible values for k there are sixteen operators to implement
in total. Only one-dimensional operators are provided, because higher-dimensional
ones can be constructed trivially as compositions of the one-dimensional ones along
each dimension of the computational domain. For example, the wedge product and
the contraction are composite operators that can be expressed in terms of P↑ and P↓
by Eq. (4.20) and Eq. (4.22), respectively.
Building Blocks
The Spex operators are defined in terms of the following basic operators which are
used as elementary building blocks:
̃ ⋆ R̄,
H=P

H−1 = P̄ ⋆−1 R,

̃R̄,
S=P

S = P̄ R,

Q = P̄ ⋆ R̄.

These operators correspond to maps between the discrete cochain spaces on the
periodic domain.
+ ̃1
(5.1)
Λ̄S0 k
−1

̃0 k

H−1

Λ̄1

In the above formulation, it is clear that the basis functions and the geometry of the
cells alone completely determine H, S, and Q through the reconstruction map R
and the reduction map P. However, any implementation that relies on the above
definitions is slow as it must construct dense matrices by evaluating the functions φ

79
on the corresponding cells σ. In this section we provide an implementation via the
discrete fast Fourier transforms (FFT) denoted F and its inverse F ∗ , as well as in
terms of other fast array operations, and use the basis formulation to validate that
the provided implementation is indeed correct.
We can efficiently express our building blocks using Fourier space:
H = F ∗ HF,

S = F ∗̂
SF,

Q = F ∗ QF,

the Fourier space operators are diagonal matrices whose entries are given by
̂ qq = 2 sin ( πk q ),
kq
iπk q
Sqq = exp (
),
̂ qq = 1 (exp ( 2iπk q ) − 1),
ik q
for q = 0, 1, . . . , N − 1, and where k q refers to the FFT sampling frequencies, i.e.,
q < (N + 1)/2
⎪q
kq = ⎨
otherwise
The time complexity of the above operators is thus the same as that of the FFT which
is O(N log N). Having defined the FFT based operators H, Q, and S, we need to
define a few more fast array operations (which we call auxiliary operators). The
auxiliary operators provide the infrastructure upon which the rest of the higher level
Spex operators are most efficiently built. Usually they act on arrays of size N and
have an output that is also an array. The auxiliary operators have time complexity
linear in the size of their input O(N). Table 5.1 lists all the operators that will be
used for defining the full hierarchy of Spex operators. A more detailed description
of each individual operator is provided in 5.A.
Periodic Domain
Derivatives:
D̄0 ( f ) = roll−1 ( f ) − f

D̄1 ( f ) = 0

̃ 0 ( f ) = f − roll+1 ( f )

̃ 1( f ) = 0

Hodge stars:
H̄0 ( f ) = H( f )

H̄1 ( f ) = H−1 ( f )

̃ 0 ( f ) = H( f )

̃ 1 ( f ) = H−1 ( f )

80
Auxilary
Operator

Description

diff

discrete difference of an array

roll

cyclically roll an array

slice

extract elements from an array

even

extract even elements from an array

odd

extract odd elements from an array

weave

merge two arrays by alternating their elements

concat

concatenate two arrays sequentially

mid

extract the middle elements of an array

rev

reverse an array in order

mirror

concatenate an array with its mirror image

Table 5.1: List of auxilary operators
Upsampling:
P̄0↑ ( f ) = weave( f , S( f ))
̃0 ( f ) = weave(̃
S( f ), f )

̃0 (H̄1 ( f ))
P̄1↑ ( f ) = P
̃1 ( f ) = P̄0 (H
̃ 1 ( f ))

Downsampling:
P̄0↓ = even( f )

P̄1↓ = evenodd(Q( f ))

̃0 = odd( f )

̃1 = evenodd(roll+1 (Q( f )))

where
evenodd( f ) = even( f ) + odd( f )
Bounded Domain
In the non-periodic Chebyshev domain, integration requires multiplication by nonuniform weights. Depending on the sampling points (primal or dual) we need to
multiply by one of the following two diagonal weight matrices:
qπ
N −1
̃ qq = sin (q + 2 )π

Wqq = sin

q = {0, 1, . . . , N − 1},
q = {0, 1, . . . , N − 2}.

81
A00 pads the array with zeros on both sides. Abb pads the array by extrapolating the
boundary values and adding them to the array.
A00 ( f ) = concat ([0], f , [0]),
Abb ( f ) = concat ([b( f , −1)], f , [b( f , +1)]),
A† ( f ) = mid( f ),
with the extrapolated boundary values given by

b( f , x) = ∑ fn φ0n (x).
n=0

A† is the left inverse of the padding operators, i.e. A† A00 =Id and A† Abb =Id.
The mirroring matrices are
M±0 ( f ) ∶= concat( f , ±mid(rev( f ))),
M±1 ( f ) ∶= concat( f , ±rev( f )),

M†0 ( f ) ∶= slice0 ∶ ⌊ 2 ⌋+1 ( f ),

M†1 ( f ) ∶= slice0 ∶ ⌊ 2 ⌋ ( f ).
M†0 and M†1 are the left inverses of the above operators and they satisfy M†0 M±0 = Id
and M†1 M±1 =Id. The mirroring operations are explicitly given in section 5.A.
Derivatives:
D̄0 ( f ) = diff(A00 ( f ))

D̄1 ( f ) = 0

̃ 0 ( f ) = diff( f )

̃ 1( f ) = 0

Hodge stars:
̃ −1 (M† (H−1 (M− ( f ))))
H̄0 ( f ) = A† (M†0 (H(M−0 (W(A00 ( f )))))) H̄1 ( f ) = W
̃ 0 ( f ) = M† (H(M− (W(
̃ f ))))

S( f ) = M†1 (S(M+0 (Abb ( f ))))
Q( f ) = M†1 (Q(M−0 (W(A00 ( f )))))

̃ 1 ( f ) = A† (W−1 (M† (H−1 (M− (A00 ( f )))))

S( f ) = A† (M†0 (̃
S(M+1 ( f ))))
̃ f ) = A† (W−1 (M† (Q(M
̃ − ( f )))))
Q(

82
Upsampling:
P̄0↑ ( f ) = weave( f , S( f ))
̃0 ( f ) = weave(̃
S( f ), f )

̃0 (H̄1 ( f ))
P̄1↑ ( f ) = P
̃1 ( f ) = P̄0 (H
̃ 1 ( f ))

Downsampling:

5.2

P̄0↓ ( f ) = odd( f )

P̄1↓ ( f ) = evenodd(Q( f ))

̃0 ( f ) = even( f )

̃1 ( f ) = evenodd(A† (Q( f )))

Convergence and Validation

The Laplace-Beltrami operator is a composition of almost all the operators that we
have implemented in the previous section. On two-dimensional domains, it acts
on 0-forms and 1-forms (2-forms are omitted because we do not use them for our
convergence and validation tests) as follows:
̃ 1D
̃ 1 H1 D0
̃ 1D
̃0
⎪H1 D1 H
∆=⎨
̃ 1D
̃ 0 H2 D1 + D0 H
̃ 2D
̃ 1 H1
1 0 ̃ 2 ̃1 ̃0 2 1 ̃ 1
⎩H D H D + D H D H

Λ̄0 → Λ̄0
̃0 → Λ
̃0
Λ̄1 → Λ̄1

̃1 → Λ
̃1

If we can show the correctness of the Laplace-Beltrami implementation, it is proof
of correctness for implementation of the operators it is composed of as well. For
that purpose, we define the 0-form f 0 and the 1-form f 1 , which will be used as test
functions that will be fed to the ∆ operator:
f 0 (x, y) = e x + e y,

f 1 (x, y) = e x dx + e y dy,

∆ f 0 (x, y) = e x + e y,

∆ f 1 (x, y) = e x dx + e y dy.

Each one will be tested on the grids generated by the diffeomorphisms ϕ ∶ D → M
that will map the reference computational grid to physical space:
ϕ0 ∶

(x, y) → (x, y),

ϕ1 ∶

(x, y) → (x + sin(πx) sin(πy), y + sin(πx) sin(πy))

= 13/100,

ϕ2 ∶

(x, y) → (x, y − sin(πx))

= 1/3,

ϕ3 ∶

(x, y) → (x 2 − y 2, 2x y).

Figure 5.1: Convergence for the Laplace operator in a square domain with rectilinear and
curvilinear Chebyshev grids for 0-forms (top) and 1-forms (bottom).

Figure 5.2: Convergence of the Laplace operator in a non-square domain.

85
In order to show that our method converges, we are going to apply the numerical
Laplacian to the functions above in each one of the grids, and compare the output
to the expected exact solution. As we increase the grid resolution we expect to see
a faster than exponential (supergeometric) decrease in the error. Refer to Figure 5.1
for convergence of the Laplace-Beltrami operator acting on f 0 and f 1 for the
grids generated by ϕ0 and ϕ1 , which have rectilinear boundaries. For curvilinear
boundaries, corresponding to the diffeomorphisms ϕ2 and ϕ3 , refer to Figure 5.2.
The figures display logarithmic plots of the L∞ error between the expected and the
actual output of the discrete operator. On these plots, spectral convergence looks
like curves that develop an increasingly negative slope relative to a straight line.
Eventually the round-off plateau is reached at which point the error flattens out and
begins to slowly increase due to the accumulation errors from the rising number
of floating point operations. Our method converges much faster on flat domains
(N ∼ 15), where no metric related manipulations are performed, and is somewhat
slower (N ∼ 60), when done on curved domains because of the extra computations
related to the transformation to and from an orthonormal frame, which contributes
to the computational error.
5.3

Applications

Having provided the detail of the implementation, and having proved its correctness,
we now proceed to show various applications to problems in physics.
For the purposes of this section, plots for 1-forms are generated using the line integral
convolution (LIC) method [16], which is a popular method to visualize vector fields.
In order to visualize 1-forms, we must first convert them to vector fields by taking
their sharp, and then use the result as input to the LIC algorithm which then outputs
a pixel map of the visualization. Taking the sharp is necessary because 1-forms are
not appropriate for generating curves along which to take convolutions, vector fields
are.
Laplace’s Equation
Laplace’s equation is ubiquitous in physics and mathematics. It is an elliptic partial
differential equation, that can be readily generalized to k-forms denoted as f using
the Laplace-Beltrami operator as
∆ f = 0.
Forms that satisfy Laplace’s equation are known as harmonic forms.

(a) f ∣B = 1

(b) ⋆d f ∣B = 1

Figure 5.3: Laplace’s equation for a 0-form f . The 0-form is represented using a greyscale
colorplot. The Dirichlet boundary condition (a) results in a uniform solution, whereas the
Neumann boundary condition (b) results in a quadratic polynomial.

As is well known from the classical theory of partial differential equations [24], for
scalar fields (i.e., 0-forms in our present framework) we distinguish between two
kinds of boundary conditions - Dirichlet and Neumann. The Dirichlet boundary
conditions fix the function at the boundary f ∣B while Neumann boundary conditions
fix the normal component of its gradient (⋆d f )∣B . We can easily reproduce these
boundary conditions in the present framework (see Figure 5.3).
Because of the generality of the Laplace-Beltrami operator and the possibility to
apply it to higher-dimensional forms, we can solve Laplace’s equation for 1-forms in
addition to scalar fields. In the case of 1-forms the treatment for boundary conditions
is richer than in the case of forms of degree 0. We can either set f ∣B and (⋆d ⋆ f )∣B
along the boundary, or we can set (⋆ f )∣B and (⋆d f )∣B . Note that in the first case we
are prescribing the tangential component and the divergence, whereas in the second
case we are setting the normal component and the curl. Figure 5.4 illustrates the
different possibilities for combining such boundary conditions. Note that the dual
solutions at the bottom row are 90○ rotations of the primal ones at the top. This is
because the two boundary conditions are related by the transformation f ↦ ⋆ f , and
therefore the solutions must also be.
For the exact details of what matrix equation is solved for Laplace’s equation with
boundary conditions and how that corresponds to the different cases described above
see Appendix 5.B.
Further illustration of the generality of our approach is that it is not limited to
flat domains. The Laplace-Beltrami operator has been implemented for non-flat
domains which induce a metric on the computational grid that may be different
from the identity. We demonstrate such domains in Figure 5.5. On the first row we
show the solution on a flat grid that is identical to the reference grid D. This solution

(a) f ∣B = 0

⋆ d ⋆ f ∣B = 1

(b) f ∣B = 1

⋆ d ⋆ f ∣B = 0

⋆ d f ∣B = 1

(d) ⋆ f ∣B = 1

⋆ d f ∣B = 0

Figure 5.4: Plots for the solutions of Laplace’s Equation on a 1-form f with different

combinations of boundary conditions. The domain is [−1, 1]2 , discretized by a 16 × 16
Chebyshev grid. On the top row we solve for f on the primal grid, whereas on the bottom
on the dual grid.

will be used for comparison with the second row, where we preserve the boundaries
exactly but internally deform the grid. Because the boundary remains the same, the
solution on the physical grid must be identical to the solution of the undeformed
domain, which, as is evident from the figure, is indeed the case. This is further
evidence for the validity of the computations and transformations that are performed
in our implementation. While the solutions on the computational reference grids
are different, their push-forwards to the physical domain match exactly. The last two
rows show the solutions for grids that have curvilinear boundaries.
So far we have only considered cases where the physical grid is embedded in 2D,
but there is no difficulty if one wishes to embed in 3D. On Figure 5.6 we illustrate
solutions with boundary conditions, which in the flat case result in a trivial uniform

Figure 5.5: Solution to Laplace’s equation for 1-forms on curved domains with boundary
conditions ⋆ f ∣B = 0 and ⋆d f ∣B = 1. The rightmost column shows the grid, the middle
column shows the solution, and the leftmost column shows the pull-back of the solution to
the reference grid.

Figure 5.6: Laplace’s equation for 1-forms on non-flat domains with boundary conditions
⋆ f ∣Bleft = 1, ⋆ f ∣Bright = −1, ⋆ f ∣Btop,bottom = 0 and ⋆d f ∣B = 0.

Figure 5.7: Laplace’s equation for 1-forms on non-flat domains with boundary conditions
⋆ f ∣B = 0 and ⋆d f ∣B = 1.

91
field. Because the Laplace-Beltrami operator is metric dependent, the deformation
modifies the solution in different ways for each surface. Figure 5.7 illustrates
solutions with boundary conditions that result in a rotational field. Again depending
on the surface, the solutions on the reference grid are modified in comparison to the
flat domain.
Poisson’s Equation
Poisson’s equation is a non-homogenuous differential equation that adds a forcing
term q to Laplace’s equation:
∆ f = q.
In this section we illustrate the usefulness of Poisson’s equation in computing
a divergence-free velocity field that matches a prescribed vorticity field. While
computing the vorticity from a velocity field is straightforward, doing the reverse
is not. Being able to compute the velocity from vorticity is important in numerical
fluid mechanics, where the vorticity field is being advected by the divergence-free
velocity field. Computing the velocity field at each time step would significantly
reduce the memory footprint of the simulation as one would need to store only the
vorticity, rather than both the vorticity and velocity.
The vorticity 2-form of a velocity field 1-form is defined as
ω = dv.
Given a vorticity 2-form, we will compute a divergence-free velocity 1-form from
it. By the Hodge decomposition theorem, the velocity 1-form can be expressed as
v = dα + ⋆d⋆ β + γ
where α ∈ Λ0 , β ∈ Λ2 , and γ ∈ Λ1 is harmonic. Because we require the velocity field
to be divergence-free, d ⋆ v = 0, we must set α = 0. For simplicity we will ignore
the harmonic term and set γ = 0. Then,
v = ⋆d⋆ β.
Applying d to both sides, we obtain
d⋆d⋆ β = ω.
This is precisely Poisson’s equation for 2-forms. Generally the right hand side of
the above equation is thought of as an input of some kind that drives the solution

Figure 5.8: Velocity reconstruction from vorticity. The velocity field is expressed as

v = ⋆d ⋆ β. The leftmost column displays the vorticity 2-form. The middle column
shows the solution for the velocity 1-form on the primal grid with a boundary condition
⋆d ⋆ β∣B = 0. The rightmost column shows the solution on the dual grid with a boundary
condition ⋆β∣B = 0.

93
of the equation. Thus, given ω we can solve for β from which we can compute the
velocity directly, or formally
v = ⋆d⋆(d⋆d⋆)−1 ω,
Figure 5.8 illustrates how the method described here can be used to reconstruct
velocity fields from vortices placed at a variety of positions and with a variety of
strengths. The vortices are modeled as Gaussian functions with widths that are
small relative to the domain size. On the last row, we place many vortices at random
locations in order to show that the method can be used with complicated inputs
rather than just nicely behaved and symmetric ones.
Diffusion Equation
The diffusion equation (also known as the heat equation)
ft = ∆ f
is a parabolic partial differential equation which describes how the density f of a
diffusing material changes with time as matter is transported from regions of high
concentration to regions of low concentration.

Figure 5.9: Diffusion equation at t = 0.0, 0.02, 0.04, 0.08, 0.16, 0.32 with periodic
boundary conditions.

We can easily simulate it using our Spex implementation. Figure 5.9 illustrates the
diffusion from a region of high concentration under periodic boundary conditions
with explicit integration in time. The concentration is shown both as a colormap
and as a heightmap. At first the colored part is separated from the clear part by
a sharp well-defined boundary, but as time passes the boundary becomes less and
less clear with the outside region becoming more intensely colored and the inside
region becoming correspondingly less intensely colored until eventually as t → ∞
(not shown on the figure) the whole domain becomes uniformly colored.

94
Wave Equation
The wave equation
ftt = ∆ f
is a hyperbolic partial differential equation that governs the evolution of a variety of
phenomena in physics such as sound waves, light waves, or water waves.
In its one-dimensional form, the wave equation describes the vibrations of a string.
When the string is fastened at each end, this corresponds to Dirichlet boundary
conditions, whereas when it is allowed to move freely we have Neumann boundary
conditions. To simulate such a string we set up a grid with periodic boundary conditions in the y-direction and periodic, Neumann, and Dirichlet boundary condtions
in the x-direction. We place a Guassian function with an initial condition set to its
derivative along the x-direction and allow it to evolve forward in time. Because of
the symmetry in the y-direction, the problem is essentially a one-dimensional one.
The results are illustrated in Figure 5.10. These solutions can be easily derived using
classical wave theory, but in this case they have been obtained by running simulations
on our Spex framework which provides further validation for its correctness.
We also provide solutions that cannot be derived analytically. Figure 5.11 illustrates
a radially symmetric Gaussian placed at the center of the computational domain
with stationary initial conditions. Figure 5.12 illustrates the same Gaussian with
initial conditions set to the spatial derivative in the x-direction.
Because we have implemented all our operators as tensor products of their onedimensional counterparts it is trivial to set up separate boundary conditions along
each dimension of the computational grid. We show such mixed boundary conditions (Dirichlet along x-direction, Neumann along y-direction) in the last column of
Figure 5.13.
When simulating the wave equation, in order to ensure numerical stability, special
care must be taken that the Courant-Friedrichs-Lewy (CFL) condition
C=

u∆t
≤ Cmax
∆x

is met. The meaning of the variable is that u is the wave velocity, ∆t is the size of
the time step, and ∆x is the size of the grid spacing. For explicit time integrators
(as we use here) Cmax is usually equal to 1. The condition dictates that the largest
stable time step ∆t is proportional to the size of the grid cell. As the computational
grid is refined to increase its spatial resolution, the impact of the CFL condition is

Figure 5.10: One-dimensional wave propagation with periodic, Neumann and Dirichlet
boundary conditions on a 64×8 grid at times t = 0.1, 0.2, 0.3, 0.4, 0.45, 0.48, 0.5, 0.52,0.55,
0.6, 0.7, 0.8, 0.9.

Figure 5.11: Wave propagation with periodic, Neumann, and Dirichlet boundary conditions on a 64 × 64 grid shown at times t = 0.0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8. The initial
condition is set to be a Gaussian function with a time derivative set to zero: ft (x, y, t)∣t=0 = 0

Figure 5.12: Wave propagation with periodic, Neumann and Dirichlet boundary conditions
on a 64 × 64 grid shown at times t = 0.0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8. The initial condition
is set to be a Gaussian function with a time derivative equal to the spacial derivative in the
x-direction: ft (x, y, t)∣t=0 = fx (x, y, 0).

Figure 5.13: Wave propagation with periodic, Neumann, Dirichlet and mixed boundary
conditions on a 64 × 64 grid at t = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 with an initial condition

such that the wave propagates radially outward ft (x, y, t)∣t=0 = [ fx (x, y, 0)2 + fy (x, y, 0)2 ] 2 .
To aid in visualizing the surface a blue light is shone from the right, and a red light is shone
from the left.

99
to decrease the time step and raise the number of iterations required to maintain a
well behaved solution.
Wave Propagation in Curved Space
Motivated by gravitational lensing and light propagation in curved space-time in
astrophysics, we now demonstrate how a wave propagates over a curved twodimensional surface. In order not to distract from the actual geometry of the
surface, we plot the wave as a color-map rather than a height-map as we have done
in the previous examples.
Figure 5.14 demonstrates the effect of a region of non-zero curvature on the propagation of a flat wave front. We let a flat wave move from left to right over a
two-dimensional domain with Neumann boundary conditions. The domain is flat
almost everywhere except a localized dip in the center with non-zero curvature.
As the wave passes over the dip it becomes distorted. The deeper the dip, the
greater the effect on the wave. Although this is a rather simple and naive prototype of gravitational lensing, in principle there should be no difficulty in extending
our Spex framework to four-dimensional space-time in order to compute the actual
propagation of a light wave around regions of high curvature.
Lie Advection
The Lie advection equation, ubiquitous in many advection phenomena in physics,
is given by
ft = −L X f ,
where f is an arbitrary differential k-form that is being advected along the vector
field X. A common context for describing advection in numerical physics is the
advection of scalar fields, whereas no unified framework exists for advecting nonscalar entities such as vorticity in fluid mechanics. Recent works [47, 29] have
shown the benefits of considering the Lie derivative in the context of discrete
exterior calculus by leveraging Cartan’s magic formula which relies on the exterior
derivative d and the contraction ⌟. Both of these operators are implemented in Spex
and therefore we can define a discrete Lie derivative with ease as their composition.
Figure 5.15 illustrates the results of advecting 1-forms along vector fields, using
the discrete Lie derivative of the Spex framework. The top row shows the vector
field X along which the 1-form f is being advected, whereas the rows below display
snapshots of the advected 1-form f at different times. In the first column, the

100

Figure 5.14: Wave equation on a 128 × 128 Chebyshev grid subject to Neumann boundary
conditions. The wave travels from left to right. Snapshots are taken at times t = 0.15,
0.30, 0.45, 0.60, 0.75, 0.93. Time integration is explicit with a time step of ∆t = 10−4 .
The 2D space in which the wave propagates is given by x(u, v) = u, y(u, v) = v, z(u, v) =
−z0 e−9(u +v ) , where z0 = 0, 31 , 12 for the first, second, and third column, respectively.

101

Figure 5.15: Lie advection ft = −LX f .

102
advection is carried along a uniform vector field which results in steady translation
along the direction of the field. In the second column, we have a narrow jet which
deflects a second vertical jet in the middle. In the last column, we have a circular
jet that deflects a vertical jet in mutually opposing directions at each end producing
a distortion that looks like the letter S.
5.A

Auxilary Operators

In this section we review the auxilary operators used in the implementation.
diff
diff computes the discrete difference of an array.
diff ∶

RN → RN−1
⎡ x ⎤
⎢ 0 ⎥ ⎡
⎢ x ⎥ ⎢⎢ x1 − x0 ⎥⎥
⎢ 1 ⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢⎢ x2 − x1 ⎥⎥
⎥ ⎢
⎢ xN−2 ⎥ ⎢⎢
⎥ ⎢
⎥ ⎣ xN−1 − xN−2 ⎥⎦
⎢ xN−1 ⎥

roll
rolln rolls an array by n steps. The order is cyclic. Elements that roll beyond the last
position are reintroduced at the beginning of the array.
rolln ∶

RN → RN
⎡ x ⎤ ⎡ x
⎢ 0 ⎥ ⎢ (0−n) mod N ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢
⎥ ⎢
⎥ ⎢
⎢ xi ⎥ ↦ ⎢ x
⎥ ⎢ (i−n) mod N ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢
⎥ ⎢
⎥ ⎢
⎢ xN−1 ⎥ ⎢ x(N−1−n) mod N ⎥
⎦ ⎣

slice
A slice extracts elements from an array based on a start index, a stop index and a
step. We specify an optional first index, an optional last index, and an optional step.
slicestart ∶ stop ∶ step
The semantics are exactly equivalent to the slice implementation for lists and tuples
a[start ∶ stop ∶ step] in the Python programming language. Slice indices have useful

103
defaults. An omitted first index defaults to zero. An omitted second index defaults
to the size of the array that is sliced. An omitted third index defaults to one, and all
elements are strided sequentially. A negative step means that the array is traversed
in the opposite direction.
even
even picks only the elements at even indices. In terms of slice notation the even
function is given by
even( f ) ∶= slice0∶∶2 ( f ).
Explicitly, the even function can be expressed as
even ∶

RN → R⌊ 2 ⌋
⎡ x ⎤ ⎡ x ⎤
⎢ 0 ⎥ ⎢ 0 ⎥
⎥ ⎢
⎢ x ⎥ ⎢ x ⎥
⎢ 1 ⎥ ⎢ 2 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎢ xN−2 ⎥ ⎢ xN−3 ⎥
⎥ ⎢
⎥ ⎢
⎢ xN−1 ⎥ ⎢ xN−1 ⎥
⎦ ⎣

odd
odd picks the elements at odd locations
odd( f ) ∶= slice1∶∶2 ( f ),
odd ∶

RN → R⌈ 2 ⌉
⎡ x ⎤ ⎡ x ⎤
⎢ 0 ⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ x ⎥ ⎢ x ⎥
⎢ 1 ⎥ ⎢ 3 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎢ xN−2 ⎥ ⎢ xN−4 ⎥
⎥ ⎢
⎥ ⎢
⎢ xN−1 ⎥ ⎢ xN−2 ⎥
⎦ ⎣

104
weave
weave mixes two arrays by alternating between their elements:
weave ∶

RN × RN → R2N ,

⎡ x ⎤
⎢ 0 ⎥
⎢ y ⎥
⎢ 0 ⎥
⎡ x ⎤ ⎡ y ⎤ ⎢ x ⎥
⎢ 0 ⎥ ⎢ 0 ⎥ ⎢ 1 ⎥
⎥ ⎢
⎥ ⎢
⎢ x ⎥ ⎢ y ⎥ ⎢ y ⎥
⎢ 1 ⎥ ⎢ 1 ⎥ ⎢ 1 ⎥
⎥ ⎢
⎥ ⎢
⎢ ⋮ ⎥×⎢ ⋮ ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎢ xN−2 ⎥ ⎢ yN−2 ⎥ ⎢ xN−2 ⎥
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎢ xN−1 ⎥ ⎢ yN−1 ⎥ ⎢ yN−2 ⎥
⎦ ⎢
⎦ ⎣
⎢ xN−1 ⎥
⎢ yN−1 ⎥

concat
Concatenate two arrays by sequentially joining them end-to-end to create an array
of twice the size:
concat ∶

RN × RN → R2N ,

⎢ x0 ⎥
⎢ x ⎥
⎢ 1 ⎥
⎢ ⋮ ⎥
⎡ x ⎤ ⎡ y ⎤ ⎢⎢
⎢ 0 ⎥ ⎢ 0 ⎥ ⎢
⎥ ⎢
⎥ ⎢ xN−2 ⎥⎥
⎢ x ⎥ ⎢ y ⎥ ⎢
⎢ 1 ⎥ ⎢ 1 ⎥ ⎢
⎥ ⎢
⎥ ⎢ xN−1 ⎥⎥
⎢ ⋮ ⎥×⎢ ⋮ ⎥↦⎢
⎥.
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎥ ⎢ y0 ⎥⎥
⎢ xN−2 ⎥ ⎢ yN−2 ⎥ ⎢
⎥ ⎢
⎥ ⎢
⎥ ⎢
⎥ ⎢ y1 ⎥⎥
⎢ xN−1 ⎥ ⎢ yN−1 ⎥ ⎢
⎦ ⎢ ⋮ ⎥⎥
⎦ ⎣
⎢ yN−2 ⎥
⎢ yN−1 ⎥

mid
Pick the middle of an array by removing its endpoints:
mid( f ) ∶= slice1∶−1∶ ( f ),

105
mid ∶

RN → RN−2
⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥ ⎡ x ⎤
⎢ 1 ⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎢ xN−2 ⎥ ⎢ xN−2 ⎥
⎥ ⎣
⎢ xN−1 ⎥

rev
Reverse an array:
rev( f ) ∶= slice∶∶−1 ( f ),
rev ∶

RN → RN
⎡ x ⎤ ⎡x ⎤
⎢ 0 ⎥ ⎢ N−1 ⎥
⎥ ⎢
⎢ x ⎥ ⎢x ⎥
⎢ 1 ⎥ ⎢ N−2 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢ ⋮ ⎥.
⎥ ⎢
⎥ ⎢
⎢ xN−2 ⎥ ⎢ x1 ⎥
⎥ ⎢
⎥ ⎢
⎢ xN−1 ⎥ ⎢ x0 ⎥
⎦ ⎣

mirror
Mirror an array so that the resultant array is (anti-)symmetric around its center. The
boundary elements can either be kept as singletons M0 , or they can be repeated
M1 :
M±0 ∶

RN → R2N−2,
⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥
⎡ x ⎤ ⎢⎢ 1 ⎥⎥
⎢ 0 ⎥ ⎢
⎢ x ⎥ ⎢⎢ ⋮ ⎥⎥
⎢ 1 ⎥ ⎢
⎢ ⋮ ⎥ ↦ ⎢⎢ xN−2 ⎥⎥ ,
⎥ ⎢
⎢ xN−2 ⎥ ⎢⎢ xN−1 ⎥⎥
⎥ ⎢
±x ⎥
⎢ xN−1 ⎥ ⎢⎢ N−2 ⎥⎥
⎦ ⎢ ⋮ ⎥
⎢ ±x1 ⎥

M±1 ∶

RN → R2N ,
⎢ x0 ⎥
⎢ x ⎥
⎢ 1 ⎥
⎢ ⋮ ⎥
⎡ x ⎤ ⎢⎢
⎢ 0 ⎥ ⎢
⎥ ⎢ xN−2 ⎥⎥
⎢ x ⎥ ⎢
⎢ 1 ⎥ ⎢
⎥ ⎢ xN−1 ⎥⎥
⎢ ⋮ ⎥↦⎢
⎥ ⎢
⎥ ⎢±xN−1 ⎥⎥
⎢ xN−2 ⎥ ⎢
⎥ ⎢
⎥ ⎢±xN−2 ⎥⎥
⎢ xN−1 ⎥ ⎢
⎦ ⎢ ⋮ ⎥⎥
⎢ ±x1 ⎥
⎢ ±x0 . ⎥

106
The left inverses of the above operators are

M†0 ∶ RN → R 2 +1,
⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥ ⎡ x ⎤
⎢ 1 ⎥ ⎢ 0 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢ x ⎥
⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ xN−1 ⎥ ↦ ⎢ ⋮ ⎥ ,
⎥ ⎢
⎥ ⎢
⎢ xN ⎥ ⎢ xN−2 ⎥
⎥ ⎢
⎥ ⎢
⎢ ⋮ ⎥ ⎢ xN−1 ⎥
⎥ ⎣
⎢ x2N−1 ⎥

5.B

M†1 ∶

RN → R 2 ,
⎡ x ⎤
⎢ 0 ⎥
⎢ x ⎥ ⎡ x ⎤
⎢ 1 ⎥ ⎢ 0 ⎥
⎥ ⎢
⎢ ⋮ ⎥ ⎢ x ⎥
⎥ ⎢ 1 ⎥
⎥ ⎢
⎢ xN−1 ⎥ ↦ ⎢ ⋮ ⎥ .
⎥ ⎢
⎥ ⎢
⎢ xN ⎥ ⎢ xN−2 ⎥
⎥ ⎢
⎥ ⎢
⎢ ⋮ ⎥ ⎢ xN−1 ⎥
⎥ ⎣
⎢ x2N ⎥

Boundary Conditions for Laplace’s Equation

In this section, we illustrate in detail how to solve Laplace’s equation with different
boundary conditions.
0-forms For 0-forms, we can pick either Dirichlet or Neumann boundary conditions by choosing where to solve the discrete equation - on the primal or the dual
grid. Note that the terms Dirichlet and Neumann are meaningful only for 0- forms.
There equivalents for 1-forms will be described later on.
Dirichlet
f ∈ Λ̄0,

b ∈ B0,

̃ 1 H D̄0 f = 0,
HD
̃ 1 H(D̄0 f + D0 b) = 0.
HD
Solve for f :
̃ 1 H D̄0 f = −H D
̃ 1 HD0 b,
HD
corresponding to the boundary condition
f ∣B = b.
Neumann

̃ 0,
f ∈Λ

β ∈ B1,

̃ 0 f = 0,
H D̄1β H D
̃ 0 f ) + D1 β] = 0.
H[D̄1 (H D

107
Solve for f :
̃ 0 f = −HD1 β,
H D̄1 H D
corresponding to the boundary condition
(⋆d f ) ∣B = β.
1-forms For 1-forms, on the other hand, we will have a richer set of boundary
conditions.
Primal
α ∈ Λ̄1,

b ∈ B0,

β ∈ B1,

̃1H + HD
̃ 0 H D̄1 )α = 0.
(D̄0b H D
Solve for α:
̃1H + HD
̃ 0 H D̄1 )α = −D0 b − H D
̃ 0 HD1 β,
(D̄0 H D
corresponding to the boundary condition

Dual

(⋆d ⋆ α) ∣B = b

α ∣B = β.

̃ 1,
α∈Λ

β ∈ B1,

b ∈ B0,

̃1 + D
̃ 0 H D̄1 H)α = 0.
(H D̄0b H D
Solve for α:
̃1 + D
̃ 0 H D̄1 H)α = −HD0 b − D
̃ 0 HD1 β,
(H D̄0 H D
corresponding to the boundary condition
(⋆dα) ∣B = b,

(⋆α) ∣B = β.

̃ with (b, β) = {(0, 1), (1, 0)}.
Figure 5.4 illustrates the four possibilities for ᾱ and α

108
Chapter 6

PROGRAMMING CONTRIBUTIONS
Over the course of this thesis a number of libraries have been developed that have
been made freely available at github.com/drufat under an open source "copyleft"
license. The objective has been to enable third parties to easily reproduce and
replicate the results in this work while avoiding unnecessary duplication of effort.
Too often nowadays, one is faced with the task of recreating existing libraries from
the ground up instead of building on top of what has already been created by others.
What follows is an overview of the various programming contributions that have
been made both as standalone packages and as additions to existing third party
libraries.
6.1

Spexy

The Python module Spexy is the reference implementation of spectral exterior calculus. It provides a number of classes for differential forms, both in their discrete
and continuous representations, with all of the associated geometric operators implemented as methods and attributes to those classes. Additionally, classes are
provided for periodic and bounded grids that discretize the abstract surfaces on
which the aforementioned forms reside.
Three differential form classes are provided in the spexy.form module
• sym.Form - forms in their symbolic representation,
• coch.Form - forms as cochains,
• comp.Form - forms as pointwise components,
along with the ability to seamlessly map between them. Each of the differential
form classes implements the standard set of exterior calculus operators, sometimes
directly and other times by transforming into a more convenient representation in
which the operator is easier to implement.
Given an instance f of a differential form class, we can access its derivative as f.D,
its Hodge star as f.H, the wedge product of two forms f0 and f1 as f0.W(f1), the
inner product of f0 and f1 as f0.I(f1), and the contraction of f with a vector v as

109
v.C(f). For convenience, the wedge product and the contraction are also provided
via the built-in binary operators ^ and | respectively. Thus, v^f is an alias for
v.W(f), and v|f for v.C(f).
The grid classes, on the other hand, contain data about the discretization of the computational domain, more specifically about its representation as a cellular complex.
They provide a complete description of the positions and relative arrangement of the
cells. The grid classes can be imported from the spexy.grid module. Currently
only one- and two-dimensional grids are supported with the prospect of extending
to higher dimensions in the future.
• Grid_1D - one-dimensional grids
• Grid_2D - two-dimensional grids
In this work we deal with two main types of one-dimensional grids – periodic and
bounded ones. Higher-dimensional grids are easily formed by taking Cartesian products of the one-dimensional ones, and different combinations along each dimension
are possible. The code below demonstrates how to instantiate grid objects:
from spexy.grid import Grid_1D, Grid_2D
g1 = Grid_1D.chebyshev(10)
g2 = Grid_2D.periodic(10, 10)
g3 = Grid_1D.periodic(10) * Grid_1D.chebyshev(12)

In the example above g1 represents a one-dimensional bounded (Chebyshev) grid
which has 10 primal edges. g2 is an instance of a two-dimensional 10 × 10 periodic
grid. g3 demonstrates the possibility of mixing different grid types by taking their
Cartesian product along each direction - periodic along the x-axis, and bounded
along the y-axis.
Every grid object stores a reference to its dual which can be accessed via g.dual.
Since the dual of a dual is the identity, g.dual.dual becomes a circular reference
to the original grid g.
Having described the grid and form classes provided by Spexy, we now turn our attention to how these two main classes interact with each other. Given a symbolic form
fsym we can project it onto a grid g via the reduction operator fcoch = fsym.P(g)
to obtain a cochain form fcoch associated with the cells of the grid g. Conversely, given a cochain form fcoch, we can reconstruct the original form via

110
the reconstruction operator fcoch.R. The upsampling operator can be accessed
via fcomp = fcoch.Pup to output a component form fcomp, whereas given a
component form fcomp one can downsample it via fcoch = fcomp.Pdown(g)
or fcoch = fcomp.Pdown(g.dual) onto a primal or a dual cochain form fcoch
respectively.
The form class encapsulates the degree of the differential form, which for a form
instance f can be accessed by f.degree. The degree is required in order to choose
the correct operator. For example, in two dimensions one can choose from multiple
Hodge star operators {H0, H1, H2 } depending on the degree of the input form. With
the degree encapsulated in the form data type, rather than the programmer having
to manually choose the operator matching the degree, f.H automatically selects
the correct operator, and assigns the correct degree to the output form. In this
way the user is spared the mundane and error prone task of degree tracking. This
convenience is provided for all of the Spex operators. The degree computations
are done in a manner that is consistent with the mathematical definition of each
operator, i.e. for any form f the assertions below are always true:
assert f.D.degree
== f.degree + 1
# Derivative
assert f.H.degree
== ndim - f.degree
# Hodge star
assert (f ^ g).degree == f.degree + g.degree # Wedge product

where ndim is the dimension of the computational domain (typically 1 or 2).
Additionally, cochain forms also encapsulate the grid with which they are associated.
Given a form f, its grid can be accessed via f.grid. With respect to the grid, the
following assertions hold:
assert f.D.grid == f.grid
assert f.H.grid == f.grid.dual

# Derivative
# Hodge star

These assertions reflect the fact that when acting on cochain forms, the derivative
leaves the grid unchanged, whereas the Hodge star flips to the dual grid.
Given these basic Spex operators, one can implement composite ones such as the
Laplacian and the Lie derivative. For example, the Laplacian can be implemented
as:
def lap(f):
return f.H.D.H.D + f.D.H.D.H

Whereas, the Lie derivative is:

111
def lie(v, f):
return v.C(f.D) + v.C(f).D

Examples
Let us illustrate what kind of computations one can do with the Spexy library. In
order to successfully execute the examples below in the Python interpreter, we need
to import a few modules and initialize a number of variables as below:
import sympy as sy
from spexy.form import sym
from spexy.grid import Grid_2D
x, y = sy.symbols('x, y')
F0, F1, F2 = sym.Form.forms(x, y)
g = Grid_2D.chebyshev(8, 8)

First, we import the symbolic form module sym. The module attribute sym.Form.forms
provides convenient initializers F0, F1, and F2 that make it easy to instantiate zero-,
one-, and two-dimensional symbolic forms respectively. Finally, on the last line we
instantiate an 8 × 8 Chebyshev grid, and we are ready to do actual computations with
forms on that grid.
Gradient First, we compute the derivative of the zero-form x + y.

fsym = F0(x + y)
f0 = fsym.P(g)
bc = fsym.Dbc(g)
f1 = f0.D + bc
plotf(f0)
plotf(f1, color='black')

We instantiate the zero-form into the variable fsym, then project it onto the primal
grid. Because the primal grid has associated boundary vertices, the derivative
operator will require an additive boundary condition term bc, which can be computed
using the attribute Dbc. The boundary condition term is obtained by sampling the
original form at the boundary vertices. On the next line, we assign the result to

112
the f1 variable. Finally, we use the plotf function to visualize the results. We
plot one-forms indirectly by converting them to vector fields via the sharp operator,
which in the case of a flat surface is trivial.
In the next example we compute the derivative of the zero-form x 2 − y 2 , but instead
of projecting onto the primal grid, we project onto the dual. Because the dual grid
has no boundary vertices, there will be no associated boundary condition term.

f0 = F0(x ** 2 - y ** 2).P(g.dual)
f1 = f0.D
plotf(f0)
plotf(f1, color='black')

The next example is identical to the previous one, except we compute the derivative
of x 2 + y 2

f0 = F0(x ** 2 + y ** 2).P(g.dual)
f1 = f0.D
plotf(f0)
plotf(f1, color='black')

As expected, the derivative of a zero-form corresponds to the gradients of the a
scalar field.
Divergence The divergence of a vector field can be computed by taking the ⋆d⋆ of
its corresponding one-form. The example below demonstrates how to compute the
divergence of the vector field corresponding to the one-form A(x dx + y dy). The
coefficient A = e−(x +y ) ensures the field is localized near the center and decays to
zero away from it.

113

A = sy.exp(-(x ** 2 + y ** 2))
f1 = F1(A*x, A*y).P(g)
f0 = f1.H.D.H
plotf(f0)
plotf(f1, color='black')

Curl The curl can be computed by applying the ⋆d operator to the corresponding
one-form, and in the example below we do the computation for A(−y dx + x dy)

A = sy.exp(-(x ** 2 + y ** 2))
f1 = F1(-A*y, A*x).P(g.dual)
f0 = f1.D.H
plotf(f0)
plotf(f1, color='black')

Hodge star On a flat 2D surface, the Hodge star of a one-form corresponds to an
anti-clockwise rotation by 90○ . The two examples below illustrate that by applying
⋆ multiple times we can achieve such rotations.

f = F1(1, 0).P(g)
plotf(f, color='blue')
plotf(f.H, color='purple')
plotf(f.H.H, color='red')
plotf(f.H.H.H, color='green')

114

f = F1(-y, x).P(g)
plotf(f, color='blue')
plotf(f.H, color='purple')
plotf(f.H.H, color='red')
plotf(f.H.H.H, color='green')

Laplacian The Laplacian is a composite operator that can be expressed as ⋆d ⋆
d + d ⋆ d⋆. In the example below, we demonstrate how to compute the Laplacian
of the one-form y 2 dx + x 2 dy. Special care must be taken to include the boundary
terms when applying the D operator on primal forms. Otherwise, one would get
incorrect results near the boundary.

fsym = F1(y**2, x**2)
= fsym.P(g)
bc0 = fsym.H.D.H.Dbc(g)
bc1 = fsym.Dbc(g)
fl = (f.H.D.H.D + bc0) + (f.D + bc1).H.D.H
plotf(f, color='blue')
plotf(fl, color='red')

Because the Laplacian is a second derivative, we expect the result of a quadratic
function to be a constant field, which is indeed the case.
Wedge product Taking the wedge product of two one-forms is like taking the
cross product of their corresponding vector fields. The resulting two-form, which
is converted to a zero-form before plotting, is proportional to the sine of the angle
between the two fields.

115

f1 = F1(x, y).P(g)
f2 = F1(y, x).P(g)
f3 = f1 ^ f2
plotf(f3.H)
plotf(f1, color='white')
plotf(f2, color='black')

Contraction Let fv be a one-form, and fk a k-form. Below we illustrates the
contraction fv♯ ⌟ fk of fk with the vector field corresponding to fv. The contraction
of a one-form f1 with fv♯ is equivalent to the inner product between f1 and fv,
and thus the resulting scalar field must be proportional to the cosine of the angle
between the two fields. The plot below demonstrates that is indeed the case.

fv = F1(x, y).P(g)
f1 = F1(y, x).P(g)
f0 = fv | f1
plotf(f0)
plotf(fv, color='white')
plotf(f1, color='black')

The next programming snippet shows the contraction of a two-form f2 with fv♯ .

fv = F1(x, y).P(g)
f2 = F2(x + y).P(g)
f1 = fv | f2
plotf(f2.H) #convert to 0-form
plotf(fv, color='white')
plotf(f1, color='black')

116
Lie derivative By relying on Cartan’s magic formula we can also compute the Lie
derivative.

fv = F1(1, 0).P(g.dual)
f1 = F1(x, x).P(g.dual)
f2 = fv.C(f1.D) + fv.C(f1).D
plotf(fv, color='blue')
plotf(f1, color='black')
plotf(f2, color='red')

The derivative of the black field is computed along the blue field and the result is
plotted as the red field. The result does indeed appear to be equal to the rate of
change of the input black field along the chosen direction.
6.2

LicPy

There are a number of techniques to image vector fields. The most common one,
known as hedgehog, is to draw an arrow at each point with a heading and a size
matching the direction and magnitude of the vector field at that particular point
(e.g. quiver in Matplotlib). Recently Cabral and Leedom proposed a new method
to image vector fields called line integral convolution (LIC) [16]. LicPy is our
implementation in Python of the LIC method.
The method requires one to take a random noisy background texture, and for each
pixel of the texture, streamlines of equal length are computed along the input vector
field. One then integrates the background pixel values along each streamline to
produce a single pixel value and a corresponding output texture. The output looks
like a blurred version of the input, where the blurring is along the direction of the
vector field that is being visualized, and is easy for the human eye to interpret.
After generating the random texture, one needs to ensure that the input vector
field and the random texture are of the same size, i.e. for each texture pixel
there must be a single corresponding vector field value. If that is not the case,
for example, if the vector field is sampled less densely, one can easily interpolate in a piecewise linear fashion to compute a value per texture pixel. Then,
given a point (marked with x in Figure 6.1) one computes a stream-line going
through that point both in the forward (as on the figure) and backward directions.

117
This is easy to do since the vector field
is piecewise constant within each pixel,
so the streamlines are straight line segments within each pixel. Computing
them is merely a matter of figuring
out where the straight line intersects
the rectangular boundaries of the pixel.
Once the streamlines are computed, the
random background texture colors are
averaged over their length, resulting in
a new color for each pixel.
Figure 6.1: Streamline
The output image produced in this fashion is more intuitive to the human eye to interpret because it looks like the smearing
of ink along a moving fluid. Figure 6.2 illustrates a number of vector fields that are
visualized using both quiver plots and LIC plots. Notice how the LIC plots remain
well-behaved even in the presence of singularities near which the quiver plots are
hard to parse.
6.3

PyBindCpp

PyBindCpp is a Python module that we have developed in order to make it easy
to interface Python with C++. In this section we describe the approach taken by
PyBindCpp, and compare it to other similar libraries.
Python has been widely successful as a programming language for scientific computations primarily due to two factors. First of all, as a high level language it is
highly readable and expressive, which allows for rapid development. Second of all,
because of its open source nature and the access to its internal CPython1 application
programming interface (API) it is possible to interface Python with the vast collection of existing scientific libraries without having to reimplement them from the
ground up. However, despite these advantages, Python has one important drawback.
Although adequate for most applications where speed is not a critical requirement,
because of its highly dynamic nature, Python tends to be slow, especially when dealing with tight for-loops, and it is not suitable for computations where such low-level
looping constructs are necessary.

The reference implementation of the Python programming language is done in C and is known
as CPython. It can be downloaded from python.org.

118

Figure 6.2: Comparison of quiver plots generated using Matplotlib with line integral
convolution plots generated by LicPy.

119
To ameliorate this problem, a number of approaches have been proposed for a hybrid
programming model where most of the high level logic is written in Python, and the
time critical parts of the application, which usually tend to be a small fraction of the
code base, are reimplemented in a low level language, significantly speeding up the
application running time.
Among the most popular approaches is the one taken by Cython[6]. In fact, Cython is
not merely a library, but an entire programming language, which is a strict superset
of the Python language. This means that all Python programs are valid Cython
programs, but the converse is not true. Cython adds type annotations to the Python
language, where each variable can be optionally assigned a static type. This extra
information allows for the source code to be translated directly into C, and compiled
and linked against Python, and then imported as a module. While the approach
taken by Cython is quite powerful, it does have its drawbacks. First, it seems to
be mainly geared towards speeding up existing Python applications by adding type
annotations to the critical parts, and not towards porting existing C/C++ libraries.
Interfacing with external libraries is quite cumbersome as one has to copy all of
the header file declarations into a form that can be parsed by Cython. This violates
the "don’t repeat yourself" (DRY) principle where every piece of knowledge must
have a single unambiguous representation within the system. In the case of Cython,
however, declarations must exist both in the C/C++ header files as well as within
Cython files. Additionally, just to bridge Python and C/C++ one has to learn a
completely new third programming language with its new syntax and semantics.
Another approach for speeding up Python is to write the low level glue code manually
in C/C++. This requires intimate knowledge of the CPython API as well as the
addition of a large amount of boilerplate code for manual data type conversion and
memory management. While this gives a lot of control to the programmer, it also
places limitations on productivity due to the time consuming nature of programming
at such a low level. Additionally, this approach intimately depends on the internals
of the original Python implementation, and is not portable to other implementations
of the same programming language such as PyPy, Jython, or IronPython.
A third approach is to use the Ctypes module that is included by default in the
Python interpreter. Ctypes is a foreign function library for Python that provides
C compatible data types that can be used to call functions directly from compiled
DLL or shared library files. Ctypes works at the application binary interface (ABI)
level rather than the application programming interface (API) level. One still needs

120
to provide descriptions in Python of the exact signatures of the functions that are
being exported. In fact one has to be quite careful, because at the ABI level if there
is any mismatch, one can get quite nasty segmentation faults and data corruption
errors, whereas the worst that can happen at the API level is for the program not to
compile. Again, like the approach taken by Cython, this violates the DRY principle
since declarations must be copied into Python.
PyBindCpp is our own proposed solution to the problem which leverages Ctypes and
C++14 in order to provide highly portable and seamless interoperability between
C/C++ and Python. Most of the logic is implemented in Python via the Ctypes
module which enables one to call functions from shared library files and do the
appropriate type conversions. The key insight is to use C++14 for type deductions,
and to feed that information back to Python in order to generate the appropriate
Ctypes wrappers. In order to deduce the types of the variables and the signatures
of the functions, we have relied heavily on recent features added to C++ such as
template metaprogramming, type traits, parameter packs, and the auto specifier.
Because of this, the minimum version of C++ necessary to compile PyBindCpp
is C++14. Without these features, the only way to obtain the type and signature
information of functions would have been to parse C++, which is a very difficult
problem given the size and complexity of the language.
In designing PyBindCpp, we have been inspired by ChaiScrit[70], and we have
borrowed heavily from its design. ChaiScript is a scripting language that targets
C++ directly and takes advantage of its many modern features. The typesystem of
ChaiScript is identical to the C++ type system with a few minor additions. Unlike
ChaiScript, PyBindCpp is not a new language, but a module of Python, and it
provides its own type conversion layer between C++ data types and Python data
types. Like ChaiScript, PyBindCCpp is headers only and does not require any
external tools or libraries to compile, just a modern C++ compiler.
PyBind11[38] is another library that is very similar in its design goals to PyBindCpp.
However, there are some significant differences. Most of the programming logic in
PyBind11 is coded in C++ and the library relies heavily on the CPython API making
it difficult to port to other implementations of the Python language. PyBindCpp,
on the other hand, uses Ctypes which is provided by most Python interpreters, and
can therefore be easily ported. A design goal of PyBindCpp is to minimize reliance
on the API to a bare minimum, perhaps even eliminate it completely, and instead
leverage the ABI via Ctypes. Additionally, PyBindCpp is coded mainly in Python,

121
and C++ is minimally used for type deduction. Many of the Python functions are
made available to the C++ level via Ctypes callbacks. This saves us the effort of
having to code them in C++, especially for code that is rarely executed (e.g., only
once at import time) where speed is not that critical. This small speed penalty is
more than offset by the simplicity of the code when written in Python.
Next we present a minimal working example of PyBindCpp usage. Consider the
following C++ code which includes variables and functions that we wish to export
to Python:
#include
constexpr double half = 0.5;
constexpr double pi = M_PI;
int
f(int N, int n, int x)
return N + n + x;
double
mycos(double x)
return cos(x);

To port the above code, we merely need to write a single function before compiling
into a shared library that can be directly imported from Python. We must define a
function which takes as its input a reference to an ExtModule object, and uses that
to register the entities that need to exported. We can use the var attribute to register
simple types, and the fun attribute to register functions. Finally, we must add the
PYMODULE_INIT preprocessor macro which takes two arguments - the name of the
module that we are creating and the name of the function that we just defined. This
is necessary in order for Python to find the correct entry point when loading the
extension module.
#include
using namespace pybindcpp;
void
module(ExtModule& m)
m.var("half", half);
m.var("pi", pi);

122
m.var("one", (long)1);
m.var("two", (unsigned long)2);
m.var("true", true);
m.var("false", false);
m.var("name", "pybindcpp");
m.fun("f", f);
m.fun("mycos", mycos);
m.fun("cos", [](double x) -> double { return cos(x); });
m.fun("sin", static_cast(sin));
PYMODULE_INIT(example, module)

After compiling the code above into a shared module, we can use the script below
to import it from Python and run some tests to ensure that everything works as
expected.
import math
import example as m
def test_example():
assert round(m.pi, 2) == 3.14
assert m.half == 0.5
assert m.one == 1
assert m.true == True
assert m.false == False
assert m.name == b'pybindcpp'
assert m.f(1, 2, 3) == 1 + 2 + 3
assert m.mycos(0.1) == math.cos(0.1)
assert m.cos(0.1) == math.cos(0.1)
assert m.sin(0.1) == math.sin(0.1)
if __name__ == '__main__':
test_example()

6.4

Third Party Modules

Additional code contributions made during the course of our graduate studies in the
form of patches and pull requests have been submitted to third party modules that
are maintained by others. The most significant additions have been made to the
Python libraries GlumPy[61] and SymPy[46].
GlumPy is a library for creating scientific visualizations. It provides an object
oriented interface to OpenGL and even a novice not acquainted with all the intricate
details of the OpenGL API can use GlumPy to generate advanced 2D and 3D
graphics and animations. The module comes packaged with many examples, along

123
with complete vertex and fragment shaders, that allow the user to quickly get up
to speed with the capabilities of the library. Our contributions have been to port
the library to Python version 3, to simplify the animation interface, to fix minor
bugs, and to add extra examples. We use GlumPy primarily for drawing scalar- and
vector-fields on non-flat surfaces.
SymPy is a library for symbolic computations and manipulations. Its capabilities
include basic arithmetic operations, expression simplification, symbolic integration
and differentiation, taking limits, and much more. An important functionality for
our purposes is the ability to "lambdify" symbolic expressions, meaning they can
be converted to programming functions that are directly callable from code to get
numerical results. Our contributions have centered around expanding the types
of expressions that can be lambdified to include sums, piecewise expressions, and
relational operators. The lambdify method is crucial for us in order to convert
symbolic differential forms to their numerical counterparts. Throughout this work
we have relied heavily on SymPy for computations with symbolic differential forms,
and to verify that basis functions associated with grid cells do indeed satisfy the
properties that are required of them.
Additionally, although no contributions have been made to them, in this work
we have extensively used NumPy[71] and Matplotlib[37]. NumPy is a Python
module that deals with multi-dimensional arrays, and provides a large number of
fast operators on those arrays, including many mathematical functions from linear
algebra. Matplotlib, on the other hand, is a plotting library which can be used for
line, contour, scatter, and vector plots, and has been helpful in generating many of
the convergence graphs and mesh visualizations in this work.
All of the third party packages mentioned so far are required dependencies of
the Spexy module, and without them the work on implementing spectral exterior
calculus would have been much more laborious and involved. The author is grateful
for the availability of the aforementioned packages as open source libraries. Access
to the source code has enabled us to learn from the work of others, to easily modify
the software when necessary, to port to newer versions of Python, and to submit
improvements to the authors for inclusion in future releases. In return, by releasing
all our modules under a permissive open source license we hope to give back to the
community, and to similarly benefit from contributions and new functionality by
others.

124
Chapter 7

CONCLUSION
In this chapter we review the main contributions of this thesis followed by a discussion of possible future avenues for extending the present work.
7.1

Review of Contributions

The main contribution of this thesis has been to describe and implement a complete
framework of spectral exterior calculus (Spex) that encompasses the entire hierarchy
of differential geometry operators (d, ⋆, ∧, ⌟, ∆, L, ⋅, ♭, ♯).
We have introduced a twin representation for discrete differential forms, both as
cochains and as pointwise components along with the ability to seamlessly transform
between them. The cochain representation is appropriate for topological operators
such as the derivative, whereas the component representation is more adapted for
operators that require multiplication, including the metric dependent ones.
The framework is general enough to work on any domain with logically rectangular
boundaries that can have any arbitrary metric. Numerous structural identities from
the continuous world are preserved in the computational realm, such as Stokes’
theorem, Cartan’s magic formula, or involutivity of the Hodge star.
Boundary conditions are straightforward to incorporate as additive terms to the
discrete derivative operator. Multiple boundary conditions as well as their combinations have been implemented and demonstrated.
Guarantees are provided for spectral convergence of all the operators, all of which
are computed efficiently and implemented in the most appropriate form using the
fast Fourier transform with a time complexity of O(N log N).
A reference implementation is provided in the Python programming language as a
self-contained open source library. The code has been extensively unit tested. It has
been initially set up in 1D, and then extended trivially to 2D through tensor products.
In principle, there should be no difficulty in extending it to even higher-dimensional
spaces like 3D space or the 4D space-time of general relativity with the Minkowski
metric.

125
7.2

Future Work

There are numerous avenues that one can pursue in order to extend the work presented
in this document.
Domain Decomposition
In this thesis we have focused on trivial topologies given by our periodic and bounded
domains. However, general manifolds can have far richer topologies, and to be able
to adequately capture their structure we need a method to combine multiple simple rectilinear domains into more complicated surfaces by patching them together.
The theory of domain decomposition [45] is already well established in numerical
analysis. It describes a way to split a computational grid into multiple domains,
that while coupled, can each be simulated separately in order to better leverage
parallel computer architectures. While we want to achieve the opposite, combine
computational domains into more complicated ones rather than split existing ones,
the basic results still apply. We can adapt domain decomposition to our present
framework and implement all the operators of Spex on a general manifold. The
individual domains that form the discrete manifold will still be fairly independent
of each other with the coupling given by the boundary term of the d operator for
each domain. For example, if A and B are neighboring domains, then the derivative
operator on A, denoted by d A, will have a boundary term due to domain B, and vice
versa for dB . All other operators will be completely independent of each other, and
information will be exchanged between the domains only via the boundary term of
the derivative operator.
Infinite and Semi-infinite Domains
In his book on numerical spectral methods [11], Boyd describes infinite and semiinfinite domains and their associated basis functions, in addition to the periodic and
bounded domains considered here. We can use those basis functions in order to
extend the Spex framework to those domains as well. Many simulations in fluid
mechanics are actually implemented in exactly such infinite domains [66].
Spectral Connection
On a Euclidian space one can compare vectors at different points directly, but on a
general manifold one needs to introduce the connection in order to compare vectors
at different tangent spaces. There have been various works proposing structurepreserving approaches to disretization of connections [19], but those have been

126
low order. Adapting the present work to a discrete connection and its associated
covariant derivative that converges spectrally may prove to be a fruitful undertaking.
It would also allow us to spectrally compute the curvature and torsion 2-forms for
any discrete computational domain.
Spex on the GPU
All of the Spex operators are ideally suited to be implemented on parallel computer
architectures. As heterogeneous computer architectures become more prevalent,
it may be useful to port the present Spex implementation to the GPU in order to
speed up its execution and enable realtime simulations with instant feedback for the
developer rather than relying on long-running computations on the CPU.

127
Appendix A

APPENDIX

A.1

Invariance of the Discrete Wedge Product

Consider a manifold M and its polyhedrization in the context of discrete exterior
calculus, and let σi denote the simplices and φi the corresponding basis functions,
such that the pairing given by integration is natural:
⟨σi, φ j ⟩ ≡
φ j = δi j .
σi

Given the reduction
P∶

Λ k → Λ̄ k

α ↦ ᾱi =

σ̄i

and reconstruction operators
R ∶ Λ̄ k → Λ k
ᾱi ↦ α = ᾱi φi
the discrete wedge product is given by
γ̄ = W(ᾱ, β̄) = P ((Rᾱ) ∧ (R β̄))
γ̄i = ∑ Wi j k ᾱ j β̄k ,
jk

where
Wi j k = [W(φ j , φ k )]i =

σi

φ j ∧ φk .

Consider now a second manifold M′ and a diffeomorphism
ϕ∶

M→M,

which induces a new polyderization with the corresponding simplices σ ′ = ϕ∗ σ and
basis functions φ′ = ϕ∗ φ. Given a simplex σ and a form φ, the following identities
hold for push-forwards and pull-backs:
ϕ∗ φ =
φ,
(A.1)
ϕ∗ σ

128
ϕ∗ (α ∧ β) = (ϕ∗ α) ∧ (ϕ∗ β).

(A.2)

The discrete wedge product on this new manifold will be elementwise equal to the
wedge product on the original manifold. To see why this is the case:
φ′j ∧ φ′k
Wi j k =
ˆ i
(ϕ∗ φ j ) ∧ (ϕ∗ φ k )
ϕ∗ σi
ϕ∗ (φ j ∧ φ k ) (eq.A.2)
ϕ∗ σi
φ j ∧ φ k (eq.A.1)
σi

= Wi j k .
So we have shown that the discrete wedge product is indeed invariant under diffeomorphisms:
Wi′ j k = Wi j k .
This result is consistent with the theory in differential geometry, where the wedge
product is thought of as a purely topological operator that is independent of the
metric.
A.2

Involutivity of the Discrete Hodge Star Operator

In the continuous case, the Hodge star is its own inverse (up to a sign):
⋆n−k ⋆ k = ⋆ k ⋆n−k = (−1) k(n−k) id,

(A.3)

where id is the identity map. We would like to mirror this property in the discrete
setting. For that, the discrete mesh (representing a discrete polyhedrization of the
continuous manifold) must meet certain requirements. We want to study under
what conditions the discrete Hodge star operator satisfies the discrete counterpart
of Eq. (A.3):
̃ = HH
̃ = (−1) k(n−k) Id,
HH
(A.4)
̃ n−k and H
̃ n−k → Λ̄ k , i.e. under what conditions is it an
̃ ∶Λ
where H ∶ Λ̄ k → Λ
exact involution. The pairing between simplices and bases functions is given by
integration
⟨σi, φ j ⟩ ≡

σi

φ j ∈ R.

129
The cardinal basis functions φ associated with the simplices σ must by definition
satisfy
⟨σi, φ j ⟩ = δi j = ⟨̃
σj , φ̃i ⟩.
The elements of the discrete Hodge star matrix acting on primal (H) and dual (H)
forms are given respectively by:
̃ ⋆ R,
H=P

̃ = P ⋆ R,

Hi j = ⟨̃
σi, ⋆φ j ⟩,

̃i j = ⟨σi, ⋆φ̃j ⟩.

These are the definitions of the discrete Hodge star operators, and in the general
̃ is an exact inverse of H as mandated by Eq. (A.4). In
case it is not necessary that H
order to see when condition Eq. (A.4) holds, let us construct the dual basis functions
as linear combinations of the primal basis functions, namely
⋆ φ̃i = ∑ Ai j φ j ,

(A.5)

where the matrix A represents the transformation between the two basis sets. The
Hodge star (⋆) is present in the equation in order for the dual basis functions to have
the correct degree. If φ is the basis function for k-forms, then φ̃ must be the basis
function for (n − k)-forms, and therefore only ⋆φ̃ can be a linear combination of φ.
It turns out that the matrix A is equal to the transpose of H:
AT = H.
To see this,
̃i j = ⟨σi, ⋆φ̃j ⟩
= ∑ A jα ⟨σi, φ̃α ⟩

= ∑ A jα δiα

= A ji .
Thus, Eq. (A.4) becomes AT H = HAT = (−1) k(n−k) Id and since H is of full rank,
this implies that A = (−1) k(n−k) H−T . Since H only depends on the the primal
basis functions φ and the dual simplices ̃
σ , then they completely determine A and
through Eq. (A.5) they completely determine the dual basis functions. In terms of
the discrete Hodge star equation Eq. (A.5) and its dual become
̃ji,
⋆φ̃i = ∑ φ j H

⋆φi = ∑ φ̃j H ji .

130
In summary, in order for the discrete Hodge star to have an exact inverse, the basis
̃ and simplices (σ and ̃
functions (φ and φ)
σ ) must satisfy the following constraints
simultaneously:
⋆ φ̃i = ∑ φ j ⟨σj , ⋆φ̃i ⟩

⋆ φi = ∑ φ̃j ⟨̃
σj , ⋆φi ⟩.

(A.6)

If the basis functions and their corresponding simplices (or cells for a regular mesh)
satisfy the constraints given by Eq. (A.6), the discrete Hodge star will be its own exact
inverse. The basis functions for the periodic grid and the bounded grid described in
this work satisfy this requirement.
A.3

Types of Hodge Star Matrices

The discrete Hodge star operator transforms k-forms on the primal mesh into (n−k)forms on the dual mesh, and vice-versa. It is completely defined by the basis
functions φ and the simplices (or cells) σ, and its matrix elements are given by
̃ij =
⋆φ̃j .
H̄i j =
⋆φ̄
σi

σ̄i

In this section we briefly review the types of Hodge star matrices categorized by
their corresponding basis functions. The reader interested in more information may
refer to, e.g.,[8, 2, 5, 67, 72].
Diagonal
The diagonal Hodge star is commonly used because of its simplicity. Its basis
functions are piecewise constant (i.e. φ ∈ C 0 ) and only have constant non-zero
values over their corresponding simplices. These basis functions are discontinuous
across cell boundaries where their value abruptly changes to zero. The Hodge star,
as implied by its name, is a diagonal matrix whose diagonal elements are given by
H̄iik =

∣̃
σin−k ∣
∣σ̄i k ∣

̃ k = ∣σ̄i ∣ .
ii
∣̃
σik ∣
n−k

From the above equations it is evident that this simple Hodge star is trivially its
own involution (see A.2). Because of the low order of its basis function, it has the
lowest convergence rate, but approaches have been proposed to reduce its error by
optimizing the underlying geometry of the simplices [47].
Whitney
Given piecewise linear basis functions associated with vertices φi , one can define
the whole set of basis functions associated with a simplex of any degree using a

131
trick first introduced by Whitney [73]:
φi0 = φi,
φi1j = φi dφ j − φ j dφi,
φi2j k = 2(φi dφ j ∧ dφ k + φ j dφ k ∧ dφi + φ k dφi ∧ dφ j ).
The basis functions above are associated with the simplices σi0 , σi1j , and σi2j k
respectively, and they result in piecewise linear functions φ ∈ C 1 . The Hodge
star matrix associated with those functions is sparse. In the 1D case it has a very
narrow band along the diagonal. However, this Hodge star matrix is not its own
involution in the exact sense, but only in the approximate sense as the mesh size h
tends to zero, h → 0.
Spectral
The spectral basis functions are smooth φ ∈ C ∞ , and they span the entire computational domain. Therefore, the Hodge star operator is no longer a sparse matrix, but
a dense one. However, it is a special kind of dense matrix. On the periodic grid H
is a circulant matrix, whereas the one on the Chebyshev grid is a centrosymmetric
matrix. For both of these cases the Hodge star is its own exact involution. Below
we explore some other properties of such matrices.
Circulant The N × N circulant matrices are defined by
C = {C ∣ Ci j = c(i− j) mod N }
or in matrix form
cN
⎛ c1
⎜ c
c1
⎜ 2
C=⎜
⎜ ⋮
⎜cN−1 cN−2
⎝ cN cN−1

c3
c4
c1
c2

c2 ⎞
c3 ⎟
⋮ ⎟
⎟.
cN ⎟
c1 ⎠

The set of matrices C forms a group. The proof is left as an exercise to the reader.
Centrosymmetric The centrosymmetric matrices are those N × N matrices that
satisfy the property
H = {H ∣ Hi j = HN+1−i,N+1− j , det H ≠ 0}

132
or in matrix forms
a2
⎛ a1
⎜b
b2
⎜ 1
H=⎜
⎜ ⋮
⎜bN bN−1
⎝aN aN−1
Let

⋯ aN−1 aN ⎞
⋯ bN−1 bN ⎟
⋮ ⎟
⎟.
⋯ b2
b1 ⎟
⋯ a2
a1 ⎠

⎡.
⎢.
J = ⎢⎢
⎢.
⎢1

. . . . . 1⎤⎥
. . . . 1 . ⎥⎥
⎪1 if i = N − j + 1
ij
otherwise
1 . . . . . ⎥⎥
. . . . . . ⎥⎦
The matrix J is sometimes referred to as the exchange matrix.
Using the exchange matrix the above set of matrices can also be conveniently defined
as
H = {A ∣ AJ = J A} .

Proposition. The set of matrices H forms a group.
Proof. It can be easily shown that the identity is in the group, i.e. I ∈ H. Also the
set is closed under multiplication, i.e. if A ∈ H and B ∈ H then AB ∈ H. To see why,
consider
ABJ = AJ B = J AB.
If A ∈ H, then AJ = J A,
(AJ)−1 = (J A)−1,
J −1 A−1 = A−1 J −1,
J A−1 = A−1 J.
Therefore, A−1 ∈ H. This completes the proof that H forms a group. QED
A.4

In Coordinates

In this section we describe some of the exterior calculus operators in coordinates.

133
Wedge Product
Consider the 2D case. The wedge product between a 0-form and a 1-form is
α0 = α,
β1 = β x dx + βy dy,
α0 ∧ β1 = αβ x dx + αβy dy.
Between two 1-forms it is
α1 = αx dx + αy dy,
β1 = β x dx + βy dy,
α1 ∧ β1 = (αx βy − β x αy ) dx ∧ dy.
In the 3D case the wedge product between a 0-form and a 1-form is
α0 = α,
β1 = β x dx + βy dy + βz dz,
α0 ∧ β1 = αβx dx + αβy dy + αβz dz.
Between two 1-forms it is
α1 =αx dx + αy dy + αz dz,
β1 =β x dx + βy dy + βz dz,
α1 ∧ β1 =(αx βy − αy β x )dx ∧ dy+
+(αy βz − αz βy )dy ∧ dz+
+(αz β x − αx βz )dz ∧ dx.
Between a 1-form and a 2-form it is
α1 = αx dx + αy dy + αz dz,
β2 = β x y dx ∧ dy + βyz dy ∧ dz + βzx dz ∧ dx,
α1 ∧ β2 = (αx βyz + αy βzx + αz β xy )dx ∧ dy ∧ dz.
Contraction
Given a manifold M, the interior product is defined as the contraction of a differential
form with a vector field:
⌟ ∶ X × Λ k → Λ k−1 .

134
The contraction of a 1-form α = αi ei with a vector field X = X i ei is given by
X ⌟ α = ⟨X, α⟩
= αi X j ⟨ei, e j ⟩
= αi X j δij
= αi X i .
The contraction of a 2-form β = βi j ei ∧ e j with the same vector field is given by
X ⌟ β = βi j X ⌟ (ei ∧ e j )
= βi j ((X ⌟ ei ) ∧ e j − ei ∧ (X ⌟ e j ))
= βi j (X i e j − X j ei ) .
The contraction of a 3-form γ = γi j k ei ∧ e j ∧ e k is going to be
X ⌟ γ = γi j k X ⌟ (ei ∧ e j ∧ e k )
= γi j k ((X ⌟ ei ) ∧ e j ∧ e k − ei ∧ (X ⌟ e j ) ∧ e k + ei ∧ e j ∧ (X ⌟ e k ))
= γi j k (X i e j ∧ e k − X j ei ∧ e k + X k ei ∧ e j ) .
In 2D the explicit formulas for the contractions are
α = αx dx + αy dy,
X ⌟ α = X x α x + X y α y,
β = β xy dx ∧ dy,
X ⌟ β = −β xy X y dx + β xy X x dy.
If 3D, on the other hand, the contractions become
α = αx dx + αy dy + αz dz
X ⌟ α = X x α x + X y α y + X z α z,
β = β xy dx ∧ dy + βyz dy ∧ dz + βzx dz ∧ dx,
X ⌟ β = (βzx X z − β xy X y ) dx+
+ (β xy X x − βyz X z ) dy+
+ (βyz X y − βzx X x ) dz,
γ = γ xyz dx ∧ dy ∧ dz,
X ⌟ γ = γ xyz (X z dx ∧ dy + X x dy ∧ dz + X y dz ∧ dx) .

135
Lie Derivative
Here we will study the Lie derivative in coordinates. Let f be a scalar (zero-form),
X a vector field, and α be a one-form, where
∂ xi
α = αi ei = αi dx i .

X = X i ei = X i

{e1, . . . , en } = { ∂∂x 1 , . . . , ∂∂x n } form the bases that span the space of vector fields,
and {e1, . . . , en } = {dx 1, . . . , dx n } are the basis elements for the space of one-forms.
Using Cartan’s magic formula the Lie derivative is given algebraically by
L X α = d(X ⌟ α) + X ⌟ dα.

(A.7)

In coordinate form the Lie derivative of a scalar, a one-form, and a vector field
become
∂f
LX f = X ⌟ d f = X i i ,
∂x
j ∂
∂Y ∂
i ∂X
L X Y = [X, Y ] = X i i
∂x ∂xj
∂ xi ∂ x j
To compute the Lie derivative of a one form we use
dα =

∂αi
dx j ∧ dx i,
∂xj

⌟ dα
∂ xk
∂αi ∂
= X k j k ⌟ (dx j ∧ dx i )
∂x ∂x
k ∂αi
(( k ⌟ dx j ) ∧ dx i − (dx j ∧ k ⌟ dx i ))
=X
∂x
∂x
∂x
∂α
= X k j (δkj dx i − δik dx j )
∂x
∂α
∂αi
= X j j dx i − X i j dx j ,
∂x
∂x

X ⌟ dα = X k

∂ xi
= X i α j i ⌟ dx j
∂x
= X α j δij

X ⌟ α = Xi ⌟

= X i αi,

136
d(X ⌟ α) = d (X i αi )
=(

∂ Xi
∂αi
αi + X i j ) dx j .
∂x
∂x

Combining the terms above we obtain the results for the Lie derivative of a 1-form
LX α = X j

∂αi i ∂ X i
dx +
αi dx j ,
∂xj
∂xj

L X α = (X j

∂αi ∂ X j
α j ) dx i .
∂ x j ∂ xi

The Lie derivative is a derivation and so it must satisfy a product rule over the
pairing of forms and vector fields:
L X (Y ⌟ α) = (L X Y ) ⌟ α + Y ⌟ (L X α).

(A.8)

To see that the above expressions for the Lie derivative do indeed satisfy Eq. (A.8):
L X (Y ⌟ α) = L X (αiY i )
= X j j (αiY i )
∂x
∂αi
∂Y i
= X j j Y i + X j αi j ,
∂x
∂x
∂ Xi
j ∂αi
(X
dx
αi dx j )
∂x
∂x
∂x
∂ Xi
∂αi
= X j j Y k ( k ⌟ dx i ) +
αiY k ( k ⌟ dx j )
∂x
∂x
∂x
∂x
∂αi
∂X
= X j j Y k δik +
αiY k δkj
∂x
∂x
∂α
= X j j Yi +
αiY j ,
∂x
∂xj

Y ⌟ L X α = (Y k

(L X Y ) ⌟ α = (X i

j ∂
∂Y j ∂
i ∂X
) ⌟ (αk dx k )
∂ xi ∂ x j
∂ xi ∂ x j

= (αk X i

∂Y j
i ∂X
) ( j ⌟ dx k )
∂x
∂x
∂x

∂Y j
i ∂X
) δ kj
∂ xi
∂ xi
∂Y j
∂X j
= α j X i i − α jY i i .
∂x
∂x

= (αk X i

137
After a simple relabeling of the dummy indices it can be shown that Eq. (A.8) is
indeed satisfied. Thus, in two dimensions it can be shown that the Lie derivative
acts on scalars, vectors, and 1-forms as follows:
∂f
∂f
+ Xy ,
∂x
∂y
L X Y = [X x
+ Xy ,Y x
+Yy ]
∂x
∂y
∂x
∂y
x ∂
x ∂
∂Y
∂Y
∂Y y ∂
∂Y y ∂
= Xx
+ Xy
+ Xx
+ Xy
∂x ∂x
∂y ∂x
∂x ∂y
∂y ∂y
∂Xx ∂
∂Xx ∂
∂Xy ∂
∂Xy ∂
−Yx
−Yy
−Yx
−Yy
∂x ∂x
∂y ∂x
∂x ∂y
∂y ∂y
∂Y
∂X
∂X
∂Y
+ Xy
−Yx
−Yy
) +
= (X x
∂x
∂y
∂x
∂y ∂x
∂Y
∂Y
Xy ∂
+ (X x
+ Xy
−Yx
−Yy
) ,
∂x
∂y
∂x
∂y ∂y
∂αy
∂αy
∂αx
∂αx
LX α = X x
dx + X y
dx + X x
dy + X y
dy+
∂x
∂y
∂x
∂y
∂Xy
∂Xx
∂Xy
∂Xx
αx dx +
αy dx +
αx dy +
αy dy
∂x
∂x
∂y
∂y
∂αx
∂αx ∂ X x
∂Xy
= (X x
+ Xy
αx +
αy ) dx+
∂x
∂y
∂x
∂x
∂αy
∂αy ∂ X x
∂Xy
+ (X x
+ Xy
αx +
αy ) dy.
∂x
∂y
∂y
∂y

LX f = X x

Given a two form
Ω = ω dx ∧ dy,
by Cartan’s magic formula
L X Ω = d(X ⌟ Ω)
= ( (ωX x ) +
(ωX y )) dx ∧ dy.
∂x
∂y

A.5

Variational Derivation of the Equation of Motion of an Ideal Fluid

For ideal fluids, the configuration space is the group of volume-preserving diffeomorphisms DiffVol (M) on the fluid container (in general the region may change
with time as a result of deformations of the underlying manifold). A particle located
at a point x0 ∈ M will travel to a point xt = ϕt (x0 ) ∈ M at time t,
ϕt ∶

M → M.

138
Since there is no potential energy, the Lagrangian is simply the kinetic energy, which
is a mapping from the tangent bundle to the real numbers:
L∶
where

TDiffVol → R,

L(ϕ, ϕ) =

∣∣ϕ ○ ϕ−1 ∣∣2 dV .

Notice that the Lagrangian satisfies the particle relabeling symmetry expressed as
invariance under right composition:
L(ϕ ○ φ, ϕ ○ φ) = L(ϕ, ϕ),
where φ, ϕ ∈ TDiffV ol . The action is given by
ˆ t1
S(ϕ(t), ϕ(t)) =
L dt.
t0

Due to the particle relabeling symmetry the system can be described in terms of a
reduced Lagrangian
l∶

diffV ol → R
ξ ↦ l(ξ) = L(e, ξ) = L(φ, ξ ○ φ),

where ξ ∈ diffV ol is in the Lie algebra of volume-preserving diffeomorphisms, and
e is the identity element of the TDiffV ol group. To obtain the reduced action:
ˆ t1
S(ϕ(t), ϕ(t)) =
L(ϕ(t), ϕ(t)) dt
t0
ˆ t1
L(e, ϕ(t) ○ ϕ(t)−1 ) dt
ˆ 0t1
l(ξ(t)) dt where ξ = ϕ ○ ϕ−1
t0

=∶ s(ξ(t)) Reduced Action Function .
Variations must satisfy the Lin Constraint:
δξ = δ(ϕϕ−1 )
= δϕϕ−1 − ϕϕ−1 δϕϕ−1
= η + [ξ, η],

where η = δϕ ○ ϕ−1 .

139
The variation in the reduced action is
ˆ ˆ
1 t2
v 2 dtdV)
δs = δ (
2 t1 F
ˆ t2 ˆ
⟨v ♭, δv⟩dtdV
ˆ 1t2 ˆ
(⟨v ♭, u⟩ + ⟨v ♭, [v, u]⟩) dt dV,
t1

⟨v ♭, [v, u]⟩ = ⟨v ♭, Lv u⟩ = Lv ⟨v ♭, u⟩ − ⟨Lv v ♭, u⟩.
For a divergence free vector field,
∇⋅v=0
Lv f =

Therefore,
δs =

ˆ t2 ˆ
t1

∂M

(v, n) f = 0.

(⟨v ♭, u⟩ − ⟨Lv v ♭, u⟩) dt dV .

Setting the variation δs = 0 and integrating by parts:
dV ⟨v ♭ + Lv v ♭, u⟩ = 0.

This implies that the integrand must be the gradient of a function
v ♭ + Lv v ♭ + dp = 0,
which is equivalent to the Euler equation in its familiar coordinate form:
v + (v ⋅ ∇)v + ∇p = 0.

140

BIBLIOGRAPHY

[1] Ralph Abraham, Jerrold E. Marsden, and Tudor Ratiu. Manifolds, Tensor
Analysis, and Applications. Applied Mathematical Sciences Vol. 75, Springer,
1988.
[2] D. N. Arnold, R. S. Falk, and R. Winther. “Finite Element Exterior Calculus,
Homological Techniques, and Applications”. In: Acta Numerica 15 (2006),
pp. 1–155.
[3] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. “Finite element
exterior calculus: from Hodge theory to numerical stability”. In: Bull. Amer.
Math. Soc. (N.S.) 47 (2010), pp. 281–354.
[4] Douglas N Arnold et al. Compatible spatial discretizations. Vol. 142. Springer
Science & Business Media, 2007.
[5] Bernhard Auchmann and Stefan Kurz. “A geometrically defined discrete
Hodge operator on simplicial cells”. In: IEEE Transactions on Magnetics
42.4 (2006), pp. 643–646.
[6] S. Behnel et al. “Cython: The Best of Both Worlds”. In: Computing in Science
Engineering 13.2 (2011), pp. 31–39. issn: 1521-9615. doi: 10.1109/MCSE.
2010.118.
[7] Pavel B Bochev and James M Hyman. “Principles of mimetic discretizations
of differential operators”. In: Compatible spatial discretizations. Springer,
2006, pp. 89–119.
[8] A. Bossavit. Computational Electromagnetism. Academic Press, Boston,
1998.
[9] Alain Bossavit. “Generating Whitney forms of polynomial degree one and
higher”. In: IEEE Transactions on Magnetics 38.2 (Mar. 2002), pp. 341–344.
doi: 10.1109/20.996092.
[10] Mick Bouman et al. “A conservative spectral element method for curvilinear domains”. In: Spectral and High Order Methods for Partial Differential
Equations. Springer, 2011, pp. 111–119.
[11] John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation,
2001.
[12] John P. Boyd and Rolfe Petschek. “The Relationships Between Chebyshev,
Legendre and Jacobi Polynomials: The Generic Superiority of Chebyshev
Polynomials and Three Important Exceptions”. In: J. Sci. Comput. 59.1 (Apr.
2014), pp. 1–27.

141
[13] F. Brezzi, A. Buffa, and G. Manzini. “Mimetic scalar products of discrete
differential forms”. In: Journal of Computational Physics 257, Part B (2014),
pp. 1228–1259.
[14] Oscar P Bruno, Youngae Han, and Matthew M Pohlman. “Accurate, highorder representation of complex three-dimensional surfaces via Fourier continuation analysis”. In: Journal of Computational Physics 227.2 (2007),
pp. 1094–1125.
[15] William L. Burke. Applied Differential Geometry. Cambridge University
Press, 1985.
[16] Brian Cabral and Leith Casey Leedom. “Imaging vector fields using line integral convolution”. In: Proceedings of the 20th annual conference on Computer
graphics and interactive techniques. ACM. 1993, pp. 263–270.
[17] C.G. Canuto et al. Spectral Methods: Fundamentals in Single Domains. Scientific Computation Series. Springer, 2006.
[18] Élie Cartan. Les Systémes Differentiels Exterieurs et leurs Applications Géometriques.
Paris: Hermann, 1945.
[19] Keenan Crane, Mathieu Desbrun, and Peter Schröder. “Trivial Connections
on Discrete Surfaces”. In: Computer Graphics Forum 29.5 (2010), pp. 1525–
1533.
[20] M. Desbrun, A. N. Hirani, and J. E. Marsden. “Discrete exterior calculus for
variational problems in computer vision and graphics”. In: Proc. CDC 42
(2003), pp. 533–538.
[21] Mathieu Desbrun, Eva Kanso, and Yiying Tong. “Discrete Differential Forms
for Computational Modeling”. In: Discrete Differential Geometry. Ed. by A.
Bobenko and P. Schröder. Springer, 2007, pp. 287–324.
[22] Mathieu Desbrun et al. “Discrete Exterior Calculus”. In: arXiv:math/0508341v2
(2005).
[23] Sharif Elcott and Peter Schroder. “Building your own DEC at home”. In:
ACM SIGGRAPH 2006 Courses. 2006, pp. 55–59.
[24] Stanley J Farlow. Partial differential equations for scientists and engineers.
Courier Corporation, 1993.
[25] Harley Flanders. Differential Forms with Applications to the Physical Sciences by Harley Flanders. Vol. 11. Elsevier, 1963.
[26] Theodore Frankel. The Geometry of Physics. Second Edition. United Kingdom: Cambridge University Press, 2004.
[27] Matteo Frigo and Steven G Johnson. “FFTW: An adaptive software architecture for the FFT”. In: Acoustics, Speech and Signal Processing, 1998.
Proceedings of the 1998 IEEE International Conference on. Vol. 3. IEEE.
1998, pp. 1381–1384.

142
[28] Marc Gerritsma. “Edge functions for spectral element methods”. In: Spectral
and High Order Methods for Partial Differential Equations. Springer, 2011,
pp. 199–207.
[29] Marc Gerritsma, Jeroen Kunnen, and Boudewijn de Heij. “Discrete Lie
Derivative”. In: Numerical Mathematics and Advanced Applications ENUMATH 2015. Springer. 2016, pp. 635–643.
[30] Fernando de Goes et al. “Subdivision Exterior Calculus for Geometry Processing”. In: ACM Trans. Graph. 35.4 (July 2016), Art. 133.
[31] Leo J. Grady and Jonathan R. Polimeni. Discrete Calculus: Applied Analysis
on Graphs for Computational Science. Springer, 2010.
[32] P. W. Gross and P. R. Kotiuga. “Electromagnetic theory and computation: a
topological approach”. In: Mathematical Sciences Research Institute Publications, Cambridge University Press 48 (2004).
[33] J. Harrison. “Ravello lecture notes on geometric calculus”. In: arXiv:mathph/0501001 (2005).
[34] R.R. Hiemstra et al. “High order geometric methods with exact conservation
properties”. In: Journal of Computational Physics 257 (2014), pp. 1444–
1471.
[35] A. N. Hirani. “Discrete Exterior Calculus”. PhD thesis. Pasadena, CA: California Institute of Technology, 2003.
[36] Michael Holst and Ari Stern. “Geometric variational crimes: Hilbert complexes, finite element exterior calculus, and problems on hypersurfaces”. In:
Found. Comput. Math. 12.3 (2012), pp. 263–293.
[37] John D Hunter. “Matplotlib: A 2D graphics environment”. In: Computing In
Science & Engineering 9.3 (2007), pp. 90–95.
[38] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan. pybind11 — Seamless
operability between C++11 and Python. https://github.com/pybind/pybind11.
2016.
[39] David A. Kopriva. “A Staggered-Grid Multidomain Spectral Method for
the Compressible Navier-Stokes Equations”. In: Journal of Computational
Physics 143.1 (1998), pp. 125–158.
[40] P Robert Kotiuga. “Theoretical limitations of discrete exterior calculus in
the context of computational electromagnetics”. In: IEEE Transactions on
Magnetics 44.6 (2008), pp. 1162–1165.
[41] E. Kreyszig. Differential Geometry. Dover Books on Mathematics. Dover
Publications, 2013.

143
[42] Mehrdad Lakestani and Mehdi Dehghan. “The use of Chebyshev cardinal
functions for the solution of a partial differential equation with an unknown
time-dependent coefficient subject to an extra measurement”. In: Journal of
Computational and Applied Mathematics 235.3 (Dec. 2010), pp. 669–678.
doi: 10.1016/j.cam.2010.06.020.
[43] Melvin Leok. “Foundations of computational geometric mechanics”. PhD
thesis. Pasadena, CA: California Institute of Technology, 2004.
[44] Nam Mai-Duy and Roger I Tanner. “A spectral collocation method based on
integrated Chebyshev polynomials for two-dimensional biharmonic boundaryvalue problems”. In: Journal of Computational and Applied Mathematics
201.1 (2007), pp. 30–47.
[45] T. Mathew. Domain Decomposition Methods for the Numerical Solution of
Partial Differential Equations. Lecture Notes in Computational Science and
Engineering. Springer Berlin Heidelberg, 2008. isbn: 9783540772057.
[46] Aaron Meurer et al. “SymPy: symbolic computing in Python”. In: PeerJ
Computer Science 3 (2017), e103.
[47] Patrick Mullen et al. “HOT: Hodge Optimized Triangulations”. In: ACM
Transactions on Graphics 30.3 (July 2011), Art. 103.
[48] P. Mullen et al. “Discrete Lie Advection of Differential Forms”. In: Foundations of Computational Mathematics 11.2 (Sept. 2010), pp. 131–149. doi:
10.1007/s10208-010-9076-y.
[49] James R. Munkres. Elements of Algebraic Topology. Menlo Park, CA: AddisonWesley, 1984.
[50] M. Nakahara. Geometry, Topology and Physics, Second Edition. Graduate
student series in physics. Taylor & Francis, 2003. isbn: 9780750306065.
[51] J.-C. Nédélec. “Mixed finite elements in R3 ”. In: Numer. Math. 35 (1980),
pp. 315–341.
[52] Roy A. Nicolaides and Xiaonan Wu. “Covolume Solutions of Three Dimensional Div-Curl Equations”. In: SIAM J. Numer. Anal. 34 (1997), p. 2195.
[53] J. O’Rourke. Computational geometry in C (second edition). Cambridge
University Press, 1998.
[54] Artur Palha and Marc Gerritsma. “Mimetic spectral element method for
Hamiltonian systems”. In: arXiv:1505.03422 (2015).
[55] Artur Palha and Marc Gerritsma. “Spectral Element Approximation of the
Hodge-⋆ Operator in Curved Elements”. In: Spectral and High Order Methods
for Partial Differential Equations. Ed. by Jan S. Hesthaven et al. Vol. 76.
Lecture Notes in Computational Science and Engineering. Springer Berlin
Heidelberg, 2011, pp. 283–291.

144
[56] Artur Palha and Marc Gerritsma. “Spectral element approximation of the
Hodge-⋆ operator in curved elements”. In: Spectral and High Order Methods
for Partial Differential Equations. Springer, 2011, pp. 283–291.
[57] Artur Palha et al. “Physics-compatible discretization techniques on single
and dual grids, with application to the Poisson equation of volume forms”.
In: Journal of Computational Physics 257 (Jan. 2014), pp. 1394–1422. doi:
10.1016/j.jcp.2013.08.005.
[58] Dmitry Pavlov et al. “Structure-preserving discretization of incompressible
fluids”. In: Physica D: Nonlinear Phenomena 240.6 (2011), pp. 443–458.
[59] Nicolas Robidoux. Polynomial Histopolation, Superconvergent Degrees Of
Freedom, And Pseudospectral Discrete Hodge Operators. Preprint. 2006.
[60] Nicolas Robidoux and Stanly Steinberg. “A discrete vector calculus in tensor
grids”. In: Comput. Methods Appl. Math. 1.2001 (2011), pp. 1–44.
[61] Nicolas P Rougier. “Glumpy”. In: EuroScipy. 2015.
[62] Dzhelil Rufat et al. “The chain collocation method: A spectrally accurate
calculus of forms”. In: Journal of Computational Physics 257 (Jan. 2014),
pp. 1352–1372. doi: 10.1016/j.jcp.2013.08.011.
[63] Jonathan Richard Shewchuk. “Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator”. In: Applied computational geometry towards
geometric engineering. Springer, 1996, pp. 203–222.
[64] Ari Stern et al. “Geometric Computational Electrodynamics with Variational
Integrators and Discrete Differential Forms”. In: Geometry, Mechanics, and
Dynamics: The Legacy of Jerry Marsden. Ed. by Dong Eui Chang et al.
Springer New York, 2015, pp. 437–475.
[65] Ari Stern et al. “Variational integrators for Maxwell’s equations with sources”.
In: PIERS Online 4.7 (2008), pp. 711–715.
[66] Kunihiko Taira and Tim Colonius. “The immersed boundary method: a
projection approach”. In: Journal of Computational Physics 225.2 (2007),
pp. 2118–2137.
[67] Timo Tarhasaari, Lauri Kettunen, and Alain Bossavit. “Some realizations of
a discrete Hodge operator: a reinterpretation of finite element techniques”.
In: IEEE Transactions on magnetics 35.3 (1999), pp. 1494–1497. doi: 10.
1109/20.767250.
[68] Yiying Tong et al. “Designing quadrangulations with discrete harmonic
forms”. In: Symposium on Geometry Processing. 2006, pp. 201–210.
[69] Lloyd N. Trefethen. Spectral Methods in MATLAB. Philadelphia: SIAM,
2000.
[70] Jonathan Turner and Jason Turner. ChaiScript: Easy to use scripting for C++.

145
[71] Stéfan van der Walt, S Chris Colbert, and Gael Varoquaux. “The NumPy
array: a structure for efficient numerical computation”. In: Computing in
Science & Engineering 13.2 (2011), pp. 22–30.
[72] Ke Wang et al. “Edge subdivision schemes and the construction of smooth
vector fields”. In: ACM Transactions on Graphics. Vol. 25(3). 2006, pp. 1041–
1048.
[73] H. Whitney. Geometric integration theory. Princeton University Press, 1957.