# A new Hodge operator in Discrete Exterior Calculus. Application to fluid mechanics

This article introduces a new and general construction of discrete Hodge operator in the context of Discrete Exterior Calculus (DEC). This discrete Hodge operator enables to circumvent the well-centeredness limitation on the mesh with the popular diagonal Hodge. It allows a dual mesh based on any interior point, such as the incenter or the barycenter. It opens the way towards mesh-optimized discrete Hodge operators. In the particular case of a well-centered triangulation, it reduces to the diagonal Hodge if the dual mesh is circumcentric. Based on an analytical development, this discrete Hodge does not make use of Whitney forms, and is exact on piecewise constant forms, whichever interior point is chosen for the construction of the dual mesh. Numerical tests oriented to the resolution of fluid mechanics problems and thermal transfer are carried out. Convergence on various types of mesh is investigated.

## Authors

• 1 publication
• 5 publications
• 2 publications
06/21/2019

### Uniform preconditioners for problems of positive order

Using the framework of operator or Caldéron preconditioning, uniform pre...
06/15/2018

### Mumford-Shah Mesh Processing using the Ambrosio-Tortorelli Functional

The Mumford-Shah functional approximates a function by a piecewise smoot...
12/03/2017

### An Introduction to Adjoints and Output Error Estimation in Computational Fluid Dynamics

In recent years, the use of adjoint vectors in Computational Fluid Dynam...
06/24/2018

### Convergence analysis of a cell centered finite volume diffusion operator on non-orthogonal polyhedral meshes

A simple but successful strategy for building a discrete diffusion opera...
03/25/2020

### Exact discrete Lagrangian mechanics for nonholonomic mechanics

We construct the exponential map associated to a nonholonomic system tha...
01/21/2020

### An arbitrary-order Cell Method with block-diagonal mass-matrices for the time-dependent 2D Maxwell equations

We introduce a new numerical method for the time-dependent Maxwell equat...
06/17/2013

### Discrete perceptrons

Perceptrons have been known for a long time as a promising tool within t...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The interest for (geometric) structure-preserving numerical integrators has grown in computing community [1, 2, 3, 4, 5, 6, 7]. This is due to the ability of these schemes in reproducing foundamental physical properties (conservation laws, …) of the equations and their robustness in long-time integration. Discrete Exterior Calculus (DEC) belongs to this family of integrators. It aims at developing a discrete version of the theory of exterior calculus, and more generally the differential geometry theory, where most equations of physics are formulated. Initially developed by Bossavit for electromagnetism in a series of papers [8, 9, 10, 11, 12, 13, 14, 15, 16], DEC was used in fluid mechanics for the resolution of Darcy’s equation and the simulation of the dynamics of ideal and viscous fluid flows in basic configurations [17, 18, 19].

The primary calculus tools of the DEC are discrete differential forms or cochains. An important advantage of DEC is that the Stokes’ theorem is naturally verified at discrete level. This is due to the construction of the discrete exterior derivative operator by duality with the boundary operator. A second advantage is that the DEC framework offers a naturally coherent discretization of derivative operators (divergence, gradient, curl in 2D and 3D) such that the usual relations

are verified at machine precision. This is due to the fact that these derivative operators are all represented by the exterior derivative in exterior calculus framework and that in DEC, the discrete obeys the relation

 ｄ2=0, (1)

just like the continuous . These properties of DEC permitted for example Elcott et al. [17] to design a circulation-preserving numerical scheme for the simulation of ideal fluid flows. These properties also avoid the apparition of spurious quantities such as an artificial mass, potential or portance in the numerical solution.

When used as discretization method, DEC can be seen as a finite volume method in exterior calculus. Indeed, in this approach, cochains are differential forms integrated over the different elements of the mesh (vertices, edges, triangles, …). Disrete exterior calculus has a finite element version, developed mainly by Arnold and his co-workers in [20, 5, 21].

A key operator in exterior calculus, needed to express constitutive laws for example, is the Hodge star operator. In a finite element approach, the discrete Hodge operator is built straightforwardly, by applying the continuous Hodge to the polynomial differential forms which constitute the basis. After a suitable inner product, this results in the mass matrix. In DEC, to which the present article is limited, defining a discrete Hodge is less natural because one has to deal with already intregrated forms. In particular, in order to ensure the bijectivity of the discrete Hodge, a dual mesh is necessary in DEC, whereas such a notion is not required in the finite element approach. The circumcentric dual is a popular choice of dual mesh, due to the orthogonal nature of simplicial cells and their circumcentric duals. This orthogonality permits a simple construction of a discrete Hodge, having a diagonal matrix representation [13, 22]. The diagonal entries of this discrete Hodge, sometimes called circumcentric or diagonal Hodge, are simply the ratios of the volumes of dual cells and primal simplices.

The diagonal Hodge was initialy designed for completely well-centered simplicial meshes, that are meshes where each simplex contains its circumcenter. In practice, it is difficult to generate this type of triangulation [23]. It can be done for some specific simple geometries in [24, 25, 26]. For more complex geometries, VanderZee et al. propose a reprocessing algorithm making a 2D mesh well-centered but it is not very practical and may not work in 3D [27]. Hirani and his co-workers extended the use of the diagonal Hodge to Delaunay and to pairwise non-Delaunay meshes, by introducing signed elementary dual volumes [28, 29]. However, this approach may still lead to zero-volume dual edges and cells when positive contributions equalize negative ones. A reprocessing is then necessary. Mullen et al. [30] proposed an optimization of the triangulation to the diagonal Hodge. To this aim, they introduce weighted duals and weighted circumcenter, and construct the mesh by minimising the error of the diagonal Hodge. This technique generates a well-shaped mesh, but not well centered. Moreover, it may become expensive, especially for problems with time-varying domains such as fluid-structure interaction. Note also that all these discrete Hodge operators require that the dual mesh is based on the circumcenter.

Another existing discrete Hodge operator, studied by Bossavit from the beginning of DEC, is the Galerkin Hodge [31, 16]. It corresponds to a mass matrix, with the Whitney forms as shape functions [32, 33]. It is built within a finite element framework. A drawback of using the Galerkin Hodge in (the finite-volume flavoured) DEC is that, not only it introduces an inconsistency but also it is not exact even for constant forms. This can be shown straighforwardly in a right but not isocele triangle. A more DEC- than finite element-flavoured discrete Hodge, but still based on Whitney forms, is proposed in [34]

. It is defined through an interpolation with Whitney forms and an integration on dual simplices and is generally called Whitney Hodge. Another discrete Hodge operator which also uses Whitney forms is proposed in

[35]. Called geometric Hodge, it does not need an interpolation nor integration concept. For symmetry reasons, it requires a dual mesh based on the barycenters. In fact, as noticed in [36] in 2D, the geometric Hodge is a symmetric approximation of the Galerkin Hodge where the integration is approximated by a one point quadrature with evaluation only at the barycenter. Like the usual Galerkin and the Whitney Hodge, the geometric Hodge is represented by a sparse but not diagonal matrix. However, numerical comparisons in [36] shows that the additional computational time induced by the lake of diagonal structure is not very high.

In the present article, we propose a new construction of a discrete Hodge operator which comply the following requirements. First, it does not need a reprocessing of the primal mesh. This condition is important since, in some situations, modifying the mesh may be harmfull. It may for instance affect the shape of surfaces. Moreover, as mentioned, the induced cost increase may become important in time-varying domains. Second, the construction must be general enough such that the centers on which the dual mesh is based can be changed at convenience, as long as the dual mesh has no degenerate element. For example, the construction of the discrete Hodge must be valid for both a barycentric and an incentric duals. Neither the diagonal nor the geometric Hodge is appropriate for this task since with these operators, the dual mesh is imposed. Lastly, the Hodge must be exact on piecewise constant differential forms, whichever choice is made on the dual mesh.

By proposing such a construction, we intend to develop a discrete Hodge which adapts to the mesh in order to minimize some quantity (the global error, …) or to comply to other physical requirements, by choosing algorithmically the most suitable center. These center may even vary from a simplex to another. This optimization process, which will be tackled in a future work, would contrast with that in [30] where it is the mesh which adapts to the (diagonal) Hodge operator.

In section 2, a brief reminder on DEC, on the continuous Hodge operator and on the diagonal Hodge is done. A deeper introduction to DEC can be found in the papers of Bossavit [8, 9, 10, 11, 12, 13, 14, 15, 16] and in [22, 37, 38] for instance. The development of the new discrete Hodge, followed by illustrative examples, a basic preliminary error analysis and a discussion, is presented in section 3. We end up with numerical experiments, carried out on Poisson equations and on isothermal and non-isothermal fluid dynamics problems, in section 4. Various types of mesh (well-centered mesh, right triangulation and non-Delaunay meshes) will be considered. Convergence analyses will be carried out. Although our goal is not to compare the new discrete Hodge to existing ones, results given by the diagonal Hodge are presented as reference, when available.

## 2 Review on DEC

DEC is a theory of exterior calculus on a discretized domain, which preserves Stokes’ theorem and the exactness of relation (1). In what follows, we make a reminder on oriented simplicial discretization of a domain. Next, since the primary tools of exterior calculus are differential forms, we briefly present how to discretize them into cochains. Lastly, the Hodge star operator and its discretization into the diagonal Hodge is reminded.

### 2.1 Discretization of the domain

Consider a system of partial differential equations defined on an

-dimensional spatial domain for some . The domain is subdivised into -dimensional simplices. Recall that a -simplex is the convex hull of points. For instance, if then is subdivised into tetrahedra. We denote the set composed of these -dimensional simplices, together with their faces, the faces of the faces, and so on untill the vertices. For example, in 3D, is composed of tetrahedra (3-simplices), triangles (2-simplices), edges (1-simplices) and vertices (0-simplices). In 2D meshes, the top-dimensional simplices are triangles. As an example, Figure 1 presents a simplicial discretization of a curved surface in . Some regularity conditions are also generally usefull [22, 37, 39].

To each element of is assigned an orientation, which can simply be defined from an ordering of its nodes. Each simplex has two possible orientations. Two adjacent top-dimensional simplices are required to have the same orientation. By contrast, the orientation of each lower dimensional simplex is arbitrary. A sample 2D oriented complex is presented in Figure 2. For simplicity, and since it is sufficient for our applications, we assume that all the vertices have the same (positive) orientation.

We denote the set of oriented -dimensional simplices of

 Kk={σ∈K, σ oriented, dimσ=k}

and

the vector space spanned by formal linear combination of elements of

 Ωk=spanKk=⎧⎨⎩c=∑σi∈Kkciσi, ci∈R⎫⎬⎭.

An element of is called a chain. In theory, the components may take any real value but, in practice, only values in have physical meaning. A value 0 means that the element does not belong to the chain, 1 means that it is present in , and when it is present but with the opposite orientation. As an example, the chain

 e1+e4+e5−e3∈Ω1 (2)

in Figure 2 constitutes a closed loop.

### 2.2 Boundary operator

The boundary of a -dimensional simplex is the sum of its oriented -dimensional faces. In this sum, each face is given a sign, depending on wether its orientation is consistent with that of the considered -simplex. For example, the boundary of the face in Figure 2 is

 ∂f1=e1+e2−e3.

The boundary of is

 ∂e1=v2−v1.

The boundary operator extends into a linear map from to by defining the boundary of a chain , for some set of indices, as follows:

 ∂⟦∑i∈Iciσi⟧=∑i∈Ici∂σi. (3)

The operator can be represented by a (sparse) matrix. For instance, the non-zero boundary operators on the mesh in Figure 2 are

 ∂|Ω2=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝101−1−100101⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,∂|Ω1=⎛⎜ ⎜ ⎜⎝−10−1001−10−10011010001−1⎞⎟ ⎟ ⎟⎠. (4)

In equation (4), is the restriction of to , which acts on -simplices.

### 2.3 Discrete differential forms and exterior derivative

A discrete -form, or a cochain, is an element of the algebraic dual of . Since the elements of form a basis of , a discrete -form is simply a map which, to each element of , assings a real number:

 ω∈Ωk(K) : Kk⟶Rσ⟼⟨ω,σ⟩. (5)

Computationally, a -cochain is represented by an array .

In DEC, a differential -form on is discretized into a -cochain by taking the pairing as an integration

 ⟨ωK,c⟩=∫cω,␣␣for all c∈Ωk

if , and as a simple evaluation

 ⟨ωK,c⟩=ω(c)

if is a point. In the sequel, the cochain will also be denoted when there is no confusion.

The discrete exterior derivative operator is defined as the dual of the boundary operator :

 ⟨ｄω,c⟩=⟨ω,∂c⟩, for all ω∈Ωk,c∈Kk+1. (6)

It is a linear map from to . It can be remarked that relation (6) would simply be the expression of Stokes’ theorem if the pairing was an integration and ω a differential form. Computationally, the matrix of the restriction of on the space of -cochains is the transpose of the matrix of :

 ｄ|Ωk(K)=∂T|Ωk+1.

Since (the boundary of a boundary is the empty set), it follows by the duality relation (6) that

 ｄ2=0. (7)

The following sequence is then, as in the continuous case, exact:

 0\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→Ω0(K)\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→Ω1(K)\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→⋯\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→Ωn−1(K)\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→Ωn(K)\makebox[14.226378pt]$ｄ$−−−−−−−−−−−−−−−→0. (8)

### 2.4 The diagonal Hodge

Having applications in fluid mechanics domain in mind, we assume that is a flat domain in , or 3, from now on, to simplify. It inherits the usual Euclidean metric and the induced right-hand oriented volume form of .

To express material laws in exterior calculus framework, the Hodge star operator is necessary. It realizes an isomorphism between and . Rather than giving its definition (which can be found for example in [40]), we simply recall that, if ,

 ⋆1=ｄx∧ｄy,⋆ｄx=ｄy,⋆ｄy=−ｄx,⋆(ｄx∧ｄy)=1. (9)

In fact, these relations completely characterize the operator . For , we have

 ⋆1=ｄx∧ｄy∧ｄz,⋆ｄx=ｄy∧ｄz,⋆ｄy=ｄz∧ｄx,⋆ｄz=ｄx∧ｄy,⋆(ｄy∧ｄz)=ｄx,⋆(ｄz∧ｄx)=ｄy,⋆(ｄx∧ｄy)=ｄz,⋆(ｄx∧ｄy∧ｄz)=1.

It can be checked that

 ⋆⋆=−Id (10)

in 2D and

 ⋆⋆=Id

in 3D, where is the identity map. Note also that the Hodge star verifies the relation

 ω(u1,…,uk)=⋆ω(uk+1,…,un) (11)

for any positively oriented orthonormal basis of the tangent space and any differential -form .

The discrete Hodge star operator should realize an isomorphism between and . This is not possible on the same grid because, generally,

 dimΩk(K)=cardKk≠dimΩn−k(K)=cardKn−k.

For example, with the mesh on Figure 2, equals 2 and if ,

 dimΩk(K)=4,dimΩn−k(K)=2.

In order to define a discrete Hodge star operator, a dual mesh is needed. This dual mesh must be constituted by cells such that there is a one-to-one correspondance between -simplices of and -dimensional cells of . One possible dual mesh is the circumcentric dual. In a 2D mesh such as in Figure 3, the circumcentric dual is such that the dual of a triangle is its circumcenter; the dual of a primal edge is the edge connecting the circumcenters of two triangles which share this primal edge; and the dual of a vertex is the 2-cell formed by connecting the circumcenters of the primal triangles which share this vertex. Figure 4 presents an example of a mesh on a square and its circumcenter-based dual. Note that the choice of circumcenters as duals of triangles implies the orthogonality between edges and their duals.

The boundary and discrete exterior derivative can be transposed to the dual mesh. The orientation on the primal mesh also induces an orientation on (see [37])222To borrow the language of Bossavit [41], the (inner-)orientation of each primal simplex, as defined in section 2.1, constitutes the outer-orientation of its dual cell.. The space of -cochains on the dual mesh is denoted .

A discrete Hodge operator, realizing an isomorphism between and , can now defined as:

 ⋆: Ωk(K)⟶Ωn−k(⋆K)

with

 ⟪⋆ω,c∗⟫|c∗|=⟪ω,c⟫|c|, for any c∈Kk, (12)

where designates the dual of the simplex . In equation (12), is a measure of , which is set to if is a point and to its Euclidean measure otherwise. Relation (12) can be understood as follows: the discrete forms and , normalized by the measure of the simplex/cell on which they are applied, represent the same density, but act in complementary-dimensional facets. Relation (12) can also be understood as a discretization of equality (11), and forming complementary geometric objects.

Definition (12) is a low-order approximation of the continuous Hodge star, which turns out to be exact in the particular case where is piecewise constant and and are flat. It induces a representation of the discrete Hodge operator in each dimension as a diagonal matrix, having ratios of cell measures as entries. For example, in a one-triangle mesh as in Figure 5, where and are the circumcenters of the triangle and the edges respectively, the Hodge matrix acting on primal 1-cochains is

 Hdiagonal=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣|e∗1||e1|000|e∗2||e2|000|e∗3||e3|⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (13)

Due to relation (10), the discrete Hodge matrix acting on dual 1-cochain is the inverse of matrix (13).

The diagonal structure of and the induced computational efficiency makes the diagonal Hodge a popular choice of discrete Hodge operator. However, to be correctly defined, it requires a well-centered mesh. This means that each -dimensional simplex of the mesh, with , must strictly contain its circumcenter. In 2D, a well-centered mesh can contain only strictly acute triangles. To figure out the problem induced by the non well-centerdness, consider a 2D domain meshed with right triangles. The mesh is not well centered. Circumcenters lie on the boundary of the triangles and some dual edges degenerate into points (see Figure 7, left, as an example), making matrix (13) singular. If the mesh contain obtuse triangles then some circumcenters are located outside the triangles. In this case, an extension, introducing a sign convention to the elementary volumes of the dual cells, has been proposed in [29] to widen the range of applicability of the diagonal Hodge, but bringing no solution for the right-triangularized mesh for example (apart from altering the mesh). Moreover, the diagonal Hodge does not allow any choice than the circumcenters as duals of top-simplices. Indeed it has been shown that if and the are not the circumcenters, the diagonal formula (13) does not lead to a convergent Hodge operator [36].

In the next section, we propose an alternative discrete Hodge operator.

## 3 Analytically constructed discrete Hodge

In the sequel, we restrict to the case of a bidimensional mesh in . We will fix a Cartesian orthonormal frame and denote and the components of a vector in this frame. In the next subsection, we present the construction of a discrete Hodge on 1-forms. We require that it is exact on locally constant forms.

### 3.1 Construction in a one-triangle mesh

To simplify, consider an oriented bidimensional mesh composed of a single triangle , as in Figure 6. In this figure, the center of the edge is an arbitrary point lying on the edge. The center of the triangle is also an arbitrary point inside . In the particular case where is the middle point of and is the barycenter (resp. circumcenter, incenter) of , the dual mesh is said barycentric (resp. circumcentric, incentric).

Let us denote and the vectors representing the oriented edges and respectively, and

 →ui=→ei∥→ei∥and→u∗i=→e∗i∥→e∗i∥

the corresponding unit vectors. These two vectors are related by a rotation:

 (u∗ixu∗iy)=(cosθi−sinθisinθicosθi)(uixuiy)

where

 θi=ˆ(→ui,→u∗i) (14)

is the angle that they form (see Figure 6).

Consider a differential 1-form

 ω=aｄx+bｄy

where and are scalar functions of the position. Its image by the continuous Hodge star is (see equation (9)):

 ⋆ω=−bｄx+aｄy.

We denote and the values of the corresponding primal and dual cochains. Assume that is constant inside the triangle. Then

 ω∗i = ∫e∗i⋆ω = ⋆ω(u∗i) |e∗i| =(−bu∗ix+au∗iy)  |e∗i| = ((auix+buiy) sinθi+(−buix+auiy) cosθi)|e∗i| = (ω(→ui) sinθi+ω(−J→ui) cosθi)|e∗i|, = (÷ω(→ei)|ei| sinθi+÷ω(−J→ei)|J→ei| cosθi)|e∗i|,

being the skew-symmetric matrix

 J=(0−110).

Knowing that , we have:

 ω∗i=|e∗i||ei|ωi sinθi  +  |e∗i||ei|ω(−J→ei) cosθi. (15)

Let’s take . As soon as the triangle is not degenerate, there exists reals and such that

 −J→e1=a21 →e2+a31 →e3.

So,

 ω∗1=|e∗1||e1| ω1 sinθ1  +  |e∗1||e1| a21ω2 cosθ1  +  |e∗1||e1| a31ω3 cosθ1. (16)

Using similar arguments for and 3, we get:

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ω∗1ω∗2ω∗3⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣|e∗1||e1|000|e∗2||e2|000|e∗3||e3|⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣sinθ1a21cosθ1a31cosθ1a12cosθ2sinθ2a32cosθ2a13cosθ3a23cosθ3sinθ3⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ω1ω2ω3⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

This relation makes the discrete Hodge star operator appear as a product of two matrices. The first one contains ratios of dual and primal volumes. In case of circumcentric dual mesh, the second matrix reduces to the identity matrix and we retrieve the usual diagonal Hodge (

13).

Of course, one can rewrite the angles and the coefficients as functions of the and . For example,

 cosθi=→ei⋅→e∗i|ei||e∗i|,sinθi=→ei˙×→e∗i|ei||e∗i|,θi∈[0,π]

where is the inner product of and . The dotted cross product is defined as the third component of the usual cross product of the extension of and to :

 →a˙×→b=axby−aybx

which can also be considered as the real-valued 2D cross product. Moreover, it is straightforward to show that, for ,

 aij=→ej⋅→ek→ei˙×→ek,

where . With these relations, the matrix of the Hodge operator acting on discrete primal 1-form is

 H=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣→e1˙×→e∗1|e1|2→e1⋅→e∗1|e1|2→e1⋅→e3→e2˙×→e3→e1⋅→e∗1|e1|→e1⋅→e2→e3˙×→e2→e2⋅→e∗2|e2|2→e2⋅→e3→e1˙×→e3→e2˙×→e∗2|e2|2→e2⋅→e∗2|e2|2→e2⋅→e1→e3˙×→e1→e3⋅→e∗3|e3|2→e3⋅→e2→e1˙×→e2→e3⋅→e∗3|e3|2→e3⋅→e1→e2˙×→e1→e3˙×→e∗3|e3|2⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

### 3.2 Examples and remarks

As illustration, consider the right triangle presented in Figure 7. The circumcentric dual (left of the figure) is singular because the dual of is empty and its measure vanishes. This leads to a singular matrix , as in the diagonal Hodge. With the barycentric dual (Figure 7, right), reduces to the following non-singular matrix:

 Hbarycenter=⎡⎢⎣÷n3m÷m6n0÷n6m÷m3n0÷n(m2−n2)6m(m2+n2)÷m(m2−n2)6n(m2+n2)÷mn3(m2+n2)⎤⎥⎦.

In an unit right triangle (), with a barycentric and an incentric dual meshes, writes:

The discrete Hodge operator for 0-forms and for 2-forms are composed of measure ratios as in the diagonal Hodge.

For a general mesh, composed of multiple triangles, an usual assembling permits to deduce the global Hodge matrix. This matrix is sparse. It takes a primal 1-cochain into a dual one. Again due to relation (10), the Hodge operator acting on dual cochains can simply be represented by the inverse of the assembled matrix. However, a full inversion of this matrix is in general expensive. To avoid this cost, an element-wise inversion will be carried out in numerical tests when needed. Moreover, the element-wise inversion maintains the sparsity of the matrix, contrarily to a full inversion. This may significantly reduce the simulation cost for realistic problems. Despite this approximation, we will see in section 4.1 that the newly defined Hodge provides competitive numerical results.

As can be remarked, the new discrete Hodge matrix is not necessarily symmetric. Yet, some authors claim that an Hogde matrix should be symmetric, since

 ⋆ω∧θ=⋆θ∧ω (17)

for two arbitrary 1-differential forms ω and θ, being the wedge product ([35, 42]). Combined to a choice of volume form, this relation represents an inner product. However, relation (17) does not imply that the discrete Hodge has to be symmetric. Likewise, a symmetry of the matrix of the Hodge operator does not help in verifying property (17). It is the combination of the discrete Hodge matrix with the discrete wedge operator between two discrete 1-forms which should be symmetric. It may also be advocated that a symmetric Hodge matrix would be more advantageous because it would result in a symmetric problem when discretizing a symmetric operator as in (17). This argument is however not convincing since the discrete wedge product may destroy the symmetry.

By definition, the symmetry of the Hodge operator matrix alone would mean, in a single-triangle mesh, that

 ⎡⎢⎣θ1θ2θ3⎤⎥⎦⋅⦅H⎡⎢⎣ω1ω2ω3⎤⎥⎦⦆=⦅H⎡⎢⎣θ1θ2θ3⎤⎥⎦⦆⋅⎡⎢⎣ω1ω2ω3⎤⎥⎦ (18)

for any pair and of primal 1-cochains, where the dot symbol designates the Euclidean scalar product of . Not only the scalar product between cochains has still to be given a sense but also there is no obvious reason to require condition (18).

In fact, the symmetry requirement may be justified in a finite element approach as proposed by Hiptmair in [42] but not in DEC. Indeed, in a finite element framework, it is the mass matrix which is called Hodge operator. Yet, this mass matrix represents, not the Hodge operator, but inner products such as (17) integrated over elements. Note also that, as recognized in [35, 43, 44], the symmetry restriction seems to be too strong in DEC context.

Which center gives more accurate results for non-constant 1-forms depends on the form. For instance, consider the differential 1-form

 ω=(x−y)(ｄx−ｄy).

Table 1 lists the error of the discrete Hodge operator on ω, being either the barycentric or the incentric Hodge. The error is defined as the -norm of the dual cochain . Two computations are carried out. The first one is on a mesh made up of the single unit triangle, having a mean edge length 1.1381. The second computation is on a square meshed with 722 right and isocele triangles (similar to Figure 8, right) where the mean edge length is . Table 1 shows that, with the considered form, the barycenter-based Hodge is more precise than the inceter-based one on the single-triangle mesh. The same conclusion can be drawn in the discretized square. However, in the latter case, the error difference is smaller. For the differential form

 ω=(x+y)(ｄx+ｄy),

the situation is reversed. Indeed, Table 1 shows that for this form, the incentric Hodge is about 1.94 times as accurate as the barycentric one, as well on the single-triangle mesh as in the right-triangularized mesh.

This very simple test shows that it is not easy to predict which center will give the best results, in terms of precision. Moreover, in more concrete problems, there may be other constraints than the accuracy to account for. The freedom to choose the center makes it conceivable to carry out an optimization of the discrete Hodge operator.

Lastly, let us remark that, in the previous development, there is no restriction on the position of the centers of the simplices, as long as the dual mesh is not degenerate. For example, the center may be outside the triangle, or even outside the domain in the case of multiple-triangle mesh.

In the following section, numerical experiments on (linear) Poisson equation and applications to (non-linear) isothermal and anisothermal fluid flow problems are presented. We will in particular focus on the convergence with the circumcenter-, the barycentrer- and the incenter-based Hodges on different types of mesh.

## 4 Numerical results

We begin with numerical tests on Poisson problem.

### 4.1 Poisson equation

In this section, we deal with Poisson equation with Dirichlet boundary conditions in a flat bidimensional domain :

 {−Δu=f in M,\paru=g\quad on ∂M, (19)

where , and the unknown are scalar functions. In exterior calculus formulation, this equation writes

 {⋆ｄ⋆ｄu=f in M,\paru=g \quad on ∂M. (20)

In the tests, is the unit square. Two meshes are considered. The first one is well-centered (see Figure 8, left) and based on a pattern proposed in [24]. It will be called “acute mesh”. On this mesh, the usual diagonal Hodge and the new Hodge operator based on barycentric and on incentric duals will be compared. The second mesh, called “right mesh” herein, is more natural and is made of right and isocele triangles (see Figure 8, right). As the diagonal Hodge cannot be defined on it, the new Hodge operator with barycentric and incentric duals will be used.

The variable is placed, as cochain, on dual vertices . This means that the unknown is the array . A consequence is that , in equation (20), is discretized into a dual 1-cochain and the right-most operator in that equation is then represeted by the inverse of the assembled Hodge matrix. As mentioned, the inversion is done element-wise numerically.

As precision indicator, we consider the relative error

 E=÷∥u−uexact∥K∗2∥uexact∥K∗2 (21)

where designates a discrete norm of a dual 0-cochain , defined as follows:

 ∥v∥2K∗2:=∑σ∗∈K∗2⟨v,σ∗⟩2|σ∗|. (22)

An object-oriented Fortran code was designed for the simulations. A direct LU solver is selected for the resolution. In a first test, the right-hand side functions and are chosen such that the exact solution is the quadratic function

 uexact=x2+y2.

Figure 9 shows the evolution of the relative error with the mean edge length . As can be seen in Figure 8(a), barycentric and incentric dual meshes provide better results than the diagonal Hodge, despite the element-wise inversion needed in the discetized Laplace operator. Note that many triangles are close to right. However, the convergence rates are almost the same with the three types of dual mesh. The left part of Table 2 reveals an almost quadratic convergence.

With the right mesh, the results are less accurate but are still very acceptable as can be observed in Figure 9b. With the barycentric dual, the relative error is about when (whereas it was with an acute mesh). The difference between the barycentric and the incentric duals are more pronounced. However, the convergence rates are similar, as can be stated in Table 2, right.

In a second test, the exact solution is

 uexact=sin(πx)sinh(πy),

such that the right-hand side vanishes. Figure 10a shows that, with the acute mesh, the three types of dual present close precisions. The barycentric solution is a slightly more accurate with the right mesh whereas the incentric solution has clearly lost precision. They have however almost the same convergence rate, as can be seen in Table 3.

From these two experiments, it is obvious that when the mesh is well-centered, one cannot predict if the diagonal Hodge is more precise than other discrete Hodge operators. It depends on the solution. By contrast, the barycentric dual gave better accuracy than the incentric dual.

In the previous tests, the equation is linear. In the next subsections, tests on non-linear systems in fluid dynamics domain are presented.

### 4.2 Incompressible fluid flow

Consider an incompressible Newtonian fluid flow of density and kinematic viscosity in a bidimensional domain. The evolution of the velocity and the pressure is governed by the Navier-Stokes equations:

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂u∂t+(u⋅∇)u+1ρgradp−νΔu=0,\pardivu=0\par (23)

To ensure a divergence-free velocity, a stream function formulation is used. Moreover, to get rid of the pressure, a curl is applied on the momentum equation. The resulting equation is rewritten in exterior calculus formulation and discretized as in [19], (see also equations (26)). The discretized stream function is placed on primal vertices and, consequently, velocity is discretized into a dual 1-cochain.

#### 4.2.1 Poiseuille flow

We consider a Dirichlet problem in the unit square with a Poiseuille flow as exact solution. The accuracy on the stream function ψ and on the velocity is measured with the relative errors defined similarly as in (21)-(22) but with the appropriate simplex/cell set in place of .

The numerical initial flow is set to zero. The simulation is run until the flow stabilizes. The convergence of the principal variable ψ with the mean edge length is presented in Figure 11. It shows a higher precision of the diagonal Hodge in the acute triangulation in the given range of , but a slightly higher speed of convergence with the new Hodge with barycentric and incentric duals. It is confirmed by the convergence rates in Table 4, second column. For the velocity however, the circumcentric dual has both a higher precision and a higher convergence speed, as can be observed in Figure 11(a) and in the third column of Table 4. Note that a first order interpolation is used for the reconstruction of the approximate velocity field from the discrete 1-form .

In the case of the right triangulation, the barycentric and the incentric duals provide similar precisions, as well for the stream function than for the velocity. For ψ, the errors are visually indistinguishable in Figure 11b. For the velocity, the convergence rate is about 1.5, which is higher than the convergence rate of the circumcentric dual in the acute triangulation case (see Table 4, last column).

#### 4.2.2 Taylor-Green vortex

We now simulate a steady Taylor-Green vortex corresponding to the vector field

 u=(−cosxsinysinxcosy) (24)

in the domain . The kinematic viscosity and the density are set to 1. As previously, the numerical initial flow is set to zero and the simulation is run until the flow stabilizes.

The error on ψ is plotted in Figure 13 against the typical edge length. It shows that the decrease of the error is the same for the three types of dual mesh, as well in the acute mesh as in the right mesh. The convergence rate is close to 2.0 in all cases, as reported in Table 5.

For the velocity, differences can be observed between the acute and the right meshes. Indeed, it can be seen in Figure 14 that the decrease is faster in the right mesh. The convergence rate is about 1.74 (for both barycentric and incentric duals) while it is close to 1.14 in the acute mesh (for the three dual types), as can be seen in Table 5.

In the next test, we deal with an anisothermal fluid flow.

### 4.3 Anisothermal fluid flow

In the case of anisothermal fluid flow with thermal expansion coefficient and thermal diffisivity , and in the limit of the Boussinesq approximation, the evolution of the velocity, pressure and temperature is governed by the equations:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂u∂t+div(u⊗u)+1ρgradp−νΔu+βgθey=0,\pardivu=0,∂θ∂t+div(uθ)−κΔθ=0 (25)

where is the gravity acceleration vector. Let . In exterior calculus framework, equation (25) writes

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂ω∂t+ιuｄω+÷1ρｄ¯¯¯p+ν⋆ｄ⋆ｄω+βgθｄy=0,⋆ｄ⋆ω=0,∂θ∂t+⋆ｄ⋆(θω)−κ⋆ｄ⋆ｄθ=0. (26)

In (26), is the interior product and is the dynamic pressure. The exterior derivative operator is applied to the first equation of (26) to get rid of the dynamic pressure. A stream-function formulation is used. This means that the equations are solved in the stream function such that instead of ω, and in the temperature θ. The time discretization scheme is similar to that in [3], page 55.

Assume that . We choose the following traveling wave solution as reference:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ux=u1eλξ/ν+κ2βabg(c+w)2(κ−ν)θ1eλξ/κuy=w−auxb,p=−ρβbgκθ1c+weλξ/κθ=θ1eλξ/κ (27)

where the traveling wave variable is

 ξ=ax+by+ct (28)

and , , , , , and are constants such that , and . The space domain is . The fluid is characterized by , , and . The gravity acceleration is set to . The values of the different constants in the reference solution are set to , , , and . The domain is discretized with a right triangulation (Figure 8, right). The mean edge length is . The numerical simulation is carried out until time such that .

Figure 15 reproduces the relative errors on the stream function along -axis. Even if this figure presents some weird plateaux due to lack of resolution (and some machine round-off), it clearly shows a good precision with both dual meshes. At worst, the error is 0.025 percent of the norm of the exact stream function along the axis. The mean error over the whole domain is still smaller. Indeed, as listed in Table 6, the overall (root mean square over the whole domain) error is only 0.008875 percent of the norm of the exact solution with the incentric dual and 3.3 times smaller with the barycentric one.

Since the velocity is computed a posteriori from the stream function, the same conclusion can be drawn on the error. Indeed, the barycentric dual is more precise along the -axis as shown in Figure 16. It is also more precise in mean in the whole domain, as can be checked in Table 6.

By contrast, for the temperature, the two dual meshes have a similar precision. The error profiles along -axis, presented in Figure 17, have the same trend. The overall errors are almost equal as can be seen in the last column of Table 6. In mean, the relative error is about 0.56 percent. These results reveals a good precision of the new discrete Hodge star on a heat transfer problem, with different choices of dual meshes. Convergence rates will be analyzed in the next section.

### 4.4 Non-structured meshes

To ease the lecture of the article, only two types of primal mesh were considered in the above tests. The first ones were well-centered meshes such that the diagonal Hodge can be used. The second type of mesh was composed of right meshes. These meshes have structured patterns. In this following section, we show that the constructed Hodge gives also satisfactory results on various non-structured and non well-centered meshes. To simplify, only results on barycentric duals are presented.

We set and choose the following exact solution:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ux=βabgθ1(c+w)(a2+b2)(x0−ξ)eλξ/κuy=w−au1b,p=−ρβbgκθ1c+weλξ/κθ=θ1eλξ/κ (29)

where, as before, is the traveling wave variable. For the simulation, the fluid characteristics are , , , . We choose , , and . The spatial domain is .

A first test is carried out on the five meshes shown in Figure 18. These meshes are increasingly fine, but are all unstructured, non well-centered and non-Delaunay. They have approximately 40 percent of non-Delaunay triangles. Here, a triangle is said non-Delaunay if its circumcircle contains other nodes than its vertices. They are darkened in the figures. Note that a non-Delaunay mesh is always non well-centered [27].

Figure 19 shows the relative overall errors on the stream function and on the temperature with the five meshes. As can be remarked, the error decreases with the mean edge length. The decrease rates, which are reported in Table 7, is about 1.3690 for and 1.1430 for .

The relative error profiles on the horizontal velocity, the stream function and the temperature along the -axis are plotted on Figure 20. For simplicity, only the results with the finest mesh, that is mesh (e) of Figure 18, are shown. The mean edge length in this mesh is . Figure 20 shows that the errors along the axis remains very small compared to the norm of the exact solutions. The overall errors are recorded in Table 8. As can be stated, the error on the velocity is about 2.316 percent of the norm of the exact solution, whereas the error on the stream function and on the temperature are less than 0.4 percent.

The meshes used for the previous results are unstructured and highly non-Delaunay. However, they are rather regular in the sense that the triangles have comparable areas and the edges have comparable lengths. In a last test, an analysis on non-Delaunay meshes with some very-badly-shaped triangles are carried out. These meshes are obtained from Delaunay meshes generated by Gmsh [45], of which some random vertices are moved, in order to obtain a prescribed ratio of non-Delaunay triangles. Contrarily to the distortion procedure in [29], the non-Delaunay triangles are not necessarily pairwise.

Three sets of four increasignly fine meshes are used. In each set, the meshes count respectively 40, 184, 676 and 2658 triangles. The four meshes in the first set comprise approximately 15 percent of non-Delaunay triangles. They are presented in Figure 21. In the second set, the ratio of non-Delaunay triangles is about 25 percent, whereas in the last set, 50 percent of the triangles are non-Delaunay. The meshes in the second set are plotted in Figure 22 and the third set can be seen in Figure 23. The mean edge lengths in the three sets range from to . It can be remarked in Figures 21 to 23 that the triangles have irregular shapes. Many of them have an angle close to 180 degrees and the area of the triangles varies very significantly. For instance, the minimum triangle area of mesh (d) in Figure 22 is whereas the maximum is , that is more than times as high.