Latest revision as of 13:17, 28 April 2020

Abstract

We revisit the theory of Discrete Exterior Calculus (DEC) in 2D for general triangulations, relying only on Vector Calculus and Matrix Algebra. We present DEC numerical solutions of the Poisson equation and compare them against those found using the Finite Element Method with linear elements (FEML).

Keywords: Discrete element method, Poisson equation

1. Introduction

The purpose of this paper is to introduce the theory of Discrete Exterior Calculus (DEC) to the widest possible audience and, therefore, we will rely mainly on Vector Calculus and Matrix Algebra. Discrete Exterior Calculus is a relatively new method for solving partial differential equations [8] based on the idea of discretizing the mathematical theory of Exterior Differential Calculus, a theory that goes back to Cartan [3] and is fundamental in the areas of Differential Geometry and Differential Topology. Although Exterior Differential Calculus is an abstract mathematical theory, it has been introduced in various fields such as in digital geometry processing [4], numerical schemes for partial differential equations [1,8], etc.

In his PhD thesis [8], Hirani laid down the fundamental concepts of Discrete Exterior Calculus (DEC), using discrete combinatorial and geometric operations on simplicial complexes (in any dimension), proposing discrete equivalents for differential forms, vector fields, differential and geometric operators, etc. Perhaps the first numerical application of DEC to PDE was given in [9] in order to solve Darcy flow and Poisson's equation. In [7], the authors develop a modification of DEC and show that in simple cases (e.g. flat geometry and regular meshes), the equations resulting from DEC are equivalent to classical numerical schemes such as finite difference or finite volume discretizations. In [11], the authors used DEC to solve the Navier-Stokes equations and, in [5] DEC was used with a discrete lattice model to simulate elasticity, plasticity and failure of isotropic materials.

In this expository paper, we review the various operators of Exterior Differential Calculus in 2D in terms of ordinary vector calculus, and introduce only the geometrical ideas that are essential to the formulation. Among those ideas is that of duality between the differentiation operator (on vector fields) and the boundary operator (on the domain) contained in Green's theorem. This duality is one of the key ideas of the method, which justifies taking the discretized derivative matrix as the transpose of the boundary operator matrix on the given mesh. Another important ingredient is the Hodge star operator, which is hidden in the notation of Vector Calculus. In order to show the necessity of the Hodge star operator, we carry out some simple calculations. In particular, we will introduce the notion of wedge product of vectors which, roughly speaking, helps us assign algebraic objects to parallelograms and carry out algebraic manipulations with them. We present DEC in the simplest terms possible using easy examples. We also review the formulation of DEC for arbitrary meshes, which was first considered in [10]. Performance of the method is tested on the Poisson equation and compared with the Finite Element Method with linear elements (FEML).

The paper is organized as follows. In Section 2, we introduce the wedge product of vectors and the geometric Hodge star operator, and rewrite Green's theorem appropriately in order to display the duality between the differentiation and the boundary operators. In Section 3, we present the operators of DEC (mesh, dual mesh, discrete derivation, discrete Hodge star operator), showing simple examples throughout. In Section 4, we present the formulation of DEC on arbitrary triangulations. In Section 5, we present the numerical solution of a Poisson equation with DEC and FEML, in order to compare their performance. In Section 6, we present our conclusions.

2. 2D Exterior Differential Calculus as Vector Calculus

In this section we introduce two geometric operators (the wedge product and the Hodge star) and explain how to use them together with the gradient operator in order to obtain the Laplacian.

2.1 Wedge product for vectors in R²

Let ${\textstyle a,b}$ be vectors in ${\textstyle \mathbb {R} ^{2}}$. We can assign to these vector the parallelogram they span, and to such a parallelogram its area. The latter is equal to the determinant of the transformation matrix sending ${\textstyle e_{1},e_{2}}$ to ${\textstyle a,b}$ respectively. In Exterior Calculus, such a parallelogram is regarded as an algebraic object ${\textstyle a\wedge b}$, a bivector, and the set of bivectors is equipped with a vector space structure to form the one-dimensional vector space

Thus, we have the following spaces

of scalars, vectors and bivectors, respectively. For instance, the square formed by ${\textstyle e_{1}}$ and ${\textstyle e_{2}}$ is represented by the symbol

which is read as "${\textstyle e_{1}}$ wedge ${\textstyle e_{2}}$" (Figure 1).

Figure 1. Wedge product of two vectors

In ${\textstyle \mathbb {R} ^{2}}$, this represents an “element” of unit area.

Note that if we list the vectors in the opposite order, we have a different orientation and, therefore, the algebraic objects must satisfy ${\textstyle e_{2}\wedge e_{1}=-e_{1}\wedge e_{2}}$ (Figure 2).

Figure 2. Change of orientation of the parallelogram implies anticommutativity in the wedge product of two vectors

More generally, given two vectors ${\textstyle a,b\in \mathbb {R} ^{2}}$, their wedge product ${\textstyle a\wedge b}$ looks as follows (Figure 3).

Figure 3. The wedge product of the vectors ${\displaystyle a}$ and ${\displaystyle b}$

The properties of the wedge product are

• it is anticommutative: ${\textstyle a\wedge b=-b\wedge a}$ (Figure 4)

Figure 4. Anticommutativity of the wedge product

• ${\textstyle a\wedge a=0}$ since it is a parallelogram with area zero (Figure 5)

Figure 5. Parallelogram with very small area, depicting what happens when ${\displaystyle b}$ tends to ${\displaystyle a}$

• it is distributive

• it is associative

For example, let

Then

i.e. the determinant

times the canonical bivector/square ${\textstyle e_{1}\wedge e_{2}}$.

2.2 Hodge star operator for vectors in R²

Now consider the following situation: given the vector ${\textstyle e_{1}}$, find another vector ${\textstyle v}$ such that the parallelogram that they form has area ${\textstyle 1}$. It is readily seen that, for instance, ${\textstyle v=e_{2},-e_{2},e_{1}+e_{2}}$ are all solutions. Requiring orthogonality and standard orientation, we see that ${\textstyle e_{2}}$ is the unique solution. This process is summarized in the Hodge star operator, which basically says that the complementary vector for ${\textstyle e_{1}}$ is ${\textstyle e_{2}}$, and the one for ${\textstyle e_{2}}$ is ${\textstyle -e_{1}}$,

In general, the equation that defines the Hodge star operator for any given vector ${\textstyle v\in \mathbb {R} ^{2}}$ is the following

for every ${\textstyle w\in \mathbb {R} ^{2}}$. In particular, if we take ${\textstyle v=w}$,

which means that ${\textstyle v}$ and ${\textstyle \star v}$ form a square of area ${\textstyle |v|^{2}}$. Thus, the Hodge star operator on a vector ${\textstyle v\in \mathbb {R} ^{2}}$ can be thought of as finding the vector that makes with ${\textstyle v}$ a positively oriented square of area ${\textstyle \vert v\vert ^{2}}$.

It is somewhat less intuitive to work out the Hodge star of a bivector. First of all, we have to treat bivectors as vectors of a different space, namely the space of bivectors ${\textstyle \textstyle \wedge ^{2}\mathbb {R} ^{2}}$. Secondly, the length of a bivector ${\textstyle v\wedge w}$ is its area

Thus, the defining equation of the Hodge star applied to ${\textstyle (v\wedge w)}$ and ${\textstyle \star (v\wedge w)}$ reads as follows

which means

Since ${\textstyle v\wedge w={Area}(v\wedge w)e_{1}\wedge e_{2}}$ is already a bivector, ${\textstyle \star (v\wedge w)}$ must be a scalar, i.e.

When ${\textstyle v=e_{1}}$ and ${\textstyle w=e_{2}}$ we have

Finally, the Hodge star of a number ${\textstyle \lambda }$ is a bivector, i.e.

2.3 The Laplacian

Let ${\textstyle f:\mathbb {R} ^{2}\rightarrow \mathbb {R} }$ and consider the gradient

Apply the Hodge star operator to it

Now take the gradient of each coefficient function together with wedge product

By taking the Hodge star of this last expression we get the ordinary Laplacian of ${\textstyle f}$

Remarks.

(i) This rather convoluted looking way of computing the Laplacian of a function is based on Exterior Differential Calculus, a theory that generalizes the operators of vector calculus (gradient, curl and divergence) to arbitrary dimensions, and is the basis for the differential topological theory of deRham cohomology.

(ii) We would like to emphasize the necessity of using the Hodge star operator ${\textstyle \star }$ in order to make the combination of differentiation and wedge product produce the correct answer.

2.4 Duality in Green's theorem

Green's theorem states that for a vector field ${\textstyle (L,M)}$ defined on a region ${\textstyle D\subset \mathbb {R} ^{2}}$,

In this section, we will explain how this identity encodes a duality between the operator of differentiation and that of taking the boundary of the domain of integration.

2.4.1 Rewriting Green's theorem

Note that if ${\textstyle F=(L,M)}$ is a vector field, we can write it as

Then, we can apply the gradient operator together with wedge product in the following fashion

Note that we have defined a new operator ${\textstyle \nabla ^{\wedge }}$ which combines differentiation and wedge product. Applying the Hodge star operator we obtain

i.e. the integrand of the left hand side of the identity of integrals in Green's theorem. Thus, as a first step, Green's theorem can be rewritten as follows:

2.4.2 The duality in Green's theorem

Let us recall the following fact from Linear Algebra. Given a linear transformation ${\textstyle A}$ of Euclidean space ${\textstyle \mathbb {R} ^{n}}$, the transpose ${\textstyle A^{T}}$ satisfies

for any two vectors ${\textstyle v,w\in \mathbb {R} ^{n}}$, where ${\textstyle \left\langle \cdot ,\cdot \right\rangle }$ denotes the standard inner/dot product. In fact, such an identity characterizes the transpose ${\textstyle A^{T}}$.

Now, let us do the following notational trick: substitute the integration symbols in Green's theorem by ${\textstyle \ll \cdot ,\cdot \gg }$ as follows:

where ${\textstyle C=\partial D}$ is the boundary of the region. Using this notational change, Green's theorem reads as follows

Roughly speaking, this means that the differential operator ${\textstyle \nabla ^{\wedge }}$ is the transpose of the boundary operator ${\textstyle \partial }$ by means of the product ${\textstyle \ll \cdot ,\cdot \gg }$.

Remark. The previous observation is fundamental in the development of DEC, since the boundary operator is well understood and easy to calculate on meshes.

3. Discrete Exterior Calculus

Now we will discretize the differentiation operator ${\textstyle \nabla ^{\wedge }}$ presented above. We will start by describing the discrete version of the boundary operator on simplices/triangles. Afterwards, we will treat the differentiation operator as the transpose of the boundary operator.

We are interested in using certain geometric subsets of a given triangular mesh of a 2D region. Such subsets include vertices/nodes, edges/sides and faces/triangles. We will describe each one of them by means of the ordered list of vertices whose convex closure constitutes the subset of interest. For instance, consider the triangular mesh of the planar hexagonal region in Figure 6,

Figure 6. Triangular mesh of a planar hexagonal region

where the shaded triangle will be denoted by ${\textstyle [v_{0},v_{1},v_{6}]}$, and its edge joining the vertices ${\textstyle v_{0}}$ and ${\textstyle v_{1}}$ will be denoted by ${\textstyle [v_{0},v_{1}]}$. For the sake of notational consistency, we will denote the vertices also enclosed in brackets, e.g. ${\textstyle [v_{0}]}$.

3.1 Boundary operator and discrete derivative

There is a well known boundary operator ${\textstyle \partial }$ for oriented triangles, edges and points:

• For points/vertices:

Figure 7. Boundary of a vertex: ${\displaystyle \partial [v_{0}]=0}$

• For sides/edges:

Figure 8. Boundary of an edge: ${\displaystyle \partial [v_{0},v_{1}]=[v_{1}]-[v_{0}]}$

• For faces/triangles:

Figure 9. Boundary of a face: ${\displaystyle \partial [v_{0},v_{1},v_{2}]=[v_{1},v_{2}]-[v_{0},v_{2}]+[v_{0},v_{1}]}$

Example. Let us consider again the mesh of the planar hexagonal (with oriented triangles) in Figure 10.

Figure 10. Oriented triangular mesh of a planar hexagonal region

We will denote a triangle by the list of its vertices listed in order according to the orientation of the tringle. Thus, we have the following ordered lists:

• list of faces

• list of edges

• list of vertices

A key idea in DEC is to consider each face as an element of a basis of a vector space. Namely, column coordinate vectors are associated to faces as follows:

Similarly, coordinate vectors are associated to the edges

Finally, we do the same with the vertices

Now, if we take the boundary of each face, we have

which, under the previous assignments of coordinate vectors, corresponds to the linear transformation given by the following matrix which will multiply the column coordinate vectors on the left

where the subindices in ${\textstyle \partial _{2,1}}$ indicate that we are taking the boundary of 2-dimensional elements and obtaining 1-dimensional ones. Similarly, taking the boundaries of all the edges gives

which, under the previous assignments of coordinate vectors, corresponds to the linear transformation given by the following matrix

Remark. Note how these matrices encode different levels of connectivity with orientations, such as who are the edges of which oriented triangle, or which are the end points of a given oriented edge.

Due to the duality between ${\textstyle \partial }$ and ${\textstyle \nabla ^{\wedge }}$, we can define the discretiztion of ${\textstyle \nabla ^{\wedge }}$ by

For instance, we see that the operator ${\textstyle \nabla _{0,1}^{\wedge }}$ reads as follows

3.2 Dual mesh

In order to discretize the Hodge star operator, we must first introduce the notion of the dual mesh of a triangular mesh.

Consider the triangular mesh ${\textstyle K}$ in Figure 11(a). The construction of the dual mesh ${\textstyle K^{*}}$ is carried out as follows:

• The vertices of the dual mesh ${\textstyle K^{*}}$ are the circumcenters of the faces/triangles of the original mesh (the blue dots in Figures 11(b) and 11(c)). For instance, the dual of the face ${\textstyle [v_{0},v_{1},v_{6}]}$ will be denoted by ${\textstyle [v_{0},v_{1},v_{6}]^{*}}$.
• The edges of ${\textstyle K^{*}}$ are the straight line segments joining the circumcenters of two adjacent triangles (those which share an edge). Note that the resulting line segments are orthogonal to one of the original edges (the blue straight line segments in in Figures 11(b) and 11(c)). For instance, the dual of the edge ${\textstyle [v_{0},v_{6}]}$ will be denoted by ${\textstyle [v_{0},v_{6}]^{*}}$.
• The faces or cells of ${\textstyle K^{*}}$ are the areas enclosed by the new polygons determined by the new edges. For instance, the dual of the vertex ${\textstyle [v_{6}]}$ is the inner blue hexagon in in Figures 11(b) and 11(c) and will be denoted by ${\textstyle [v_{6}]^{*}}$.

(a)	(b)	(c)
Figure 11. Dual mesh construction. (a) Triangular mesh ${\displaystyle K}$. (b) Dual mesh ${\displaystyle K^{*}}$ superimposed on the mesh ${\displaystyle K}$. (c) Dual mesh ${\displaystyle K^{*}}$

The orientation of the dual edges is given by the following recipe. If we have two adjacent triangles oriented as in Figure 12(a), the dual edge crossing the edge of adjacency is oriented as in Figure 12(b).

(a)	(b)
Figure 12. (a) Two adjacent oriented triangles. (b) Compatibly oriented dual edge

3.3 Boundary operator on the dual mesh

Consider the dual mesh in Figure 13 with the given labels and orientations.

Figure 13. Oriented dual mesh

The boundary operator is applied in a similar fashion as it was applied to triangles. In this case we have

If we assign coordinate vectors to the dual faces and the dual edges as before

and