PyDEC: Software and Algorithms for Discretization of Exterior Calculus

PyDEC: Software and Algorithms for
Discretization of Exterior Calculus

Nathan Bell NVIDIA Corporation, Anil N. Hirani Department of Computer Science, University of Illinois at Urbana-Champaign,

This paper describes the algorithms, features and implementation of PyDEC, a Python library for computations related to the discretization of exterior calculus. PyDEC facilitates inquiry into both physical problems on manifolds as well as purely topological problems on abstract complexes. We describe efficient algorithms for constructing the operators and objects that arise in discrete exterior calculus, lowest order finite element exterior calculus and in related topological problems. Our algorithms are formulated in terms of high-level matrix operations which extend to arbitrary dimension. As a result, our implementations map well to the facilities of numerical libraries such as NumPy and SciPy. The availability of such libraries makes Python suitable for prototyping numerical methods. We demonstrate how PyDEC is used to solve physical and topological problems through several concise examples.

Categories and Subject Descriptors:

G.4 [Mathematical Software]; G.1.8 [Numerical Analysis]: Partial Differential Equations – Finite element methods, Finite volume methods, Discrete exterior calculus, Finite element exterior calculus; I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling – Geometric algorithms, languages, and systems, Computational topology

1 Introduction

Geometry and topology play an increasing role in the modern language used to describe physical problems [1, 25]. A large part of this language is exterior calculus which generalizes vector calculus to smooth manifolds of arbitrary dimensions. The main objects are differential forms (which are anti-symmetric tensor fields), general tensor fields, and vector fields defined on a manifold. In addition to physical applications, differential forms are also used in cohomology theory in topology [14].

Once the domain of interest is discretized, it may not be smooth and so the objects and operators of exterior calculus have to be reinterpreted in this context. For example, a surface in may be discretized as a two dimensional simplicial complex embedded in , i.e., as a triangle mesh. Even when the domain is a simple domain in space, such as an open subset of the plane or space the discretization is usually in the form of some mesh. The various objects of the problem then become defined in a piecewise varying fashion over such a mesh and so a discrete calculus is required there as well. After discretizing the domains, objects, and operators, one can compute numerical solutions of partial differential equations (PDEs), and compute some topological invariants using the same discretizations. Both these classes of applications are considered here.

There have been several recent developments in the discretization of exterior calculus and in the clarification of the role of algebraic topology in computations. These go by various names, such as covolume methods [46], support operator methods [50], mimetic discretization [36, 10], discrete exterior calculus (DEC) [31, 19, 27], compatible discretization, finite element exterior calculus [2, 4, 35], edge and face elements or Whitney forms [13, 30], and so on. PyDEC provides an implementation of discrete exterior calculus and lowest order finite element exterior calculus using Whitney forms.

Within pure mathematics itself, ideas for discretizing exterior calculus have a long history. For example the de Rham map that is commonly used for discretizing differential forms goes back at least to [16]. The reverse operation of interpolating discrete differential forms via the Whitney map appears in [52]. A combinatorial (discrete) Hodge decomposition theorem was proved in [21] and the idea of a combinatorial Hodge decomposition dates to [23]. More recent work on discretization of Hodge star and wedge product is in [53, 54]. Discretizations on other types of complexes have been developed as well [28, 49].

1.1 Main contributions

In this paper we describe the algorithms and design of PyDEC, a Python software library implementing various complexes and operators for discretization of exterior calculus, and the algorithms and data structures for those. In PyDEC all the discrete operators are implemented as sparse matrices and we often reduce algorithms to a sequence of standard high-level operations, such as sparse matrix-matrix multiplication [6], as opposed to more specialized techniques and ad hoc data structures. Since these high-level operations are ultimately carried out by efficient, natively-compiled routines (e.g. C or Fortran implementations) the need for further algorithmic optimization is generally unnecessary.

As is commonly done, in PyDEC we implement discrete differential forms as real valued cochains which will be defined in Section 2. PyDEC has been used in a thesis [8], in classes taught at University of Illinois, in experimental parts of some computational topology papers [20, 22, 32], in Darcy flow [34], and in least squares ranking on graphs [33]. The PyDEC source code and examples are publicly available [7]. We summarize here our contributions grouped into four areas.

Basic objects and methods: (1) Data structures for : simplicial complexes of dimension embedded in , ; abstract simplicial complexes; Vietoris-Rips complexes for points in any dimension; and regular cube complexes of dimension embedded in ; (2) Cochain objects for the above complexes; (3) Discrete exterior derivative as a coboundary operator, implemented as a method for cochains on various complexes.

Finite element exterior calculus: (1) Fast algorithm to construct sparse mass matrices for Whitney forms by eliminating repeated computations; and (2) Assembly of stiffness matrices for Whitney forms from mass matrices by using products of boundary and stiffness matrices. Note that only the lowest order () elements of finite element exterior calculus are implemented in PyDEC.

Discrete exterior calculus: (1) Diagonal sparse matrix discrete Hodge star for well-centered (circumcenters inside simplices) and Delaunay simplicial complexes (with an additional boundary condition); (2) Circumcenter calculation for -simplex in an -dimensional simplicial complex embedded in using a linear system in barycentric coordinates; and (3) Volume calculations for primal simplices and circumcentric dual cells.

Examples: (1) Resonant cavity curl-curl problem; (2) Flow in porous medium modeled as Darcy flow, i.e., Poisson’s equation in first order (mixed) form; (3) Cohomology basis calculation for a simplicial mesh, using harmonic cochain computation using Hodge decomposition; (4) Finding sensor network coverage holes by modeling an abstract, idealized sensor network as a Rips complex; and (5) Least squares ranking on graphs using Hodge decomposition of partial pairwise comparison data.

2 Overview of PyDEC

One common type of discrete domain used in scientific computing is triangle or tetrahedral mesh. These and their higher dimensional analogues are implemented as -dimensional simplicial complexes embedded in , . Simplicial complexes are useful even without an embedding and even when they don’t represent a manifold, for example in topology and ranking problems. Such abstract simplicial complexes without any embedding for vertices are also implemented in PyDEC. The other complexes implemented are regular cube complexes and Rips complexes. Regular cubical meshes are useful since it is easy to construct domains even in high dimensions whereas simplicial meshing is hard enough in 3 dimensions and rarely done in 4 or larger dimensions. Rips complexes are useful in applications such as topological calculations of sensor network coverage analysis [17]. The representations used for these four types of complexes are described in Section 3-6. A complex that is a manifold (i.e., locally Euclidean) will be referred to as mesh.

The definitions here are given for simplicial complexes and generalize to the other types of complexes implemented in PyDEC. In PyDEC we only consider integer valued chains and real-valued cochains. Also, we are only interested in finite complexes, that is, ones with a finite number of cells. Let be a finite simplicial complex and denote its underlying space by . Give the subspace topology as a subspace of (a set in is open iff is open in ). For a finite complex this is the same as the standard way of defining topology for [45, pages 8-9] and is a closed subspace of .

An oriented simplex with vertices will be written as and given names like with the superscript denoting the dimension and subscript denoting its place in some ordering of -simplices. Sometimes the dimensional superscript and/or the indexing subscript will be dropped. The orientation of a simplex is one of two equivalence classes of vertex orderings. Two orderings are equivalent if one is an even permutation of the other. For example and denote the same oriented triangle while is the oppositely oriented one.

A -chain of is a function from oriented -simplices of to the set of integers , such that where is the simplex oriented in the opposite way. Two chains are added by adding their values. Thus -chains are formal linear combinations (with integer coefficients) of oriented -dimensional simplices. The space of -chains is denoted and it is a free abelian group. See [45, page 21]. Free abelian groups have a basis and one does not need to impose a vector space structure. For example, a basis for is the set of integer valued functions that are 1 on a -simplex and 0 on the rest, with one such basis element corresponding to each -simplex. These are called elementary chains and the one corresponding to a -simplex will also be referred to as . The existence of this basis and the addition and negation of chains is the only aspect that is important for this paper. The intuitive way to think of chains is that they play a role similar to that played by the domains of integration in the smooth theory. The negative sign allows one to talk about orientation reversal and the integer coefficient allows one to say how many times integration is to be done on that domain.

Sometimes we will need to refer to a dual mesh which will in general be a cell complex obtained from a subdivision of the given complex . We’ll refer to the dual complex as . For a discrete Hodge star diagonal matrix of DEC, the dual mesh is the one obtained from circumcentric subdivision of a well-centered or Delaunay simplicial complex and such a Hodge star is described in Section 10.

Homomorphisms from the -chain group to are called -cochains of and denoted . This set is an abelian group and also a vector space over . Similarly the dual -cochains are denoted or . The discretization map from space of smooth -forms to -cochains is called the de Rham map or . See [16, 21]. For a smooth -form , the de Rham map is defined as for any chain . We will denote the evaluation of the cochain on a chain as . A basis for is the set of elementary cochains. The elementary cochain is the one that takes value 1 on elementary chain and 0 on the other elementary chains. Thus the vector space dimension of is the number of -simplices in . We’ll denote this number by . Thus will be the number of vertices, the number of edges, the number of triangles and so on.

Like most of the numerical analysis literature mentioned in Section 1 we assume that the smooth forms are either defined in the embedding space of the simplicial complex, or on the complex itself, or can be obtained by pullback from the manifold that the complex approximates. In contrast, most mathematics literature quoted including [52, 21] uses simplicial complex defined on the smooth manifold as a “curvilinear” triangulation. In the applied literature, the complex approximates the manifold. Many finite element papers deal with open subsets of the plane or so they are working with triangulations of a manifold with piecewise smooth boundaries. Surface finite element methods have been studied outside of exterior calculus [18]. A variational crimes methodology is used for finite element exterior calculus on simplicial approximations of manifolds in [35]. In the computer graphics literature, piecewise-linear triangle mesh surfaces embedded in are common and convergence questions for operators on such surfaces have been studied [29]. In light of all of these, PyDEC’s framework of using simplicial or other approximations of manifolds is appropriate.

Operators such as the discrete exterior derivative () and Hodge star () can be implemented as sparse matrices. At each dimension, the exterior derivative can be easily determined by the incidence structure of the given simplicial mesh. For DEC the Hodge star is a diagonal matrix whose entries are determined by the ratios of primal and dual volumes. Care is needed for dual volume calculation when the mesh is not well-centered. For finite element exterior calculus we implement Whitney forms. The corresponding Hodge star is the mass matrix which is sparse but not diagonal. One of the stiffness matrices can be obtained from it by combining it with the exterior derivative.

Once the matrices implementing the basic operators have been determined, they can be composed together to obtain other operators such as the codifferential () and Laplace-deRham (). While this composition could be performed manually, i.e. the appropriate set of matrices combined to form the desired operation, it is prone to error. In PyDEC this composition is handled automatically. For example, the function which implements the exterior derivative, looks at the dimension of its argument to determine the appropriate matrix to apply. The same method can be applied to the codifferential function , which then makes their composition work automatically. This automation eliminates a common source of error and makes explicit which operators are being used throughout the program.

PyDEC is intended to be fast, flexible, and robust. As an interpreted language, Python by itself is not well-suited for high-performance computing. However, combined with numerical libraries such as NumPy and SciPy one can achieve fast execution with only a small amount of overhead. The NumPy array structure, used extensively throughout SciPy and PyDEC, provides an economical way of storing N-dimensional arrays (comparable to C or Fortran) and exposes a C API for interfacing Python with other, potentially lower-level libraries [51]. In this way, Python can be used to efficiently “glue” different highly-optimized libraries together with greater ease than a purely C, C++, or Fortran implementation would permit [47]. Indeed, PyDEC also makes extensive use of the sparse module in SciPy which relies on natively-compiled C++ routines for all performance-sensitive operations, such as sparse matrix-vector and matrix-matrix multiplication. PyDEC is therefore scalable to large data sets and capable of solving problems with millions of elements [8].

Even large-scale, high-performance libraries such as Trilinos provide Python bindings showing that Python is useful beyond the prototyping stage. We also make extensive use of Python’s built-in unit testing framework to ensure PyDEC’s robustness. For each non-trivial component of PyDEC, a number of examples with known results are used to check for consistency.

2.1 Previous work

Discrete differential forms now appear in several finite element packages such as FEMSTER [15], DOLFIN [43] and deal.II [5]. These libraries support arbitrary order conforming finite element spaces in two and three dimensions. In contrast, for finite elements PyDEC supports simplicial and cubical meshes of arbitrary dimension, albeit with lowest order elements. In addition, PyDEC also supports the operators of discrete exterior calculus and complexes needed in topology. We note that Exterior [41], an experimental library within the FEniCS [42] project, realizes the framework developed by Arnold et al. [3] which generalizes to arbitrary order and dimension. Exterior uses symbolic methods and supports integration of forms on the standard simplex. PyDEC supports mass and stiffness matrices on simplicial and cubical complexes. The discovery of lower dimensional faces in a complex and the computation of all the boundary matrices is also implemented in PyDEC.

The other domain where PyDEC is useful is in computational topology. There are several packages in this domain as well, and again PyDEC has a different set of features and aims from these. In [39] efficient techniques are developed for finding meaningful topological structures in cubical complexes, such as digital images. In addition to simplicial and cubical manifolds, PyDEC also provides support for abstract simplicial complexes such as the Rips complex of a point set. The Applied and Computational Topology group at Stanford University has been the source for several packages for computational topology. These include various versions of PLEX such as JPlex and javaPlex which are designed for persistent homology calculations. Another package from the group is Dionysus, a C++ library that implements persistent homology and cohomology [24, 55] and other interesting topological and geometric algorithms. In contrast, we view the role of PyDEC in computational topology as providing a tool to specify and represent different types of complexes, compute their boundary matrices, and compute cohomology representatives with or without geometric information.

3 Simplicial Complex Representation

Before detailing the algorithms used to implement discretizations of exterior calculus, we discuss the representation of various complexes, starting in this section with simplicial complexes. Consider the triangle mesh shown in Figure 1 with vertices and faces enumerated as shown. This example mesh is represented by arrays

where the subscript denotes the dimension of the simplices. The -th row of contains the spatial coordinates of the -th vertex. Likewise the -th row of simplex array contains the indicies of the vertices that form the -th triangle. The indices of each simplex in in this example are ordered in a manner that implies a counter-clockwise orientation for each. For an -dimensional discrete manifold, or mesh, arrays and suffice to describe the computational domain.

Figure 1: Simplicial mesh with enumerated vertices and simplices.

In addition to and , an -dimensional simplicial complex is comprised by its -dimensional faces, . In the case of Figure 1, these are

which correspond to the vertices (0-simplices) and oriented edges (1-simplices) of the complex. A graphical representation of this simplicial complex is shown in Figure 2. Since the orientation of the lower () dimensional faces is arbitrary, we use the convention that face indices will be in sorted order. Furthermore, we require the rows of to be sorted in lexicographical order. As pointed out in Section 7 and 9, these conventions facilitate efficient construction of differential operators and stiffness matrices.

Figure 2: Simplicial complex with oriented edges and triangles.

4 Regular Cube Complex Representation

PyDEC provides a regular cube complex of dimension embedded in for any . As mentioned earlier, in dimension higher than 3, constructing simplicial manifold complexes is hard. In fact, even construction of good tetrahedral meshes is still an active area in computational geometry. This is one reason for using regular cube complexes in high dimensions. Moreover, for some applications, like topological image analysis or analysis of voxel data, the regular cube complex is a very convenient framework [39].

A regular cube complex can be easily specified by an -dimensional array of binary values (bitmap) and a regular -dimensional cube is placed where the bit is on. For example the cube complex shown in Figure 3 can be created by specifying the bitmap array

A bitmap suffices to describe the top level cubes, but a cube array (like simplex array) is used during construction of differential operators and for computing faces. In this paper we describe the construction of exterior derivative, Hodge star and Whitney forms on simplicial complexes. For cube complexes we describe only the construction of exterior derivative and lower dimensional faces. However, the other operators and objects are also implemented in PyDEC for such complexes. For example, Whitney-like elements for hexahedral grids are described in [11] and are implemented in PyDEC.

Figure 3: Regular cube mesh with enumerated vertices and faces.
Figure 4: Regular cube complex with oriented edges and faces.

Converting a bitmap representation of a mesh into a cube array representation is straightforward. For example, the cube array representation of the mesh in Figure 3 is

As with the simplex arrays, the rows of correspond to individual two-dimensional cubes in the mesh. The two left-most columns of encode the origins of each two-dimensional cube, namely , and . The remaining two columns encode the coordinate directions spanned by the cube. Since represents the top-level elements, all cubes span both the (coordinate ) and (coordinate ) dimensions. In general, the first columns of encode the origin or corner of a cube while the remaining columns identify the coordinate directions swept out by the cube. We note that the cube array representation is similar to the cubical notation used by Sen [49].

The edges of the mesh in Figure 3 are represented by the cube array

where again the first two columns encode the origin of each edge and the last column indicates whether the edge points in the or direction. For example, the row corresponds to edge in Figure 4 which begins at and extends one unit in the direction. Similarly the row encodes an edge starting at extending one unit in the direction. Since zero-dimensional cubes (points) have no spatial extent their cube array representation

contains only their coordinate locations.

The cube array provides a convenient representation for regular cube complexes. While a bitmap representation of the top-level cubes is generally more compact, the cube array representation generalizes naturally to lower-dimensional faces and is straightforward to manipulate.

5 Rips Complex Representation

The Rips complex, or Vietoris-Rips complex of a point set is defined by forming a simplex for every subset of vertices with diameter less than or equal to a given distance . For example, if pair of vertices are no more than distance apart, then the Rips complex contains an edge (1-simplex) between the vertices. In general, a set of vertices forms a -simplex when all pairs of vertices in the set are separated by at most .

Figure 5: Broadcast radii and Rips complex for a sensor network.

In recent work, certain sensor network coverage problems have been shown to reduce to finding topological properties of the network’s Rips complex, at least for an abstract model of sensor networks [17]. Such coordinate-free methods rely only on pairwise communication between nodes and do not require the use of positioning devices. These traits are especially important in the context of ad-hoc wireless networks with limited per-node resources. Figure 5 depicts a planar sensor network and its associated Rips complex.

In this section we describe an efficient method for computing the Rips complex for a set of points. Although we consider only the case of points embedded in Euclidean space, our methodology applies to more general metric spaces. Indeed, only the construction of the -skeleton of the Rips complex requires metric information. The higher-dimensional simplices are constructed directly from the -skeleton.

We compute the -skeleton of the Rips complex with a kD-Tree data structure. Specifically, for each vertex we compute the set of neighboring vertices . The hierarchical structure of the kD-Tree allows such queries to be computed efficiently.

The -skeleton of the Rips complex is stored in an array , using the convention discussed in Section 3. Additionally, the (oriented) edges of the -skeleton are used to define , a directed graph stored in a sparse matrix format. Specifically, if is an edge of the Rips complex, and zero otherwise. For the Rips complex depicted in Figure 6,

are the corresponding simplex array and directed graph respectively.

Figure 6: Five directed edges form the -skeleton of the Rips complex.

The arrays of higher dimensional simplices can be computed as follows. Let denote the (sparse) matrix whose rows are identified with the -simplices as specified by . Each row of encodes the vertices which form the corresponding simplex. Specifically, takes the value if the -th simplex contains vertex and zero otherwise. For the example shown in Figure 6,

encodes the edges stored in . Once is constructed we compute the sparse matrix-matrix product . For our example the result is

Like , the product is a matrix that relates the -simplices to the vertices: the matrix entry of counts the number of directed edges that exist from the vertices of simplex to vertex . When the value of entry of is equal to , we form a -simplex of the Rips complex by concatenating simplex with vertex . In the example, matrix entries and of are equal to which implies that the -skeleton of the Rips complex contains two simplices, formed by appending vertex to the -simplices and , or

in array format. This process may be applied recursively to develop higher dimensional simplices as required by the application. Thus our algorithm computes simplices of the Rips complex with a handful of sparse and dense matrix operations.

6 Abstract Simplicial Complex Representation

In Section 5 we saw an example of a simplicial complex which was not a manifold complex (Figure 5). Rips complexes described in Section 5 demonstrate one way to construct such complexes in PyDEC, starting from locations of vertices. There are other applications, for example in topology, where we would like to create a simplicial complex which is not necessarily a manifold. In addition we would like to do this without requiring that the location of vertices be given. For example, in topology, surfaces are often represented as a polygon with certain sides identified. One way to describe such an object is as an abstract simplicial complex [45, Section 3]. This is a collection of finite nonempty sets such that if a set is in the collection, so are all the nonempty subsets of it. Figure 7 shows two examples of abstract simplicial complexes created in PyDEC.

Figure 7: Examples of abstract simplicial complexes. The one of the left represents the triangulation of a Möbius strip and the one on the right that of a projective plane.

In PyDEC, abstract simplicial complexes are created by specifying a list of arrays. Each array contains simplices of a single dimension, specified as an array of vertex numbers. Lower dimensional faces of a simplex need not be specified explicitly. For example the Möbius strip triangulation shown in Figure 7 can be created by giving the array

as input to PyDEC. Abstract simplicial complexes need not be a triangulation of a manifold. For example one consisting of 2 triangles with an extra edge attached and a standalone vertex may be created using a list consisting of the arrays

as input.

The boundary matrices of a simplicial complex encode the connectivity information and can be computed from a purely combinatorial description of a simplicial complex. The locations of the vertices are not required. Thus the abstract simplicial complex structure is all that is required to compute these matrices as will be described in the next section.

7 Discrete Exterior Derivative

Given a manifold , the exterior derivative which acts on differential -forms, generalizes the derivative operator of calculus. When combined with metric dependent operators Hodge star, sharp, and flat appropriately, the vector calculus operators div, grad and curl can be generated from . But itself is a metric independent operators whose definition does not require any Riemannian metric on the manifold. See [1] for details. The discrete exterior derivative (which we will also denote as ) in PyDEC is defined as is usual in the literature, as the coboundary operator of algebraic topology [45]. Thus

for arbitrary -cochain and -chain . Recall that the boundary operator on cochains, is defined by extension of its definition on an oriented simplex. The boundary operator on a -simplex is given in terms of its -dimensional faces ( in number) as


where means that is omitted. Therefore, given an -dimensional simplicial complex represented by , for the discrete exterior derivative, it suffices to compute . As is usual in algebraic topology, in PyDEC we compute matrix representations of these in the elementary chain basis. Boundary matrices are useful in finite elements since their transposes are the coboundary operators. They are also useful in computational topology since homology and cohomology groups are the quotient groups of kernel and image of boundary matrices [45].

For the complex pictured in Figure 2 the boundary operators are

In the following we describe an algorithm that takes as input and computes both and . This procedure is applied recursively to produce all faces of the complex and their boundary operator at that dimension.

The first step of the algorithm converts a simplex array into a canonical format. In the canonical format each simplex (row of ) is replaced by the simplex with sorted indices and the relative parity of the original ordering to the sorted order. For instance, the simplex becomes since an odd number of transpositions (namely one) are needed to transform into . Similarly, the canonical format for simplex is since an even number of transpositions are required to sort the simplex indices. Since the complex dimension is typically small (i.e. ), a simple sorting algorithm such as insertion sort is employed at this stage. We denote the aforementioned process canonical_format where the rightmost column of contains the simplex parity. Applying canonical_format to in our example yields

Once a simplex array has been transformed into canonical format, the -dimensional faces and boundary operator are readily obtained. We denote this process

In order to establish the correspondence between the -dimensional simplices and their faces, we first enumerate the simplices by adding another column to to form . For example,

The formula (7.1) is applied to in a columnwise fashion by excluding the -th column of simplex indices, multiplying the parity column by , and carrying the last column over unchanged. For example,

The resultant array is then sorted by the first columns in lexicographical order, allowing the unique faces to then be extracted.

Furthermore, a Compressed Sparse Row (CSR) [48] sparse matrix representation of as

is obtained from the sorted matrix. For these are


where indices and data correspond to the fourth and third rows of the sorted matrix.

This process is then applied to and so on down the dimension. Since the lower dimensional simplex array rows are already sorted, those arrays are already in canonical format. Thus a single algorithm generates the lower dimensional faces as well as the boundary matrices at all the dimensions. The boundary matrices, and hence the coboundary operators, are generated in a convenient sparse matrix format.

7.1 Generalization to Abstract Complexes

The boundary operators and faces of an abstract simplicial complex are computed with a straightforward extension of the boundary_faces algorithm. Recall from Section 6 that an abstract simplicial complex is specified by a list of simplex arrays of different dimensions, where the lower-dimensional simplex arrays represent simplices that are not a face of any higher-dimensional simplex. Generalizing the previous scheme to the case of abstract simplicial complexes is accomplished by (1) augmenting the set of computed faces with the user-specified simplices and (2) modifying the computed boundary operator accordingly.

Consider the abstract simplicial complex represented by the simplex arrays

which consists of two triangles, an edge, and an isolated vertex. Applying boundary_faces to produces an array of face edges and corresponding boundary operator

which includes all but the user-specified edge . User-specified simplices are then incorporated into the simplex array in a three-stage process: (1) user-specified simplices are concatenated to the computed face array; (2) the rows of the combined simplex array are sorted lexicographically; (3) redundant simplices (if any) are removed from the sorted array. Upon completion, the augmented simplex array contains the union of the face simplices and the user-specified simplices. Continuing the example, the edge is incorporated into as follows

In the final stage of the procedure, the computed boundary operator ( in the example) is updated to reflect the newly incorporated simplices. Since the new simplices do not lie in the boundary of any higher-dimensional simplex, we may simply add empty rows into the sparse matrix representation of the boundary operator for each newly added simplex. Therefore, the boundary operator update procedure amounts to a simple remapping of row indices. In the example, the addition of the edge into the fifth row of the simplex array requires the addition of an empty row into the boundary operator at the corresponding position,

The Rips complex of Section 5 does have the location information for the vertices. However, ignoring those, such a complex is an abstract simplicial complex. Thus the boundary matrices for a Rips complex can be computed as described above. In practice some efficiency can be obtained by ordering the computation differently, so that the matrices are built as the complex is being built from the edge skeleton of the Rips complex. That is how it is implemented in PyDEC.

7.2 Boundary operators and faces for cubical complexes

The algorithm used to compute the faces and boundary operator of a given cube array () is closely related to the procedure discussed in Section 7 for simplex arrays. Consider a general -cube in -dimensions, denoted by the pair where are the coordinates of cube’s origin and are the directions which the -cube spans. Note that the values correspond exactly to a row of the cube array representation introduced in Section 4. Using this notation, the boundary of a -cube is given by the expression


where denotes the omission of the -th spanning direction and is the corresponding coordinate. For example, the boundary of a square centered at the location is


The canonical format for a -cube is the one where the spanning directions are specified in ascending order. For instance, the -cube is in the canonical format because . As with simplices, each cube has a unique canonical format, through which duplicates are easily identified. Since the top-level cube array is generated from a bitmap it is already in the canonical format and no reordering of indices or parity tracking is necessary.

Applying Equation 7.2 to a -cube array with members generates a collection oriented faces. In the mesh illustrated in Figure 3 the three squares in are initially expanded into

where the fourth column of encodes the orientation of the face and the fifth column records the -cube to which each face belongs. Sorting the rows of in lexicographical order

allows the unique faces to be extracted. Lastly, a sparse matrix representation of the boundary operator is obtained from the sorted cube array in the same manner as for simplices.

8 Review of Whitney Map and Whitney Forms

In this section we review and collect some material, most of which is well-known in DEC and finite element exterior calculus research communities. It is included here partly to fix notation. In this section we also give the monomials based definition of inner product of differential forms. This is not the way inner product of forms is usually defined in most textbooks, [44] being one exception we know of. The monomial form leads to an efficient algorithm for computation of stiffness and mass matrices for Whitney forms given in Section 9.

The basic function spaces that are useful with exterior calculus are the space of square integrable -forms on a manifold and Sobolev spaces derived from that. Let be a Riemannian manifold, a manifold on which a smoothly varying inner product is defined on the tangent space at each point. Let be its metric, a smooth tensor field that defines the inner product on the tangent space at each point on .

For differential forms on such a manifold , the space of square integrable forms is denoted . One can then define the spaces which generalize the spaces and used in mixed finite element methods [4]. To define one has to define an inner product on the space of forms which is our starting point for this section. All these function spaces have been discussed in [2] and [4]. The definitions and properties of Whitney map and Whitney forms is in [21], the geometric analysis background is in [38] and the basic definition of inner products on forms is in [1, page 411].

8.1 Inner product of forms

To define the spaces and more precisely, we recall the definitions related to inner products of forms. We will need the exterior calculus operators wedge product and Hodge star which we recall first. For a manifold the wedge product is an operator for building -forms from -forms and -forms. It is defined as the skew-symmetrization of the tensor product of the two forms involved. For a Riemannian manifold of dimension , the Hodge star operator is an isomorphism between the spaces of and -forms. For more details, see [1, page 394] for wedge products and [1, page 411] for Hodge star.

Definition 8.1.

Given two smooth -forms , on a Riemannian manifold , their pointwise inner product at point is defined by


where is the volume form associated with the metric induced by the inner product on .

The pointwise inner products of forms can be defined in another way, which will be more useful to us in computations. The second definition given below in Definition 8.2 is equivalent to the one given above in Definition 8.1. The operator (the sharp operator) used below is an isomorphism between 1-forms and vector fields and is defined by for given 1-form and all vector fields . See [1] for details.

Definition 8.2.

Let and be 1-forms on a Riemannian manifold . By analogy with polynomials we’ll call -forms of the type and monomial -forms. Define the following operator at a point :


where is the matrix obtained by taking . In the equation above, all the 1-forms are evaluated at the point . Extend the operation in (8.2) bilinearly pointwise to the space of all -forms. It can be shown that this defines a pointwise inner product of -forms equivalent to the one defined in (8.1). Note that if for all , the expression on the right in (8.2) is the Gram determinant.

Remark 8.3.

Note that unlike (8.1) the definition in (8.2) does not involve wedge product and Hodge star explicitly. This is an advantage of the latter form since a discrete wedge product is not available in PyDEC. The RHS of (8.2) does involve the sharp operator, but as we will see in the next section, this is easy to interpret for the purpose of discretization in this context.

Definition 8.4.

The pointwise innerproduct in (8.1), or equivalently in (8.2), induces an inner product on as


The space of -forms obtained by completion of under this inner product is the Hilbert space of square integrable -forms . The other useful space mentioned at the beginning of this section is the Sobolev space .

8.2 Whitney map and Whitney forms

Let be an -dimensional manifold simplicial complex embedded in and the underlying space. The metric on the interiors of the top dimensional simplices of will be the one induced from the embedding Euclidean space . As is usual in finite element methods, finite dimensional subspaces of the function spaces described in the previous paragraph are used in the numerical solution of PDEs. The finite dimensional spaces can be obtained by “embedding” the space of cochains into these spaces by using an interpolation. For example, to embed into one can use the Whitney map , which will be reviewed in this subsection. The image is a linear vector subspace of and is the space of Whitney -forms [52, 21, 13] and is denoted in [3, 4]. (We use instead of used in [4].) The embedding of cochains is analogous to how scalar values at discrete sample points would be interpolated to get a piecewise affine function. In the scalar case also, the space of such functions is a vector subspace of square integrable functions. In fact, , the space of Whitney -forms is the space of continuous piecewise affine functions on . The Whitney map for is actually built from barycentric coordinates which are the building blocks of piecewise linear interpolation. Thus the embedding of into for can be considered to be a generalization of the embedding of into . Thus Whitney forms enable only low order methods. However, arbitrary degree polynomial spaces suitable for use in finite element exterior calculus have been discovered [2, 3]. These however are not yet a part of PyDEC.

The space of Whitney -forms is the space of piecewise smooth differential -forms obtained by applying the Whitney map to -cochains. It can be thought of as a method for interpolating values given on -simplices of a simplicial complex. For example, inside a tetrahedron Whitney forms allow the interpolation of numbers on edges or faces to a smooth 1-form or 2-form respectively. As mentioned above, for 0-cochains, i.e. scalar functions sampled at vertices, the interpolation is the one obtained using the standard scalar piecewise affine basis functions on each simplex, that is the barycentric coordinates corresponding to each vertex of the simplex. We recall the definition of barycentric coordinates followed by the definition of the Whitney map.

Definition 8.5.

Let be a -simplex embedded in . The affine functions , , which when restricted to take the value 1 on vertex and 0 on the others, are called the barycentric coordinates in .

Definition 8.6.

Let be an oriented -simplex in an -dimensional manifold complex , and the corresponding elementary -cochain. We define


where is the barycentric coordinate function with respect to vertex and the notation indicates that the term is omitted from the wedge product. The Whitney map is the above map extended to all of by requiring that be a linear map. is called the Whitney form corresponding to , and for a general cochain , is called the Whitney form corresponding to .

For example, the Whitney form corresponding to the edge is and the Whitney form corresponding to the triangle in a tetrahedron is

Remark 8.7.

If we were using local coordinate charts on a manifold then at any point a -form would be a linear combination of monomials. Note from (8.4) that the Whitney form is a sum of monomials with coefficients. Thus Whitney forms allow us to treat forms at a point as a linear combination of monomials even though we are not using local coordinate charts.

We emphasize again that this section was a review of known material. We have tried to present this material in a manner which makes it easier to explain the examples of Section 11 and the construction of mass matrix for Whitney forms described in the next section.

9 Whitney Inner Product of Cochains

Given a manifold simplicial complex , an inner product between two -cochains and can be defined by first embedding these cochains into using Whitney map and then taking the inner product of the resulting Whitney forms [21].

Definition 9.1.

Given two -cochains , their Whitney inner product is defined by


using the inner product on forms given in (8.3). The matrix for Whitney inner product of -forms in the elementary -cochain basis will be denoted . That is, is a square matrix of order (the number of -simplices in ) such that the entry in row and column is , where and are the elementary -cochains corresponding to the -simplices and with index number and respectively.

The integral in (9.1) is the sum of integrals over each top dimensional simplex in . Inside each such simplex the inner product of smooth forms applies since the Whitney form in each simplex is smooth all the way up to the boundary of the simplex. The interior of each top dimensional simplex is given an inner product that is induced from the standard inner product of the embedding space .

Remark 9.2.

Given cochains we will refer to their representations in the elementary cochain basis also as and . Then the matrix representation of the Whitney inner product of and is .

The inner product of cochains defined in this way is a key concept that connects exterior calculus to finite element methods and different choices of the inner product lead to different discretizations of exterior calculus. This is because the inner product matrix is the mass matrix of finite element methods based on Whitney forms. The details of the efficient computation of for any and will be given in Section 9.2 and 9.3.

Recall that for a Riemannian manifold , if is the codifferential, then the Laplace-deRham operator on -forms is . For a boundaryless the codifferential is the adjoint of the exterior derivative . In case has a boundary, we have instead that


See [1, Exercise 7.5E] for a derivation of the above. Now consider Poisson’s equation on -forms defined on a -dimensional simplicial manifold complex . For simplicity, we’ll consider the weak form of this using smooth forms rather than Sobolev spaces of forms. See [2, 4] for a proper functional analytic treatment. We will also assume that the correct boundary conditions are satisfied, so that the boundary term in (9.2) is 0. In our simple treatment, the weak form of the Poisson’s equation is to find a such that . Thus it is clear that a Galerkin formulation using Whitney forms will require the computation of a term like for cochains and . By the commuting property of Whitney forms (where the second is the coboundary operator on cochains) we have that the above inner product is equal to . (See [21] for a proof of the commuting property.) Now, by definition of the Whitney inner product of cochains in (9.1) this is equal to in the inner product on -cochains. The matrix form of this inner product can be obtained from the mass matrix as . This is what we mean when we say that the stiffness matrix can be computed easily from the mass matrix. The term on the right in the weak form will use the mass matrix . Since codifferential of Whitney forms is 0, the first term in the weak form has to be handled in another way, as described in [4].

Remark 9.3.

Exploiting the aforementioned commutativity of Whitney forms to compute the stiffness matrix represents a significant simplification to our software implementation. While computing the stiffness matrix directly from the definition is possible, it is a complex operation which requires considerable programmer effort, especially if the performance of the implementation is important. In contrast, our formulation requires no additional effort and has the performance of the underlying sparse matrix-matrix multiplication implementation, an optimized and natively-compiled routine. All of the complex indexing, considerations of relative orientation, mappings between faces and indices, etc. is reduced to a simple linear algebra expression. The lower dimensional faces are oriented lexicographically and the orientation information required in stiffness matrix assembly is implicit in the boundary matrices.

9.1 Computing barycentric differentials

Given that the Whitney form in (8.4) is built using wedges of differentials of barycentric coordinates, it is clear that the algorithm for computing an inner product of Whitney forms involves computation of the gradients or differentials of the barycentric coordinates. The following lemma shows how these are computed using simple linear algebra operations.

Lemma 9.4.

Let be a -simplex embedded in , where the vertices are given in some basis for . Let be a matrix whose -th column consists of t