{'fraction_non_alphanumeric': 0.042996924797316184, 'fraction_numerical': 0.01181623334265213, 'mean_word_length': 3.9819251005880534, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 64, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite element exterior calculus: from Hodge theory to numerical stability

This article reports on the confluence of two streams of research, one emanating from the fields of numerical analysis and scientific computation, the other from topology and geometry. In it we consider the numerical discretization of partial differential equations that are related to differential complexes so that de Rham cohomology and Hodge theory are key tools for the continuous problem. After a brief introduction to finite element methods, the discretization methods we consider, we develop an abstract Hilbert space framework for analyzing stability and convergence. In this framework, the differential complex is represented by a complex of Hilbert spaces and stability is obtained by transferring Hodge theoretic structures from the continuous level to the discrete. We show stable discretization is achieved if the finite element spaces satisfy two hypotheses: they form a subcomplex and there exists a bounded cochain projection from the full complex to the subcomplex. Next, we consider the most canonical example of the abstract theory, in which the Hilbert complex is the de Rham complex of a domain in Euclidean space. We use the Koszul complex to construct two families of finite element differential forms, show that these can be arranged in subcomplexes of the de Rham complex in numerous ways, and for each construct a bounded cochain projection. The abstract theory therefore applies to give the stability and convergence of finite element approximations of the Hodge Laplacian. Other applications are considered as well, especially to the equations of elasticity. Background material is included to make the presentation self-contained for a variety of readers.


Introduction
Numerical algorithms for the solution of partial differential equations are an essential tool of the modern world. They are applied in countless ways every day in problems as varied as the design of aircraft, prediction of climate, development of cardiac devices, and modeling of the financial system. Science, engineering, and technology depend not only on efficient and accurate algorithms for approximating the solutions of a vast and diverse array of differential equations which arise in applications, but also on mathematical analysis of the behavior of computed solutions, in order to validate the computations, determine the ranges of applicability, compare different algorithms, and point the way to improved numerical methods. Given a partial differential equation (PDE) problem, a numerical algorithm approximates the solution by the solution of a finite dimensional problem which can be implemented and solved on a computer. This discretization depends on a parameter-representing, for example, a grid spacing, mesh size, or time stepwhich can be adjusted to obtain a more accurate approximate solution, at the cost of a larger finite dimensional problem. The mathematical analysis of the algorithm aims to describe the relationship between the true solution and the numerical solution as the discretization parameter is varied. For example, at its most basic, the theory attempts to establish convergence of the discrete solution to the true solution in an appropriate sense as the discretization parameter tends to zero.
Underlying the analysis of numerical methods for PDEs is the realization that convergence depends on consistency and stability of the discretization. Consistency, whose precise definition depends on the particular PDE problem and the type of numerical method, aims to capture the idea that the operators and data defining the discrete problem are appropriately close to those of the true problem for small values of the discretization parameter. The essence of stability is that the discrete problem is well-posed, uniformly with respect to the discretization parameter. Like well-posedness of the PDE problem itself, stability can be very elusive. One might think that well-posedness of the PDE problem, which means invertibility of the operator, together with consistency of the discretization, would imply invertibility of the discrete operator, since invertible operators between a pair of Banach spaces form an open set in the norm topology. But this reasoning is bogus. Consistency is not and cannot be defined to mean norm convergence of the discrete operators to the PDE operator, since the PDE operator, being an invertible operator between infinite dimensional spaces, is not compact and so is not the norm limit of finite dimensional operators. In fact, in the first part of the preceding century, a fundamental, and initially unexpected, realization was made: that a consistent discretization of a well-posed problem need not be stable [36,86,29]. Only for very special classes of problems and algorithms does well-posedness at the continuous level transfer to stability at the discrete level. In other situations, the development and analysis of stable, consistent algorithms can be a challenging problem, to which a wide array of mathematical techniques have been applied, just as for establishing the well-posedness of PDEs.
In this paper we will consider PDEs that are related to differential complexes, for which de Rham cohomology and Hodge theory are key tools for exploring the well-posedness of the continuous problem. These are linear elliptic PDEs, but they are a fundamental component of problems arising in many mathematical models, including parabolic, hyperbolic, and nonlinear problems. The finite element exterior calculus, which we develop here, is a theory which was developed to capture the key structures of de Rham cohomology and Hodge theory at the discrete level and to relate the discrete and continuous structures, in order to obtain stable finite element discretizations.
1.1. The finite element method. The finite element method, whose development as an approach to the computer solution of PDEs began over 50 years ago and is still flourishing today, has proven to be one of the most important technologies in numerical PDEs. Finite elements not only provide a methodology to develop numerical algorithms for many different problems, but also a mathematical framework in which to explore their behavior. They are based on a weak or variational form of the PDE problem, and fall into the class of Galerkin methods, in which the trial space for the weak formulation is replaced by a finite dimensional subspace to obtain the discrete problem. For a finite element method, this subspace is a space of piecewise polynomials defined by specifying three things: a simplicial decomposition of the domain, and, for each simplex, a space of polynomials called the shape functions, and a set of degrees of freedom for the shape functions, i.e., a basis for their dual space, with each degree of freedom associated to a face of some dimension of the simplex. This allows for efficient implementation of the global piecewise polynomial subspace, with the degrees of freedom determining the degree of interelement continuity.
Weak solutions to this problem are sought in the Sobolev space H 1 (−1, 1) consisting of functions in L 2 (−1, 1) whose first derivative also belong to L 2 (−1, 1). Indeed, the solution u can be characterized as the minimizer of the energy functional over the spaceH 1 (−1, 1) (which consists of H 1 (−1, 1) functions vanishing at ±1), or, equivalently, as the solution of the weak problem: Find u ∈H 1 (−1, 1) such that It is easily seen by integrating by parts that a smooth solution of the boundary value problem satisfies the weak formulation and that a solution of the weak formulation which posesses appropriate smoothness will be a solution of the boundary value problem.
Letting V h denote a finite dimensional subspace ofH 1 (−1, 1), called the trial space, we may define an approximate solution u h ∈ V h as the minimizer of the functional J over the trial space (the classical Ritz method), or, equivalently by Galerkin's method, in which u h ∈ V h is determined by requiring that the variation given in the weak problem hold only for functions in V h , i.e., by the equations By choosing a basis for the trial space V h , the Galerkin method reduces to a linear system of algebraic equations for the coefficients of the expansion of u h in terms of the basis functions. More specifically, if we write u h = M j=1 c j φ j , where the functions φ j form a basis for the trial space, then the Galerkin equations hold if and only if Ac = b where the coefficient matrix of the linear system is given by Since this is a square system of linear equations, it is nonsingular if and only if the only solution for f = 0 is u h = 0. This follows immediately by choosing v = u h .
The simplest finite element method is obtained applying Galerkin's method with the trial space V h consisting of all element ofH 1 (−1, 1) which are linear on each subinterval of some chosen mesh −1 = x 0 < x 1 < · · · < x N = 1 of the domain (−1, 1). Figure 1.1 compares the exact solution u = cos(πx/2) and the finite element solution u h in the case of a uniform mesh with N = 14 subintervals. The derivatives are compared as well. For this simple problem, u h is simply the orthogonal projection of u into V h with respect to the inner product defined by the left hand side of (3), and the finite element method gives good approximation even with a fairly coarse mesh. Higher accuracy can easily be obtained by using a finer mesh or piecewise polynomials of higher degree. The weak formulation (3) associated to minimization of the functional (2) is not the only variational formulation that can be used for discretization, and in more complicated situations other formulations may bring important advantages. In this simple situation, we may, for example, start by writing the differential equation This is called a mixed formulation of the boundary value problem. Note that for the mixed formulation, the Dirichlet boundary condition is implied by the formulation, as can be seen by integration by parts. Note also that in this case the solution is a saddle point of the functional I, not an extremum: I(σ, v) ≤ I(σ, u) ≤ I(τ, u) for τ ∈ H 1 (−1, 1), v ∈ L 2 (−1, 1). Although the mixed formulation is perfectly well-posed, it may easily lead to a discretization which is not. If we apply Galerkin's method with a simple choice of trial subspaces Σ h ⊂ H 1 (−1, 1) and V h ⊂ L 2 (−1, 1), we obtain a finite-dimensional linear system, which, however, may be singular, or may become increasingly unstable as the mesh is refined. This concept will be formalized and explored in the next section, but the result of such instability is clearly visible in simple computations. For example, the choice of continuous piecewise linear functions for both Σ h and V h leads to a singular linear system. The choice of continuous piecewise linear functions for Σ h and piecewise constants for V h leads to stable discretization and good accuracy. However choosing piecewise quadratics for Σ h and piecewise constants for V h gives a nonsingular system but unstable approximation (see [25] for further discussion of this example). The dramatic difference between the stable and unstable methods can be seen in Figure 1.2.
In one dimension, finding stable pairs of finite dimensional subspaces for the mixed formulation of the two-point boundary value problem is easy. For any integer r ≥ 1, the combination of continuous piecewise polynomials of degree at most r for σ and arbitrary piecewise polynomials of degree at most r−1 for u is stable as can be verified via elementary means (and which can be viewed as a very simple application of the theory presented in this paper). In higher dimensions, the problem of finding stable combinations of elements is considerably more complicated. This is discussed in Section 2.3.1 below. In particular, we shall see that the choice of continuous piecewise linear functions for σ and piecewise constant functions for u is not stable in more than one dimension. However stable element choices are known for this problem and again may be viewed as a simple application of the finite element exterior calculus developed in this paper.
1.2. The contents of this paper. The brief introduction to the finite element method just given will be continued in Section 2. In particular, there we formalize the notions of consistency and stability and establish their relation to convergence. Approximation of the mixed formulation for −u = f in one dimension with two choices of elements, piecewise constants for u and piecewise linears for σ (a stable method, shown in green), or piecewise constants for u and piecewise quadratics for σ (unstable, shown in red). The left plot shows u and the right plot shows σ, with the exact solution in blue. (In the right plot, the blue curve essentially coincides with the green curve and hence is not visible.) We shall also give several more computational examples. While seemingly simple, some of these examples may be surprising even to specialists, and they illustrate the difficulty in obtaining good methods and the need for a theoretical framework in which to understand such behaviors.
Like the theory of weak solutions of PDEs, the theory of finite element methods is based on functional analysis and takes its most natural form in a Hilbert space setting. In Section 3 of this paper, we develop an abstract Hilbert space framework which captures key elements of Hodge theory, and can be used to explore the stability of finite element methods. The most basic object in this framework is a cochain complex of Hilbert spaces, referred to as a Hilbert complex. Function spaces of such complexes will occur in the weak formulations of the PDE problems we consider, and the differentials will be differential operators entering into the PDE problem. The most canonical example of a Hilbert complex is the L 2 de Rham complex of a Riemannian manifold, but it is a far more general object with other important realizations. For example, it allows the definition of spaces of harmonic forms and the proof that they are isomorphic to the cohomology groups. A Hilbert complex includes enough structure to define an abstract Hodge Laplacian, defined from a variational problem with a saddle point structure. However, for these problems to be well-posed, we need the additional property of a closed Hilbert complex, i.e., that the range of the differentials are closed.
In this framework, the finite element spaces used to compute approximate solutions are represented by finite dimensional subspaces of the spaces in the closed Hilbert complex. We identify two key properties of these subspaces: first, they should combine to form a subcomplex of the Hilbert complex, and, second, there should exist a bounded cochain projection from the Hilbert complex to this subcomplex. Under these hypotheses and a minor consistency condition, it is easy to show that the subcomplex inherits the cohomology of the true complex, i.e., that the cochain projections induce an isomorphism from the space of harmonic forms to the space of discrete harmonic forms, and to get an error estimate on the difference between a harmonic form and its discrete counterpart. In the applications, this will be crucial for stable approximation of the PDEs. In fact, a main theme of finite element exterior calculus is that the same two assumptions, the subcomplex property and the existence of a bounded cochain projection, are the natural hypotheses to establish the stability of the corresponding discrete Hodge Laplacian, defined by the Galerkin method.
In Section 4 we look in more depth at the canonical example of the de Rham complex for a bounded domain in Euclidean space, beginning with a brief summary of exterior calculus. We interpret the de Rham complex as a Hilbert complex and discuss the PDEs most closely associated with it. This brings us to the topic of Section 5, the construction of finite element de Rham subcomplexes, which is the heart of finite element exterior calculus and the reason for its name. In this section, we construct finite element spaces of differential forms-piecewise polynomial spaces defined via a simplicial decomposition and specification of shape functions and degrees of freedom-which combine to form a subcomplex of the L 2 de Rham complex admitting a bounded cochain projection. First we construct the spaces of polynomial differential forms used for shape functions, relying heavily on the Koszul complex and its properties, and then we construct the degrees of freedom. We next show that the resulting finite element spaces can be efficiently implemented, have good approximation properties, and can be combined into de Rham subcomplexes. Finally, we construct bounded cochain projections, and, having verified the hypotheses of the abstract theory, draw conclusions for the finite element approximation of the Hodge Laplacian.
In the final two sections of the paper, we make other applications of the abstract framework. In the last section, we study a differential complex we call the elasticity complex, which is quite different from the de Rham complex. In particular, one of its differentials is a partial differential operator of second order. Via the finite element exterior calculus of the elasticity complex, we have obtained the first stable mixed finite elements using polynomial shape functions for the equations of elasticity, with important applications in solid mechanics.
1.3. Antecedents and related approaches. We now discuss some of the antecedents of finite element exterior calculus and some related approaches. While the first comprehensive view of finite element exterior calculus, and the first use of that phrase, was in the 2006 paper [8], this was certainly not the first intersection of finite element theory and Hodge theory. In 1957, Whitney [88] published his complex of Whitney forms, which is, in our terminology, a finite element de Rham subcomplex. Whitney's goals were geometric. For example, he used these forms in a proof of de Rham's theorem identifying the cohomology of a manifold defined via differential forms (de Rham cohomology) with that defined via a triangulation and cochains (simplicial cohomology). With the benefit of hindsight, we may view this, at least in principle, as an early application of finite elements to reduce the computation of a quantity of interest defined via infinite dimensional function spaces and operators, to a finite dimensional computation using piecewise polynomials on a triangulation. The computed quantities are the Betti numbers of the manifold, i.e., the dimensions of the de Rham cohomology spaces. For these integer quantities, issues of approximation and convergence do not play much of a role. The situation is different in the 1976 work of Dodziuk [39] and Dodziuk and Patodi [40], who considered the approximation of the Hodge Laplacian on a Riemannian manifold by a combinatorial Hodge Laplacian, a sort of finite difference approximation defined on cochains with respect to a triangulation. The combinatorial Hodge Laplacian was defined in [39] using the Whitney forms, thus realizing the finite difference operator as a sort of finite element approximation. A key result in [39] was a proof of some convergence properties of the Whitney forms. In [40] the authors applied them to show that the eigenvalues of the combinatorial Hodge Laplacian converge to those of the true Hodge Laplacian. This is indeed a finite element convergence result, as the authors remark. In 1978, Müller [71] further developed this work and used it to prove the Ray-Singer conjecture. This conjecture equates a topological invariant defined in terms of the Riemannian structure with one defined in terms of a triangulation, and was the original goal of [39,40]. (Cheeger [30] gave a different, independent proof of the Ray-Singer conjecture at about the same time.) Other spaces of finite element differential forms have appeared in geometry as well, especially the differential graded algebra of piecewise polynomial forms on a simplicial complex introduced by Sullivan [83,84]. Baker [15] calls these Sullivan-Whitney forms, and, in an early paper bringing finite element analysis techniques to bear on geometry, gives a numerical analysis of their accuracy for approximating the eigenvalues of the Hodge Laplacian.
Independently of the work of the geometers, during the 1970s and 1980s numerical analysts and computational engineers reinvented various special cases of the Whitney forms and developed new variants of them to use for the solution of partial differential equations on two and three-dimensional domains. In this work, naturally, implementational issues, rates of convergence, and sharp estimates played a more prominent role than in the geometry literature. The pioneering paper of Raviart and Thomas [76], presented at a finite element conference in 1975, proposed the first stable finite elements for solving the scalar Laplacian in two dimensions using the mixed formulation. The mixed formulation involves two unknown fields, the scalar-valued solution, and an additional vector-valued variable representing its gradient. Raviart and Thomas proposed a family of pairs of finite element spaces, one for each polynomial degree. As was realized later, in the lowest degree case the space they constructed for the vector-valued variable is just the space of Whitney 1-forms, while they used piecewise constants, which are Whitney 2-forms, for the scalar variable. For higher degrees, their elements are the higher order Whitney forms. In three-dimensions, the introduction of Whitney 1-and 2-forms for finite element computations and their higher degree analogues was made by Nédélec [72] in 1980, while the polynomial mixed elements which can be viewed as Sullivan-Whitney forms were introduced as finite elements by Brezzi, Douglas, and Marini [26] in 1985 in two dimensions, and then by Nédélec [72] in 1986 in three dimensions.
In 1988 Bossavit, in a paper in the IEE Transactions on Magnetics [21], made the connection between Whitney's forms used by geometers and some of the mixed finite element spaces that had been proposed for electromagnetics, inspired in part by Kotiuga's Ph.D. thesis in electrical engineering [66]. Maxwell's equations are naturally formulated in terms of differential forms, and the computational electromagnetics community developed the connection between mixed finite elements and Hodge theory in a number of directions. See, in particular, [17,37,57,58,59,70].
The methods we derive here are examples of compatible discretization methods, which means that at the discrete level they reproduce, rather than merely approximate, certain essential structures of the continuous problem. Other examples of compatible discretization methods for elliptic PDEs are mimetic finite difference methods [16,27] including covolume methods [74] and the discrete exterior calculus [38]. In these methods, the fundamental object used to discretize a differential k-form is typically a simplicial cochain, i.e., a number is assigned to each k-dimensional face of the mesh representing the integral of the k-form over the face. This is more of a finite difference, rather than finite element, point of view, recalling the early work of Dodziuk on combinatorial Hodge theory. Since the space of k-dimensional simplicial cochains is isomorphic to the space of Whitney k-forms, there is a close relationship between these methods and the simplest methods of the finite element exterior calculus. In some simple cases, the methods even coincide. In contrast to the finite element approach, these cochain-based approaches do not naturally generalize to higher order methods. Discretizations of exterior calculus and Hodge theory have also been used for purposes other than solving partial differential equations. For example, discrete forms which are identical or closely related to cochains or the corresponding Whitney forms play an important role in geometric modeling, parameterization, and computer graphics. See for example [50,54,56,87].

1.4.
Highlights of the finite element exterior calculus. We close this introduction by highlighting some of the features that are unique or much more prominent in the finite element exterior calculus than in earlier works.
• We work in an abstract Hilbert space setting that captures the fundamental structures of the Hodge theory of Riemannian manifolds, but applies more generally. In fact, the paper proceeds in two parts, first the abstract theory for Hilbert complexes, and then the application to the de Rham complex and Hodge theory and other applications. • Mixed formulations based on saddle point variational principles play a prominent role. In particular, the algorithms we use to approximate the Hodge Laplacian are based on a mixed formulation, as is the analysis of the algorithms. This is in contrast to the approach in the geometry literature, where the underlying variational principle is a minimization principle. In the case of the simplest elements, the Whitney elements, the two methods are equivalent. That is, the discrete solution obtained by the mixed finite element method using Whitney forms, is the same as obtained by Dodziuk's combinatorial Laplacian. However the different viewpoint leads naturally to different approaches to the analysis. The use of Whitney forms for the mixed formulation is obviously a consistent discretization, and the key to the analysis is to establish stability (see the next section for the terminology). However, for the minimization principle, it is unclear whether Whitney forms provide a consistent approximation, because they do not belong to the domain of the exterior coderivative, and, as remarked in [40], this greatly complicates the analysis. The results we obtain are both more easily proven and sharper. • Our analysis is based on two main properties of the subspaces used to discretize the Hilbert complex. First, they can be formed into subcomplexes, which is a key assumption in much of the work we have discussed. Second, there exist a bounded cochain projection from the Hilbert complex to the subcomplex. This is a new feature. In previous work, a cochain projection often played a major role, but it was not bounded, and the existence of bounded cochain projections was not realized. In fact, they exist quite generally (see Theorem 3.7), and we review the construction for the de Rham complex in Section 5.5. • Since we are interested in actual numerical computations, it is important that our spaces be efficiently implementable. This is not true for all piecewise polynomial spaces. As explained in the next section, finite element spaces are a class of piecewise polynomial spaces that can be implemented efficiently by local computations thanks to the existence of degrees of freedom, and the construction of degrees of freedom and local bases is an important part of the finite element exterior calculus. • For the same reason, high order piecewise polynomials are of great importance, and all the constructions and analysis of finite element exterior calculus carries through for arbitrary polynomial degree. • A prominent aspect of the finite element exterior calculus is the role of two families of spaces of polynomial differentials form, P r Λ k and P − r Λ k , where the index r ≥ 1 denotes the polynomial degree and k ≥ 0 the form degree. These are the shape functions for corresponding finite element spaces of differential k-forms which include, as special cases, the Lagrange finite element family, and most of the stable finite element spaces that have been used to define mixed formulations of the Poisson or Maxwell's equations. The space P − 1 Λ k is the classical space of Whitney k-forms. The finite element spaces based on P r Λ k are the spaces of Sullivan-Whitney forms. We show that for each polynomial degree r, there are 2 n−1 ways to form these spaces in de Rham subcomplexes for a domain in n dimensions. The unified treatment of the spaces P r Λ k and P − r Λ k , particularly their connections via the Koszul complex, is new to the finite element exterior calculus.
The finite element exterior calculus unifies many of the finite element methods that have been developed for solving PDEs arising in fluid and solid mechanics, electromagnetics, and other areas. Consequently, the methods developed here have been widely implemented and applied in scientific and commercial programs such as GetDP [42], FEniCS [44], EMSolve [45], deal.II [46], Diffpack [61], Getfem++ [77], and NGSolve [78]. We also note that, as part of a recent programming effort connected with the FEniCS project, Logg and Mardal [69] have implemented the full set of finite element spaces developed in this paper strictly following the finite element exterior framework as laid out here and [8].

Finite element discretizations
In this section we continue the introduction to the finite element method begun above. We move beyond the case of one dimension, and consider not only the formulation of the method, but also its analysis. To motivate the theory developed later in this paper, we present further examples that illustrate how for some problems, even rather simple ones, deriving accurate finite element methods is not a straightforward process.
2.1. Galerkin methods and finite elements. We consider first a simple problem, which can be discretized in a straightforward way, namely the Dirichlet problem for Poisson's equation in a polyhedral domain Ω ⊂ R n : This is the generalization to n dimensions of the problem (1) discussed in the introduction, and the solution may again be characterized as the minimizer of an energy functional analogous to (2) or as the solution of a weak problem analogous to (3). This leads to discretization just as for the one-dimensional case, by choosing a trial space V h ⊂H 1 (Ω) and defining the approximate solution u h ∈ V h by Galerkin's method: As in one dimension, the simplest finite element method is obtained by using the trial space consisting of all piecewise linear functions with respect to a given simplicial triangulation of the domain Ω, which are continuous and vanish on ∂Ω.
A key to the efficacy of this finite element method is the existence of a basis for the trial space consisting of functions which are locally supported, i.e., vanish on all but a small number of the elements of the triangulation. See Figure 2.1. Because of this, the coefficient matrix of the linear system is easily computed and is sparse, and so the system can be solved efficiently. More generally, a finite element method is a Galerkin method for which the trial space V h is a space of piecewise polynomial functions which can be obtained by what is called the finite element assembly process. This means that the space can be defined by specifying the triangulation T h and, for each element T ∈ T h , a space of polynomial functions on T called the shape functions, and a set of degrees of freedom. By degrees of freedom on T , we mean a set of functionals on the space of shape functions, which can be assigned values arbitrarily to determine a unique shape function. In other words, the degrees of freedom form a basis for the dual space of the space of shape functions. In the case of piecewise linear finite elements, the shape functions are of course the linear polynomials on T , a space of dimension n + 1, and the degrees of freedom are the n + 1 evaluation functionals p → p(x), where x varies over the vertices of T . For the finite element assembly process, we also require that each degree of freedom be associated to a face of some dimension of the simplex T . For example, in the case of piecewise linear finite elements, the degree of freedom p → p(x) is associated to the vertex x. Given the triangulation, shape functions, and degrees of freedom, the finite element space V h is defined as the set of functions on Ω (possibly multivalued on the element boundaries) whose restriction to any T ∈ T h belongs to the given space of shape functions on T , and for which the degrees of freedom are single-valued in the sense that when two elements share a common face, the corresponding degrees of freedom take on the same value. Returning again to the example of piecewise linear functions, V h is the set of functions which are linear polynomials on each element, and which are single-valued at the vertices. It is easy to see that this is precisely the space of continuous piecewise linear functions, which is a subspace of H 1 (Ω). As another example, we could take the shape functions on T to be the polynomials of degree at most 2, and take as degrees of freedom the functions p → p(x), x a vertex of T , and p → e p ds, e an edge of T . The resulting assembled finite element space is the space of all continuous piecewise quadratics. The finite element assembly process insures the existence of a computable locally supported basis, which is important for efficient implementation.
2.2. Consistency, stability, and convergence. We now turn to the important problem of analyzing the error in the finite element method. To understand when a Galerkin method will produce a good approximation to the true solution, we introduce the standard abstract framework. Let V be a Hilbert space, B : V × V → R a bounded bilinear form, and F : V → R a bounded linear form. We assume the problem to be solved can be stated in the form: Find u ∈ V such that This problem is called well-posed if for each F ∈ V * , there exists a unique solution u ∈ V and the mapping F → u is bounded, or, equivalently, if the operator L : V → V * given by Lu, v = B(u, v) is an isomorphism. For the Dirichlet problem for Poisson's equation, A generalized Galerkin method for the abstract problem begins with a finitedimensional normed space V h (not necessarily a subspace of V ), a bilinear form A Galerkin method is the special case of a generalized Galerkin method for which V h is a subspace of V and the forms B h and F h are simply the restrictions of the forms B and F to the subspace. The more general setting is important since it allows the incorporation of additional approximations, such as numerical integration to evaluate the integrals, and also allows for situations in which V h is not a subspace of V . Although we do not treat approximations such as numerical integration in this paper, for the fundamental discretization method we study, namely the mixed method for the abstract Hodge Laplacian introduced in Section 3.4, the trial space V h is not a subspace of V , since it involves discrete harmonic forms which will not, in general, belong to the space of harmonic forms. The generalized Galerkin method (6) may be written If the finite-dimensional problem is nonsingular, then we define the norm of the discrete solution operator, , as the stability constant of the method.
Of course, in approximating the original problem determined by V , B, and F , by the generalized Galerkin method given by V h , B h , and F h , we intend that the space V h in some sense approximates V and that the discrete forms B h and F h in some sense approximate B and F . This is the essence of consistency. Our goal is to prove that the discrete solution u h approximates u in an appropriate sense (convergence). In order to make these notions precise, we need to compare a function in V to a function in V h . To this end, we suppose that there is a restriction operator π h : V → V h , so that π h u is thought to be close to u. Then the consistency error is simply L h π h u − F h and the error in the generalized Galerkin method which we wish to control is π h u − u h . We immediately get a relation between the error and the consistency error , and so the norm of the error is bounded by the product of the stability constant and the norm of the consistency error: Stated in terms of the bilinear form B h , the norm of the consistency error can be written As for stability, the finite dimensional problem is nonsingular if and only if and the stability constant is then given by γ −1 h . Often we consider a sequence of spaces V h and forms B h and F h where we think of h > 0 as an index accumulating at 0. The corresponding generalized Galerkin method is consistent if the V h norm of the consistency error tends to zero with h and it is stable if the stability constant γ −1 h is uniformly bounded. For a consistent, stable generalized Galerkin method, π h u − u h V h tends to zero, i.e., the method is convergent.
In the special case of a Galerkin method, we can bound the consistency error In this case it is natural to choose the restriction π h to be the orthogonal projection onto V h , and so the consistency error is bounded by the norm of the bilinear form times the error in the best approximation of the solution. Thus we obtain Combining this with the triangle inequality, we obtain the basic error estimate for Galerkin methods (In fact, in this Hilbert space setting, the quantity in parentheses can be replaced with γ −1 h B , see [89].) Note that a Galerkin method is consistent as long as the sequence of subspaces V h is approximating in V in the sense that A consistent, stable Galerkin method converges, and the approximation given by the method is quasioptimal, i.e., up to multiplication by a constant, it is as good as the best approximation in the subspace. In practice, it can be quite difficult to show that the finite dimensional problem is nonsingular and to bound the stability constant, but there is one important case in which it is easy. Namely when the form B is coercive, i.e., there is a positive constant α for which , v ∈ V, and so γ h ≥ α. The bilinear form (5) for Poisson's equation is coercive, as follows from Poincaré's inequality. This explains, and can be used to prove, the good convergence behavior of the method depicted in Figure 1.1.

Computational examples.
2.3.1. Mixed formulation of the Laplacian. For an example of a problem that fits in the standard abstract framework with a noncoercive bilinear form, we take the mixed formulation of the Dirichlet problem for Poisson's equation, already introduced in one dimension in Section 1.1. Just as there, we begin by writing Poisson's equation as the first order system The pair (σ, u) can again be characterized variationally as the unique critical point (a saddle point) of the functional over H(div; Ω) × L 2 (Ω) where H(div; Ω) = {σ ∈ L 2 (Ω) : div σ ∈ L 2 (Ω)}. Equivalently, it solves the weak problem: Find σ ∈ H(div; Ω), u ∈ L 2 (Ω) satisfying This mixed formulation of Poisson's equation fits in the abstract framework if we define V = H(div; Ω) × L 2 (Ω), In this case the bilinear form B is not coercive, and so the choice of subspaces and the analysis is not so simple as for the standard finite element method for Poisson's equation, a point we already illustrated in the one-dimensional case.
Finite element discretizations based on such saddle point variational principles are called mixed finite element methods. Thus a mixed finite element for Poisson's equation is obtained by choosing subspaces Σ h ⊂ H(div; Ω) and V h ⊂ L 2 (Ω) and seeking a critical point of I over Σ h × V h . The resulting Galerkin method has the form: This again reduces to a linear system of algebraic equations.
Since the bilinear form is not coercive, it is not automatic that the linear system is nonsingular, i.e., that for f = 0, the only solution is σ h = 0, u h = 0. Choosing τ = σ h and v = u h and adding the discretized variational equations, it follows immediately that when f = 0, σ h = 0. However, u h need not vanish unless the condition that Ω u h div τ dx = 0 for all τ ∈ Σ h implies that u h = 0. In particular, this requires that dim(div Σ h ) ≥ dim V h . Thus, even nonsingularity of the approximate problem depends on a relationship between the two finite dimensional spaces. Even if the linear system is nonsingular, there remains the issue of stability, i.e., a uniform bound on the inverse operator.
As mentioned earlier, the combination of continous piecewise linear elements for σ and piecewise constants for u is not stable in two dimensions. The simplest stable elements use the piecewise constants for u, and the lowest order Raviart-Thomas elements for σ. These are finite elements defined with respect to a triangular mesh by shape functions of the form (a + bx 1 , c + bx 2 ) and one degree of freedom for each edge e, namely σ → e σ · n ds. We show in Fig. 2.2 below two numerical computations that demonstrate the difference between an unstable and stable choice of elements for this problem. The stable method accurately approximates the true solution u = x(1 − x)y(1 − y) on (0, 1) × (0, 1) with a piecewise constant, while the unstable method is wildly oscillatory.  This problem is a special case of the Hodge Laplacian with k = n as discussed briefly in Section 4.2; see especially Section 4.2.4. The error analysis for a variety of finite element methods for this problem, including the Raviart-Thomas elements, is thus a special case of the general theory of this paper, yielding the error estimates in Section 5.6.

2.3.2.
The vector Laplacian on a nonconvex polygon. Given the subtlety of finding stable pairs of finite element spaces for the mixed variational formulation of Poisson's equation, we might choose to avoid this formulation, in favor of the standard formulation, which leads to a coercive bilinear form. However, while the standard formulation is easy to discretize for Poisson's equation, additional issues arise already if we try to discretize the vector Poisson equation. For a domain Ω in R 3 with unit outward normal n, this is the problem (10) − grad div u + curl curl u = f, in Ω, u · n = 0, curl u × n = 0, on ∂Ω.
The solution of this problem can again be characterized as the minimizer of an appropriate energy functional, but this time over the space H(curl; Ω) ∩H(div; Ω), where H(curl; Ω) = {u ∈ L 2 (Ω) | curl u ∈ L 2 (Ω)} andH(div; Ω) = {u ∈ H(div; Ω) | u · n = 0 on ∂Ω} with H(div; Ω) defined above. This problem is associated to a coercive bilinear form, but a standard finite element method based on a trial subspace of the energy space H(curl; Ω) ∩H(div; Ω), e.g., using continuous piecewise linear vector functions, is very problematic. In fact, as we shall illustrate shortly, if the domain Ω is a nonconvex polyhedron, for almost all f the Galerkin method solution will converge to a function that is not the true solution of the problem! The essence of this unfortunate situation is that any piecewise polynomial subspace of H(curl; Ω) ∩ H(div; Ω) is a subspace of H 1 (Ω) ∩H(div; Ω), and this space is a closed subspace of H(curl; Ω) ∩H(div; Ω). For a nonconvex polyhedron, it is a proper closed subspace and generally the true solution will not belong to it, due to a singularity at the reentrant corner. Thus the method, while stable, is inconsistent. For more on this example, see [35]. An accurate approximation of the vector Poisson equation can be obtained from a mixed finite element formulation, based on the system: σ = − div u, grad σ + curl curl u = f in Ω, u · n = 0, curl u × n = 0 on ∂Ω.
Writing this system in weak form, we obtain the mixed formulation of the problem: Find σ ∈ H 1 (Ω), u ∈ H(curl; Ω) satisfying In contrast to a finite element method based on minimizing the energy (11), a finite element approximation based on the mixed formulation uses separate trial subspaces of H 1 (Ω) and H(curl; Ω), rather than a single subspace of the intersection H(curl; Ω) ∩H(div; Ω).
We now illustrate the nonconvergence of a Galerkin method based on energy minimization and the convergence of one based on the mixed formulation, via computations in two space dimensions (so now the curl of a vector u is the scalar ∂u 2 /∂x 1 − ∂u 1 /∂x 2 ). For the trial subspaces we make the simplest choices: for the former method we use continuous piecewise linear functions and for the mixed method we use continuous piecewise linear functions to approximate σ ∈ H 1 (Ω) and a variant of the lowest order Raviart-Thomas elements, for which the shape functions are the infinitesimal rigid motions (a − bx 2 , c + bx 1 ) and the degrees of freedom are the tangential moments u → e u · s ds for e an edge. The discrete solutions obtained by the two methods for the problem when f = (−1, 0) are shown in Figure 2.3. As we shall show later in this paper, the mixed formulation gives an approximation that provably converges to the true solution, while, as can be seen from comparing the two plots, the first approximation scheme gives a completely different (and therefore inaccurate) result. This problem is again a special case of the Hodge Laplacian, now with k = 1. See Section 4.2.2. The error analysis thus falls within the theory of this paper, yielding estimates as in Section 5.6.

2.3.3.
The vector Laplacian on an annulus. In the example just considered, the failure of a standard Galerkin method based on energy minimization to solve the vector Poisson equation was related to the reentrant corner of the domain and the resulting singular behavior of the solution. A quite different failure mode for this method occurs if we take a domain which is smoothly bounded, but not simply connected, e.g., an annulus. In that case, as discussed below in Section 3.2, the boundary value problem (10) is not well-posed except for special values of the forcing function f . In order to obtain a well-posed problem, the differential equation should be solved only modulo the space of harmonic vector fields (or harmonic 1forms)-which is a one-dimensional space for the annulus-and the solution should be rendered unique by enforcing orthogonality to the harmonic vector fields. If we choose the annulus with radii 1/2 and 1, and forcing function f = (0, x), the resulting solution, which can be computed accurately with a mixed formulation falling within the theory of this paper, is displayed on the right in Figure 2  The true solution shown here on the right in an (accurate) approximation by a mixed method. It is orthogonal to the harmonic fields, and satisfies the differential equation only modulo harmonic fields. The standard Galerkin solution using continuous piecewise linear vector fields shown on the left, is totally different.

2.3.4.
The Maxwell eigenvalue problem. Another situation where a standard finite element method gives unacceptable results, but a mixed method succeeds, arises in the approximation of elliptic eigenvalue problems related to the vector Laplacian or Maxwell's equation. This will be analyzed in detail later in this paper, and here we only present a simple but striking computational example. Consider the eigenvalue problem for the vector Laplacian discussed above, which we write in mixed form as: Find nonzero (σ, u) ∈ H 1 (Ω) × H(curl; Ω) and λ ∈ R satisfying As explained in Section 3.6.1, this problem can be split into two subproblems. In particular, if 0 = u ∈ H(curl; Ω) and λ ∈ R solves the eigenvalue problem and λ is not equal to zero, then (σ, u), λ is an eigenpair for (12) with σ = 0.
We now consider the solution of the eigenvalue problem (13), with two different choices of trial subspaces in H(curl; Ω). Again, to make our point it is enough to consider a two dimensional case, and we consider the solution of (13) with Ω a square of side length π. For this domain, the positive eigenvalues can be computed by separation of variables. They are of the form m 2 + n 2 with m and n integers: 1, 1, 2, 4, 4, 5, 5, 8, . . .. If we approximate (13) using the space of continuous piecewise linear vector fields as the trial subspace of H(curl; Ω), the approximation fails badly. This is shown for an unstructured mesh in Figure 2.5 and for a structured crisscross mesh in Figure 2.6, where the nonzero discrete eigenvalues are plotted. Note the very different mode of failure for the two mesh types. For more discussion of the spurious eigenvalues arising using continuous piecewise linear vector fields on a crisscross mesh see [20]. By contrast, if we use the lowest order Raviart-Thomas approximation of H(curl; Ω), as shown on the right of Figure 2.5, we obtain a provably good approximation for any mesh. This is a very simple case of the general eigenvalue approximation theory presented in Section 3.6 below.  (13) on an unstructured mesh of the square (left) using continuous piecewise linear finite elements (middle) and Raviart-Thomas elements (right). For the former, the discrete spectrum looks nothing like the true spectrum, while for the later it is very accurate. Approximation of the nonzero eigenvalues of (13) using continuous piecewise linear elements on the structured mesh shown. The first seven discrete nonzero eigenvalues converge to true eigenvalues, but the eighth converges to a spurious value.

Hilbert complexes and their approximation
In this section, we construct a Hilbert space framework for finite element exterior calculus. The most basic object in this framework is a Hilbert complex, which extracts essential features of the L 2 de Rham complex. Just as the Hodge Laplacian is naturally associated with the de Rham complex, there is a system of variational problems, which we call the abstract Hodge Laplacian, associated to any Hilbert complex. Using a mixed formulation we prove that these abstract Hodge Laplacian problems are well-posed. We next consider approximation of Hilbert complexes using finite dimensional subspaces. Our approach emphasizes two key properties, the subcomplex property and the existence of bounded cochain projections. These same properties prove to be precisely what is needed both to show that the approximate Hilbert complex accurately reproduces geometrical quantities associated to the complex, like cohomology spaces, and also to obtain error estimates for the approximation of the abstract Hodge Laplacian source and eigenvalue problems, which is our main goal in this section. In the following section of the paper we will derive finite element subspaces in the concrete case of the de Rham complex and verify the hypotheses needed to apply the results of this section.
Although the L 2 de Rham complex is the canonical example of a Hilbert complex, there are many others. In this paper, in Section 6 we consider some variations of de Rham complex that allow us to treat more general PDEs and boundary value problems. In the final section we briefly discuss the equations of elasticity, for which a very different Hilbert complex, in which one of the differentials is a second order PDE, is needed. Another useful feature of Hilbert complexes, is that a subcomplex of Hilbert complex is again such, and so the properties we establish for them apply not only at the continuous, but also at the discrete level.
3.1. Basic definitions. We begin by recalling some basic definitions of homological algebra and functional analysis and establishing some notation.
3.1.1. Cochain complexes. Consider a cochain complex (V, d) of vector spaces, i.e., a sequence of vector spaces V k and linear maps d k , called the differentials: Equivalently, we may think of such a complex as the graded vector space V = V k , equipped with a graded linear operator d : V → V of degree +1 satisfying d • d = 0. A chain complex differs from a cochain complex only in that the indices decrease. All the complexes we consider are nonnegative and finite, meaning that V k = 0 whenever k is negative or sufficiently large.
Given a cochain complex (V, d), the elements of the null space Z k = Z k (V, d) of d k are called the k-cocycles and the elements of the range B k = B k (V, d) of d k−1 the k-coboundaries. The kth cohomology group is defined to be the quotient space Given two cochain complexes (V, d) and (V , d ), a set of linear maps f k : is called a cochain map. When f is a cochain map, f k maps k-cocycles to k-cocycles and k-coboundaries to k-coboundaries, and hence induces a mapf on the cohomology spaces. This map is functorial, i.e., it respects composition. Let Then the inclusion i h : V h → V is a cochain map and so induces a map of cohomology. If there exists a cochain projection of and similarly for B k h ). We conclude that in this caseī h is injective andπ h is surjective. In particular, the dimension of the cohomology spaces of the subcomplex is at most that of the supercomplex.

3.1.2.
Closed operators on Hilbert space. This material can be found in many places, e.g., [ By an operator T from a Hilbert space X to a Hilbert space Y , we mean a linear operator from a subspace V of X, called the domain of T , into Y . The operator T is not necessarily bounded and the domain V is not necessarily closed. We say that the operator T is closed if its graph It is easy to check that this makes V a Hilbert space (i.e., complete), if and only if T is closed, and moreover, that T is a bounded operator from V to Y . Of course, the null space of T is the set of elements of its domain which it maps to 0, and the range of T is T (V ). The null space of a closed operator from X to Y is a closed subspace of X, but its range need not be closed in Y (even if the operator is defined on all of X and is bounded).
The operator T is said to be densely defined if its domain V is dense in X. In this case the adjoint operator T * from Y to X is defined to be the operator whose domain consists of all y ∈ Y for which there exists x ∈ X with in which case T * y = x (well-defined since V is dense). If T is closed and densely defined, then T * is as well and T * * = T . Moreover, the null space of T * is the orthogonal complement of the range of T in Y . Finally, by the closed range theorem, the range of T is closed in Y if and only if the range of T * is closed in X.
If the range of a closed linear operator is of finite codimension, i.e., dim Y /T (V ) < ∞, then the range is closed [60, Lemma 19.1.1].

Hilbert complexes.
A Hilbert complex is a sequence of Hilbert spaces W k and closed, densely-defined linear operators d k from W k to W k+1 such that the range of d k is contained in the domain of d k+1 and d k+1 • d k = 0. A Hilbert complex is bounded if, for each k, d k is a bounded linear operator from W k to W k+1 . In other words, a bounded Hilbert complex is a cochain complex in the category of Hilbert spaces. A Hilbert complex is closed if for each k, the range of d k is closed in W k+1 . A Fredholm complex is a Hilbert complex for which the range of d k is finite codimensional in the null space of d k+1 (and so is closed). Hilbert and Fredholm complexes have been discussed by various authors working in geometry and topology. Brüning and Lesch [28] have advocated for them as an abstraction of elliptic complexes on manifolds and applied them to spectral geometry on singular spaces. Glotko [53] used them to define a generalization of Sobolev spaces on Riemannian manifolds and study their compactness properties and Gromov and Shubin [55] to define topological invariants of manifolds.
Associated to any Hilbert complex (W, d) is a bounded Hilbert complex (V, d), called the domain complex, for which the space V k is the domain of d k , endowed with the inner product associated to the graph norm: Then d k is a bounded linear operator from V k to V k+1 , and so (V, d) is indeed a bounded Hilbert complex. The domain complex is closed or Fredholm if and only if the original complex (W, d) is.
Of course, for a Hilbert complex (W, d), we have the null spaces and ranges Z k and B k . Utilizing the inner product, we define the space of harmonic forms H k = Z k ∩ B k⊥ , the orthogonal complement of B k in Z k . It is isomorphic to the reduced cohomology space Z k /B k or, for a closed complex, to the cohomology space Z k /B k . For a closed Hilbert complex, we immediately obtain the Hodge decomposition For the domain complex (V, d), the null space, range, and harmonic forms are the same spaces as for the original complex, and the Hodge decomposition is The third summand, Continuing to use the Hilbert space structure, we define the dual complex (W, d * ), which is a Hilbert chain complex rather than cochain complex. The dual complex uses the same spaces W i , with the differential d * k being the adjoint of d k−1 , so d * k is a closed, densely-defined operator from W k to W k−1 , whose domain we denote by V * k . The dual complex is closed or bounded if and only if the original complex is.
k is the space of harmonic forms both for the original complex and the dual complex. Since Z k⊥ W = B * k , the Hodge decomposition (14) can be written Then d k is a bounded bijection from Z k⊥ to the Hilbert space B k+1 and hence, by Banach's bounded inverse theorem, there exists a constant c P such that which we refer to as a Poincaré inequality. We remark that the condition that B k is closed is not only sufficient to obtain (16), but also necessary.
, and is continuously included in W k . We say that the Hilbert complex (W, d) has the compactness property if V k ∩ V * k is dense in W k and the inclusion is a compact operator. Restricted to the space H k of harmonic forms, the V k ∩ V * k norm is equal to the W k norm (times √ 2). Therefore the compactness property implies that the inclusion of H k into itself is compact, so H k is finite dimensional. In summary, for Hilbert complexes, compactness property =⇒ Fredholm =⇒ closed.
3.2. The abstract Hodge Laplacian and the mixed formulation. Given a Hilbert complex (W, d), the operator L = dd * + d * d is an unbounded operator W k → W k called, in the case of the de Rham complex, the Hodge Laplacian. We refer to it as the abstract Hodge Laplacian in the general situation. Its domain is Note that, in this equation, and henceforth, we use · , · and · without subscripts, meaning the inner product and norm in the appropriate W k space.
The harmonic functions measure the extent to which the Hodge Laplace problem (17) is well-posed. The solutions to the homogeneous problem (f = 0) are precisely the functions in H k . Moreover, a necessary condition for a solution to exist for For computational purposes, a formulation of the Hodge Laplacian based on (17) may be problematic, even when there are no harmonic forms, because it may not be possible to construct an efficient finite element approximation for the space V k ∩V * k . We have already seen an example of this in the discussion of the approximation of a boundary value problem for the vector Laplacian in Section 2.3.2. Instead we introduce another formulation, which is a generalization of the mixed formulation discussed in Section 2 and which, simultaneously, accounts for the nonuniqueness associated with harmonic forms. With (W, d) a Hilbert complex, (V, d) the associated domain complex, and f ∈ W k given, we define the mixed formulation of the abstract Hodge Laplacian as the problem of finding (σ, u, p) Remark. The equations (18) are the Euler-Lagrange equations associated to a variational problem. Namely, if we define the quadratic functional I : Thus the critical point is a saddle point.
An important result is that if the Hilbert complex is closed, then the mixed formulation is well-posed. The requirement that the complex is closed is crucial, since we rely on the Poincaré inequality.
The mixed formulation of the abstract Hodge Laplacian is well-posed. That is, for We shall prove Theorem 3.1 in Section 3.2.2. First, we interpret the mixed formulation.
3.2.1. Interpretation of the mixed formulation. The first equation states that u belongs to the domain of d * and d * u = σ ∈ V k−1 . The second equation similarly states that du belongs to the domain of d * and d * du = f − p − dσ. Thus u belongs to the domain D L of L and solves the abstract Hodge Laplacian equation The harmonic form p is simply the orthogonal projection P H f of f onto H k , required for existence of a solution. Finally the third equation fixes a particular solution, through the condition u ⊥ H k . Thus Theorem 3.1 establishes that for any f ∈ W k there is a unique u ∈ D L such that Lu = f − P H f and u ⊥ H k . We define Kf = u, so the solution operator K : W k → W k is a bounded linear operator mapping into D L . The solution to the mixed formulation is The mixed formulation (18) is also intimately connected to the Hodge decomposition (15). Since dσ ∈ B k , p ∈ H k , and d * du ∈ B * k , the expression f = dσ+p+d * du is precisely the Hodge decomposition of f . In other words where P B and P B * are the W k -orthogonal projections onto B k and B * k , respectively. We also note that K commutes with d and d * in the sense that and both du and df are orthogonal to harmonic forms. This establishes that du = Kdf , i.e., dKf = Kdf . The second equation is established similarly.
If we restrict the data f in the abstract Hodge Laplacian to an element of B * k or of B k , we get two other problems which are also of great use in applications.
The solution u can be characterized as the unique element of Z k⊥ such that and any solution to this problem is a solution of Lu = f , and so is uniquely determined.
and any solution to this problem is a solution of Lu = f , σ = d * u, and so is uniquely determined.

3.2.2.
Well-posedness of the mixed formulation. We now turn to the proof of Theorem 3.1. Let B : X × X → R be a symmetric bounded bilinear form on a Hilbert space X which satisfies the inf-sup condition Then the problem of finding x ∈ X such that B(x, y) = F (y) for all y ∈ X is well-posed: it has a unique solution x for each F ∈ X * , and x X ≤ γ −1 F X * [13]. The abstract Hodge Laplacian problem (18) is of this form, where B : The following theorem establishes the inf-sup condition, and so implies Theorem 3.1.
There exists a constant γ > 0, depending only on the constant c P in the Poincaré inequality (16), such that for any (σ, u, p) Proof. By the Hodge decomposition, where c P ≥ 1 is the constant in Poincareś inequality. Let From (21) and the orthogonality of the Hodge decomposition, we have We also get, from a simple computation using (21) and (22), that The theorem easily follows from this bound and (23).
We close this section with two remarks. First we note that in fact Theorem 3.2 establishes more than the well-posedness of the problem (18) stated in Theorem 3.1. It establishes that, for any G ∈ (V k−1 ) * , F ∈ (V k ) * , and R ∈ (H k ) * (these are the dual spaces furnished with the dual norms), there exists a unique (σ, u, p) and moreover the correspondence (σ, u, p) Second, we note that the above result bears some relation to a fundamental result in the theory of mixed finite element methods, due to Brezzi [24], which we state here. Theorem 3.3. Let X and Y be Hilbert spaces and a : and suppose that there exists positive constants α and γ such that with the constant c depending only on α, γ, and the norms of the bilinear forms a and b.
If we make the additional assumption (usually satisfied in applications of this theorem), that the bilinear form a is symmetric and satisfies a(x, x) > 0 for all 0 = x ∈ X, then this theorem can be viewed as a special case of Theorem 3.2. In fact, we define W 0 as the completion of X in the inner product given by a, let W 1 = Y , and define d as the closed linear operator from W 0 to W 1 with domain X given by dx, y Y = b(x, y), x ∈ X, y ∈ Y. In this way we obtain a Hilbert complex (with just two spaces W 0 and W 1 ). The inf-sup condition of Theorem 3.3 implies that d has closed range, so it is a closed Hilbert complex, and so Poincaré's inequality holds. The associated abstract Hodge Laplacian is just the system (24), and from Theorem 3.2 and the first remark above, we have . But, using the coercivity in the kernel, the decomposition X = Z + Z ⊥ , and Poincaré's inequality, we get u X ≤ C[ a(u, u)+ du Y ], which gives the estimate from Brezzi's theorem. Finally, we mention that we could dispense with the extra assumption about the symmetry and positivity of the bilinear form a, but this would require a slight generalization of Theorem 3.2 which we do not consider here.

3.3.
Approximation of Hilbert complexes. The remainder of this section will be devoted to approximation of quantities associated to a Hilbert complex, such as the cohomology spaces, harmonic forms, and solutions to the Hodge Laplacian, by quantities associated to a subcomplex.
Let (W, d) be a Hilbert complex with domain complex (V, d), and suppose we We also take W k h to be the same subspace V k h but endowed with the norm of W k . In this way we obtain a closed (even bounded) Hilbert complex (W h , d) with domain complex (V h , d), and all the results of Sections 3.1.3 and 3.2 apply to this subcomplex. Although the differential for the subcomplex is just the restriction of d, and so does not need a new notation, with the obvious meanings. We use the term discrete when we wish to emphasize quantities associated to the subcomplex. For example, we refer to H k h as the space of discrete harmonic k-forms, and the discrete Hodge decomposition is In order that the subspaces V k h can be used effectively to approximate quantities associated to the original complex, we require not only that they form a subcomplex, but of course need to know something about the approximation of V k by V k h , i.e., an assumption that inf v∈V k h u − v V is sufficiently small for some or all u ∈ V k . A third assumption, which plays an essential role in our analysis, is that there exists a bounded cochain projection π h from the complex ( In other words, we have the following commuting diagram relating the complex (V, d) to the subcomplex (V h , d): Note that a bounded projection gives quasioptimal approximation: We now present two results indicating that, under these assumptions, the space H k h of discrete harmonic forms provides a faithful approximation of H k . In the first result we show that a bounded cochain projection into a subcomplex of a bounded closed Hilbert complex which satisfies a rather weak approximability assumption (namely (25) below), induces, not only a surjection, but an isomorphism on cohomology.
Hilbert subcomplex, and π h a bounded cochain projection. Suppose that for all k Then the induced map on cohomology is an isomorphism.
Proof. We already know that the induced map is a surjection, so it is sufficient to prove that it is an injection. Thus, given z ∈ Z k with π h z ∈ B k h , we must prove that z ∈ B k . By the Hodge decomposition, z = q + b with q ∈ H k and b ∈ B k . We have that π h z ∈ B k h by assumption and π h b ∈ B k h since b ∈ B k and π h is a cochain map. Thus π h q = π h z − π h b ∈ B k h ⊂ B k , and so π h q ⊥ q. In view of (25), this implies that q = 0, and so z ∈ B k as desired.
Remark. In applications, the space of harmonic forms, H k , is a finite-dimensional space of smooth functions, and π h is a projection operator associated to a triangulation with mesh size h. The estimate (25) will then be satisfied for h sufficiently small. However in the most important application, in which (V, d) is the de Rham complex and (V h , d) is a finite element discretization, π h induces an isomorphism on cohomology not only for h sufficiently small, but in fact for all h. See Section 5.6.
The second result relating H k and H k h is quantitative in nature, bounding the distance, or gap, between these two spaces. Recall that the gap between two subspaces E and F of a Hilbert space V is defined [64, Chapter IV, § 2.1] Hilbert subcomplex, and π h a bounded cochain projection. Then , and so is orthogonal to both the discrete harmonic form q and the harmonic form P H q. Therefore (29) is an immediate consequence of (27) and (28).
Next we deduce another important property of a Hilbert subcomplex with a bounded cochain projection. Since the complex (V, d) is closed and bounded, the Poincaré inequality (16) holds (with the W and V norms coinciding). Now we obtain the Poincaré inequality for the subcomplex with a constant that depends only on the Poincaré constant for the supercomplex and the norm of the cochain projection. In the applications, we will have a sequence of such subcomplexes related to a decreasing mesh size parameter and this theorem will imply that the discrete Poincaré inequality is uniform with respect to the mesh parameter, an essential step in proving stability for numerical methods.
and the result follows.
We have established several important properties possessed by a subcomplex of a bounded closed Hilbert complex with bounded cochain projection. We also remark that from (28) and the triangle inequality, we have (30) q V ≤ c P H q V , q ∈ H k h . We close this section by presenting a converse result. Namely we show that if the discrete Poincaré inequality and the bound (30) hold, then a bounded cochain projection exists.
for some constants c 1 and c 2 . Then there exists a bounded cochain projection π h from (V, d) to (V h , d), and the V operator norm π h can be bounded in terms of c 1 and c 2 .
Proof. As a first step of the proof we define an operator Q h : This operator is bounded in V k and invariant on V k h , since the three terms correspond exactly to the discrete Hodge decomposition in this case. Furthermore, 3.4. Stability and convergence of the mixed method. Next we consider a closed Hilbert complex (W, d) and the Galerkin discretization of its Hodge Laplacian using finite dimensional subspaces V k h of the domain spaces V k . Our main assumptions are those of Section 3.3: first, that dV k h ⊂ V k+1 h , so that we obtain a subcomplex Let f ∈ W k . In view of the mixed formulation (18), we take as an approximation scheme: (Recall that we use · , · and · without subscripts for the W inner product and norm.) The discretization (31) is a generalized Galerkin method as discussed in Section 2.2. In the case that there are no harmonic forms (and therefore no discrete harmonic forms), it is a Galerkin method, but in general not, since H k h is not in general a subspace of H k . We may write the solution of (31), which always exists and is unique in view of the results of Section 3.2.2, as h is the W k -orthogonal projection. As in Section 2.2, we will bound the error in terms of the stability of the discretization and the consistency error. We start by establishing a lower bound on the inf-sup constant, i.e., an upper bound on the stability constant.
Proof. This is just Theorem 3.2 applied to the Hilbert complex (V h , d), combined with the fact that the constant in the Poincaré inequality for V k h is c P π h by Theorem 3.6.
From this stability result, we obtain the following error estimate.
Proof. First observe that (σ, u, p) satisfies Theorem 3.8 then gives Using (27) and the boundedness of the projection π h we have (33) p

Next we show that
By Theorem 3.5, there exists r ∈ H k (and so r ⊥ B k ) with r ≤ 1 and since π h is a bounded projection. The theorem follows from (32)-(34) and the triangle inequality.
To implement the discrete problem, we need to be able to compute the discrete harmonic forms. The following lemma shows one way to do this; namely it shows that the discrete harmonic forms can be computed as the elements of the null space of a matrix. For finite element approximations of the de Rham sequence, which are the most canonical example of this theory and which will be discussed below, it is often possible to compute the discrete harmonic forms more directly. See, for example, [2].
On the other hand, if (σ h , u h ) is a solution, by taking τ = σ h , v = u h , and combining the two equations, we find that σ h 2 + du h 2 = 0, so that σ h = 0 and du h = 0. Then the first equation implies that dτ, u h = 0 for all τ ∈ V k−1 h , so indeed u h ∈ H k h . 3.5. Improved error estimates. Suppose we have a family of subcomplexes (V h , d) of the domain complex (V, d) of a closed Hilbert complex, parametrized by h with uniformly bounded cochain projections π h . Assuming also that the subspaces V k h are approximating in V k in the sense of (8), we can conclude from Theorem 3.9 that σ h → σ, u h → u, and p h → p as h → 0 (in the norms of V k−1 and V k ). In other words, the Galerkin method for the Hodge Laplacian is convergent.
The rate of convergence will depend on the approximation properties of the subspaces V k h , the particular component considered (σ h , u h , or p h ), the norm in which we measure the error (e.g., W or V ), as well as properties of the data f and the corresponding solution. For example, in Section 5 we will consider approximation of the de Rham complex using various subcomplexes for which the spaces V k h consist of piecewise polynomial differential forms with respect to a triangulation T h of the domain with mesh size h. The space W k is the space of L 2 differential k-forms in this case. One possibility we consider for the solution of the Hodge Laplacian for k-forms using the mixed formulation is to take subspaces P r+1 Λ k−1 (T h ) and P r Λ k (T h ). The space P r Λ k (T h ), which is defined in Section 5.2, consists piecewise of all k-forms of polynomial degree at most r. Assuming that the solution u to the Hodge Laplacian is sufficiently smooth, application of Theorem 3.9 will give, in this case, . Approximation theory tells us that this rate is the best possible for u − u h V , but we might hope for a faster rate for u − u h and for σ − σ h V and σ − σ h .
In order to obtain improved error estimates, we make two addition assumptions, first that the complex (W, d) satisfies the compactness property introduced at the end of Section 3.1, and second, that the cochain projection is bounded not only in V but in W : • The intersection V k ∩ V * k is a dense subset of W k with compact inclusion. • The cochain projections π k h are bounded in L(W k , W k ) uniformly with respect to h. The second property implies that π h extends to a bounded linear operator W k → V k h . Since the subspaces V k h are approximating in W k as well as in V k (by density), it follows that π h converges pointwise to the identity in W k . Finally, note that if we are given a W -bounded cochain projection mapping W k → V k h , the restrictions to V k define a V -bounded cochain projection.
Next, we note that on the Hilbert space V k ∩ V * k , the inner product given by u, v V ∩V * := d * u, d * v + du, dv + P H u, P H v is equivalent to the usual intersection inner product which is the sum of the inner products for V k and V * k . This can be seen by Hodge decomposing u as P B u + P H u + P B * u, and using the Poincaré bound P B * u ≤ c dP B * u = c du and the analogous bound is the compact inclusion operator. Hence K is a compact operator W k → V k ∩ V * k and, a fortiori, compact as an operator from W k to itself. As an operator on W k , K is also self-adjoint, since for all f, g ∈ W k . Furthermore, if we follow K by one of the bounded operators d : V k → W k+1 or d * : V * k → W k−1 , the compositions dK and d * K are also compact operators from W k to itself. Since we have assumed the compactness property, dim H k < ∞, and so P H k is also a compact operator on W k . Define (Note that µ already appeared in Theorem 3.9.) Recalling that composition on the right with a compact operator converts pointwise convergence to norm convergence, we see that η, δ, µ → 0 as h → 0. In the applications in Sections 5 and 6, the spaces Λ k h will consist of piecewise polynomials. We will then have η = O(h), δ = O(h min(2,r+1) ), and µ = O(h r+1 ) where r denotes the largest degree of complete polynomials in the space Λ k h . In Theorem 3.9, the error estimates were given in terms of the best approximation error afforded by the subspaces in the V norm. The improved error estimates will be in terms of the best approximation error in the W norm, for which we introduce the notation The following theorem gives improved the improved error estimates. Its proof incorporates a variety of techniques developed in the numerical analysis literature in the last decades, e.g., [48,41,8].  (31). Then for some constant C independent of h and (σ, u, p), we have We now develop the proof of Theorem 3.11 in a series of lemmas.
Proof. The second estimate is just (34). Using the Hodge decomposition of f and the fact that by (27). Applying (28) we get Combining these results completes the proof of the lemma. Proof. From (18) and (31), (41) and set v = d * Kdv h = P B * v h . Using Lemmas 3.12 and (27), we get The result follows by the triangle inequality.
Then using Lemmas 3.12 and 3.15, we get We next estimate the term v, u − u h . Since dv = dv h and Kv ∈ Z k⊥ , we get Combining these results, and using the previous lemmas, we obtain The final result of the lemma follows immediately.
It is now an easy matter to prove (39) and so complete the proof Theorem 3.11. We write , and so (39) follows from Lemmas 3.13, 3.14, and 3.16.
To get a feeling for these results, we return to the example mentioned earlier, where V k−1 h = P r+1 Λ k−1 (T h ) and V k h = P r Λ k (T h ) are used to approximate the k-form Hodge Laplacian with some r ≥ 1. If the domain is convex, we may apply elliptic regularity to see that Kf belongs to the Sobolev space H 2 Λ k for f ∈ L 2 Λ k , so dKf ∈ H 1 Λ k+1 and d * Kf ∈ H 1 Λ k−1 , and then standard approximation theory shows that η = O(h), δ = O(h 2 ), and µ = O(h 2 ). From Theorem 3.11, we then obtain that assuming the solution u is sufficiently smooth. That is, all components converge with the optimal order possible given the degree of polynomial approximation.
Finally we note a corollary of Theorem 3.11, which will be useful in the analysis of the eigenvalue problem.
Corollary 3.17. Under the assumptions of Theorem 3.11, there exists a constant C such that Therefore all three operator norms converge to zero with h.
Proof. Let f ∈ W k and set u = Kf , σ = d * Kf , p = P H f , and We end this section with an estimate of the difference between the true Hodge decomposition and the discrete Hodge decomposition of an element of V k h . While this result is not needed in our approach, such an estimate was central to the estimation of the eigenvalue error for the Hodge Laplacian using Whitney forms made in [40].
h u denote the Hodge and discrete Hodge decompositions of u ∈ V k h . Then The first estimate follows immediately from Lemma 3.12. To obtain the second estimate, we write . The second estimate now follows directly from (28) and the first estimate and the final estimate follows from the first two by the triangle inequality. 3.6. The eigenvalue problem. The purpose of this section is to study the eigenvalue problem associated to the abstract Hodge Laplacian (18). As in the previous section, we will assume (W, d) is a Hilbert complex satisfying the compactness property and that the cochain projections π k h are bounded in L(W k , W k ), uniformly in h. A pair (λ, u) ∈ R × V k , where u = 0, is referred to as an eigenvalue/eigenvector of the problem (18) if there exists (σ, p) ∈ V k−1 × H k such that In operator terms, u = λKu, σ = d * u, p = 0. Note that it follows from this system that λ u 2 = du 2 + d * u 2 > 0, so Ku = λ −1 u. Since the operator K ∈ L(W k , W k ) is compact and self-adjoint, we can conclude that the problem (42) has at most a countable set of eigenvalues, each of finite multiplicity. We denote these by where each eigenvalue is repeated according to its multiplicity. Furthermore, when W k is infinite dimensional, we have lim j→∞ λ j = ∞. We denote by {v i } a corresponding orthonormal basis of eigenvectors for W k . The corresponding discrete eigenvalue problem takes the form As above, we can conclude that p h = 0, λ h > 0, and that λ −1 h is an eigenvalue for the operator (43), repeated according to multiplicity, and by {v i,h } a corresponding W k -orthonormal eigenbasis for V k h . Next, we will study how the discrete eigenvalue problem (43) converges to the eigenvalue problem (42), i.e., how the eigenvalues and eigenvectors of the operator K h converge to those of K. We let E i and E i,h denote the 1-dimensional spaces spanned by v i and v i,h . For every positive integer j, let m(j) denote the number of eigenvalues less than or equal to the jth distinct eigenvalue of the Hodge-Laplace problem (42). Thus m(j) i=1 E i is the space spanned by the eigenspaces associated to the first j distinct eigenvalues and does not depend on the choice of the eigenbasis.
The discrete eigenvalue problem (43) is said to converge to the exact eigenvalue problem (42) if, for any > 0 and integer j > 0, there exists a mesh parameter h 0 > 0 such that, for all h ≤ h 0 , we have (44) max where the gap between two subspaces E and F of a Hilbert space is defined by (26). The motivation for this rather strict concept of convergence is that it not only implies that each eigenvalue is approximated by the appropriate number of discrete eigenvalues, counting multiplicities, and that the eigenspace is well-approximated by the corresponding discrete eigenspaces, but it also rules out spurious discrete eigenvalues and eigenvectors. In particular, this rules out the behavior exemplified in Figure 2.6 in Section 2.
A key result in the perturbation theory of linear operators is that, for eigenvalue problems of the form we consider, corresponding to the bounded compact self-adjoint Hilbert space operators K and K h , convergence of the eigenvalue approximation holds if the operators K h P h converge to K in L(W k , W k ). This result, widely used in the theory of mixed finite element eigenvalue approximation, essentially follows from the contour integral representation of the spectral projection, and can be extracted from [ [20], this operator norm convergence is sufficient as well as necessary for obtaining convergence of the eigenvalue approximations in the sense above. As a consequence of Corollary 3.17, we therefore obtain the following theorem. It is also possible to use the theory developed above to obtain rates of convergence for the approximation of eigenvalues, based on the following result (Theorem 7.3 of [14]). Here, we consider only the case of a simple eigenvalue, but with a small modification the results extend to eigenvalues of positive multiplicity. Lemma 3.20. If λ is a simple eigenvalue and u a normalized eigenvector, then Proof. To estimate the first term in Lemma 3.20, we write where we have used the orthogonality condition in the second last line above. The theorem then follows from this estimate and the fact that Order of convergence estimates now follow directly from Theorem 3.11, with f replaced by u, u replaced by w, and σ replaced by φ. In the example following Theorem 3.11, i.e., V k−1 as long as the domain is convex and the solution w = Ku sufficiently smooth. Inserting these results into Theorem 3.21, we find that the eigenvalue error |λ−λ h | = O(h 2r ), double the rate achieved for the source problem. As another example, one can check that for the Whitney forms, V k−1 , improving on the O(h| log h|) estimate of [40].
Remark. It has long been observed that for mixed finite element approximation of eigenvalue problems, stability and approximability alone, while sufficient for convergence of approximations of the source problem, is not sufficient for convergence of the eigenvalue problem. An extensive literature developed in order to obtain eigenvalue convergence, and a wide variety of additional properties of the finite element spaces have been defined and hypothesized. In particular, the first convergence results for the the important case of electromagnetic eigenvalue problems, were obtained by Kikuchi based on the discrete compactness property [65,19], and the more recent approach by Boffi and collaborators emphasize the Fortid property [18]. In our context, the discrete compactness property says that whenever v h ∈ Z k⊥ h is uniformly bounded in V k , then there exists a sequence h i converging to zero such that v hi converges in W k , and the Fortid property says that lim h→0 I − π k h L(V k ∩V * k ,W k ) = 0. These additional properties play no role in the theory as presented here. We prove eigenvalue convergence under the same sort of assumptions we use to establish stability and convergence for the source problem, namely subcomplexes, bounded cochain projections, and approximability, but we require boundedness in W for eigenvalue convergence, while for stability we only require boundedness in V . We believe that subcomplexes with bounded cochain projections provide an appropriate framework for the analysis of the eigenvalue problem, since, as far as we know, these properties hold in all examples where eigenvalue convergence has been obtained by other methods. Moreover, it is easy to show that the discrete compactness property and Fortid property hold whenever there exist W -bounded cochain projections.
3.6.1. Related eigenvalue problems. Recall that the source problem for the Hodge Laplacian, given by (18), can be decomposed into the problems (19) and (20), referred to as the B * and the B problem, respectively. More precisely, these problems arise if the right hand side f of (18) is in B * k or B k , respectively. In a similar way the eigenvalue problem (42) can be decomposed into a B * and the B problem. To see this, note if (λ, u) ∈ R × V k is an eigenvalue/eigenvector for (42), then u ∈ H k⊥ , and therefore u = u B + u ⊥ , where u B ∈ B k and u ⊥ ∈ B * k ∩ V k = Z k⊥ . It is straightforward to check that u ⊥ satisfies the B * k eigenvalue problem given by (45) du Here u ⊥ can be equal to zero even if u = 0 is an eigenvector. On the other hand, if the pair (λ, u ⊥ ) ∈ R × Z k⊥ is an eigenvalue/eigenvector for (45), then (λ, u ⊥ ) is also an eigenvalue/eigenvector for (42), where σ and p are both zero. Furthermore, for the extended eigenvalue problem where u ⊥ is sought in V k , the nonzero eigenvalues correspond precisely to the eigenvalues of (45), while the eigenvalue λ = 0 has the eigenspace Z k . The pair (λ, u B ) satisfies the corresponding B k eigenvalue problem given by Furthermore, any solution (λ, u B ) of (47), where u B = 0, corresponds to an eigenvalue/eigenvector of the full Hodge Laplacian (42). In fact, any eigenvalue of (42) is an eigenvalue of either the B * or B problem, or both, and all eigenvalues of the B * and B problems correspond to an eigenvalue of (42). In short, the eigenvalue problem (42) can be decomposed into the two problems (45) and (47).
Also, note that if (λ, u B ) is an eigenvalue/eigenvector for the B k problem (47), and σ = d * u B ∈ Z (k−1)⊥ , then σ = 0, and by taking v = dτ in (47), we obtain (48) dσ, dτ = λ σ, τ , τ ∈ Z (k−1)⊥ , which is a B * k−1 problem. Hence, we conclude that any eigenvalue of the B k problem is an eigenvalue of the B * k−1 problem. The converse is also true. To see this, let (λ, σ) ∈ R × Z (k−1)⊥ , σ = 0, be a solution of (48), and define u B ∈ B k such that d * u B = σ. Since d * : B k → Z (k−1)⊥ is a bijection, this determines u B uniquely, and (λ, u B ) is an eigenvalue/eigenvector for the B k problem (47). Therefore, the two problems (47) and (48) are equivalent, with the same eigenvalues, and with eigenvectors related by d * u B = σ. In particular, the eigenvalues of the full Hodge Laplace problem (42) correspond precisely to the eigenvalues of the B * k−1 and B * k problems.
In the special case of the de Rham complex, the B * 1 problem is closely related to Maxwell's equations. Discretizations of this problem have therefore been intensively studied in the literature (e.g., see [19] and the references therein). However, it is usually not straightforward to compute the eigenvalues and eigenvectors of the discrete B * k problem from the formulation (45) since we do not, a priori, have a basis for the corresponding space Z k⊥ h ⊂ V k h available. Usually, we will have only have a basis for the space V k h at our disposal, and a direct computation of a basis for Z k⊥ h from this basis is costly. A better alternative is therefore to solve the discrete version of the extended eigenvalue problem (46), and observe that the positive eigenvalues, and the corresponding eigenvectors, are precisely the solutions of the corresponding problem (45). Alternatively, we can solve the full Hodge Laplace problem (43), and then throw away all eigenvalues corresponding to σ h = 0.

Exterior calculus and the de Rham complex
We next turn to the most important example of the preceding theory, in which the Hilbert complex is the de Rham complex associated to a bounded domain Ω in R n . We begin by a quick review of basic notions from exterior calculus. Details can be found in many references, e.g., [8,12,22,49,62,67,85].

4.1.
Basic notions from exterior calculus. For a vector space V and a nonnegative integer k, we denote by Alt k V the space of real-valued k-linear forms on V . If dim V = n, then dim Alt k V = n k . The wedge product ω ∧ η ∈ Alt j+k V of ω ∈ Alt j V and η ∈ Alt k V is given by where the sum is over all permutations σ of {1, . . . , j + k}, for which σ(1) < σ(2) < · · · < σ(j) and σ(j + 1) < σ(j + 2) < · · · < σ(j + k). An inner product on V induces an inner product on Alt k V : ω, η = σ ω(e σ(1) , . . . , e σ(k) )η(e σ(1) , . . . , e σ(k) ), ω, η ∈ Alt k V, where the sum is over increasing sequences σ : {1, . . . , k} → {1, . . . , n} and e 1 , . . ., e n is any orthonormal basis (the right-hand side being independent of the choice of orthonormal basis). If we orient V by assigning positive orientation to some particular ordered basis (thereby assigning a positive or negative orientation to all ordered bases, according to the determinant of the change of basis transformation), then we may define a unique volume form vol in Alt n V , n = dim V , characterized by vol(e 1 , . . . , e n ) = 1 for any positively oriented ordered orthonormal basis e 1 , . . . , e n . The Hodge star operator is an isometry of Alt k V onto Alt n−k V given by Given a smooth manifold Ω, possibly with boundary, a differential k-form ω is a section of the k-alternating bundle, i.e., a map which assigns to each x ∈ Ω an element ω x ∈ Alt k T x Ω where T x Ω denotes the tangent space to Ω at x. We write C m Λ k (Ω) for the space of m times continuously differentiable k-forms, i.e., forms for which x → ω x (v 1 (x), . . . , v n (x)) belongs to C m (Ω) for any smooth vector fields v i . Similarly, we may define C ∞ Λ k , C ∞ c Λ k (smooth forms with compact support contained in the interior of Ω), etc. If Ω is a Riemannian manifold, and so has a measure defined on it, we may similarly define Lebesgue spaces L p Λ k and Sobolev spaces W m p Λ k and H m Λ k = W m 2 Λ k . The spaces H m Λ k are Hilbert spaces. In particular, the inner product in L 2 Λ k = H 0 Λ k is given by We also write Λ k or Λ k (Ω) for the space of all smooth differential k-forms, or at least sufficiently smooth, as demanded by the context.
For any smooth manifold Ω and ω ∈ Λ k (Ω), the exterior derivative dω is a (k + 1)-form, which itself has vanishing exterior derivative: d(dω) = 0. On an oriented n-dimensional piecewise smooth manifold, a differential n-form (with, e.g., compact support) can be integrated, without recourse to a measure or metric.
A smooth map φ : Ω → Ω between manifolds, induces a pullback of differential forms from Ω to Ω . Namely, if ω ∈ Λ k (Ω), the pullback φ * ω ∈ Λ k (Ω ) is defined by The pullback respects exterior products and exterior derivatives: If φ is an orientation-preserving diffeomorphism of oriented n-dimensional manifolds, and ω is an n-form on Ω, then If Ω is a submanifold of Ω, then the pullback of the inclusion map is the trace map tr Ω,Ω , written simply tr Ω or tr when the manifolds can be inferred from context. We recall the trace theorem which states that if Ω is a submanifold of codimension 1, then the trace map extends to a bounded operator from H 1 Λ k (Ω) to L 2 Λ k (Ω ), or, more precisely, to a bounded surjection of H 1 Λ k (Ω) onto H 1/2 Λ k (Ω ). A particularly important situation is when Ω = ∂Ω, in which case Stokes's theorem relates the integrals of the exterior derivative and trace of an (n − 1)-form ω on an oriented n-dimensional manifold-with-boundary Ω: (A common abuse of notation is to write ∂Ω ω for the ∂Ω tr ω.) Applying Stokes's theorem to the differential form ω ∧ η with ω ∈ Λ k−1 , η ∈ Λ n−k , and using the Leibniz rule d(ω ∧ η) = dω ∧ η + (−1) k−1 ω ∧ dη, we obtain the integration-by-parts formula for differential forms: On an oriented n-dimensional Riemannian manifold, there is a volume form vol ∈ Λ n (Ω) which at each point x of the manifold is equal to the volume form on T x Ω. Consequently, the Hodge star operation takes ω ∈ Λ k (Ω) to ω ∈ Λ n−k (Ω) satisfying for all µ ∈ Λ n−k (Ω). Introducing the coderivative operator δ : Λ k → Λ k−1 defined by (50) δω = (−1) k d ω, and setting η = µ, the integration-by-parts formula becomes (51) dω, µ = ω, δµ + ∂Ω tr ω ∧ tr µ, ω ∈ Λ k−1 , µ ∈ Λ k .
Since the Hodge star operator is smooth and an isometry at every point, every property of k-forms yields a corresponding property for n − k forms. For example, the Hodge star operator maps the spaces C m Λ k , W m p Λ k , etc., isometrically onto C m Λ n−k , W m p Λ n−k , etc. By definition, δ : Λ n−k → Λ n−k−1 corresponds to d : Λ k → Λ k+1 via the Hodge star isomorphism. Thus for every property of d there is a corresponding property of δ.
In case Ω is a domain in R n , we may (but usually will not) use the standard coordinates of R n to write a k-form as where the a σ ∈ L 2 (Ω) and dx j : R n → R is the linear form which associates to a vector its jth coordinate. In this case, the exterior derivative is given by the simple formula extended by linearity to a sum of such terms. If Ω is a domain in R n , then Ω a dx 1 ∧ · · · ∧ dx σn ), has the value suggested by the notation. In the case of a domain in R n , the volume form is simply dx 1 ∧ · · · ∧ dx n . We remark that the exterior derivative and integral of differential forms can be computed on arbitrary manifolds from the formulas on subdomains on R n and pullbacks through charts.

4.2.
The de Rham complex as a Hilbert complex. Henceforth we restrict attention to the case that Ω is a bounded domain in R n with a piecewise smooth, Lipschitz boundary. In this section we show that the de Rham complex is a Hilbert complex which satisfies the compactness property, and so the abstract theory of Section 3 applies. We then interpret the results in the case of the de Rham complex. We begin by indicating how the exterior derivative d can be viewed as a closed densely-defined operator from W k = L 2 Λ k (Ω) to W k+1 = L 2 Λ k+1 (Ω). Let ω ∈ L 2 Λ k (Ω). In view of (51), we say that η ∈ L 2 Λ k+1 (Ω) is the weak exterior derivative of ω if ω, δµ = η, µ , µ ∈ C ∞ c Λ k+1 . The weak exterior derivative of ω, if one exists, is unique and we denote it by dω.
In analogy with the definition of Sobolev spaces (cf., e.g., [47, Section 5.2.2]), we define HΛ k to be the space of forms in L 2 Λ k with a weak derivative in L 2 Λ k+1 . With the inner product this is easily seen to be a Hilbert space and clearly d is a bounded map from HΛ k to L 2 Λ k+1 . A standard smoothing argument, as in [47,Section 5.3], implies that C ∞ Λ k (Ω) is dense in HΛ k . We take HΛ k (Ω) as the domain V k of the exterior derivative, which is thus densely defined in L 2 Λ k (Ω). Since HΛ k is complete, d is a closed operator.
Thus the spaces L 2 Λ k (Ω) and the exterior derivative operators d form a Hilbert complex with the associated domain complex (52) 0 → HΛ 0 (Ω) This is the L 2 de Rham complex. We shall see below that it satisfies the compactness property.
Note that H * Λ k consists of those differential form in L 2 Λ k for which a weak coderivative exists in L 2 Λ k−1 , where the weak exterior coderivative is defined in exact analogy to the weak exterior derivative. Its inner product is The space H * Λ k is isometric to HΛ n−k via the Hodge star, but is quite different from HΛ k .
Theorem 4.1. Let d be the exterior derivative viewed as an unbounded operator L 2 Λ k−1 → L 2 Λ k with domain HΛ k . Then the adjoint d * , as an unbounded operator L 2 Λ k → L 2 Λ k−1 , has as its domainH * Λ k , and coincides with the operator δ defined in (50).
As a corollary, we obtain a concrete characterization of the harmonic forms: In other words, a k-form is harmonic if it satisfies the differential equations dω = 0 and δω = 0 together with the boundary conditions tr ω = 0. If the boundary of Ω is smooth, then HΛ k ∩H * Λ k is contained in H 1 Λ k [52], and hence, by the Rellich theorem, we obtain the compactness property of Section 3.1.3. For a general Lipschitz boundary, e.g., for a polygonal domain, the inclusion of HΛ k ∩H * Λ k in H 1 Λ k need not hold, but the compactness property remains valid [75]. Thus all the results of Section 3 apply to the de Rham complex. In particular, we have the Hodge decomposition of L 2 Λ k and of HΛ k , the Poincaré inequality, well-posedness of the mixed formulation of the Hodge Laplacian, and all the approximation results established in Sections 3.3-3.6. We now interpret these results a bit more concretely in the present setting.
First of all, the cohomology groups associated to the complex (52) are the de Rham cohomology groups, whose dimensions are the Betti numbers of the domain. Turning next to the Hodge Laplacian problem given in the abstract case by (18), we get that (σ, u, p) ∈ HΛ k−1 × HΛ k × H k is a solution if and only if tr u = 0, tr du = 0 on ∂Ω, The first differential equation and the first boundary condition are implied by the first equation in (18), and the second differential equation and boundary condition by the second, while the third equation in (18) is simply the side condition u ⊥ H k . Note that both boundary conditions are natural in this variational formulation: they are implied but not imposed in the spaces where the solution is sought. Essential boundary conditions could be imposed instead. We discuss this in Section 6.2.
To make things more concrete, we now restrict to a domain Ω ⊂ R 3 , and consider the Hodge Laplacian for k-forms, k = 0, 1, 2, and 3. We also discuss the B * and B problems given by (19) and (20) for each k. We shall encounter many of the most important partial differential equations of mathematical physics: the Laplacian, the vector Laplacian, div-curl problems, and curl-curl problems. These PDEs arise in manifold applications in electromagnetism, solid mechanics, fluid mechanics, and many other fields.
We begin by noting that, on any oriented Riemannian manifold of dimension n, we have a natural way to view 0-forms and n-forms as real-valued functions, and 1-forms and (n − 1)-forms as vector fields. In fact, 0-forms are real-valued functions and 1-forms are covector fields, which can be identified with vector fields via the inner product. The Hodge star operation then carries these identifications to nforms and (n − 1)-forms. In the case of a 3-dimensional domain in R 3 , via these identifications all k-forms can be viewed as either scalar or vector fields (sometimes called proxy fields). With these identifications, the Hodge star operation becomes trivial in the sense if a certain vector field is the proxy for, e.g., a 1-form ω, the exact same vector field is the proxy for the 2-form ω. Via proxy fields, the exterior derivatives coincide with standard differential operators of calculus: d 0 = grad, d 1 = curl, d 2 = div, and the de Rham complex (52) is realized as The exterior coderivatives δ become, of course, − div, curl, and − grad, when acting on 1-forms, 2-forms, and 3-forms, respectively. The trace operation on 0-forms is just the restriction to the boundary, and the trace operator on 3-forms vanishes (since there are no nonzero 3-forms on ∂Ω). The trace operator from 1-forms on Ω to 1-forms on the boundary takes a vector field u on Ω to a tangential vector field on the boundary, namely at each boundary point x, (tr u) x is the tangential projection of u x . For a 2-form u, the trace corresponds to the scalar u · n (with n the unit normal) at each boundary point. Similarly, the first boundary condition in (60) vanishes, while the second is the Neumann condition grad u · n = 0 on ∂Ω.
The side condition (61) specifies a unique solution by requiring its average value to be zero. Nothing additional is obtained by considering the split into the B * and B subproblems, since the latter is trivial. Furthermore, the eigenvalue problem (42) is precisely the corresponding eigenvalue problem for the scalar Laplacian, with the 0 eigenspace H 0 filtered out.

4.2.2.
The Hodge Laplacian for k = 1. In this case the differential equations and boundary conditions are which is a formulation of a boundary value problem for the vector Laplacian curl curl − grad div. (Here we have used the fact that the vanishing of the tangential component of a vector is equivalent to the vanishing of its cross product with the normal.) The solution is determined uniquely by the additional condition that u be orthogonal to H 1 , which in this case consists of those vector fields satisfying curl p = 0, div p = 0 in Ω, p · n = 0, on ∂Ω.
The dimension of H 1 is equal to the first Betti number, i.e., the number of handles, of the domain (so H 1 = 0 if the domain is simply-connected). The B * 1 problem (19) is defined for L 2 vector fields f which are orthogonal to both gradients and the vector fields in H 1 . In that case, the solution to (62) has σ = 0 and p = 0, while u satisfies curl curl u = f, div u = 0 in Ω, u · n = 0, curl u × n = 0 on ∂Ω.
The orthogonality condition u ⊥ H 1 again determines the solution uniquely.
Next we turn to the B 1 problem. For source functions of the form f = grad F for some F ∈ H 1 , which we normalize so that Ω F = 0, (62) reduces to the problem of finding σ ∈ H 1 and u ∈ B 1 = grad H 1 such that: The differential equations may be simplified to − grad div u = f , and the condition that u ∈ B 1 can be replaced by the differential equation curl u = 0, together with orthogonality to H 1 . Now grad(σ − F ) = 0 and Ω σ = − ∂Ω u · n = 0 = Ω F , so σ = F , and we may rewrite the system as − div u = F, curl u = 0 in Ω, u · n = 0, on ∂Ω, which, again, has a unique solution subject to orthogonality to H 1 .
The eigenvalue problem (42) is the corresponding eigenvalue problem for the vector Laplacian with the boundary conditions of (62), and with the eigenspace H 1 of the eigenvalue λ = 0 filtered out. As mentioned above, in this case the B * 1 eigenvalue problem, given by (45), is important for models based on Maxwell's equations. This problem takes the form (63) curl curl u = λu, div u = 0 in Ω, u · n = 0, curl u × n = 0 on ∂Ω, u ⊥ H 1 .

4.2.3.
The Hodge Laplacian for k = 2. The differential equations and boundary conditions are This is again a formulation of a boundary value problem for the vector Laplacian curl curl − grad div, but with different boundary conditions than for (62), and this time stated in terms of two vector variables, rather than one vector and one scalar. This time uniqueness is obtained by imposing orthogonality to H 2 , the space of vector fields satisfying curl p = 0, div p = 0 in Ω, p × n = 0, on ∂Ω, which has dimension equal to the second Betti number, i.e., the number of voids in the domain. The B * 2 problem arises for source functions of the form f = grad F for some F ∈H 1 (H 1 is the function space corresponding toH * Λ 3 ). We find σ = 0, and u solves − div u = F, curl u = 0 in Ω, u × n = 0, on ∂Ω, i.e., the same differential equation as for B 1 , but with different boundary conditions, an extra assumption on F , and, of course, now uniqueness is determined by orthogonality to H 2 . If div f = 0 and f ⊥ H 2 , we get the B 2 problem for which the differential equations are σ = curl u, curl σ = f and the condition div u = 0 arising from the membership of u in B 2 . Thus u solves curl curl u = f, div u = 0 in Ω, u × n = 0, on ∂Ω, the same differential equation as for B * 1 , but with different boundary conditions. The eigenvalue problem (42) is the corresponding eigenvalue problem for the vector Laplacian with the boundary conditions of (64), and with the eigenspace H 2 of the eigenvalue λ = 0 filtered out, while the corresponding B 2 problem, of the form (47), takes the form curl curl u = λu, div u = 0 in Ω, u × n = 0, on ∂Ω, u ⊥ H 2 .
Note that this is the same problem as the B * 1 eigenvalue problem(63), but with different boundary conditions. However, if we define σ = curl u, then it is straightforward to check that the pair (λ, σ) will indeed solve the B * 1 problem (63). This is an instance of the general equivalence between the B k problem, given by (47), and the corresponding B * k−1 problem, which was pointed out in Section 3.6.1.

The Hodge
Laplacian for k = 3. In this case the Hodge Laplacian problem, which coincides with the B 3 problem, is which is the Dirichlet problem for Poisson's equation. There are no non-zero harmonic forms, and the problem has a unique solution. Furthermore, the eigenvalue problem (42) is the corresponding eigenvalue problem for the scalar Laplacian with Dirichlet boundary conditions.

Finite element approximation of the de Rham complex
Our goal in this section is to discretize the de Rham complex so that we may apply the abstract results on approximation of Hilbert complexes from Section 3. Hence we need to construct finite dimensional subspaces Λ k h of HΛ k (Ω). As we saw in Section 3.4, the key properties these spaces must possess is, first, that dΛ k h ⊂ Λ k+1 h so they form a subcomplex (Λ h , d) of the de Rham complex, second, that there exist uniformly bounded cochain projections π h from (L 2 Λ, d) to (Λ h , d), and third, good approximation properties. We may then use these spaces in a Galerkin method based on the mixed formulation, as described in Section 3.4, and the error estimates given in Theorems 3.9 and 3.11 bound the error in the Galerkin solution in terms of the approximation error afforded by the subspaces. In this section we will construct the spaces Λ k h as spaces of finite element differential forms and show that they satisfy all these requirements, and can be efficiently implemented.
As discussed in Section 2, a finite element space is a space of piecewise polynomials which is specified by the finite element assembly process, i.e., by giving a triangulation of the domain, a finite dimensional space of polynomial functions (or, in our case, differential forms) on each element of the triangulation, called the shape functions, and a set of degrees of freedom for the shape functions associated with the faces of various dimensions of the elements, which will be used to determine the degree of interelement continuity. Spaces constructed in this way can be implemented very efficiently, since they admit basis functions with small support and so lead to sparse algebraic systems.
In Section 5.1, we will discuss the spaces of polynomial differential forms which we will use as shape functions. In Section 5.2, we specify degrees of freedom for these spaces on simplices, and study the resulting finite element spaces of differential forms. In particular, we show that they can be combined in a variety of ways to form subcomplexes of the de Rham complex. In Section 5.3, we briefly describe recent work related to the implementation of such finite elements with explicit local bases. In Sections 5.4 and 5.5, we construct L 2 -bounded cochain projections into these spaces and obtain error estimates for them. Finally Section 5.6 is simply a matter of collecting our results in order to obtain error estimates for the resulting approximations of the Hodge Laplacian. Many of the results of this section have appeared previously, primarily in [8], and therefore many proofs are omitted. 5.1. Polynomial differential forms and the Koszul complex. In this section, we consider spaces of polynomial differential forms, which lead to a variety of subcomplexes of the de Rham complex. These will be used in the next section to construct finite element spaces of differential forms. The simplest spaces of polynomial differential k-forms are the spaces P r Λ k (R n ) consisting of all differential k-forms on R n whose coefficients are polynomials of degree at most r. In addition to these spaces, we will use another family of polynomial form spaces, denoted P − r Λ k (R n ), which will be constructed and analyzed using the Koszul differential and the associated Koszul complex. Spaces taken from these two families can be combined into polynomial subcomplexes of the de Rham complex in numerous ways (there are essentially 2 n−1 such subcomplexes associated to each polynomial degree). These will lead to finite element de Rham subcomplexes, presented in Section 5.2. Some of these have appeared in the literature previously, with the systematic derivation of all of them first appearing in [6]. 5.1.1. Polynomial differential forms. Let P r (R n ) and H r (R n ) denote the spaces of polynomials in n variables of degree at most r and of homogeneous polynomial functions of degree r, respectively. We interpret these spaces to be the zero space if r < 0. We can then define spaces of polynomial differential forms, P r Λ k (R n ), H r Λ k (R n ), etc., as those differential forms which, when applied to a constant vector field, have the indicated polynomial dependence. For brevity, we will at times suppress R n from the notation and write simply P r , H r , P r Λ k , etc.
The dimensions of these spaces are easily calculated: and dim H r Λ k (R n ) = dim P r Λ k (R n−1 ). For each polynomial degree r ≥ 0, we get a homogeneous polynomial subcomplex of the de Rham complex: We shall verify below the exactness of this sequence. More precisely, the cohomology vanishes if r > 0 and also for r = 0 except in the lowest degree, where the cohomology space is R (reflecting the fact that the constants are killed by the gradient).
Taking the direct sum of the homogeneous polynomial de Rham complexes over all polynomial degrees gives the polynomial de Rham complex: for which the cohomology space is R in the lowest degree, and vanishes otherwise.

5.1.2.
The Koszul complex. Let x ∈ R n . Since there is a natural identification of R n with the tangent space T 0 R n at the origin, there is a vector in T 0 R n corresponding to x. (The origin is chosen for convenience here, but we could use any other point instead.) Then the translation map y → y + x induces an isomorphism from T 0 R n to T x R n , and so there is an element X(x) ∈ T x R n corresponding to x. (Essentially X(x) is the vector based at x which points opposite to the origin and whose length is |x|.) Contraction with the vector field X defines a map κ from Λ k (R n ) to Λ k−1 (R n ) called the Koszul differential : It is easy to see that κ is a graded differential, i.e., κ • κ = 0 and κ(ω ∧ η) = (κω) ∧ η + (−1) k ω ∧ (κη), ω ∈ Λ k , η ∈ Λ l . In terms of coordinates, if ω x = a(x) dx σ1 ∧ · · · ∧ dx σ k , then where the notation dx σi means that the term is omitted. Note that κ maps H r Λ k to H r+1 Λ k−1 , i.e., κ increases polynomial degree and decreases form degree, the exact opposite of the exterior derivative d.
The Koszul differential gives rise to the homogeneous Koszul complex [68, Chapter 3.4.6], (67) 0 We show below that this complex is exact for r > 0. Adding over polynomial degrees, we obtain the Koszul complex (for any r ≥ 0), for which all the cohomology spaces vanish, except the rightmost, which is equal to R.
To prove the exactness of the homogeneous polynomial de Rham and Koszul complexes, we note a key connection between the exterior derivative and the Koszul differential. In the language of homological algebra, this says that the Koszul differential is a contracting homotopy for the homogeneous polynomial de Rham complex.
Proof. This result can be established by a direct computation or by using the homotopy formula of differential geometry (Cartan's magic formula). See [8] for the details.
As a simple consequence of (68), we prove the injectivity of d on the range of κ and vice versa.
Another easy application of (68) is to establish the claimed cohomology of the Koszul complex and polynomial de Rham complex. Suppose that ω ∈ H r Λ k for some r, k ≥ 0 with r + k > 0, and that κω = 0. From (68), we see that ω = κη with η = dω/(r + k) ∈ H r−1 Λ k+1 . This establishes the exactness of the homogeneous Koszul complex (67) (except when r = 0 and the sequence reduces to 0 → R → 0). A similar argument establishes the exactness of (65).
Remark. We defined the Koszul differential as contraction with the vector field X, where X(x) is the translation to x of the vector pointing from the origin in R n to x. The choice of the origin as a base point is arbitrary-any point in R n could be used. That is, if y ∈ R n , we can define a vector field X y by assigning to each point x the translation to x of the vector pointing from y to x, and then define a Koszul differential κ y by contraction with X y . It is easy to check that for ω ∈ P r−1 Λ k+1 and any two points y, y ∈ R n , the difference κ y ω − κ y ω ∈ P r−1 Λ k . Hence the space P − r Λ k = P r−1 Λ k + κ y P r−1 Λ k+1 does not depend on the particular choice of the point y. This observation is important, because it allows us to define P − r Λ k (V ) for any affine subspace V of R n . We simply set where y is any point of V . Note that if ω ∈ P − r Λ k (R n ), then the trace of ω on V belongs to P − r Λ k (V ).
Remark. The spaces P − r Λ k (R n ) are affine-invariant, i.e., if φ : R n → R n is an affine map, then the pullback φ * maps this space into itself. Of course, the full polynomial space P r Λ k (R n ) is affine-invariant as well. In [8, Section 3.4], all the finite dimensional affine-invariant spaces of polynomial differential forms are determined. These are precisely the spaces in the P and P − families together with one further family of spaces which is of less interest to us.

5.1.4.
Exact sequences of polynomial differential forms. We have seen the polynomial de Rham complex (66) is a subcomplex of the de Rham complex on R n for which cohomology vanishes except for the constants at the lowest order. In other words, for any r ≥ 0 the sequence is exact, i.e., is a resolution of R.
As we shall soon verify, the complex is another resolution of R, for any r > 0. Note that in this complex, involving the P − r Λ k spaces, the polynomial degree r is held fixed, while in (71), the polynomial degree decreases as the form order increases. Recall that the 0th order spaces in these complexes, P r Λ 0 and P − r Λ 0 , coincide. In fact, the complex (71) is a subcomplex of (72), and these two are the extreme cases of a set of 2 n−1 different resolutions of R, each a subcomplex of the next, and all of which have the space P r Λ 0 in the 0th order. To exhibit them, we begin with the inclusion R → P r Λ 0 . We may continue the complex with either the map d : P r Λ 0 → P r−1 Λ 1 or d : P r Λ 0 → P − r Λ 1 , the former being a subcomplex of the latter. With either choice, the cohomology vanishes at the first position. Next, if we made the first choice, we can continue the complex with either d : P r−1 Λ 1 → P r−2 Λ 2 or d : P r−1 Λ 1 → P − r−1 Λ 2 . Or, if we made the second choice, we can continue with either d : In the first case, we may use the exactness of (71) to see that the second cohomology space vanishes. In the second case, this follows from Lemma 5.4. Continuing in this way at each order, k = 1, . . . , n − 1, we have two choices for the space of k-forms (but only one choice for k = n, since P r−1 Λ n coincides with P − r Λ n ), and so we obtain 2 n−1 complexes. These form a totally ordered set with respect to subcomplex inclusion. For r ≥ n these are all distinct (but for small r some coincide because the later spaces vanish).
In the case n = 3, the four complexes so obtained are:

5.2.
Degrees of freedom and finite element differential forms. Having introduced the spaces of polynomial differential forms P r Λ k (R n ) and P − r Λ k (R n ), we now wish to create finite element spaces of differential forms. We begin with notation for spaces of polynomial differential forms on simplices. If f is a simplex (of any dimension) in R n , we define The spaces P − r Λ k (f ) andP − r Λ k (f ) are defined similarly. Now let Ω be a bounded polyhedral domain which is triangulated, i.e., partitioned into a finite set T of n-simplices determining a simplicial decomposition of Ω. This means that the union of the elements of T is the closure of Ω, and the intersection of any two is either empty or a common subsimplex of each. By way of notation, for any simplex T we denote by ∆ d (T ) the set of subsimplices of T of dimension d, and by ∆(T ) the set of all subsimplices of T . We also set ∆ d (T ) := T ∈T ∆ d (T ) and ∆(T ) := T ∈T ∆(T ).
Corresponding to the P and P − families of spaces of polynomial differential forms, we will define two families of spaces of finite element differential forms with respect to the triangulation T , denoted P r Λ k (T ) and P − r Λ k (T ). We shall show that these are subspaces of HΛ k (Ω) and can be collected, in various ways, into subcomplexes of the de Rham complex.
The spaces P r Λ k (T ) and P − r Λ k (T ) will be obtained by the finite element assembly process. For each T ∈ T , we will choose the corresponding polynomial space P r Λ k (T ) or P − r Λ k (T ) to be used as shape functions. The other ingredient needed to define the finite element space is a set of degrees of freedom for the shape function spaces, that is, a basis for the dual space, in which each degree of freedom is associated with a particular subsimplex. When a subsimplex is shared by more than one simplex in the triangulation, we will insist that the degrees of freedom associated with that subsimplex be single-valued in a sense made precise below, and this will determine the interelement continuity.
The degrees of freedom for P r Λ k (T ) and P − r Λ k (T ) which we shall associate to a d-dimensional subsimplex f of T will be of the following form: for some (d − k)form η on f , the functional will be ω → f tr T,f ω ∧ η. The span of all the degrees of freedom associated to f is a subspace of the dual space of the shape function space, and so we obtain a decomposition of the dual space of the shape functions on T into a direct sum of subspaces indexed by the subsimplices of T . It is really this geometric decomposition of the dual space that determines the interelement continuity rather than the particular choice of degrees of freedom, since we may choose any convenient basis for each space in the decomposition and obtain the same assembled finite element space.
The geometric decompositions of the dual spaces of P r Λ k (T ) or P − r Λ k (T ) are given specifically in the following theorem, which is proven in [8, Sections 4.5 and 4.6].
Theorem 5.5. Let r, k, and n be integers with 0 ≤ k ≤ n and r > 0, and let T be an n-simplex in R n .
1. To each f ∈ ∆(T ), associate a space W k r (T, f ) ⊂ P r Λ k (T ) * : via the obvious correspondence, and Note that the spaces W k r (T, f ) and W k− r (T, f ) vanish if dim f < k. Note also that the dual space of P r Λ k (T ) is expressed in terms of spaces in the P − family, and vice versa. This intimate connection between the P and P − families of spaces of polynomial differential forms is most clearly seen in the following algebraic proposition, which is closely related to Theorem 5.5, and is also proved in [8, Sections 4.5 and 4.6].
Lemma 5.6. With r, k, n as above, With the decompositions in Theorem 5.5, we can define the finite element spaces. Thus P r Λ k (T ) consists of all forms ω ∈ L 2 Λ k (Ω) such that ω| T belongs to the shape function space P r Λ k (T ) for all T ∈ T , and for which the quantities f tr f ω ∧ η are single-valued for all f ∈ ∆(T ) and all η ∈ P − r+k−dim f Λ dim f −k (f ). More precisely, this means that if f is a common face of T 1 , T 2 ∈ T , then for all such f and η. The P − family of spaces is defined analogously.
The degrees of freedom determine the amount of interelement continuity enforced on the finite element space. Of course we need to know that the assembled spaces belong to HΛ k (Ω). In fact, the degrees of freedom we imposed enforce exactly the continuity needed, as shown in the following theorem. This is proved in [8, Section 5.1], where the equations below are taken as definitions, and it is shown that the assembly process leads to the same spaces.
Next we note that these spaces of finite element differential forms can be collected in subcomplexes of the de Rham complex. In view of Theorem 5.7, we have dP r Λ k (T ) ⊂ P r−1 Λ k+1 (T ) and dP − r Λ k (T ) ⊂ P − r Λ k+1 (T ). Corresponding to the resolutions defined in Section 5.1.4, we obtain 2 n−1 de Rham subcomplexes for each value of r and each mesh T . Each complex begins with the space P r Λ 0 (T ) = P − r Λ 0 (T ). The maps making up the subcomplexes are all one of the following types: for some s and some k. In Section 5.5, we shall show that all the subcomplexes admit bounded cochain projections. These 2 n−1 complexes, which are distinct for r ≥ n, are linearly ordered by inclusion. The maximal complex is in this complex are referred to as the higher order Whitney forms, since for r = 1, this is exactly the complex introduced by Whitney [88]. The minimal complex with this starting space is the complex of polynomial differential forms 0 → P r Λ 0 (T ) This complex was used extensively by Sullivan [83,84] and is sometimes referred to as the complex of Sullivan-Whitney forms [15]. It was introduced into the finite element literature in [37]. The intermediate complexes involve both higher order Whitney spaces and full polynomial spaces. Finally, we note that the degrees of freedom used to define the space determine a canonical projection I h : CΛ k (Ω) → P r Λ k (T ). Namely, I h ω ∈ P r Λ k (T ) is determined by Similar considerations apply to P − r Λ k (T ). The canonical projection can be viewed as a map from the (sufficiently) smooth de Rham complex to any one of the complexes built from the maps (73). The degrees of freedom were chosen exactly so that the canonical projection is a cochain map, i.e., commutes with d. This is verified using Stokes's theorem. See [8, Theorem 5.2].
Remark. Given a simplicial triangulation T , we have defined, for every form degree k and every polynomial degree r, two finite element spaces of differential k-forms of degree at most r: P r Λ k (T ) and P − r Λ k (T ). In n = 2 and n = 3 dimensions, we may use proxy fields to identify these spaces of finite element differential forms with finite element spaces of scalar and vector functions. The space P r Λ 0 (T ) corresponds to the Lagrange elements [33], and the spaces P r Λ n (T ) and P − r Λ n (T ) correspond to the space of discontinuous piecewise polynomials of degree ≤ r and ≤ r − 1, respectively. When n = 2, the spaces P r Λ 1 (T ) and P − r Λ 1 (T ) correspond to the Brezzi-Douglas-Marini H(div) elements of degree ≤ r, introduced in [26] and the Raviart-Thomas H(div) elements of degree ≤ r − 1 introduced in [76]. These spaces were generalized to three dimensions by Nédélec [72], [73]. The spaces P r Λ 1 (T ), P r Λ 2 (T ), P − r Λ 1 (T ), and P − r Λ 2 (T ) then correspond to the Nédélec 2nd kind H(curl) and H(div) elements of degree ≤ r and the Nédélec 1st kind H(curl) and H(div) elements of degree ≤ r − 1, respectively.

Computational bases.
This subsection relates to the implementation of the P and P − families of finite element differential forms. It is not essential to the rest of the paper. Because the spaces P − r Λ k (T ) and P r Λ k (T ) were constructed through a finite element assembly procedure (shape functions and degrees of freedom), we have at hand a basis for their dual spaces. Consider the space P − r Λ k (T ), for example, with reference to the decomposition of the dual space of the shape function space P − r Λ k (T ) given in Theorem 5.5. Choose any f ∈ ∆(T ) of dimension d ≥ k, and choose any convenient basis for P r+k−d Λ d−k (f ). For each of the basis functions η, we obtain an element of P − r Λ k (T ) * : (This is meaningful since the integral is, by construction, single-valued.) Taking the union over f ∈ ∆(T ) of the sets of elements of P − r Λ k (T ) * obtained in this way, gives a basis for that space. An interesting case is that of the Whitney forms P − 1 Λ k (T ). In this case there is exactly one dual basis function for each f ∈ ∆ k (T ), namely ω → f tr f ω.
For computation with finite elements, we also need a basis for the finite element space itself, not only for the dual space. One possibility which is commonly used is to use the dual basis to the basis for the dual space just discussed. Given an element of the basis of the dual space, associated to some f ∈ ∆(T ), it is easy to check that the corresponding basis function of the finite element space vanishes on all simplices T ∈ T which do not contain f . Thus we have a local basis, which is very efficient for computation.
In the case of the Whitney space P − 1 Λ k (T ), the dual basis can be written down very easily. We begin with the standard dual basis for P 1 Λ 0 (T h ): the piecewise linear function λ i associated to the vertex x i is determined by λ i (x j ) = δ ij (after picking some numbering x 1 , . . . , x N of the vertices). Then given a k-face f with vertices x σ(0) , . . . , x σ(k) we define the Whitney form [88, p. 229] For higher degree finite element spaces it does not seem possible to write down the dual basis explicitly, and it must be computed. The computation comes down to inverting a matrix of size d × d where d is the dimension of the space of shape functions. This can be carried out once on a single reference simplex and the result transferred to any simplex via an affine transformation.
An alternative to the dual basis, which is often preferred, is to use a basis for the finite element space which can be written explicitly in terms of barycentric coordinates. In particular, for the Lagrange finite element space P r Λ 0 (T ) consisting of continuous piecewise polynomials of degree at most r, one often uses the Bernstein basis, defined piecewise by monomials in the barycentric coordinates, instead of the dual or Lagrange basis. It turns out that explicit bases analogous to the Bernstein basis can be given for all the finite element spaces in the P and P − families, as was shown in [10]. Here we content ourselves with displaying a few typical cases.
Bases for the spaces P − r Λ 1 (T ) and P r Λ 2 (T ) are summarized in Tables 5.1 and 5.2, respectively, for n = 3 dimensions and polynomial degrees r = 1, 2, and 3. To explain the presentation, we interpret the second line of Table 5.1. We are assuming that T is a triangulation of a 3-dimensional polyhedron. The table indicates that for the space P − 2 Λ 1 (T ), there are two basis functions associated to each edge of T and two basis functions associated to each 2-dimensional face of T . If an edge has vertices x i and x j with i < j, then on an simplex T containing the edge, the corresponding basis functions are given by λ i φ ij and λ j φ ij , where λ i and λ j are the barycentric coordinates functions on T equal to 1 at the vertices x i and x j , respectively, and φ ij = λ i dλ j − λ j dλ i is the Whitney form associated to the edge. Similarly, if T contains a face with vertices x i , x j , x k , i < j < k, then on T the two basis functions associated to the face are given by λ k φ ij and λ j φ ik .  r

Approximation properties.
To apply the general theory for approximation of Hilbert complexes to the finite element exterior calculus, outlined above, we need to construct bounded cochain projections from L 2 Λ k (Ω) onto the various finite element spaces. The error estimates given in Theorems 3.9 and 3.11 bound the error in the finite element solution in terms of the approximation error afforded by the subspaces. In this subsection, we show that the spaces P r Λ k (T ) and P − r Λ k (T ) provide optimal order approximation of differential k-forms as the mesh size tends to zero, where optimal order means that the rate of convergence obtained is the highest possible given the degree of the piecewise polynomials and the smoothness of the form being approximated. In the next subsection, we shall construct L 2bounded cochain projections which attain the same accuracy.
Let {T h } be a family of simplicial triangulations of Ω ⊂ R n , indexed by decreasing values of the mesh parameter h given as h = max T ∈T h diam T . We will assume throughout that the family {T h } is shape regular, i.e., that the ratio of the volume of the circumscribed to the inscribed ball associated to any element T is bounded uniformly for all the simplices in all the triangulations of the family. Parts of the construction below simplify if we assume, in addition, that the family {T h } is quasiuniform, i.e., that the ratio h/ diam T is bounded for all T ∈ T h , uniformly over the family. However, we do not require quasi-uniformity, only shape-regularity.
Throughout this and the following subsection, Λ k−1 h will denote a subspace of HΛ k−1 (Ω) and Λ k h a subspace of HΛ k (Ω). Motivated by Section 5.1.4, for each r = 0, 1, . . . , we consider the following possible pairs of spaces: Note that Λ k h contains all polynomials of degree r, but not r + 1 in each case. As we have seen in Section 5.2, the degrees of freedom define canonical projection operators I h = I k h mapping CΛ k (Ω) boundedly to Λ k h and commuting with the exterior derivative. However, the operators I h do not extend boundedly to all of L 2 Λ k (Ω), or even to HΛ k (Ω) and so do not fulfill the requirements of bounded cochain projections of the abstract theory developed in Section 3. In this subsection, we will describe the Clément interpolant, which is a modification of the canonical projection which is bounded on L 2 Λ k and can be used to establish the approximation properties of the finite element space. However, the Clément interpolant is not a projection and does not commute with the exterior derivative. So in the next subsection we construct a further modification which regains these properties and maintains the approximation properties of the Clément interpolant.
The difficulty with the canonical interpolation operators I h is that they make use of traces onto lower dimensional simplexes, and as a consequence, they cannot be extended boundedly to L 2 Λ k (Ω) (except for k = n). The Clément interpolant [34] is a classical tool of finite element theory developed to overcome this problem in the case of 0-forms. To define this operator (for general k), we need some additional notation. For any f ∈ ∆(T h ), we let Ω f ⊂ Ω be the union of elements containing f : and let P f : L 2 Λ k (Ω f ) → P r Λ k (Ω f ) be the L 2 projection onto polynomial k-forms of degree at most r. For ω ∈ L 2 Λ k (Ω), we determine the Clément interpolant The Clément interpolant is local in the sense that for any In fact, by scaling, as in [8, Section 5.3], for example, it can be seen that where the constants c 0 , c 1 may depend on the polynomial degree r and the dimension n, but are independent of T ∈ T h and h, thanks to the shape regularity assumption. Therefore, the Clément interpolant is uniformly bounded in (with the constant independent of h).
Another key property of the Clément interpolant is that it preserves polynomials locally in the sense thatĨ h ω| T = ω| T if ω ∈ P r Λ k (T * ). It is then standard, using the Bramble-Hilbert lemma [23] and scaling, to obtain the following error estimate.
. Then there is a constant c, independent of h, such that the Clément interpolant Note that the estimate implies that any sufficiently smooth k-form is approximated by elements Λ k h with order O(h r+1 ) in L 2 . Since the polynomial spaces used to construct Λ k h contain P r Λ k but not P r+1 Λ k , this is the optimal order of approximation.

5.5.
Bounded cochain projections. The Clément interpolantĨ h is both uniformly bounded in L 2 and gives optimal error bounds for smooth functions. However, it is not a bounded cochain projection in the sense of the theory of Section 3. Indeed, it is neither a projection operator-it does not leave Λ k h invariant-nor does it commute with the exterior derivative. Therefore, to construct bounded cochain projections, we consider another modification of the canonical projection I h , in which the operator I h is combined with a smoothing operator.
This construction-key ingredients of which were contributed by Schöberl [79] and Christiansen [31]-was discussed in detail in [8] (where it was called a smoothed projection), under the additional assumption that the family of triangulations {T h } is quasi-uniform, and then in [32] in the general shape regular case. Therefore, we will just give a brief outline of this construction here.
Let Λ k h be one of the spaces P − r+1 Λ k (T h ) or, if r ≥ 1, P r Λ k (T h ), and let I h : CΛ k (Ω) → Λ k h be the corresponding canonical interpolant. To define an appropriate smoothing operator, we let ρ : R n → R be a nonnegative smooth function with support in the unit ball and with integral equal to one. In the quasi-uniform case, we utilize a standard convolution operator of the form ω → |y|≤1 ρ(y)ω(x + δy) dy, where the domainΩ ⊃ Ω is sufficiently large so that the convolution operator is well defined. We set the smoothing parameter δ = h, where the proportionality constant > 0 is a parameter to be chosen. In the more general shape regular case, we need to generalize this operator slightly. As in [32], we introduce a Lipschitz continuous function g h : Ω → R + , with Lipschitz constant bounded uniformly in h, such that g h (x) is uniformly equivalent to diam T for T ∈ T h and x ∈ T . Hence, g h | T approximates diam T . The appropriate smoothing operator in the general case, mapping L 2 Λ k (Ω) into CΛ k (Ω), is now given by where the map Φ y h : R n → R n is defined by Φ y h (x) = x + g h (x)y. Since this smoothing operator is defined as an average of pullbacks, it will indeed commute with the exterior derivative. By combining it with an appropriate extension operator E : L 2 Λ k (Ω) → L 2 Λ k (Ω), constructed such that it commutes with the exterior derivative, we obtain an operator R h : L 2 Λ k (Ω) → CΛ k (Ω). The operator Q h = I h • R h maps L 2 Λ k (Ω) into Λ k h , and commutes with the exterior derivative. Furthermore, for each > 0, the operators Q h are bounded in L 2 , uniformly in h.
However, the operator Q h is not invariant on the subspace Λ k h . The remedy to fix this is to establish that the operators Q h | Λ k h : Λ k h → Λ k h converge to the identity in the L 2 operator norm as tends to zero, uniformly in h. Therefore, for a fixed > 0, taken sufficiently small, this operator has an inverse J h : Λ k h → Λ k h . The desired projection operator, π h = π k h , is now given as π h = J h • Q h . The operators π h are uniformly bounded with respect to h as operators in L(L 2 Λ k (Ω), L 2 Λ k (Ω)). Furthermore, since they are projections onto Λ k h , we obtain (75) ω Based on these considerations, we obtain the the following theorem. Theorem 5.9. 1. Let Λ k h be one of the spaces P − r+1 Λ k (T h ) or, if r ≥ 1, P r Λ k (T h ) and π k h : L 2 Λ k (Ω) → Λ k h the smoothed projection operator constructed above. Then π k h is a projection onto Λ k h and satisfies 3. Let Λ k−1 h and Λ k h be taken as in (74) and π k−1 h and π k h the corresponding smoothed projections. Then dπ k−1 h = π k h d.
Proof. We have already established that the π k h are uniformly bounded projections which commute with d. The error estimate in the first statement of the theorem then follows from (75). Statement 2 follows from 1 and 3. 5.6. Approximation of the de Rham complex and the Hodge Laplacian. Let Λ n h be the space of piecewise polynomial n-forms of degree at most s for some nonnegative integer s. We showed in Section 5.2 that there are 2 n−1 distinct discrete de Rham complexes Thus the projections give a cochain projection from the de Rham complex to the discrete de Rham complex, and so induce a surjection on cohomology. Since lim h→0 ω −π h w HΛ k (Ω) = 0, we can use Theorem 3.4 (with V = HΛ and V h = Λ h ) to see immediately that for h sufficiently small, this is in each case an isomorphism on cohomology. The simplest finite element de Rham complex is the complex of Whitney forms, Via the basic Whitney forms φ f , the space P − 1 Λ k (T h ) is isomorphic to the space of simplicial cochains of dimension k associated to the triangulation T h . In this way, the complex (77) of Whitney forms is isomorphic, as a cochain complex, to the simplicial cochain complex, and its cohomology is isomorphic to the simplicial cohomology of Ω. Thus Theorem 3.4, when applied to the Whitney forms, gives an isomorphism between the de Rham cohomology and simplicial cohomology for a sufficiently fine triangulation. Since the simplicial cohomology is independent of the triangulation (being isomorphic to the singular cohomology), for the complex of Whitney forms, the isomorphism on cohomology given by Theorem 3.4 holds for any triangulation, not just one sufficiently fine. This proves de Rham's theorem on the equality of de Rham and simplicial cohomology.
For all of the discrete complexes (76), as for the Whitney forms, the cohomology is independent of the triangulation (and equal to the de Rham cohomology). To see this, following [31], we consider the Whitney forms complex (77) as a subcomplex of (76). The canonical projections I h define cochain projections. Note that the I h are defined on the finite element spaces Λ k h , because all of the trace moments they require are single-valued on Λ k h . From the commuting diagram → 0, we conclude that the cohomology of the top row, which we have already seen to be an image of the de Rham cohomology, maps onto the cohomology of the bottom row, which is isomorphic to the de Rham cohomology. Hence the dimension of all the corresponding cohomology groups are equal and both cochain projections induce an isomorphism on cohomology.
Consider now the numerical solution of the Hodge Laplacian problem (59). We approximate this by a Galerkin method: For the finite element spaces, we may choose Λ k−1 h and Λ k h to be any of the pairs given in (74). We have verified that d maps Λ k−1 h into Λ k h and this map can be extended to a full subcomplex of the de Rham complex admitting a bounded cochain projection. Therefore we may combine the error estimates derived in the abstract setting in Sections 3.5 and 3.6 and approximation estimates from Section 5.5 to obtain convergence and rates of convergence. The rates of convergence are limited by three factors: (1) the smoothness of the data f , (2) the amount of elliptic regularity (determined by the smoothness of the domain), and (3) the degree of complete polynomials contained in the finite element shape functions. Specifically, from Theorem 3.11 we get the following result. Assume the regularity estimate holds for 0 ≤ t ≤ t max . Then for 0 ≤ s ≤ t max the following estimates hold: Thus, the error in each case is the optimal order allowed by the subspace if we have sufficient elliptic regularity.
Concerning the eigenvalue problem, applying Theorem 3.21, we immediately obtain the following error estimate: This estimate does not apply in case Λ k h = P − 1 Λ k (T ), i.e., the Whitney forms. In that case, we get the estimate:

Variations of the de Rham complex
In Sections 4 and 5, we applied the abstract theory of Section 3 to the Hilbert complex obtained by choosing W k to be L 2 Λ k (Ω) and d k to be the exterior derivative with domain V k = HΛ k (Ω). This led to finite element approximations of certain boundary value problems for the Hodge Laplacian, and a variety of related problems (B and B * problems, eigenvalue problems). In this section, we show how some small variations of these choices lead to methods for important related problems, namely generalizations of the Hodge Laplacian with variable coefficients and problems with essential boundary conditions. 6.1. Variable coefficients. In this section, d k is again taken to be the exterior derivative with domain V k = HΛ k (Ω) and W k is again taken to be the space L 2 Λ k , but W k is furnished with an inner product which is equivalent to, but not equal to, the standard L 2 Λ k inner product. Specifically, we let a k : W k → W k be a bounded, symmetric, positive definite operator and choose With this choice, we can check that d * k = (d k−1 ) * is given by d * k u = (a k−1 ) −1 δ k a k u, with the domain of d * k given by . The harmonic forms are then As an example of the utility of this generalization, we consider the B * 1 eigenvalue problem in 3 dimensions. We take a 1 = ε and a 2 = µ −1 where ε and µ are symmetric positive definite 3 × 3 matrix fields. Then the B * 1 eigenvalue problem is curl µ −1 curl u = λεu. This is the standard Maxwell eigenvalue problem where ε is the dielectric tensor of the material and µ is its magnetic permeability. These are scalar for an isotropic material, but given by matrices in general. They are constant for a homogeneous material, but vary if the material is not constant across the domain.
We notice that the Hilbert complex in this case is exactly the de Rham complex (52), except that the spaces are furnished with different, but equivalent, inner products. Thus the finite element de Rham subcomplexes we constructed in Section 5 apply for this problem, and satisfy the subcomplex and bounded cochain projection properties. In other words, the same finite element spaces which were developed for the constant coefficient Hodge Laplacian, can be applied equally well for variable coefficient problems.
6.2. The de Rham complex with boundary conditions. For another important application of the abstract theory of Section 3, we again choose W k = L 2 Λ k (Ω) and d k to be the exterior derivative, but this time we take its domain V k to be the spaceHΛ k (Ω). This space is dense in W k , and is a closed subspace of HΛ k , so is complete, and hence d is again a closed operator with this domain. Thus we obtain the following complex of Hilbert spaces, called the de Rham complex with boundary conditions: In Section 4.2, we used the fact that HΛ k ∩H * Λ k is compactly included in L 2 Λ k . Replacing k by n − k and applying the Hodge star operator, we conclude the analogous result, thatHΛ k ∩ H * Λ k is compactly included in L 2 Λ k . Then, just as in Section 4.2, all the results of Section 3 apply to the de Rham complex with boundary conditions. In particular, we have the Hodge decomposition of L 2 Λ k and ofHΛ k , the Poincaré inequality, well-posedness of the mixed formulation of the Hodge Laplacian, and all the approximation results established in Sections 3.3-3.6. We denote the spaces of k-cocycles, k-coboundaries, and harmonic k-forms as Now, from (53), we have that (79) dω, µ = ω, δµ , ω ∈HΛ k−1 , µ ∈ H * Λ k . Therefore This gives a concrete characterization of the harmonic forms analogous to (58): Now, from (54), (55), and (50), the Hodge star operator maps the set of ω ∈ HΛ n−k with dω = 0 isomorphically onto the set of µ ∈ H * Λ k with δµ = 0, and maps the set of ω ∈H * Λ n−k with δω = 0 isomorphically onto the set of µ ∈HΛ k with dµ = 0. Comparing (58) and (80), we conclude that the Hodge star operator maps H n−k isomorphically ontoH k . This isomorphism is called Poincaré duality. In particular, dimH k is the (n − k)th Betti number of Ω. Finally, from (79), we can characterize d * , the adjoint of the exterior derivative d viewed as an unbounded operator L 2 Λ k−1 → L 2 Λ k with domainHΛ k . Namely, in this case d * has domain H * Λ k and coincides with the operator δ defined in (50). By contrast, when we took the domain of d to be all of HΛ k , the domain of d * turned out to be the smaller spaceH * Λ k .
We now consider the Hodge Laplacian problem given in the abstract case by (18). In this case, we get that (σ, u, p) ∈HΛ k−1 ×HΛ k ×H k is a solution if and only if tr σ = 0, tr u = 0 on ∂Ω, Note that now both boundary conditions are essential in the sense that they are implied by membership in the spacesHΛ k−1 andHΛ k where the solution is sought, rather than by the variational formulation.
To make things more concrete, we restrict to a domain Ω ⊂ R 3 , and consider the Hodge Laplacian for k-forms, k = 0, 1, 2, and 3. We also discuss the B * and B problems given by (19) and (20) for each k. Again we have d 0 = grad, d 1 = curl, d 2 = div, and the de Rham complex (78) is now realized as 6.2.1. The Hodge Laplacian for k = 0. For k = 0, the boundary value problem (81)-(83) is the Dirichlet problem for the ordinary scalar Laplacian. The space HΛ −1 is understood to be 0, so σ vanishes. The spaceH 0 of harmonic 0-forms vanishes. The first differential equation of (81) vanishes, and the second gives Poisson's equation − div grad u = f in Ω. Similarly, the first boundary condition in (82) vanishes, while the second is the Dirichlet condition u = 0 on ∂Ω. The side condition (83) is automatically satisfied. Nothing additional is obtained by considering the split into the B * and B subproblems, since the latter is trivial. 6.2.2. The Hodge Laplacian for k = 1. In this case the differential equations and boundary conditions are (84) σ = − div u, grad σ + curl curl u = f − p in Ω, σ = 0, u × n = 0 on ∂Ω, which is a formulation of a boundary value problem for the vector Laplacian curl curl − grad div. The solution is determined uniquely by the additional condition that u be orthogonal toH 1 = H 2 , i.e., those vector fields satisfying curl p = 0, div p = 0 in Ω, p × n = 0, on ∂Ω.
The B * 1 problem (19) is defined for L 2 vector fields f which are orthogonal to both gradients of functions inH 1 (Ω) and the vector fields inH 1 . In that case, the solution to (84) has σ = 0 and p = 0, while u satisfies curl curl u = f, div u = 0 in Ω, u × n = 0 on ∂Ω.
The orthogonality condition u ⊥H 1 again determines the solution uniquely.
Next we turn to the B 1 problem. For source functions of the form f = grad F for some F ∈H 1 , (84) reduces to the problem of finding σ ∈H 1 and u ∈ B 1 = gradH 1 such that: σ = − div u, grad σ = f in Ω.
The differential equations may be simplified to − grad div u = f and the condition that u ∈ B 1 can be replaced by the differential equation curl u = 0 and the boundary condition u × n = 0, together with orthogonality toH 1 . Now grad(σ − F ) = 0 and σ − F = 0 on ∂Ω, so σ = F . Hence, we may rewrite the system as − div u = F, curl u = 0 in Ω, u × n = 0, on ∂Ω, which, again, has a unique solution subject to orthogonality toH 1 . 6.2.3. The Hodge Laplacian for k = 2. The differential equations and boundary conditions are σ = curl u, curl σ − grad div u = f − p in Ω, σ × n = 0, u · n = 0 on ∂Ω.
This is again a formulation of a boundary value problem for the vector Laplacian curl curl − grad div, but with different boundary conditions than for (84), and this time stated in terms of two vector variables, rather than one vector and one scalar. This time uniqueness is obtained by imposing orthogonality toH 2 = H 1 , the space of vector fields satisfying curl p = 0, div p = 0 in Ω, p · n = 0, on ∂Ω.
The B * 2 problem arises for source functions of the form f = grad F for some F ∈ H 1 which we may take to have integral 0. We find σ = 0, and u solves − div u = F, curl u = 0 in Ω, u · n = 0, on ∂Ω, i.e., the same differential equation as for B 1 , but with different boundary conditions and, of course, now uniqueness is determined by orthogonality toH 2 .
If div f = 0, f · n = 0 on ∂Ω, and f ⊥H 2 , we get the B 2 problem for which the differential equations are σ = curl u, curl σ = f and the condition div u = 0 arising from the membership of u in B 2 . Thus u solves curl curl u = f, div u = 0 in Ω, u · n = 0, on ∂Ω, the same differential equation as for B * 1 , but with different boundary conditions. 6.2.4. The Hodge Laplacian for k = 3. The spaceH 3 consists of the constants (or, if Ω is not connected, of functions which are constant on each connected component). The Hodge Laplacian problem is σ = − grad u, div σ = f − p in Ω, σ · n = 0 on ∂Ω, Ω u q dx = 0, q ∈H 3 , which is the Neumann problem for Poisson's equation, where the (piecewise) constant p is required for there to exist a solution and the final condition fixes a unique solution.

The elasticity complex
In this section, we present another complex, which we call the elasticity complex, and apply it to the development of numerical methods for linear elasticity. In contrast to the other examples presented above, the elasticity complex is not a simple variant of the de Rham complex, and, in particular, it involves a differential operator of second order. However, the two complexes are related in a subtle manner, via a construction known as the Bernstein-Bernstein-Gelfand resolution, explained for example in [8,9,43] and references given there.
The equations of elasticity are of great importance in modeling solid structures, and the need to solve them in engineering applications was probably the primary reason for the development of the finite element method in the 1960s. The simplest finite element methods for elasticity are based on displacement approaches, in which the displacement vector field is characterized as the minimizer of an elastic energy functional. The design and analysis of displacement finite element methods is standard and discussed in many textbooks, (e.g., [33]). However, for more general material models, arising, for example, for incompressible materials or some viscoelastic or plastic materials, the displacement method is either not feasible or performs poorly, and a mixed approach, in which the elastic stress and the displacement are taken as unknowns, is the natural alternative. In fact, mixed finite element methods for elasticity were proposed in the earliest days of finite elements [51], and stable finite elements for the mixed formulation of elasticity have been sought for over four decades. These proved very elusive. Indeed, one of the motivations of the pioneering work of Raviart and Thomas [76] on mixed finite elements for the Laplacian, was the hope that the solution to this easier problem would pave the way to such elements for elasticity, and there were many attempts to generalize their elements to the elasticity system [4,5,63,80,81,82]. However, it was not until 2002 that the first stable mixed finite elements for elasticity using polynomial shape functions were discovered by two of the present authors [11], based on techniques and insights from the then nascent finite element exterior calculus. See also [7]. They developed and analyzed a family of methods for plane (2-dimensional) elasticity. The basic elements in the family involve a space of shape functions intermediate between quadratics and cubics for the stress, with 24 degrees of freedom, while the displacement field consists of the complete space of piecewise linear vector fields, without continuity constraints. It is these elements we will describe here.
The plane elasticity elements of [11] were generalized to three dimensions in [3]. See also [1]. These elements are, however, quite complicated (162 degrees of freedom per tetrahedron for the stress). A much more promising approach involves a variant Hilbert complex for elasticity, the elasticity complex with weak symmetry. In [8] and [9], stable families of finite elements are derived for this formulation in n-dimensions, again using techniques from the finite element exterior calculus. For the lowest order stable elements, the stress space is discretized by piecewise linears with just 12 degrees of freedom per triangle in two dimensions, 36 per tetrahedron in three dimensions, while the other variables are approximated by piecewise constants. 7.1. The elasticity system. The equations of linear elasticity arise as a system consisting of a constitutive equation and an equilibrium equation: (85) Aσ = u, div σ = f in Ω.
Here the unknowns σ and u denote the stress and displacement fields engendered by a body force f acting on a linearly elastic body which occupies a region Ω ⊂ R n . The constitutive equation posits a linear relationship between the linearized deformation or strains due to the displacement and the stresses. The equilibrium equation states that the stresses, which measure the internal forces in the body, are in equilibrium with the externally applied force. The stress field σ takes values in the space S := R n×n sym of symmetric matrices and the displacement field u takes values in V := R n . The differential operator is the symmetric part of the gradient, the div operator is applied row-wise to a matrix, and the compliance tensor A = A(x) : S → S is a bounded and symmetric, uniformly positive definite operator reflecting the material properties of the body. If the body is clamped on the boundary ∂Ω of Ω, then the appropriate boundary condition for the system (85) is u = 0 on ∂Ω. For simplicity, this boundary condition will be assumed here.
The pair (σ, u) can alternatively be characterized as the unique critical point of the Hellinger-Reissner functional The critical point is sought among τ ∈ H(div; Ω; S), the space of square-integrable symmetric matrix fields with square-integrable divergence, and v ∈ L 2 (Ω; V), the space of square-integrable vector fields. Equivalently, the pair (σ, u) ∈ H(div; Ω; S) ×L 2 (Ω; V) is the unique solution to the weak formulation of the system (85): (86) Ω (Aσ : τ + div τ · u) dx = 0, τ ∈ H(div; Ω; S), 7.2. The elasticity complex and its discretization. The discretization of the problem (86) is closely tied to discretization of the elasticity complex. For a two dimensional domain Ω, the elasticity complex takes the form The operator J is a second-order differential operator mapping scalar fields into symmetric matrix fields, namely the rotated Hessian We now show that this complex is closed, i.e., that the range of the two operators J and div are closed subspaces. The null space of J is just the space P 1 of linear polynomials and we have q H 2 ≤ c Jq L 2 for all q ∈ H 2 (Ω) with q ⊥ P 1 , which implies that the range of J is closed. We now show that div maps H(div; Ω; S) onto L 2 (Ω; R 2 ), using the fact that the scalar-valued divergence maps H 1 (Ω; R 2 ) onto L 2 (Ω). Given f ∈ L 2 (Ω; R 2 ) the latter result implies the existence of τ ∈ H 1 (Ω; R 2×2 ) with div τ = f . This τ may not be symmetric. However, for any v ∈ H 1 (Ω; R 2 ), we let curl v denote the matrix with ith row equal to curl v i = (−∂v i /∂x 2 , ∂v i /∂x 1 ). Then div(τ + curl v) = f and τ + curl v ∈ H(div; Ω; S) if div v = τ 21 − τ 12 . Thus div : H(div; Ω; S) → L 2 (Ω; R 2 ) is indeed surjective. Setting W 0 = L 2 (Ω), W 1 = L 2 (Ω; S), and W 2 = L 2 (Ω; R 2 ), we obtain a closed Hilbert complex with domain complex (87). Now we describe only the simplest discretization from [11]. The space V 0 h ⊂ V 0 = H 2 (Ω) is the finite element space whose shape functions on each triangle T are the quintic polynomials P 5 (T ) and whose degrees of freedom are the values, first, and second derivatives at each vertex, together with the mean value of the normal derivatives along each edge, which gives 21 = dim P 5 (T ) degrees of freedom per triangle. The resulting finite element space V 0 h is referred to in the finite element literature as the Argyris space, and is easily seen to be contained in C 1 (Ω), and hence in H 2 (Ω). In fact, the Argyris space is in a certain sense the simplest finite element space contained in C 1 (Ω). However the Argyris space does not coincide with the space of all C 1 piecewise quintics (which is not a finite element space, as it cannot be defined via degrees of freedom). Elements of the Argyris space have extra smoothness at the vertices, where they are C 2 .
The space V 1 h ⊂ V 1 = H(div; Ω; S) consists of those continuous piecewise cubic symmetric matrix fields whose divergence is piecewise linear. The degrees of freedom determining an element σ ∈ V 1 h are • the value of σ at each vertex, • the moments of degree 0 and 1 of the two normal components of σ on each edge, and • the moments of degree 0 of σ on each triangle. Hence, the restriction of σ to a triangle T ∈ T h is uniquely determined by 24 degrees of freedom. Finally, the space V 2 h ⊂ V 2 = L 2 (Ω; R 2 ) consists of all piecewise linear vector fields with respect to T h , without any imposed continuity. For degrees of freedom, we choose the moments of degrees 0 and 1 for each component on each triangle T . Using the degrees of freedom, we also define projections I k h v ∈ V k h defined for sufficiently smooth v ∈ V k by φ(I k h v) = φ(v) for all degrees of freedom φ.
It is straightforward to check that we have the subcomplex property JV 0 h ⊂ V 1 h and div V 1 h ⊂ V 2 h and that the interpolation operators I h := I k h give a commuting diagram (88) h → 0 However, since I 0 h involves point values of the second derivative and I 1 h involves point values, these operators do not extend to bounded operators on the Hilbert spaces V 0 and V 1 .
In order to apply the general theory of approximation of Hilbert complexes developed in Section 3, we will modify the interpolation operators I k h in (88) to obtain bounded cochain projections π k h : W k → V k h for k = 0, 1, 2. As in Section 5.5, the π k h will be obtained from I k h via smoothing. 7.3. Bounded cochain projections for the elasticity complex. First, we define an appropriate pullback of an affine map, and show that it defines a cochain map for the elasticity sequence. Let F : Ω → Ω ⊂ R 2 be an affine map, so that DF = DF (x) is a constant 2 × 2 matrix. Let B be the matrix B = det(DF )(DF ) −1 = O(DF ) T O T and define pullbacks F * = F * k : W k (Ω ) → W k (Ω) by It is straightforward to check that, as long as F is affine, F * commutes with the differential operators J and div. In other words, we obtain the commuting diagram To define a smoothing operator which maps W k into itself, we need to make the assumption that the domain Ω is star-shaped with respect to some point in its interior (which we assume, without further loss of generality, to be the origin). We then let a = max x∈∂Ω |x| −1 and define the dilation F δ : R 2 → R 2 by F δ (x) = x/(1 + aδ), δ > 0. Then F δ maps the δ-neighborhood of Ω, Ω δ , into Ω and so F * δ maps L 2 (Ω) to L 2 (Ω δ ). Composing with a standard mollification v → |z|<1 ρ(z)v(x + δz) dz, we obtain a smoothing operator R δ mapping W k into W k . Here ρ : R 2 → R is a smooth, nonnegative function supported in the unit ball and with integral equal to one. This smoothing commutes with any differential operator with constant coefficients, in particular with J and div.
The construction now proceeds as in Section 5.5, and we just give an outline of it. For simplicity, we assume a quasi-uniform family of meshes, although extension to general shape regular meshes could be made as in Section 5.5. Let > 0 be fixed.
In fact, it can be seen by a scaling argument, that the operators Q k h are bounded as operators in L(W k , W k ), uniformly in h. Furthermore, Q k h commutes with J and div, i.e., However, the operators Q k h are not projections, since they do not restrict to the identity on V k h . Therefore, we define the desired projections π k h : W k → V k h as π k h = (Q k h | V k h ) −1 Q k h . To justify this definition, we need to argue that (Q k h | V k h ) −1 exists. In fact, for sufficiently small, but fixed, (Q k h | V k h ) −1 is even uniformly bounded in L(W k h , W k h ), where W k h = V k h , but equipped with the norm of W k . This follows from the fact that the operators Q k h | V k h converge to the identity with in L(W k h , W k h ), uniformly in h, established by a scaling argument. We conclude that if the parameter > 0 is fixed, but sufficiently small, then the operators π k h are well-defined and uniformly bounded with respect to h in L(W k , W k ). Furthermore, the diagram h → 0, commutes, as a consequence of (89).
Thus, under the assumption of a star-shaped domain, we have established the existence of bounded cochain projections and so the results of Section 3.4 apply to the elasticity complex. This means that the corresponding Hodge Laplacians are well-posed and that the discretizations are stable. In particular, using the spaces V 1 h and V 2 h , we obtain a stable, convergent mixed discretization of the elasticity equations.
Remark. Stability and convergence of the mixed finite element method for elasticity using the spaces V 1 h and V 2 h described above is proven in [11] without restricting to