Discrete Exterior Calculus

This work presents a theory of discrete exterior calculus (DEC) motivated by potential applications in computational methods for field theories such as elasticity, fluids, and electromagnetism. In addition, it provides much needed mathematical machinery to enable a systematic development of numerical schemes that mirror the approach of geometric mechanics.

This theory has a long history that we shall outline below in §2, but we aim at a comprehensive, systematic, as well as useful, treatment. Many previous works, as we shall review, are incomplete both in terms of the objects that they treat as well as the types of meshes that they allow.

Our vision of this theory is that it should proceed ab initio as a discrete theory that parallels the continuous one. General views of the subject area of DEC are common in the literature (see, for instance, Mattiussi [2000]), but they usually stress the process of discretizing a continuous theory and the overall approach is tied to this goal. However, if one takes the point of view that the discrete theory can, and indeed should, stand in its own right, then the range of application areas naturally is enriched and increases.

Convergence and consistency considerations alone are inadequate to discriminate between the various choices of discretization available to the numerical analyst, and only by requiring, when appropriate, that the discretization exhibits discrete analogues of continuous properties of interest can we begin to address the question of what makes a discrete theory a canonical discretization of a continuous one.

Applications to Variational Problems. One of the major application areas we envision is to variational problems, be they in mechanics or optimal control. One of the key ingredients in this direction that we imagine will play a key role in the future is that of AVI's (asynchronous variational integrators) designed for the numerical integration of mechanical systems, as in Lew et al. [2003]. These are integration algorithms that respect some of the key features of the continuous theory, such as their multi-symplectic nature and exact conservation laws. They do so by discretizing the underlying variational principles of mechanics rather than discretizing the equations. It is well-known (see the reference just mentioned for some of the literature) that variational problems come equipped with a rich exterior calculus structure and so on the discrete level, such structures will be enhanced by the availability of a discrete exterior calculus. One of the objectives of this chapter is to fill this gap.

Structured Constraints. There are many constraints in numerical algorithms that naturally involve differential forms, such as the divergence constraint for incompressibility of fluids, as well as the fact that differential forms are naturally the fields in electromagnetism, and some of Maxwell's equations are expressed in terms of the divergence and curl operations on these fields. Preserving, as in the mimetic differencing literature, such features directly on the discrete level is another one of the goals, overlapping with our goals for variational problems.

Lattice Theories. Periodic crystalline lattices are of important practical interest in material science, and the anisotropic nature of the material properties arises from the geometry and connectivity of the intermolecular bonds in the lattice. It is natural to model these lattices as inherently discrete objects, and an understanding of discrete curvature that arises from DEC is particularly relevant, since part of the potential energy arises from stretched bonds that can be associated with discrete curvature in the underlying relaxed configuration of the lattice. In particular, this could yield a more detailed geometric understanding of what happens at grain boundaries. Lattice defects can also be associated with discrete curvature when appropriately interpreted. The introduction of a discrete notion of curvature will lay the foundations for a better understanding of the role of geometry in the material properties of solids.

Some of the Key Theoretical Accomplishments. Our development of discrete exterior calculus includes discrete differential forms, the Hodge star operator, the wedge product, the exterior derivative, as well as contraction and the Lie derivative. For example, this approach leads to the proper definition of discrete divergence and curl operators and has already resulted in applications like a discrete Hodge type decomposition of 3D vector fields on irregular grids—see Tong et al. [2003].

Context. We present the theory and some applications of DEC in the context of simplicial complexes of arbitrary finite dimension.

Methodology. We believe that the correct way to proceed with this program is to develop, as we have already stressed, ab initio, a calculus on discrete manifolds which parallels the calculus on smooth manifolds of arbitrary finite dimension. Chapters 6 and 7 of Abraham et al. [1988] are a good source for the concepts and definitions in the smooth case. However we have tried to make this chapter as self-contained as possible. Indeed, one advantage of developing a calculus on discrete manifolds, as we do here, is pedagogical. By using concrete examples of discrete two- and three-dimensional spaces one can explain most of calculus on manifolds at least formally as we will do using the examples in this chapter. The machinery of Riemannian manifolds and general manifold theory from the smooth case is, strictly speaking, not required in the discrete world. The technical terms that are used in this introduction will be defined in subsequent sections, but they should be already familiar to someone who knows the usual exterior calculus on smooth manifolds.

The Objects in DEC. To develop a discrete theory, one must define discrete differential forms along with vector fields and operators involving these. Once discrete forms and vector fields are defined, a calculus can be developed by defining the discrete exterior derivative ($\backslash mathbf\{d\}$), codifferential ($\backslash delta$) and Hodge star ($*$) for operating on forms, discrete wedge product ($\backslash wedge$) for combining forms, discrete flat ($\backslash flat$) and sharp ($\backslash sharp$) operators for going between vector fields and 1-forms and discrete contraction operator ($\backslash mathbf\{i\}\_X$) for combining forms and vector fields. Once these are done, one can then define other useful operators. For example, a discrete Lie derivative ($\backslash mathcal\{L\}\_X$) can be defined by requiring that the Cartan magic (or homotopy) formula hold. A discrete divergence in any dimension can be defined. A discrete Laplace–deRham operator ($\backslash Delta$) can be defined using the usual definition of $\backslash mathbf\{d\}\backslash delta\; +\; \backslash delta\backslash mathbf\{d\}$. When applied to functions, this is the same as the discrete Laplace–Beltrami operator ($\backslash nabla^2$), which is defined as $\backslash text\{div\}\; \backslash circ\; \backslash text\{curl\}$. We define all these operators in this chapter.

The discrete manifolds we work with are simplicial complexes. We will recall the standard formal definitions in §3 but familiar examples of simplicial complexes are meshes of triangles embedded in $\backslash mathbb\{R\}^3$ and meshes made up of tetrahedra occupying a portion of $\backslash mathbb\{R\}^3$. We will assume that the angles and lengths on such discrete manifolds are computed in the embedding space $\backslash mathbb\{R\}^N$ using the standard metric of that space. In other words, in this chapter we do not address the issue of how to discretize a given smooth Riemannian manifold, and how to embed it in $\backslash mathbb\{R\}^N$, since there may be many ways to do this. For example, $\backslash text\{SO\}(3)$ can be embedded in $\backslash mathbb\{R\}^9$ with a constraint, or as the unit quaternions in $\backslash mathbb\{R\}^4$. Another potentially important consideration in discretizing the manifold is that the topology of the simplicial complex should be the same as the manifold to be discretized. This can be verified using the methods of computational homology (see, for example, Kaczynski et al. [2004]), or discrete Morse theory (see, for example, Forman [2002], Wood [2003]). For the purposes of discrete exterior calculus, only local metric information is required, and we will comment towards the end of §3 how to address the issue of embedding in a local fashion, as well as the criterion for a good global embedding.

Our development in this chapter is for the most part formal in that we choose appropriate geometric definitions of the various objects and quantities involved. For the most part, we do not prove that these definitions converge to the smooth counterparts. The definitions are chosen so as to make some important theorems like the generalized Stokes’ theorem true by definition. Moreover, in the cases where previous results are available, we have checked that the operators we obtain match the ones obtained by other means, such as variational derivations.

The use of simplicial chains and cochains as the basic building blocks for a discrete exterior calculus has appeared in several papers. See, for instance, Sen et al. [2000], Adams [1996], Bossavit [2002c], and references therein. These authors view forms as linearly interpolated versions of smooth differential forms, a viewpoint originating from Whitney [1957], who introduced the Whitney and deRham maps that establish an isomorphism between simplicial cochains and Lipschitz differential forms.

We will, however, view discrete forms as real-valued linear functions on the space of chains. These are inherently discrete objects that can be paired with chains of oriented simplices, or their geometric duals, by the bilinear pairing of evaluation. In the next chapter, where we consider applications involving the curvature of a discrete space, we will relax the condition that discrete forms are real-valued, and consider group-valued forms.

Intuitively, this natural pairing of evaluation can be thought of as integration of the discrete form over the chain. This difference from the work of Sen et al. [2000] and Adams [1996] is apparent in the definitions of operations like the wedge product as well.

There is also much interest in a discrete exterior calculus in the computational electromagnetism community, as represented by Bossavit [2001, 2002a,b,c], Gross and Kotiuga [2001], Hiptmair [1999, 2001a,b, 2002], Mattiussi [1997, 2000], Nicolaides and Wang [1998], Teixeira [2001], and Tonti [2002].

Many of the authors cited above, for example, Bossavit [2002c], Sen et al. [2000], and Hiptmair [2002], also introduce the notions of dual complexes in order to construct the Hodge star operator. With the exception of Hiptmair, they use barycentric duals. This works if one develops a theory of discrete forms and does not introduce discrete vector fields. We show later that to introduce discrete vector fields into the theory the notion of circumcentric duals seems to be important.

Other authors, such as Moritz [2000], Moritz and Schwalm [2001], Schwalm et al. [1999], have incorporated vector fields into the cochain based approach to exterior calculus by identifying vector fields with cochains, and having them supported on the same mesh. This is ultimately an unsatisfactory approach, since dual meshes are essential as a means of encoding physically relevant phenomena such as fluxes across boundaries.

The use of primal and dual meshes arises most often as staggered meshes in finite volume and finite difference methods. In fluid computations, for example, the density is often a cell-centered quantity, which can either be represented as a primal object by being associated with the 3-cell, or as a dual object associated with the 0-cell at the center of the 3-cell. Similarly, the flux across boundaries can be associated with the 2-cells that make up the boundary, or the 1-cell which is normal to the boundary.

Another approach to a discrete exterior calculus is presented in Dezin [1995]. He defines a one-dimensional discretization of the real line in much the same way we would. However, to generalize to higher dimensions he introduces a tensor product of this space. This results in logically rectangular meshes. Our calculus, however, is defined over simplicial meshes. A further difference is that like other authors in this field, Dezin [1995] does not introduce vector fields into his theory.

A related effort for three-dimensional domains with logically rectangular meshes is that of Mansfield and Hydon [2001], who established a variational complex for difference equations by constructing a discrete homotopy operator. We construct an analogous homotopy operator for simplicial meshes in proving the discrete Poincaré lemma.

In constructing the discretization of a continuous problem in the context of our formulation of discrete exterior calculus, we first discretize the manifold of interest as a simplicial complex. While this is typically in the form of a simplicial complex that is embedded into Euclidean space, it is only necessary to have an abstract simplicial complex, along with a local metric defined on adjacent vertices. This abstract setting will be addressed further toward the end of this section.

We will now recall some basic definitions of simplices and simplicial complexes, which are standard from simplicial algebraic topology. A more extensive treatment can be found in Munkres [1984].

Definition 3.1. A k-simplex is the convex span of $k\; +\; 1$ geometrically independent points,

The points $v\_0,\; \backslash dots,\; v\_k$ are called the vertices of the simplex, and the number $k$ is called the dimension of the simplex. Any simplex spanned by a (proper) subset of $\backslash \{v\_0,\; \backslash dots,\; v\_k\backslash \}$ is called a (proper) face of $\backslash sigma^k$. If $\backslash sigma^l$ is a proper face of $\backslash sigma^k$, we denote this by $\backslash sigma^l\; \backslash prec\; \backslash sigma^k$.

Example 3.1. Consider 3 non-collinear points $v\_0,\; v\_1$ and $v\_2$ in $\backslash mathbb\{R\}^3$. Then, these three points individually are examples of 0-simplices, to which an orientation is assigned through the choice of a sign. Examples of 1-simplices are the oriented line segments $[v\_0,\; v\_1]$, $[v\_1,\; v\_2]$ and $[v\_0,\; v\_2]$. By writing the vertices in that order we have given orientations to these 1-simplices, i.e., $[v\_0,\; v\_1]$ is oriented from $v\_0$ to $v\_1$. The triangle $[v\_0,\; v\_1,\; v\_2]$ is a 2-simplex oriented in counterclockwise direction. Note that the orientation of $[v\_0,\; v\_2]$ does not agree with that of the triangle.

Definition 3.2. A simplicial complex $K$ in $\backslash mathbb\{R\}^N$ is a collection of simplices in $\backslash mathbb\{R\}^N$, such that,

Definition 3.3. A simplicial triangulation of a polytope $|K|$ is a simplicial complex $K$ such that the union of the simplices of $K$ recovers the polytope $|K|$.

Definition 3.4. If $L$ is a subcollection of $K$ that contains all faces of its elements, then $L$ is a simplicial complex in its own right, and it is called a subcomplex of $K$. One subcomplex of $K$ is the collection of all simplices of $K$ of dimension at most $k$, which is called the k-skeleton of $K$, and is denoted $K^\{(k)\}$.

Circumcentric Subdivision. We will also use the notion of a circumcentric dual or Voronoi mesh of the given primal mesh. We will point to the importance of this choice later on in §7 and 9. We call the Voronoi dual a circumcentric dual since the dual of a simplex is its circumcenter (equidistant from all vertices of the simplex).

Definition 3.5. The circumcenter of a $k$-simplex $\backslash sigma^k$ is given by the center of the $k$-circumsphere, where the $k$-circumsphere is the unique $k$-sphere that has all $k\; +\; 1$ vertices of $\backslash sigma^k$ on its surface. Equivalently, the circumcenter is the unique point in the $k$-dimensional affine space that contains the $k$-simplex that is equidistant from all the $k\; +\; 1$ nodes of the simplex. We will denote the circumcenter of a simplex $\backslash sigma^k$ by $c(\backslash sigma^k)$.

The circumcenter of a simplex $\backslash sigma^k$ can be obtained by taking the intersection of the normals to the $(k\; -\; 1)$-dimensional faces of the simplex, where the normals are emanating from the circumcenter of the face. This allows us to recursively compute the circumcenter.

If we are given the nodes which describe the primal mesh, we can construct a simplicial triangulation by using the Delaunay triangulation, since this ensures that the circumcenter of a simplex is always a point within the simplex. Otherwise we assume that a nice mesh has been given to us, i.e., it is such that the circumcenters lie within the simplices. While this is not essential for our theory it makes some proofs simpler. For some computations the Delaunay triangulation is desirable in that it reduces the maximum aspect ratio of the mesh, which is a factor in determining the rate at which the corresponding numerical scheme converges. But in practice there are many problems for which Delaunay triangulations are a bad idea. See, for example, Schewchuck [2002]. We will address such computational issues in a separate work.

Definition 3.6. The circumcentric subdivision of a simplicial complex is given by the collection of all simplices of the form

where $\backslash sigma\_0\; \backslash prec\; \backslash sigma\_1\; \backslash prec\; \backslash dots\; \backslash prec\; \backslash sigma\_k$, or equivalently, that $\backslash sigma\_i$ is a proper face of $\backslash sigma\_j$ for all .

Circumcentric Dual. We construct a circumcentric dual to a $k$-simplex using the circumcentric duality operator, which is introduced below.

Definition 3.7. The circumcentric duality operator is given by

where the $\backslash epsilon\_\{\backslash sigma^k,\; \backslash dots,\; \backslash sigma^n\}$ coefficient ensures that the orientation of $[c(\backslash sigma^k),\; c(\backslash sigma^\{k+1\}),\; \backslash dots,\; c(\backslash sigma^n)]$ is consistent with the orientation of the primal simplex, and the ambient volume-form.

Orienting $\backslash sigma^k$ is equivalent to choosing a ordered basis, which we shall denote by $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^k$. Similarly, $[c(\backslash sigma^k),\; c(\backslash sigma^\{k+1\}),\; \backslash dots,\; c(\backslash sigma^n)]$ has an orientation denoted by $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$. If the orientation corresponding to $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$ is consistent with the volume-form on the manifold, then $\backslash epsilon\_\{\backslash sigma^k,\; \backslash dots,\; \backslash sigma^n\}\; =\; 1$, otherwise it takes the value $-1$.

We immediately see from the construction of the circumcentric duality operator that the dual elements can be realized as a submesh of the first circumcentric subdivision, since it consists of elements of the form $[c(\backslash sigma\_0),\; \backslash dots,\; c(\backslash sigma\_k)]$, which are, by definition, part of the first circumcentric subdivision.

Example 3.2. The circumcentric duality operator maps a 0-simplex into the convex hull generated by the circumcenters of $n$-simplices that contain the 0-simplex,

and the circumcentric duality operator maps a $n$-simplex into the circumcenter of the $n$-simplex,

This is more clearly illustrated in Figure 1, where the primal and dual elements are color coded to represent the dual relationship between the elements in the primal and dual mesh.

The choice of a circumcentric dual is significant, since it allows us to recover geometrically important objects such as normals to $(n-1)$-dimensional faces, which are obtained by taking their circumcentric dual, whereas, if we were to use a barycentric dual, the dual to a $(n-1)$-dimensional face would not be normal to it.

Orientation of the Dual Cell. Notice that given an oriented simplex $\backslash sigma^k$, which is represented by $[v\_0,\; \backslash dots,\; v\_k]$, the orientation is equivalently represented by $(v\_1\; -\; v\_0)\; \backslash wedge\; (v\_2\; -\; v\_1)\; \backslash wedge\; \backslash dots\; \backslash wedge\; (v\_k\; -\; v\_\{k-1\})$, which we denote by,

which is an equivalence at the level of orientation. It would be nice to express our criterion for determining the orientation of the dual cell in terms of the $(k+1)$-vertex representation.

To determine the orientation of the $(n-k)$-simplex given by $[c(\backslash sigma^k),\; c(\backslash sigma^\{k+1\}),\; \backslash dots,\; c(\backslash sigma^n)]$, or equivalently, $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$, we consider the $n$-simplex given by $[c(\backslash sigma^0),\; \backslash dots,\; c(\backslash sigma^n)]$, where $\backslash sigma^0\; \backslash prec\; \backslash dots\; \backslash prec\; \backslash sigma^k$. This is related to the expression $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$, up to a sign determined by the relative orientation of $[c(\backslash sigma^0),\; \backslash dots,\; c(\backslash sigma^k)]$ and $\backslash sigma^k$. Thus, we have that

Then, we need to check that $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$ is consistent with the volume-form on the manifold, which is represented by the orientation of $\backslash sigma^n$. Thus, we have that the correct orientation for the $[c(\backslash sigma^k),\; c(\backslash sigma^\{k+1\}),\; \backslash dots,\; c(\backslash sigma^n)]$ term is given by,

These two representations of the choice of orientation for the dual cells are equivalent, but the combinatorial definition above might be preferable for the purposes of implementation.

Example 3.3. We would like to compute the orientation of the dual of a 1-simplex, in two dimensions, given the orientation of the two neighboring 2-simplices.

Given a simplicial complex, as shown in Figure 2(a), we consider a 2-simplex of the form $[c(\backslash sigma^0),\; c(\backslash sigma^1),\; c(\backslash sigma^2)]$, which is illustrated in Figure 2(b).

Notice that the orientation is consistent with the given orientation of the 2-simplex, but it is not consistent with the orientation of the primal 1-simplex, so the orientation should be reversed, to give the dual cell illustrated in Figure 2(c).

We summarize the results for the induced orientation of dual cells for the other 2-simplices of the form $[c(\backslash sigma^0),\; c(\backslash sigma^1),\; c(\backslash sigma^2)]$, in Table 1.

Orientation of the Dual of a Dual Cell. While the circumcentric duality operator is a map from the primal simplicial complex to the dual cell complex, we can formally extend the circumcentric duality operator

to a map from the dual cell complex to the primal simplicial complex. However, we need to be slightly careful about the orientation of primal simplex we recover from applying the circumcentric duality operator twice.

We have that, $\backslash star\backslash star(\backslash sigma^k)\; =\; \backslash pm\backslash sigma^k$, where the sign is chosen to ensure the appropriate choice of orientation. If, as before, $\backslash sigma^k$ has an orientation represented by $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^k$, and $\backslash star\backslash sigma^k$ has an orientation represented by $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$, then the orientation of $\backslash star\backslash star(\backslash sigma^k)$ is chosen so that $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n\; \backslash wedge\; dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^k$ is consistent with the ambient volume-form. Since, by construction, $\backslash star(\backslash sigma^k)$, $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$ has an orientation consistent with the ambient volume-form, we need only compare $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n\; \backslash wedge\; dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^k$ with $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$. Notice that it takes $n-k$ transpositions to get the $dx^1$ term in front of the $dx^\{k+1\}\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^n$ terms, and we need to do this $k$ times for each term of $dx^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; dx^k$, so it follows that the sign is simply given by $(-1)^\{k(n-k)\}$, or equivalently,

A similar relationship holds if we use a dual cell instead of the primal simplex $\backslash sigma^k$.

Support Volume of a Primal Simplex and Its Dual Cell. We can think of a cochain as being constructed out of a basis consisting of cosimplices or cocells with value 1 on a single simplex or cell, and 0 otherwise. The way to visualize this cochain is that it is associated with a differential form that has support on what we will refer to as the support volume associated with a given simplex or cell.

Definition 3.8. The support volume of a simplex $\backslash sigma^k$ is a $n$-volume given by the convex hull of the geometric union of the simplex and its circumcentric dual. This is given by

The intersection with $|K|$ is necessary to ensure that the support volume does not extend beyond the polytope $|K|$ which would otherwise occur if $|K|$ is nonconvex.

We extend the notion of a support volume to a dual cell $\backslash star\backslash sigma^k$ by similarly defining

To clarify this definition, we will consider some examples of simplices, their dual cells, and their corresponding support volumes. For two-dimensional simplicial complexes, this is illustrated in Table 2.

The support volume has the nice property that at each dimension, it partitions the polytope $|K|$ into distinct non-intersecting regions associated with each individual $k$-simplex. For any two distinct $k$-simplices, the intersection of their corresponding support volumes have measure zero, and the union of the support volumes of all $k$-simplices recovers the original polytope $|K|$.

Notice, from our construction, that the support volume of a simplex and its dual cell are the same, which suggests that there is an identification between cochains on $k$-simplices and cochains on $(n\; -\; k)$-cells. This is indeed the case, and is a concept associated with the Hodge star for differential forms.

Examples of simplices, their dual cells, and the corresponding support volumes in three dimensions are given in Table 3.

In our subsequent discussion, we will assume that we are given a simplicial complex $K$ of dimension $n$ in $\backslash mathbb\{R\}^N$. Thus, the highest-dimensional simplex in the complex is of dimension $n$ and each 0-simplex (vertex) is in $\backslash mathbb\{R\}^N$. One can obtain this, for example, by starting from 0-simplices, i.e., vertices, and then constructing a Delaunay triangulation, using the vertices as sites. Often, our examples will be for two-dimensional discrete surfaces in $\backslash mathbb\{R\}^3$ made up of triangles (here $n\; =\; 2$ and $N\; =\; 3$) or three-dimensional manifolds made of tetrahedra, possibly embedded in a higher-dimensional space.

Cell Complexes. The circumcentric dual of a primal simplicial complex is an example of a cell complex. The definition of a cell complex follows.

Definition 3.9. A cell complex $\backslash star\; K$ in $\backslash mathbb\{R\}^N$ is a collection of cells in $\backslash mathbb\{R\}^N$ such that,

We will see in the next section that the notion of boundary in the circumcentric dual has to be modified slightly from the geometric notion of a boundary in order for the circumcentric dual to be made into a cell complex.

While it is computationally more convenient to have a global embedding of the simplicial complex into a higher-dimensional ambient space to account for non-flat manifolds it suffices to have an abstract simplicial complex along with a local metric on vertices. The metric is local in the sense that distances between two vertices are only defined if they are part of a common $n$-simplex in the abstract simplicial complex. Then, the local metric is a map $d\; :\; \backslash \{(v\_0,\; v\_1)\; \backslash mid\; v\_0,\; v\_1\; \backslash in\; K^\{(0)\},\; [v\_0,\; v\_1]\; \backslash prec\; \backslash sigma^n\; \backslash in\; K\backslash \}\; \backslash rightarrow\; \backslash mathbb\{R\}$.

The axioms for a local metric are as follows,

Positive.: $d(v\_0,\; v\_1)\; \backslash geq\; 0$, and $d(v\_0,\; v\_0)\; =\; 0,\; \backslash forall\; [v\_0,\; v\_1]\; \backslash prec\; \backslash sigma^n\; \backslash in\; K$.

Strictly Positive.: If $d(v\_0,\; v\_1)\; =\; 0$, then $v\_0\; =\; v\_1,\; \backslash forall\; [v\_0,\; v\_1]\; \backslash prec\; \backslash sigma^n\; \backslash in\; K$.

Symmetry.: $d(v\_0,\; v\_1)\; =\; d(v\_1,\; v\_0),\; \backslash forall\; [v\_0,\; v\_1]\; \backslash prec\; \backslash sigma^n\; \backslash in\; K$.

Triangle Inequality.: $d(v\_0,\; v\_2)\; \backslash leq\; d(v\_0,\; v\_1)\; +\; d(v\_1,\; v\_2),\; \backslash forall\; [v\_0,\; v\_1,\; v\_2]\; \backslash prec\; \backslash sigma^n\; \backslash in\; K$.

This allows us to embed each $n$-simplex locally into $\backslash mathbb\{R\}^n$, and thereby compute all the necessary metric dependent quantities in our formulation. For example, the volume of a $k$-dual cell will be computed as the sum of the $k$-volumes of the dual cell restricted to each $n$-simplex in its local embedding into $\backslash mathbb\{R\}^n$.

This notion of local metrics and local embeddings is consistent with the point of view that exterior calculus is a local theory with operators that operate on objects in the tangent and cotangent space of a fixed point. The issue of comparing objects in different tangent spaces is addressed in the discrete theory of connections on principal bundles in Leok et al. [2003].

This also provides us with a criterion for evaluating a global embedding. The embedding should be such that the metric of the ambient space $\backslash mathbb\{R\}^N$ restricted to the vertices of the complex, thought of as points in $\backslash mathbb\{R\}^N$, agrees with the local metric imposed on the abstract simplicial complex. A global embedding that satisfies this condition will produce the same numerical results in discrete exterior calculus as that obtained using the local embedding method.

It is essential that the metric condition we impose is local, since the notion of distances between points in a manifold which are far away is not a well-defined concept, nor is it particularly useful for embeddings. As the simple example below illustrates, there may not exist any global embeddings into Euclidean space that satisfies a metric constraint imposed for all possible pairs of vertices.

Example 4.1. Consider a circle, with the distance between two points given by the minimal arc length. Consider a discretization given by 4 equidistant points on the circle, labelled $v\_0,\; \backslash dots,\; v\_3$, with the metric distances as follows,

where the indices are evaluated modulo 4, and this distance function is extended to a metric on all pairs of vertices by symmetry. It is easy to verify that this distance function is indeed a metric on vertices.

If we only use the local metric constraint, then we only require that adjacent vertices are separated by 1, and the following is an embedding of the simplicial complex into $\backslash mathbb\{R\}^2$,

If, however, we use the metric defined on all possible pairs of vertices, by considering $v\_0,\; v\_1,\; v\_2$, we have that $d(v\_0,\; v\_1)\; +\; d(v\_1,\; v\_2)\; =\; d(v\_0,\; v\_2)$. Since we are embedding these points into a Euclidean space, it follows that $v\_0,\; v\_1,\; v\_2$ are collinear.

Similarly, by considering $v\_0,\; v\_2,\; v\_3$, we conclude that they are collinear as well, and that $v\_1,\; v\_3$ are coincident, which contradicts $d(v\_1,\; v\_3)\; =\; 2$. Thus, we find that there does not exist a global embedding of the circle into Euclidean space if we require that the embedding is consistent with the metric on vertices defined for all possible pairs of vertices.

We will now define discrete differential forms. We will use some terms (which we will define) from algebraic topology, but it will become clear by looking at the examples that one can gain a clear and working notion of what a discrete form is without any algebraic topology. We start with a few definitions for which more details can be found on page 26 and 27 of Munkres [1984].

Definition 5.1. Let $K$ be a simplicial complex. We denote the free abelian group generated by a basis consisting of oriented $k$-simplices by $C\_k(K;\; \backslash mathbb\{Z\})$. This is the space of finite formal sums of the $k$-simplices, with coefficients in $\backslash mathbb\{Z\}$. Elements of $C\_k(K;\; \backslash mathbb\{Z\})$ are called $k$-chains.

Example 5.1. Figure 3 shows examples of 1-chains and 2-chains.

We view discrete $k$-forms as maps from the space of $k$-chains to $\backslash mathbb\{R\}$. Recalling that the space of $k$-chains is a group, we require that these maps be homomorphisms into the additive group $\backslash mathbb\{R\}$. Thus, discrete forms are what are called cochains in algebraic topology. We will define cochains below in the definition of forms but for more context and more details readers can refer to any algebraic topology text, for example, page 251 of Munkres [1984].

This point of view of forms as cochains is not new. The idea of defining forms as cochains appears, for example, in the works of Adams [1996], Dezin [1995], Hiptmair [1999], and Sen et al. [2000]. Our point of departure is that the other authors go on to develop a theory of discrete exterior calculus of forms only by introducing interpolation of forms, which we will be able to avoid. The formal definition of discrete forms follows.

Definition 5.2. A primal discrete $k$-form $\backslash alpha$ is a homomorphism from the chain group $C\_k(K;\; \backslash mathbb\{Z\})$ to the additive group $\backslash mathbb\{R\}$. Thus, a discrete $k$-form is an element of $\backslash text\{Hom\}(C\_k(K),\; \backslash mathbb\{R\})$, the space of cochains. This space becomes an abelian group if we add two homomorphisms by adding their values in $\backslash mathbb\{R\}$. The standard notation for $\backslash text\{Hom\}(C\_k(K),\; \backslash mathbb\{R\})$ in algebraic topology is $C^k(K;\; \backslash mathbb\{R\})$. But we will often use the notation $\backslash Omega\_d^k(K)$ for this space as a reminder that this is the space of discrete (hence the $d$ subscript) $k$-forms on the simplicial complex $K$. Thus,

Note that, by the above definition, given a $k$-chain $\backslash sum\_i\; a\_i\; c\_i^k$ (where $a\_i\; \backslash in\; \backslash mathbb\{Z\}$) and a discrete $k$-form $\backslash alpha$, we have that

and for two discrete $k$-forms $\backslash alpha,\; \backslash beta\; \backslash in\; \backslash Omega\_d^k(K)$ and a $k$-chain $c\; \backslash in\; C\_k(K;\; \backslash mathbb\{Z\})$,

In the usual exterior calculus on smooth manifolds integration of $k$-forms on a $k$-dimensional manifold is defined in terms of the familiar integration in $\backslash mathbb\{R\}^k$. This is done roughly speaking by doing the integration in local coordinates, and showing that the value is independent of the choice of coordinates, due to the change of variables theorem in $\backslash mathbb\{R\}^k$. For details on this, see the first few pages of Chapter 7 of Abraham et al. [1988]. We will not try to introduce the notion of integration of discrete forms on a simplicial complex. Instead the fundamental quantity that we will work with is the natural bilinear pairing of cochains and chains, defined by evaluation. More formally, we have the following definition.

Definition 5.3. The natural pairing of a $k$-form $\backslash alpha$ and a $k$-chain $c$ is defined as the bilinear pairing

As mentioned above, in discrete exterior calculus, this natural pairing plays the role that integration of forms on chains plays in the usual exterior calculus on smooth manifolds. The two are related by a procedure done at the time of discretization. Indeed, consider a simplicial triangulation $K$ of a polyhedron in $\backslash mathbb\{R\}^n$, i.e., consider a “flat” discrete manifold. If we are discretizing a continuous problem, we will have some smooth forms defined in the space $|K|\; \backslash subset\; \backslash mathbb\{R\}^n$. Consider such a smooth $k$-form $\backslash alpha^k$. In order to define the discrete form $\backslash alpha\_d^k$ corresponding to $\backslash alpha^k$, one would integrate $\backslash alpha^k$ on all the $k$-simplices in $K$. Then, the evaluation of $\backslash alpha\_d^k$ on a $k$-simplex $\backslash sigma^k$ is defined by $\backslash alpha\_d^k(\backslash sigma^k)\; :=\; \backslash int\_\{\backslash sigma^k\}\; \backslash alpha^k$. Thus, discretization is the only place where integration plays a role in our discrete exterior calculus.

In the case of a non-flat manifold, the situation is somewhat complicated by the fact that the smooth manifold, and the simplicial complex, as geometric sets embedded in the ambient space do not coincide. A smooth differential form on the manifold can be discretized into the cochain representation by identifying the vertices of the simplicial complex with points on the manifold, and then using a local chart to identify $k$-simplices with $k$-volumes on the manifold.

There is the possibility of $k$-volumes overlapping even when their corresponding $k$-simplices do not intersect, and this introduces a discretization error that scales like the mesh size. One can alternatively construct geodesic boundary surfaces in an inductive fashion, which yields a partition of the manifold, but this can be computationally prohibitive to compute.

Now we can define the discrete exterior derivative which we will call $\backslash mathbf\{d\}$, as in the usual exterior calculus. The discrete exterior derivative will be defined as the dual, with respect to the natural pairing defined above, of the boundary operator, which is defined below.

Definition 5.4. The boundary operator $\backslash partial\_k\; :\; C\_k(K;\; \backslash mathbb\{Z\})\; \backslash rightarrow\; C\_\{k-1\}(K;\; \backslash mathbb\{Z\})$ is a homomorphism defined by its action on a simplex $\backslash sigma^k\; =\; [v\_0,\; \backslash dots,\; v\_k]$,

where $[v\_0,\; \backslash dots,\; \backslash hat\{v\}\_i,\; \backslash dots,\; v\_k]$ is the $(k-1)$-simplex obtained by omitting the vertex $v\_i$. Note that $\backslash partial\_k\; \backslash circ\; \backslash partial\_\{k+1\}\; =\; 0$.

Example 5.2. Given an oriented triangle $[v\_0,\; v\_1,\; v\_2]$ the boundary, by the above definition, is the chain $[v\_1,\; v\_2]\; -\; [v\_0,\; v\_2]\; +\; [v\_0,\; v\_1]$, which are the three boundary edges of the triangle.

Definition 5.5. On a simplicial complex of dimension $n$, a chain complex is a collection of chain groups and homomorphisms $\backslash partial\_k$, such that,

and $\backslash partial\_k\; \backslash circ\; \backslash partial\_\{k+1\}\; =\; 0$.

Definition 5.6. The coboundary operator, $\backslash delta^k\; :\; C^k(K)\; \backslash rightarrow\; C^\{k+1\}(K)$, is defined by duality to the boundary operator, with respect to the natural bilinear pairing between discrete forms and chains. Specifically, for a discrete form $\backslash alpha^k\; \backslash in\; \backslash Omega\_d^k(K)$, and a chain $c\_\{k+1\}\; \backslash in\; C\_\{k+1\}(K;\; \backslash mathbb\{Z\})$, we define $\backslash delta^k$ by

That is to say

This definition of the coboundary operator induces the cochain complex,

where it is easy to see that $\backslash delta^\{k+1\}\; \backslash circ\; \backslash delta^k\; =\; 0$.

Definition 5.7. The discrete exterior derivative denoted by $\backslash mathbf\{d\}\; :\; \backslash Omega\_d^k(K)\; \backslash rightarrow\; \backslash Omega\_d^\{k+1\}(K)$ is defined to be the coboundary operator $\backslash delta^k$.

Remark 5.1. With the above definition of the exterior derivative, $\backslash mathbf\{d\}\; :\; \backslash Omega\_d^k(K)\; \backslash rightarrow\; \backslash Omega\_d^\{k+1\}(K)$, and the relationship between the natural pairing and integration, one can regard equation 5.1 as a discrete generalized Stokes' theorem. Thus, given a $k$-chain $c$, and a discrete $k$-form $\backslash alpha$, the discrete Stokes' theorem, which is true by definition, states that

Furthermore, it also follows immediately that $\backslash mathbf\{d\}^\{k+1\}\; \backslash mathbf\{d\}^k\; =\; 0$.

Dual Discrete Forms. Everything we have said above in terms of simplices and the simplicial complex $K$ can be said in terms of the cells that are duals of simplices and elements of the dual complex $\backslash star\; K$. One just has to be a little more careful in the definition of the boundary operator, and the definition we construct below is well-defined on the dual cell complex. This gives us the notion of cochains of cells in the dual complex and these are the dual discrete forms.

Definition 5.8. The dual boundary operator, $\backslash partial\_k\; :\; C\_k(\backslash star\; K;\; \backslash mathbb\{Z\})\; \backslash rightarrow\; C\_\{k-1\}(\backslash star\; K;\; \backslash mathbb\{Z\})$, is a homomorphism defined by its action on a dual cell $\backslash hat\{\backslash sigma\}^k\; =\; \backslash star\; \backslash sigma^\{n-k\}\; =\; \backslash star[v\_0,\; \backslash dots,\; v\_\{n-k\}]$,

where $\backslash sigma^\{n-k+1\}$ is oriented so that it is consistent with the induced orientation on $\backslash sigma^\{n-k\}$.

Just as discrete forms come in two flavors, primal and dual (being linear functionals on primal chains or chains made up of dual cells), discrete vector fields also come in two flavors. Before formally defining primal and dual discrete vector fields, consider the examples illustrated in Figure 4. The distinction lies in the choice of basepoints, be they primal or dual vertices, to which we assign vectors.

Definition 7.1. Let $K$ be a flat simplicial complex, that is, the dimension of $K$ is the same as that of the embedding space. A primal discrete vector field $X$ on a flat simplicial complex $K$ is a map from the zero-dimensional primal subcomplex $K^\{(0)\}$ (i.e., the primal vertices) to $\backslash mathbb\{R\}^N$. We will denote the space of such vector fields by $\backslash mathfrak\{X\}\_d(K)$. The value of such a vector field is piecewise constant on the dual $n$-cells of $\backslash star\; K$. Thus, we could just as well have called such vector fields dual and defined them as functions on the $n$-cells of $\backslash star\; K$.

Definition 7.2. A dual discrete vector field $X$ on a simplicial complex $K$ is a map from the zero-dimensional dual subcomplex $(\backslash star\; K)^\{(0)\}$ (i.e, the circumcenters of the primal $n$ simplices) to $\backslash mathbb\{R\}^N$ such that its value on each dual vertex is tangential to the corresponding primal $n$-simplex. We will denote the space of such vector fields by $\backslash mathfrak\{X\}\_d(\backslash star\; K)$. The value of such a vector field is piecewise constant on the $n$-simplices of $K$. Thus, we could just as well have called such vector fields primal and defined them as functions on the $n$-simplices of $K$.

Remark 7.1. In this paper we have defined the primal vector fields only for flat meshes. We will address the issue of non-flat meshes in separate work.

As in the smooth exterior calculus, we want to define the flat ($\backslash flat$) and sharp ($\backslash sharp$) operators that relate forms to vector fields. This allows one to write various vector calculus identities in terms of exterior calculus.

Definition 7.3. Given a simplicial complex $K$ of dimension $n$, the discrete flat operator on a dual vector field, $\backslash flat\; :\; \backslash mathfrak\{X\}\_d(\backslash star\; K)\; \backslash rightarrow\; \backslash Omega^d(K)$, is defined by its evaluation on a primal 1 simplex $\backslash sigma^1$,

where $X\; \backslash cdot\; \backslash bar\{\backslash sigma\}^1$ is the usual dot product of vectors in $\backslash mathbb\{R\}^N$, and $\backslash bar\{\backslash sigma\}^1$ stands for the vector corresponding to $\backslash sigma^1$, and with the same orientation. The sum is over all $\backslash sigma^n$ containing the edge $\backslash sigma^1$. The volume factors are in dimension $n$.

Definition 7.4. Given a simplicial complex $K$ of dimension $n$, the discrete sharp operator on a primal 1-form, $\backslash sharp\; :\; \backslash Omega^d(K)\; \backslash rightarrow\; \backslash mathfrak\{X\}\_d(\backslash star\; K)$, is defined by its evaluation on a given vertex $v$,

where the outer sum is over all 1-simplices containing the vertex $v$, and the inner sum is over all $n$-simplices containing the 1-simplex $[v,\; \backslash sigma^0]$. The volume factors are in dimension $n$, and the vector $\backslash hat\{n\}\_\{[v,\; \backslash sigma^0]\}$ is the normal vector to the simplex $[v,\; \backslash sigma^0]$, pointing into the $n$-simplex $\backslash sigma^n$.

For a discussion of the proliferation of discrete sharp and flat operators that arise from considering the interpolation of differential forms and vector fields, please see Hirani [2003].

Lemma 8.1. The discrete wedge product, $\backslash wedge\; :\; C^k(K)\; \backslash times\; C^l(K)\; \backslash rightarrow\; C^\{k+l\}(K)$, is anti-commutative, i.e.,

Proof. We first rewrite the expression for the discrete wedge product using the following computation,

Here, we used the elementary fact, from permutation group theory, that a $k+1$ cycle can be written as the product of $k$ transpositions, which accounts for the $(-1)^k$ factors. Also, $\backslash rho$ is a transposition of 0 and $k$. Then, the discrete wedge product can be rewritten as

For ease of notation, we denote $[v\_0,\; \backslash dots,\; v\_k]$ by $\backslash sigma^k$, and $[v\_0,\; v\_\{k+1\},\; \backslash dots,\; v\_\{k+l\}]$ by $\backslash sigma^l$. Then, we have

Furthermore, we denote $[v\_0,\; v\_\{l+1\},\; \backslash dots,\; v\_\{k+l\}]$ by $\backslash bar\{\backslash sigma\}^k$, and $[v\_0,\; v\_1,\; \backslash dots,\; v\_l]$ by $\backslash bar\{\backslash sigma\}^l$. Then,

Consider the permutation $\backslash theta\; \backslash in\; S\_\{k+l+1\}$, given by

which has the property that

Then, we have

By making the substitution, $\backslash tau\; =\; \backslash bar\{\backslash tau\}\backslash theta$, and noting that $S\_\{k+l+1\}\backslash theta\; =\; S\_\{k+l+1\}$, we obtain

To obtain the desired result, we simply need to compute the sign of $\backslash theta$, which is given by

This follows from the observation that in order to move each of the last $l$ vertices of $\backslash sigma^\{k+l\}$ forward, we require $k$ transpositions with $v\_1,\; \backslash dots,\; v\_k$. Therefore, we obtain

and

Lemma 8.2. The discrete wedge product satisfies the Leibniz rule,

Proof. The proof of the Leibniz rule for discrete wedge products is directly analogous to the proof of the coboundary formula for the simplicial cup product on cochains, which can be found on page 206 of Hatcher [2001]. This is because the discrete exterior derivative is precisely the boundary operator, and the wedge product is constructed out of weighted sums of cup products.

The cup product satisfies the Leibniz rule for an given partial ordering of the vertices, and the permutations in the signed sum in the discrete wedge product correspond to different choices of partial ordering. We then obtain the Leibniz rule for the discrete wedge product by applying it term-wise for each choice of permutation.

Consider

and

The last set of terms, $i\; =\; k+1$, of the first expression cancels the first set of terms, $i\; =\; k$, of the second expression, and what remains is simply $\backslash langle\; \backslash alpha^k\; \backslash wedge\; \backslash beta^l,\; \backslash partial\; \backslash sigma^\{k+l+1\}\; \backslash rangle$. Therefore, we can conclude that

or simply that the Leibniz rule for discrete differential forms holds,

Associativity for the Wedge Product. The discrete wedge product which we have introduced is not associative in general. This is a consequence of the fact that the stencil for the two possible triple wedge products are not the same. In the expression for $\backslash langle\; \backslash alpha^k\; \backslash wedge\; (\backslash beta^l\; \backslash wedge\; \backslash gamma^m),\; \backslash sigma^\{k+l+m\}\; \backslash rangle$, each term in the double summation consists of a geometric factor multiplied by $\backslash langle\; \backslash alpha^k,\; \backslash sigma^k\; \backslash rangle\; \backslash langle\; \backslash beta^l,\; \backslash sigma^l\; \backslash rangle\; \backslash langle\; \backslash gamma^m,\; \backslash sigma^m\; \backslash rangle$ for some $k,\; l,\; m$ simplices $\backslash sigma^k,\; \backslash sigma^l,\; \backslash sigma^m$.

Since $\backslash beta^l$ and $\backslash gamma^m$ are wedged together first, $\backslash sigma^l$ and $\backslash sigma^m$ will always share a common vertex, but $\backslash sigma^k$ could have a vertex in common with only $\backslash sigma^l$, or only $\backslash sigma^m$, or both. We can represent this in a graph, where the nodes denote the three simplices, which are connected by an edge if, and only if, they share a common vertex. The graphical representation of the terms which arise in the two possible triple wedge products are given in Figure 6.

For the wedge product to be associative for all forms, the two stencils must agree. Since the stencils for the two possible triple wedge products differ, the wedge product is not associative in general. However, in the case of closed forms, we can rewrite the terms in the sum so that all the discrete forms are evaluated on triples of simplices that share a common vertex. This is illustrated graphically in Figure 7.

This result is proved rigorously in the follow lemma.

Lemma 8.3. The discrete wedge product is associative for closed forms. That is to say, for $\backslash alpha^k\; \backslash in\; C^k(K)$, $\backslash beta^l\; \backslash in\; C^l(K)$, $\backslash gamma^m\; \backslash in\; C^m(K)$, such that $\backslash mathbf\{d\}\backslash alpha^k\; =\; 0$, $\backslash mathbf\{d\}\backslash beta^l\; =\; 0$, $\backslash mathbf\{d\}\backslash gamma^m\; =\; 0$, we have that

Proof.

Here, either $\backslash rho\backslash tau(k)\; =\; \backslash tau(k+l)$, in which case all three permuted simplices share $v\_\{\backslash tau(k+l)\}$ as a common vertex, or we need to rewrite either $\backslash langle\; \backslash alpha^k,\; \backslash rho\backslash tau[v\_0,\; \backslash dots,\; v\_k]\; \backslash rangle$ or $\backslash langle\; \backslash beta^l,\; \backslash rho\backslash tau[v\_k,\; \backslash dots,\; v\_\{k+l\}]\; \backslash rangle$, using the fact that $\backslash alpha^k$ and $\backslash beta^l$ are closed forms.

If $v\_\{\backslash tau(k+l)\}\; \backslash notin\; \backslash rho\backslash tau[v\_0,\; \backslash dots,\; v\_k]$, then we need to rewrite $\backslash langle\; \backslash alpha^k,\; \backslash rho\backslash tau[v\_0,\; \backslash dots,\; v\_k]\; \backslash rangle$ by considering the simplex obtained by adding the vertex $v\_\{\backslash tau(k+l)\}$ to $\backslash rho\backslash tau[v\_0,\; \backslash dots,\; v\_k]$, which is $[v\_\{\backslash tau(k+l)\},\; v\_\{\backslash rho\backslash tau(0)\},\; \backslash dots,\; v\_\{\backslash rho\backslash tau(k)\}]$. Then, since $\backslash alpha^k$ is closed, we have that

or equivalently,

Notice that all the simplices in the sum, with the exception of the last one, will share two vertices, $v\_\{\backslash tau(k+l)\}$ and $v\_\{\backslash rho\backslash tau(k)\}$ with $\backslash rho\backslash tau[v\_k,\; \backslash dots,\; v\_\{k+l\}]$, and so their contribution in the triple wedge product will vanish due to the anti-symmetrized sum.

Similarly, if $v\_\{\backslash tau(k+l)\}\; \backslash notin\; \backslash rho\backslash tau[v\_k,\; \backslash dots,\; v\_\{k+l\}]$, using the fact that $\backslash beta^l$ is closed yields

As before, all the simplices in the sum, with the exception of the last one, will share two vertices, $v\_\{\backslash tau(k+l)\}$ and $v\_\{\backslash rho\backslash tau(k)\}$ with $\backslash rho\backslash tau[v\_0,\; \backslash dots,\; v\_k]$, and so their contribution in the triple wedge product will vanish due to the anti-symmetrized sum.

This allows us to rewrite the triple wedge product in the case of closed forms as

where $\backslash tau\; \backslash in\; S\_\{k+l+m\}$ is thought of as acting on the set $\backslash \{1,\; \backslash dots,\; k+l+m\backslash \}$, and $\backslash rho\_i$ is a transposition of 0 and $i$. A similar argument allows us to write $\backslash alpha^k\; \backslash wedge\; (\backslash beta^l\; \backslash wedge\; \backslash gamma^m)$ in the same form, and therefore, the wedge product is associative for closed forms. $\backslash square$

Remark 8.1. This lemma is significant, since if we think of a constant smooth differential form, and discretize it to obtain a discrete differential form, this discrete form will be closed. As such, this lemma states that in the infinitesimal limit, the discrete wedge product we have defined will be associative.

In practice, if we have a mesh with characteristic length $\backslash Delta\; x$, then we will have that

which is to say that the average of the associativity defect is of the order of the mesh size, and therefore vanishes in the infinitesimal limit.

In this section, we will illustrate the application of some of the DEC operations we have previously defined to the construction of new discrete operators such as the divergence and Laplace–Beltrami operators.

Divergence. The divergence of a vector field is given in terms of the Lie derivative of the volume-form, by the expression, $(\backslash operatorname\{div\}(\backslash mathbf\{X\})\backslash mu\; =\; \backslash mathcal\{L\}\_\{\backslash mathbf\{X\}\}\backslash mu$. Physically, this corresponds to the net flow per unit volume of an infinitesimal volume about a point.

We will define the discrete divergence by using the formulas defining them in the smooth exterior calculus. The divergence definition will be valid for arbitrary dimensions. The resulting expressions involve operators that we have already defined and so we can actually perform some calculations to express these quantities in terms of geometric quantities. We will show that the resulting expression in terms of geometric quantities is the same as that derived by variational means in Tong et al. [2003].

Definition 9.1. For a discrete dual vector field $X$ the divergence $\backslash operatorname\{div\}(X)$ is defined to be

Remark 9.1. The above definition is a theorem in smooth exterior calculus. See, for example, page 458 of Abraham et al. [1988].

As an example, we will now compute the divergence of a discrete dual vector field on a two-dimensional simplicial complex $K$, as illustrated in Figure 8.

A similar derivation works in higher dimensions, where one needs to be mindful of the sign that arises from applying the Hodge star twice, $**\backslash alpha^k\; =\; (-1)^\{k(n-k)\}\backslash alpha^k$. Since $\backslash operatorname\{div\}(X)\; =\; -\backslash delta\; X^b$, it follows that $\backslash operatorname\{div\}(X)\; =\; *\backslash mathbf\{d\}\; *\; X^b$. Since this is a primal 0-form it can be evaluated on a 0-simplex $\backslash sigma^0$, and we have that

Using the definition of discrete Hodge star, and the discrete generalized Stokes' theorem, we get

The second equality is obtained by applying the definition of the Hodge star, and the last equality is obtained by applying the discrete generalized Stokes' theorem. But,

as given by the expression for the boundary of a dual cell in Equation 5.8. Thus,

This expression has the nice property that the divergence theorem holds on any dual $n$-chain, which, as a set, is a simply connected subset of $|K|$. Furthermore, the coefficients we computed for the discrete divergence operator are the unique ones for which a discrete divergence theorem holds.

Laplace–Beltrami. The Laplace–Beltrami operator is the generalization of the Laplacian to curved spaces. In the smooth case the Laplace–Beltrami operator on smooth functions is defined to be $\backslash nabla^2\; =\; \backslash operatorname\{div\}\; \backslash circ\; \backslash operatorname\{curl\}\; =\; \backslash delta\; d$. See, for example, page 459 of Abraham et al. [1988]. Thus, in the smooth case, the Laplace–Beltrami on functions is a special case of the more general Laplace–deRham operator, $\backslash Delta\; :\; \backslash Omega^k(M)\; \backslash rightarrow\; \backslash Omega^k(M)$, defined by $\backslash Delta\; =\; \backslash delta\; d\; +\; d\backslash delta$.

As an example, we compute $\backslash Delta\; f$ on a primal vertex $\backslash sigma^0$, where $f\; \backslash in\; \backslash Omega\_d^0(K)$, and $K$ is a (not necessarily flat) triangle mesh in $\backslash mathbb\{R\}^3$, as illustrated in Figure 9.

This calculation is done below.

where $\backslash partial\; \backslash sigma^1\; =\; v\; -\; \backslash sigma^0$. But, the above is the same as the formula involving cotangents found by Meyer et al. [2002] without using discrete exterior calculus.

Another interesting aspect, which will be discussed in §12, is that the characterization of harmonic functions as those functions which vanish when the Laplace–Beltrami operator is applied is equivalent to that obtained from a discrete variational principle using DEC as the means of discretizing the Lagrangian.

In this section we will discuss some more operators that involve vector fields, namely contraction, and Lie derivatives.

For contraction, we will first define the usual smooth contraction algebraically, by relating it to Hodge star and wedge products. This yields one potential approach to defining discrete contraction. However, since in the discrete theory we are only concerned with integrals of forms, we can use the interesting notion of extrusion of a manifold by the flow of a vector field to define the integral of a contracted discrete differential form.

We learned about this definition of contraction via extrusion from Bossavit [2002b], who goes on to define discrete extrusion in his paper. Thus, he is able to obtain a definition of discrete contraction. Extrusion turns out to be a very nice way to define integrals of operators involving vector fields, and we will show how to define integrals of Lie derivatives via extrusion, which will yield discrete Lie derivatives.

Definition 10.1. Given a manifold $M$, and $S$, a $k$-dimensional submanifold of $M$, and a vector field $X\; \backslash in\; \backslash mathfrak\{X\}(M)$, we call the manifold obtained by sweeping $S$ along the flow of $X$ for time $t$ as the extrusion of $S$ by $X$ for time $t$, and denote it by $E\_X^t(S)$. The manifold $S$ carried by the flow for time $t$ will be denoted $\backslash varphi\_X^t(S)$.

Example 10.1. Figure 10 illustrates the 2-simplex that arises from the extrusion of a 1-simplex by a discrete vector field that is interpolated using a linear shape function.

Contraction (Extrusion). We first establish an integral property of the contraction operator.

Lemma 10.1.

Proof. Prove instead that

Then, by first fundamental theorem of calculus, the desired result will follow. To prove the above, simply take coordinates on $S$ and carry them along with the flow and define the transversal coordinate to be the flow of $X$. This proof is sketched in Bossavit [2002b]. $\backslash square$

This lemma allows us to interpret contraction as being the dual, under the integration pairing between $k$-forms and $k$-volumes, to the geometric operation of extrusion. The discrete contraction operator is then given by

where the evaluation of the RHS will typically require that the discrete differential form and the discrete vector field are appropriately interpolated.

Remark 10.1. Since the dynamic definition of the contraction operator only depends on the derivative of pairing of the differential form with the extruded region, it will only depend on the vector field in the region $S$, and not on its extension into the rest of the domain.

In addition, if the interpolation for the discrete vector field satisfies a superposition principle, then the discrete contraction operator will satisfy a corresponding superposition principle.

Contraction (Algebraic). Contraction is an operator that allows one to combine vector fields and forms. For a smooth manifold $M$, the contraction of a vector field $X\; \backslash in\; \backslash mathfrak\{X\}(M)$ with a $(k+1)$-form $\backslash alpha\; \backslash in\; \backslash Omega^\{k+1\}(M)$ is written as $\backslash mathbf\{i\}\_X\; \backslash alpha$, and for vector fields $X\_1,\; \backslash dots,\; X\_k\; \backslash in\; \backslash mathfrak\{X\}(M)$, the contraction in smooth exterior calculus is defined by

We define contraction by using an identity that is true in smooth exterior calculus. This identity originally appeared in Hirani [2003], and we state it here with proof.

Lemma 10.2 (Hirani [2003]). Given a smooth manifold $M$ of dimension $n$, a vector field $X\; \backslash in\; \backslash mathfrak\{X\}(M)$, and a $k$-form $\backslash alpha\; \backslash in\; \backslash Omega^k(M)$, we have that

Proof. Recall that for a smooth function $f\; \backslash in\; \backslash Omega^0(M)$, we have that $\backslash mathbf\{i\}\_X\; \backslash alpha\; =\; f\; \backslash mathbf\{i\}\_X\; \backslash alpha$. This, and the multilinearity of $\backslash alpha$, implies that it is enough to show the result in terms of basis elements. In particular, let $\backslash tau\; \backslash in\; S\_n$ be a permutation of the numbers $1,\; \backslash dots,\; n$, such that , and . Let $X\; =\; e\_\{\backslash tau(j)\}$, for some $j\; \backslash in\; 1,\; \backslash dots,\; n$. Then, we have to show that

It is easy to see that the LHS is 0 if $j\; >\; k$, and it is

otherwise, where $\backslash widehat\{e^\{\backslash tau(j)\}\}$ means that $e^\{\backslash tau(j)\}$ is omitted from the wedge product. Now, on the RHS of Equation 10.2, we have that

Thus, the RHS is equal to

which is 0 as required if $j\; >\; k$. So, assume that $1\; \backslash leq\; j\; \backslash leq\; k$. We need to compute

which is given by

where the sign $s\; =\; \backslash pm\; 1$, such that the equation,

holds for the standard volume-form, $\backslash mu\; =\; e^1\; \backslash wedge\; \backslash dots\; \backslash wedge\; e^n$. This implies that

Then, $\backslash text\{RHS\}\; =\; \backslash text\{LHS\}$ as required. $\backslash square$

Since we have expressions for the discrete Hodge star $(*)$, wedge product $(\backslash wedge)$, and flat $(\backslash flat)$, we have the necessary ingredients to use the algebraic expression proved in the above lemma to construct a discrete contraction operator.

One has to note, however, that the wedge product is only associative for closed forms, and as a consequence, the Leibniz rule for the resulting contraction operator will only hold for closed forms as well. This is, however, sufficient to establish that the Leibniz rule for the discrete contraction will hold in the limit as the mesh is refined.

Lie Derivative (Extrusion). As was the case with contraction, we will establish an integral identity that allows the Lie derivative to be interpreted as the dual of a geometric operation on a volume. This involves the flow of a volume by a vector field, and it is illustrated in the following example.

Example 10.2. Figure 11 illustrates the flow of a 1-simplex by a discrete vector field interpolated using a linear shape function.

Lemma 10.3.

Proof.

□

This lemma allows us to define a discrete Lie derivative as follows,

where, as before, evaluating the RHS will require the discrete differential form and discrete vector field to be appropriately interpolated.

Lie Derivative (Algebraic). Alternatively, as we have expressions for the discrete contraction operator ($\backslash mathbf\{i\}\_X$), and exterior derivative ($\backslash mathbf\{d\}$), we can construct a discrete Lie derivative using the Cartan magic formula,

As is the case with the algebraic definition of the discrete contraction, the discrete Lie derivative will only satisfy a Leibniz rule for closed forms. As before, this is sufficient to establish that the Leibniz rule will hold in the limit as the mesh is refined.

In this section, we will prove the discrete Poincaré lemma by constructing a homotopy operator through a generalized cocone construction. This section is based on the work in Desbrun et al. [2003].

The standard cocone construction fails at the discrete level, since the cone of a simplex is not, in general, expressible as a chain in the simplicial complex. As such, the standard cocone does not necessarily map $k$-cochains to $(k-1)$-cochains.

An example of how the standard cone construction fails to map chains to chains is illustrated in Figure 12. Given the simplicial complex on the left, consisting of triangles, edges and nodes, we wish, in the center figure, to consider the cone of the bold edge with respect to the top most node. Clearly, the resulting cone in the right figure, which is shaded grey, cannot be expressed as a combination of the triangles in the original complex.

In this subsection, a generalized cone operator that is valid for chains is developed which has the essential homotopy properties to yield a discrete analogue of the Poincaré lemma.

We will first consider the case of trivially star-shaped complexes, followed by logically star-shaped complexes, before generalizing the result to contractible complexes.

Definition 11.1. Given a $k$-simplex $\backslash sigma^k\; =\; [v\_0,\; \backslash dots,\; v\_k]$ we construct the cone with vertex $w$ and base $\backslash sigma^k$, as follows,

Lemma 11.1. The geometric cone operator satisfies the following property,

Proof. This is a standard result from simplicial algebraic topology. □

We recall that discrete variational mechanics is based on a discrete analogue of Hamilton's principle, and they yield the discrete Euler–Lagrange equations. A particularly interesting property of DEC arises when it is used to construct the discrete Lagrangian for harmonic functions, and Maxwell's equations.

In particular, for these examples, the following diagram commutes,

Which is to say that directly discretizing the differential equations for harmonic functions, and Maxwell's equations using DEC results in the same expressions as the discrete Euler–Lagrange equations associated with a discrete Lagrangian which is discretized from the corresponding continuous Lagrangian by using DEC as the discretization scheme.

This is significant, since it implies that when DEC is used to discretize these equations, the corresponding numerical scheme which is obtained is variational, and consequently exhibits excellent structure-preserving properties.

In the variational principles for both harmonic functions and Maxwell's equations, we require the $L^2$ norm obtained from the $L^2$ inner product on $\backslash Omega^k(M)$, which is given by

The discrete analogue of this requires a primal-dual wedge product, which is given below for forms of complementary dimension.

Definition 12.1. Given a primal discrete $k$-form $\backslash alpha^k\; \backslash in\; \backslash Omega\_d^k(K)$, and a dual discrete $(n-k)$-form $\backslash hat\{\backslash beta\}^\{n-k\}\; \backslash in\; \backslash Omega\_d^\{n-k\}(*K)$, the discrete primal-dual wedge product is defined as follows,

where $V\_\{\backslash sigma^k\}$ is the $n$-dimensional support volume obtained by taking the convex hull of the simplex $\backslash sigma^k$ and its dual cell $\backslash star\; \backslash sigma^k$.

The corresponding $L^2$ inner product is as follows.

Definition 12.2. Given two primal discrete $k$-forms, $\backslash alpha^k,\; \backslash beta^k\; \backslash in\; \backslash Omega\_d^k(K)$, their discrete $L^2$ inner product, $\backslash langle\; \backslash alpha^k,\; \backslash beta^k\; \backslash rangle\_d$, is given by

Remark 12.1. Notice that it would have been quite natural from the smooth theory to propose the following metric tensor $\backslash langle\backslash langle\; \backslash cdot,\; \backslash cdot\; \backslash rangle\backslash rangle$ for differential forms,

where the $|V\_\{\backslash sigma^k\}|$ is the factor arising from integrating the volume-form over $V\_\{\backslash sigma^k\}$, and

is what we would expect for $\backslash langle\backslash langle\; \backslash alpha^k,\; \backslash beta^k\; \backslash rangle\backslash rangle$, if the forms $\backslash alpha^k$ and $\backslash beta^k$ were constant on $\backslash sigma^k$, which is the product of the average values of $\backslash alpha^k$, and $\backslash beta^k$.

If we adopt this as our definition of the metric tensor for forms, we can recover the definition we obtained in §6 for the Hodge star operator. Starting from the definition from the smooth theory,

and expanding this in terms of the metric tensor for discrete forms, and the primal-dual wedge operator, we obtain

When we eliminate common factors from both sides, we obtain the expression,

which is the expression we previously obtained in Definition 6.1 of §6.

The $L^2$ norm for discrete differential forms is given below.

Definition 12.3. Given a primal discrete $k$-form $\backslash alpha^k\; \backslash in\; \backslash Omega\_d^k(K)$, its discrete $L^2$ norm is given by

Given these definitions, we can now reproduce some computations that were originally shown in Castrillón-López [2003].

Harmonic Functions. Harmonic functions $\backslash phi\; :\; M\; \backslash rightarrow\; \backslash mathbb\{R\}$ can be characterized in a variational fashion as extremals of the following action functional,

where $\backslash mathbf\{v\}$ is a Riemannian volume-form in $M$. The corresponding Euler–Lagrange equation is given by

which is the familiar characterization of harmonic functions in terms of the Laplace–Beltrami operator.

The discrete action functional can be expressed in terms of the $L^2$ norm we introduced above for discrete forms,

The basic variations needed for the determination of the discrete Euler–Lagrange operator are obtained from variations that vary the value of the function $\backslash phi$ at a given vertex $v\_0$, leaving the other values fixed. These variations have the form,

where $\backslash tilde\{\backslash eta\}\; \backslash in\; \backslash Omega^0(M;\; \backslash mathbb\{R\})$ is such that $\backslash langle\; \backslash tilde\{\backslash eta\},\; v\_0\; \backslash rangle\; =\; 1$, and $\backslash langle\; \backslash tilde\{\backslash eta\},\; v\; \backslash rangle\; =\; 0$, for any $v\; \backslash in\; K^\{(0)\}\; -\; \backslash \{v\_0\backslash \}$. This family of variations is enough to establish the variational principle. That is, we have

where $\backslash operatorname\{sgn\}(\backslash sigma^1;\; v)$ stands for the sign of $\backslash sigma^1$ with respect to $v$. Which is to say, $\backslash operatorname\{sgn\}(\backslash sigma^1;\; v)\; =\; 1$ if $\backslash sigma^1\; =\; [v\text{'},\; v]$, and $\backslash operatorname\{sgn\}(\backslash sigma^1;\; v)\; =\; -1$ if $\backslash sigma^1\; =\; [v,\; v\text{'}]$. On the other hand,

where in the second to last equality, one has to note that the border of the dual cell of a vertex $v\_0$ consists, up to orientation, in the dual of all the 1-simplices starting from $v\_0$. This is illustrated in Figure 21, and follows from a general expression for the boundary of a dual cell that was given in Definition 5.8.

The sign factor comes from the relation between the orientation of the dual of the 1-simplices and that of $\backslash partial\; *\; v\_0$. From this, we conclude that the variational discrete equation, given in Equation 12.1, is equivalent to the vanishing of the discrete Laplace–Beltrami operator,

Maxwell Equations. We can formulate the Maxwell equations of electromagnetism in a covariant fashion by considering the 1-form $A$ (the potential) as our fundamental variable in a Lorentzian manifold $X$. The action functional for a Lagrangian formulation of electromagnetism is given by,

where $\backslash |\backslash cdot\backslash |$ is the norm on forms induced by the Lorentzian metric on $X$, and $\backslash mathbf\{v\}$ is the pseudo-Riemannian volume-form. The 1-form $A$ is related to the 4-vector potential encountered in the relativistic formulation of electromagnetism (see, for example Jackson [1998]).

The Euler–Lagrange equation corresponding to this action functional is given by

In terms of the field strength, $F\; =\; \backslash mathbf\{d\}A$, the last equation is usually rewritten as

which is the geometric formulation of the Maxwell equations.

For the purposes of simplicity of exposition, we consider the special case where the Lorentzian manifold decomposes into $X\; =\; M\; \backslash times\; \backslash mathbb\{R\}$, where $(M,\; g)$ is a compact Riemannian 3-manifold. In formulating the discrete version of this variational problem, we need to generalize the notion of a discrete Hodge dual to take into account the pseudo-Riemannian metric structure. This can be subtle in practice, and to overcome this, we consider a special family of complexes instead.

Let $K\text{'}$ be a simplicial complex modelling $M$. For the sake of simplicity we consider $M\; =\; \backslash mathbb\{R\}^3$ although this is not strictly necessary. We now consider a discretization $\backslash \{t\_n\backslash \}\_\{n\; \backslash in\; \backslash mathbb\{Z\}\}$ of $\backslash mathbb\{R\}$. We define the complex $K$, modelling $X\; =\; \backslash mathbb\{R\}^4$, the cells of which are the sets $\backslash sigma\; =\; \backslash sigma\text{'}\; \backslash times\; \backslash \{t\_n\backslash \}\; \backslash subset\; \backslash mathbb\{R\}^3\; \backslash times\; \backslash mathbb\{R\}$, and $\backslash sigma\; =\; \backslash sigma\text{'}\; \backslash times\; (t\_n,\; t\_\{n+1\})\; \backslash subset\; \backslash mathbb\{R\}^3\; \backslash times\; \backslash mathbb\{R\}$ for any $\backslash sigma\text{'}\; \backslash in\; K\text{'}$ and $n\; \backslash in\; \backslash mathbb\{Z\}$. Of course, this is not a simplicial complex but rather a “prismal” complex, as shown in Figure 22.

The advantage of these cell complexes is the existence of the Voronoi dual. More precisely, given any prismal cells $\backslash sigma\text{'}\; \backslash times\; \backslash \{t\_n\backslash \}\; \backslash in\; K$ and $\backslash sigma\text{'}\; \backslash times\; (t\_n,\; t\_\{n+1\})\; \backslash in\; K$, the Lorentz orthonormal to any of its edges coincide with the Euclidean one in $\backslash mathbb\{R\}^4$ and the existence of the circumcenter is thus guaranteed. In other words, the Lorentz circumcentric dual $\backslash star\; K$ to $K$ is the same as the Euclidean one in $\backslash mathbb\{R\}^4$.

Remark 12.2. Much of the construction above can be carried out more generally by considering arbitrary cell complexes in $\backslash mathbb\{R\}^4$ that are not necessarily prismal, as long as none of its 1-cells are lightlike. This causality condition is necessary to ensure that the circumcentric dual complex is well-behaved. However, it is sufficient for computational purposes that the complex is well-centered, in the sense that the Lorentzian circumcenter of each cell is contained inside the cell. These issues will be addressed in future work.

Recall that the Hodge star $\backslash star$ is uniquely defined by satisfying the following expression,

for all $\backslash alpha,\; \backslash beta\; \backslash in\; \backslash Omega^k(X)$. The upshot of this is that the Hodge star operator depends on the metric, and since we have a pseudo-Riemannian metric, there is a sign that is introduced in our expression for the discrete Hodge star (Definition 6.1) that depends on whether the cell it is applied to is either spacelike or timelike. The discrete Hodge star for prismal complexes in Lorentzian space is given below.

Definition 12.4. The discrete Hodge star for prismal complexes in Lorentzian space is a map $\backslash star\; :\; \backslash Omega\_d^k(K)\; \backslash rightarrow\; \backslash Omega\_d^k(\backslash star\; K)$ defined by giving its action on cells in a prismal complex as follows,

where $|\backslash cdot|$ stands for the volume and the causality sign $\backslash kappa(\backslash sigma^k)$ is defined to be +1 if all the edges of $\backslash sigma^k$ are spacelike, and -1 otherwise.

The causality sign of 2-cells in a $(2+1)$-space-time is summarized in Table 4. We should note that the causality sign for a 0-simplex, $\backslash kappa(\backslash sigma^0)$, is always 1. This is because a 0-simplex has no edges, and as such the statement that all of its edges are spacelike is trivially true.

This causality term in the discrete Hodge star has consequences for the expression for the discrete norm (Definition 12.3), which is now given.

Definition 12.5. Given a primal discrete $k$-form $\backslash alpha^k\; \backslash in\; \backslash Omega\_d^k(K)$, its discrete $L^2$ Lorentzian norm is given by,

Having defined the discrete Lorentzian norm, we can express the discrete action as

The basic variations needed to determine the discrete Euler–Lagrange operator are obtained from variations that vary the value of the 1-form $A$ at a given 1-simplex $\backslash sigma\_0^1$, leaving the other values fixed. These variations have the form,

where $\backslash tilde\{\backslash eta\}\; \backslash in\; \backslash Omega\_d^1(K)$ is given by $\backslash langle\; \backslash tilde\{\backslash eta\},\; \backslash sigma\_0^1\; \backslash rangle\; =\; 1$ for a fixed interior $\backslash sigma\_0^1\; \backslash in\; K$ and $\backslash langle\; \backslash tilde\{\backslash eta\},\; \backslash sigma^1\; \backslash rangle\; =\; 0$ for $\backslash sigma^1\; \backslash neq\; \backslash sigma\_0^1$. The derivation of the variational principle gives

which vanishes for all the basic variations above. On the other hand, we now expand the discrete 1-form $\backslash star\; \backslash mathbf\{d\}\; \backslash star\; \backslash mathbf\{d\}A$. For any $\backslash sigma\_0^1\; \backslash in\; K$, we have that

where the sign $\backslash text\{sgn\}(\backslash sigma^2,\; \backslash sigma^1)$ stands for the relative orientation between $\backslash sigma^2$ and $\backslash sigma^1$. Which is to say, $\backslash text\{sgn\}(\backslash sigma^2,\; \backslash sigma^1)\; =\; 1$ if the orientation induced by $\backslash sigma^2$ on $\backslash sigma^1$ coincides with the orientation of $\backslash sigma^1$, and $\backslash text\{sgn\}(\backslash sigma^2,\; \backslash sigma^1)\; =\; -1$ otherwise. For the second to last equality, one has to note that the border of the dual cell of an edge $\backslash sigma\_0^1$ consists, conveniently oriented with the $\backslash text\{sgn\}$ operator, of the union of the duals of all the 2-simplices containing $\backslash sigma\_0^1$. This statement is the content of Definition 5.8, which gives the expression for the boundary of a dual cell, and was illustrated in Figure 21 for the case of $n$-dimensional dual cells.

By comparing the two computations, we find that for an arbitrary choice of $\backslash sigma\_0^1\; \backslash in\; K$, $\backslash langle\; \backslash star\; \backslash mathbf\{d\}\; \backslash star\; \backslash mathbf\{d\}A,\; \backslash sigma\_0^1\; \backslash rangle$ is equal (up to a non-zero constant) to $\backslash delta\; \backslash mathcal\{S\}\_d(A)$, which always vanishes. It follows that the variational discrete equations obtained above is equivalent to the discrete Maxwell equations,

It is desirable to leverage the exactness properties of the operators of discrete exterior calculus to construct numerical algorithms with discrete conservation properties. For these purposes, it is appropriate to extend the scope of DEC to incorporate dynamical behavior, by addressing the issue of discrete diffeomorphisms and flows.

As discussed in the previous section, DEC and discrete mechanics have interesting synergistic properties, and in this section we will explore a groupoid interpretation of discrete mechanics that is particularly appropriate to formulating the notion of pull-back and push-forward of discrete differential forms.

13.1. Groupoid Interpretation of Discrete Variational Mechanics. The groupoid formulation of discrete mechanics is particularly fruitful and natural, and it serves as a unifying tool for understanding the variational formulation of discrete Lagrangian mechanics, and discrete Euler–Poincaré reduction, as discussed in the work of Weinstein [1996] and Marsden et al. [1999, 2000].

The groupoid interpretation of discrete mechanics is most clearly illustrated if we consider the discretization of trajectories on $TQ$ in two stages. Given a curve $\backslash gamma\; :\; \backslash mathbb\{R\}^+\; \backslash rightarrow\; TQ$, we consider a discrete sampling given by

We then approximate $TQ$ by $Q\; \backslash times\; Q$, and associate to $g\_i$ two elements in $Q$. We denote this by

Or equivalently, in the language of groupoids, see Cannas da Silva and Weinstein [1999], Weinstein [2001], we have

where $\backslash alpha$ is the source map, and $\backslash beta$ is the target map. Then,

This can be visualized as

A product $\backslash cdot\; :\; G^\{(2)\}\; \backslash rightarrow\; G$ is defined on the set of composable pairs,

The groupoid composition $g\; \backslash cdot\; h$ is defined by

This can be represented graphically as follows,

The set of composable pairs is the discrete analogue of the set of second-order curves on $TQ$. A curve $\backslash gamma\; :\; \backslash mathbb\{R\}^+\; \backslash rightarrow\; TQ$ is said to be second-order if there exists a curve $q\; :\; \backslash mathbb\{R\}^+\; \backslash rightarrow\; Q$, such that,

The corresponding condition for discrete curves is that given a sequence of points in $Q\; \backslash times\; Q$, $(q\_1^0,\; q\_1^1),\; \backslash dots,\; (q\_p^0,\; q\_p^1)$, we require that

This implies that the discrete curve on $Q\; \backslash times\; Q$ is derived from a $(p+1)$-pointed curve $(q\_0,\; \backslash dots,\; q\_p)$ on $Q$, where

This condition has a direct equivalent in groupoids,

Which is to say that the sequence of points in $Q\; \backslash times\; Q$ are composable. In general, this hierarchy of sets is denoted by

where $G^\{(0)\}\; \backslash simeq\; Q$.

In addition, the groupoid inverse is defined by the following,

This is represented as follows,

Visualizing Groupoids. In summary, composition of groupoid elements, and the inverse of groupoid elements can be illustrated by Figure 23. As we will see in the next subsection, representing discrete

diffeomorphisms as pair groupoids is the natural method of ensuring that the mesh remains nondegenerate.

13.2. Discrete Diffeomorphisms and Discrete Flows. We will adopt the point of view of representing a discrete diffeomorphism as a groupoid, which was first introduced in Pekarsky and West [2003], and appropriately modify it to reflect the simplicial nature of our mesh. In addition, we will address the induced action of a discrete diffeomorphism on the dual mesh.

Definition 13.1. Given a complex $K$ embedded in $V$, and its corresponding abstract simplicial complex $M$, a discrete diffeomorphism, $\backslash varphi\; \backslash in\; \backslash text\{Diff\}\_d(M)$, is a pair of simplicial complexes $K\_1,\; K\_2$, which are realizations of $M$ in the ambient space $V$. This is denoted by $\backslash varphi(M)\; =\; (K\_1,\; K\_2)$.

Definition 13.2. A one-parameter family of discrete diffeomorphisms is a map $\backslash varphi\; :\; I\; \backslash rightarrow\; \backslash text\{Diff\}\_d(M)$, such that,

Since we are concerned with evolving equations represented by these discrete diffeomorphisms, and mesh degeneracy causes the numerics to fail, we introduce the notion of non-degenerate discrete diffeomorphisms,

Definition 13.3. A non-degenerate discrete diffeomorphism $\backslash varphi\; =\; (K\_1,\; K\_2)$ is such that $K\_1$ and $K\_2$ are non-degenerate realizations of the abstract simplicial complex $M$ in the ambient space $V$.

Notice that it is sufficient to define the discrete diffeomorphism on the vertices of the abstract complex $M$, since we can extend it to the entire complex by the relation

If $X\; \backslash in\; K^\{(0)\}$ is a material vertex of the manifold, corresponding to the abstract vertex $w$, that is to say, $\backslash pi\_1\backslash varphi\_t(w)\; =\; X,\; \backslash forall\; t\; \backslash in\; I$, the corresponding trajectory followed by $X$ in space is $x\; =\; \backslash pi\_2\backslash varphi\_t(w)$. Then, the material velocity $V(X,\; t)$ is given by

and the spatial velocity $v(x,\; t)$ is given by

The distinction between the spatial and material representation is illustrated in Figure 24.

The material velocity field can be thought of as a discrete vector field with the vectors based at the vertices of $K$, which is to say that $T\backslash varphi\_t\; \backslash in\; \backslash mathfrak\{X\}\_d(K)$, is a discrete primal vector field. Notice that $\backslash varphi\_t$ on $K$ induces a map $\backslash star\backslash varphi\_t$ on the vertices of the dual $\backslash star\; K$, by the following,

Similarly then, $T\backslash star\backslash varphi\_t\; \backslash in\; \backslash mathfrak\{X\}\_d(\backslash star\; K)$ is a discrete dual vector field.

Comparison with Interpolatory Methods. At first glance, the groupoid formulation seems like a cumbersome way to define a one-parameter family of discrete diffeomorphisms, and one may be tempted to think of extending $\backslash varphi\_t$ to the ambient space. We would then be thinking of $\backslash varphi\_t\; :\; V\; \backslash rightarrow\; V$. This is undesirable since given $\backslash varphi\_t$ and $\backslash psi\_s$ which are non-degenerate flows, their composition $\backslash varphi\_t\; \backslash circ\; \backslash psi\_s$, which is defined, may result in a degenerate mesh when applied to $K$. Thus, non-degenerate flows are not closed under this notion of composition.

If we adopt groupoid composition instead at the level of vertices, we can always be sure that if we compose two nondegenerate discrete diffeomorphisms, they will remain a nondegenerate discrete diffeomorphism.

Discrete Diffeomorphisms as Pair Groupoids. The space of discrete diffeomorphisms naturally has the structure of a pair groupoid. The discrete analogue of $T\backslash text\{Diff\}(M)$ from the point of view of temporal discretization is the pair groupoid $\backslash text\{Diff\}(M)\; \backslash times\; \backslash text\{Diff\}(M)$. In addition, we discretize $\backslash text\{Diff\}(M)$ using $\backslash text\{Diff\}\_d(M)$, which is in turn a pair groupoid involving realizations of an abstract simplicial complex in an ambient space.

13.3. Push-Forward and Pull-Back of Discrete Vector Fields and Discrete Forms. For us to construct a discrete theory of exterior calculus that admits dynamic problems, it is critical that we introduce the notion of push-forward and pull-back of discrete vector fields and discrete forms under a discrete flow.

It is sometimes desirable, particularly in the context of multigrid, multiscale, and multiresolution computations, to be able to represent a discrete differential form which is given as a cochain on a prescribed mesh, as one which is supported on a new mesh. Given a differential form $\backslash omega^k\; \backslash in\; \backslash Omega^k(K)$, and a new mesh $M$ such that $|K|\; =\; |M|$, we can define it at the level of cosimplices,

and extend this by linearity to cochains. Here, $\backslash text\{sgn\}(\backslash tau^k,\; \backslash sigma^k)$ is $+1$ if the orientation of $\backslash tau^k$ and $\backslash sigma^k$ are consistent, and $-1$ otherwise. Since $k$-skeletons of meshes that are not related by subdivision may not have nontrivial intersections, intersections of support volumes are used in the remeshing formula, as opposed to intersections of the $k$-simplices.

We denote this transformation at the level of cochains as, $T\_\{K,M\}\; :\; C^k(K)\; \backslash rightarrow\; C^k(M)$. This has the natural property that if we have a $k$-volume $U^k$ that can be represented as a chain in either the complex $K$ or the complex $M$, that is to say, $U^k\; =\; \backslash sigma\_1^k\; +\; \backslash dots\; +\; \backslash sigma\_m^k\; =\; \backslash tau\_1^k\; +\; \backslash dots\; +\; \backslash tau\_l^k$, then we have

Which is to say that the integral of the differential form over $U^k$ is well-defined, and independent of the representation of the differential form.

Note that, in particular, if we choose to coarsen the mesh, the value the form takes on a cell in the coarser mesh is simply the sum of the values the form takes on the old cells of the fine mesh which make up the new cell in the coarser mesh.

Non-Flat Manifolds. The case of non-flat manifolds presents a challenge in remeshing akin to that encountered in the discretization of differential forms. In particular, if the two meshes represent different discretizations of a non-flat manifold, they will in general correspond to different polyhedral regions in the embedding space, and not have the same support region.

We assume that our discretization of the manifold is sufficiently fine that for every simplex, all its vertices are contained in some chart. Then, by using these local charts, we can identify support volumes in the computational domain with $n$-volumes in the manifold, and thereby make sense of the remeshing formula.

We have presented a framework for discrete exterior calculus using the cochain representation of discrete differential forms, and introduced combinatorial representations of discrete analogues of differential operators on discrete forms and discrete vector fields. The role of primal and dual cell complexes in the theory are developed in detail. In addition, extensions to dynamic problems and multi-resolution computations are discussed.

In the next few paragraphs, we will describe some of the future directions that emanate from the current work on discrete exterior calculus.

Relation to Computational Algebraic Topology Since we have introduced a discrete Laplace-deRham operator, one can hope to develop a discrete Hodge-deRham theory, and relate the deRham cohomology of a simplicial complex to its simplicial cohomology.

Extensions to Non-Flat Manifolds. The intrinsic notion of what constitutes the discrete tangent space to a node on a non-flat mesh remains an open question. It is possible that this notion is related to a choice of discrete connection on the mesh, and it is an issue that deserves further exploration.

Generalization to Arbitrary Tensors. The discretization of differential forms as cochains is particularly natural, due to the pairing between forms and volumes by integration. When attempting to discretize an arbitrary tensor, the natural discrete analogue is unclear. In particular, while it is possible to expand an arbitrary tensor using the tensor product of covariant and contravariant one-tensors, this would be cumbersome to represent on a mesh. In Leok et al. [2003], which is on discrete connections, we will see Lie group-valued discrete 1-forms, and one possible method of discretizing a $(p,\; q)$-tensor that is alternating in the contravariant indices, is to consider it as a $(0,\; q)$-tensor-valued discrete $p$-form.

It would be particularly interesting to explore this in the context of the elasticity complex (see, for example, Arnold [2002]),

where $\backslash mathbb\{S\}$ is the space of $3\; \backslash times\; 3$ symmetric matrices. One approach to discretize this was suggested in Arnold [2002], which cites the use of the Bernstein–Gelfand–Gelfand resolution in Eastwood [2000] to derive the elasticity complex from the deRham complex. Alternatively, it might be appropriate in the context of the elasticity complex to consider Lie algebra-valued discrete differential forms.

Convergence and Higher-Order Theories. The natural question from the point of view of numerical analysis would be to carefully analyze the convergence properties of these discrete differential geometric operators. In addition, higher-order analogues of the discrete theory of exterior calculus are desirable from the point of view of computational efficiency, but the cochain representation is attractive due to its conceptual simplicity and the elegance of representing discrete operators as combinatorial operations on the mesh.

It would therefore be desirable to reconcile the two, by ensuring that high-order interpolation and combinatorial operations are consistent. As a low-order example, Whitney forms, which are used to interpolate differential forms on a simplicial mesh, have the nice property that taking the Whitney form associated with the coboundary of a simplicial cochain is equal to taking the exterior derivative of the Whitney form associated with the simplicial cochain. As such, the coboundary operation, which is a combinatorial operation akin to finite differences, is an exact discretization of the exterior derivative, when applied to the degrees of freedom associated to the finite-dimensional function space of Whitney forms.

It would be interesting to apply subdivision surface techniques to construct interpolatory spaces that are compatible with differential geometric operations that are combinatorial operations on the degrees of freedom. This will result in a massively simplified approach to higher-order theories of discrete exterior calculus, by avoiding the use of symbolic computation, which would otherwise be necessary to compute the action of continuous exterior differential operators on the polynomial expansions for differential forms.

158-79, COMPUTER SCIENCE, CALTECH, PASADENA, CA 91125.

E-mail address: mathieu@caltech.edu

DEPARTMENT OF COMPUTER SCIENCE, UNIVERSITY OF ILLINOIS, URBANA, IL 61801.

E-mail address: hirani@cs.uiuc.edu

DEPARTMENT OF MATHEMATICS, UNIVERSITY OF MICHIGAN, ANN ARBOR, MI 48109.

E-mail address: mleok@umich.edu

107-81, CONTROL AND DYNAMICAL SYSTEMS, CALTECH, PASADENA, CA 91125.

E-mail address: marsden@cds.caltech.edu