Anisotropy in 2D Discrete Exterior Calculus

TL;DR

This paper introduces a local formulation for 2D Discrete Exterior Calculus that effectively handles material heterogeneity and anisotropic fluxes, demonstrating comparable or improved accuracy over FEML in solving anisotropic Poisson equations.

Contribution

It presents a novel local formulation of 2D DEC that incorporates anisotropy and heterogeneity, aligning its computational cost with FEML while improving solution accuracy.

Findings

DEC solutions converge numerically for anisotropic Poisson problems.

DEC solutions are comparable to FEML on fine meshes.

DEC performs slightly better than FEML on coarse meshes.

Abstract

We present a local formulation for 2D Discrete Exterior Calculus (DEC) similar to that of the Finite Element Method (FEM), which allows a natural treatment of material heterogeneity (element by element). It also allows us to deduce, in a robust manner, anisotropic fluxes and the DEC discretization of the pullback of 1-forms by the anisotropy tensor, i.e. we deduce how the anisotropy tensor acts on primal 1-forms. Due to the local formulation, the computational cost of DEC is similar to that of the Finite Element Method with Linear interpolations functions (FEML). The numerical DEC solutions to the anisotropic Poisson equation show numerical convergence, are very close to those of FEML on fine meshes and are slightly better than those of FEML on coarse meshes.

Figures

Equations

$$\partial [v_1,v_2,v_3]=[v_2,v_3]-[v_1,v_3]+[v_1,v_2],$$

$$\partial [v_1,v_2]$$

$$\partial [v_2,v_3]$$

$$\partial [v_3,v_1]$$

$$\left(\begin{array}[]{r}1\\ 1\\ 1\end{array}\right),$$

$$\left(\begin{array}[]{rrr}-1&0&1\\ 1&-1&0\\ 0&1&-1\end{array}\right).$$

$$f\sim\left(\begin{array}[]{l}f_{1}\\ f_{2}\\ f_{3}\end{array}\right).$$

$$\displaystyle\left(\begin{array}[]{rrr}-1&0&1\\ 1&-1&0\\ 0&1&-1\end{array}\right)^{T}\left(\begin{array}[]{l}f_{1}\\ f_{2}\\ f_{3}\end{array}\right)$$

$$f_2-f_1\approx d{f}_{v_1}(v_2-v_1).$$

$$D_{0}:=\left(\begin{array}[]{rrr}-1&1&0\\ 0&-1&1\\ 1&0&-1\end{array}\right).$$

$$\frac{\mathrm{length}[p_1,c]}{\mathrm{length}[v_1,v_2]}(f_2-f_1)$$

$$\frac{\mathrm{length}[p_2,c]}{\mathrm{length}[v_2,v_3]}(f_3-f_2)$$

$$\frac{\mathrm{length}[p_3,c]}{\mathrm{length}[v_3,v_1]}(f_1-f_3)$$

$$M_{1}=\left(\begin{array}[]{ccc}{{\rm length}[p_{1},c]\over{\rm length}[v_{1},v_{2}]}&0&0\\ 0&{{\rm length}[p_{2},c]\over{\rm length}[v_{2},v_{3}]}&0\\ 0&0&{{\rm length}[p_{3},c]\over{\rm length}[v_{3},v_{1}]}\end{array}\right).$$

$$\mathrm{Area}[v_1,p_1,c,p_3]f_1$$

$$\mathrm{Area}[v_2,p_2,c,p_1]f_2$$

$$\mathrm{Area}[v_3,p_3,c,p_2]f_2$$

$$M_{0}=\left(\begin{array}[]{ccc}{\rm Area}[v_{1},p_{1},c,p_{3}]&0&0\\ 0&{\rm Area}[v_{2},p_{2},c,p_{1}]&0\\ 0&0&{\rm Area}[v_{3},p_{3},c,p_{2}]\end{array}\right).$$

$$d{f}_p(w)$$

$$d{f}_{v_1}(v_2-v_1)$$

$$d{f}_{v_2}(v_3-v_2)$$

$$d{f}_{v_3}(v_1-v_3)$$

$$d{f}_{v_1}(v_2-v_1)$$

$$d{f}_{v_2}(v_3-v_2)$$

$$d{f}_{v_3}(v_1-v_3)$$

$$W_1\cdot (v_2-v_1)$$

$$W_1\cdot (v_3-v_1)$$

$$v_1$$

$$v_2$$

$$v_3$$

$$W_1$$

$$W_2\cdot (v_1-v_2)$$

$$W_2\cdot (v_3-v_2)$$

$$f_3-f_1$$

$$W_1\cdot (v_3-v_2)=f_3-f_2.$$

Full text

1 Introduction

The theory of Discrete Exterior Calculus (DEC) is a relatively recent discretization [7] of the classical theory of Exterior Differential Calculus, a theory developed by E. Cartan [2] which has been a fundamental tool in Differential Geometry and Topology for over a century. The aim of DEC is to solve partial differential equations preserving their geometrical and physical features as much as possible. There are only a few papers about implementions of DEC to solve certain PDEs, such as the Darcy flow and Poisson’s equation [8], the Navier-Stokes equations [9], the simulation of elasticity, plasticity and failure of isotropic materials [4], some comparisons with the finite differences and finite volume methods on regular flat meshes [6], as well as applications in digital geometry processing [3].

In this paper, we describe a local formulation of DEC which is reminiscent of that of the Finite Element Method (FEM) since, once the local systems of equations have been established, they can be assembled into a global linear system. This local formulation is also efficient and helpful in understanding various features of DEC that can otherwise remain unclear while dealing with an entire mesh. We will, therefore, take a local approach when recalling all the objects required by DEC [5]. Our main results are the following:

• •

We develop a local formulation of DEC analogous to that of FEM, which allows a natural treatment of heterogeneous material properties assigned to subdomains (element by element) and eliminates the need of dealing with it through ad hoc modifications of the global discrete Hodge star operator.

• •

Guided by the local formulation, we also deduce a natural way to approximate the flux/gradient-vector of a discretized function, as well as the anisotropic flux vector. We carry out a comparison of the formulas defining the flux in both DEC and Finite Element Method with linear interpolation functions (FEML).

• •

From the local formulation, we deduce the local DEC-discretization of the anisotropic Poisson equation. More precisely, in Exterior Differential Calculus the anisotropy tensor acts by pullback on the differential of the unknown function. Here, we deduce how the anisotropy tensor acts on primal 1-forms. We also carry out an algebraic comparison of the DEC and FEML local formulations of the anisotropic Poisson equation.

• •

We present three numerical examples of the approximate solutions to the stationary anisotropic Poisson equation on different domains using DEC and FEML. The numerical examples show numerical convergence and a competitive performance of DEC, as well as a computational cost similar to that of FEML. In fact, the numerical solutions with both methods on fine meshes are identical, and DEC shows a slightly better performance than FEML on coarse meshes.

The paper is organized as follows. In Section 2, we describe the local versions of the discrete derivative operator, the dual mesh and the discrete Hodge star operator. In Section 3, we deduce the natural way of computing flux vectors in DEC (which turns out to be equivalent to the FEML procedure), as well as the anisotropic flux vectors. In Section 4, we present the local DEC formulation of the 2D anisotropic Poisson equation and compare it with the local system of FEML, proving that the diffusion terms are identical while the source terms are discretized differently due to a different area-weight assignment for the nodes. In Section 5, we re-examine some of the local DEC quantities. In Section 6, we present and compare numerical examples of DEC and FEML approximate solutions to the 2D anisotropic Poisson equation on different domains with meshes of various resolutions. In Section 7 we summarize the contributions of this paper.

Acknowledgements. The second named author was partially supported by a CONACyT grant, and would like to thank the International Centre for Numerical Methods in Engineering (CIMNE) and the University of Swansea for their hospitality. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research.

2 Preliminaries on DEC from a local viewpoint

Let us consider a primal mesh made up of a single (positively oriented) triangle.

2.1 Boundary operator

There is a well known boundary operator

[TABLE]

which describes the boundary of the triangle as an alternated sum of its oriented edges $[v_{1},v_{2}]$, $[v_{2},v_{3}]$ and $[v_{3},v_{1}]$. Similarly, one can compute the boundary of each edge

[TABLE]

If we consider

• •

the symbol $[v_{1},v_{2},v_{3}]$ as a basis vector of a 1-dimensional vector space,

• •

the symbols $[v_{1},v_{2}]$, $[v_{2},v_{3}]$, $[v_{3},v_{1}]$ as an ordered basis of a 3-dimensional vector space,

• •

the symbols $[v_{1}]$, $[v_{2}]$, $[v_{3}]$ as an ordered basis of a 3-dimensional vector space,

then the map (1), which sends the oriented triangle to a sum of its oriented edges, is represented by the matrix

[TABLE]

while the map (2), which sends the oriented edges to sums of their oriented vertices, is represented by the matrix

[TABLE]

2.2 Discrete derivative

It has been argued that the DEC discretization of the differential of a function is given by the transpose of the matrix of the boundary operator on edges (see [7, 5]). More precisely, suppose we have a function discretized by its values at the vertices

[TABLE]

Its discrete derivative, according to DEC, is

[TABLE]

Indeed, such differences are rough approximations of the directional derivatives of $f$. For instance, $f_{2}-f_{1}$ is a rough approximation of the directional derivative of $f$ at $v_{1}$ in the direction of the vector $v_{2}-v_{1}$, i.e.

[TABLE]

It is precisely in this sense that, according to DEC,

• •

the value $f_{2}-f_{1}$ is assigned to the edge $[v_{1},v_{2}]$,

• •

the value $f_{3}-f_{2}$ is assigned to the edge $[v_{2},v_{3}]$,

• •

and the value $f_{1}-f_{3}$ is assigned to the edge $[v_{3},v_{1}]$.

Let

[TABLE]

2.3 Dual mesh

The dual mesh of the primal mesh consisting of a single triangle is constructured as follows:

• •

To the 2-dimensional triangular face $[v_{1},v_{2},v_{3}]$ will correspond the 0-dimensional point given by the circumcenter $c$ of the triangle.

• •

To the 1-dimensional edge $[v_{1},v_{2}]$ will correspond the 1-dimensional straight line segment $[p_{1},c]$ joining the midpoint $p_{1}$ of the edge $[v_{1},v_{2}]$ to the circumcenter $c$. Similarly for the other edges.

• •

To the 0-dimensional vertex/node $[v_{1}]$ will correspond the 2-dimensional quadrilateral $[v_{1},p_{1},c,p_{3}]$.

2.4 Discrete Hodge star

For the Poisson equation in 2D, we need two matrices: one relating original edges to dual edges, and another relating vertices to dual cells.

• •

The discrete Hodge star map $M_{1}$ applied to the discrete differential of a discretized function $f\sim(f_{1},f_{2},f_{3})$ is given as follows:

• –

the value $f_{2}-f_{1}$ assigned to the edge $[v_{1},v_{2}]$ is changed to the new value

[TABLE]

assigned to the segment $[p_{1},c]$;

• –

the value $f_{3}-f_{2}$ assigned to the edge $[v_{2},v_{3}]$ is changed to the new value

[TABLE]

assigned to the segment $[p_{2},c]$;

• –

the value $f_{1}-f_{3}$ assigned to the edge $[v_{3},v_{1}]$ is changed to the new value

[TABLE]

assigned to the edge $[p_{3},c]$.

In other words,

[TABLE]

• •

Similarly, the discrete Hodge star map $M_{0}$ on values on vertices is given as follows

• –

the value $f_{1}$ assigned to the vertex $[v_{1}]$ is changed to the new value

[TABLE]

assigned to the quadrilateral $[v_{1},p_{1},c,p_{3}]$;

• –

the value $f_{2}$ assigned to the vertex $[v_{2}]$ is changed to the new value

[TABLE]

assigned to the quadrilateral $[v_{2},p_{2},c,p_{1}]$;

• –

the value $f_{3}$ assigned to the vertex $[v_{3}]$ is changed to the new value

[TABLE]

assigned to the quadrilateral $[v_{3},p_{3},c,p_{2}]$.

In other words,

[TABLE]

3 Flux and anisotropy

In this section, we deduce the DEC formulae for the local flux, the local anisotropic flux and the local anisotropy operator for primal 1-forms.

3.1 The flux in local DEC

We wish to find a natural construction for the discrete flux (discrete gradient vector) of a discrete function. Recall from Vector Calculus that the directional derivative of a differentiable function $f:\mathbb{R}^{2}\longrightarrow\mathbb{R}$ at a point $p\in\mathbb{R}^{2}$ in the direction of $w\in\mathbb{R}^{2}$ is defined by

[TABLE]

Thus, we have three Vector Calculus identities

[TABLE]

As in subsection 2.2, the rough approximations to directional derivatives of a function $f$ in the directions of the (oriented) edges are given as follows

[TABLE]

Thus, if we want to find a discrete gradient vector $W_{1}$ of $f$ at the point $v_{1}$, we need to solve the equations of approximations

[TABLE]

If

[TABLE]

then

[TABLE]

Now, if we were to find a discrete gradient vector $W_{2}$ of $f$ at the point $v_{2}$, we need to solve the equations

[TABLE]

The vectors $W_{2}$ solving these equations is actually equal to $W_{1}$. Indeed, consider

[TABLE]

so that

[TABLE]

Thus, adding up (5) and (7) we get

[TABLE]

Subtracting (6) from (8) we get

[TABLE]

Since $v_{2}-v_{1}$ and $v_{3}-v_{2}$ are linearly independent and the two inner products in (9) and (10) vanish,

[TABLE]

Analogously, the corresponding gradient vector $W_{3}$ of $f$ at the vertex $v_{3}$ is equal to $W_{1}$. This means that the three approximate gradient vectors at the three vertices coincide. Let us call this unique vector $W$. Note that discrete flux $W$ satisfies

[TABLE]

This means that the primal 1-form discretizing $df$ can be obtained by the dot products of the discrete flux $W$ with the vectors of the triangle’s edges.

Remark. More generally, we can see that any vector which is constant on the triangle, naturally gives a primal 1-form on the edges of the triangle by means of its dot products with the triangle’s edge-vectors.

3.1.1 Comparison of DEC and FEML local fluxes

The local flux (gradient) of $f$ in FEML is given by

[TABLE]

where

[TABLE]

and

[TABLE]

is the area of the triangle. Explicitly

[TABLE]

so that the FEML flux is given by

[TABLE]

and we can see that its formula coincides with that of the DEC flux.

3.2 The anisotropic flux vector in local DEC

We will now discuss how to discretize anisotropy in 2D DEC. Let $K$ denote the symmetric anisotropy tensor

[TABLE]

and recall the anisotropic Poisson equation

[TABLE]

As in Subsection 3.1, we wish to find a vector $W^{\prime}$ which will play the role of a discrete version of the anisotropic flux vector $K\nabla f$.

First observe that, since $K$ is symmetric, for any $w\in\mathbb{R}^{2}$

[TABLE]

where $K^{*}df_{p}$ is called the pullback of $df_{p}$ by $K$. These identities mean that in order to discretize the anisotropic flux we need to understand the discretization of the linear functional $df_{p}\circ K$. Let us suppose that $K$* is constant on our triangle*. As before, we have three natural vectors on the triangle,

[TABLE]

Given the vector $Kw_{1}$, we have the option to write it down as a linear combination of two of the three aforementioned vectors. Since $w_{1}$ is being used already, we use the other two vectors, i.e.

[TABLE]

for some $\lambda_{1},\mu_{1}\in\mathbb{R}$. Similarly,

[TABLE]

for some $\lambda_{2},\mu_{2},\lambda_{3},\mu_{3}\in\mathbb{R}$. These equations can be solved for $\lambda_{1},\lambda_{2},\lambda_{3},\mu_{1},\mu_{2},\mu_{3}$. Now

[TABLE]

Similarly, for the vectors $w_{2}$ and $w_{3}$ we have the identities

[TABLE]

These equations lead to the three equations of approximations

[TABLE]

where $W^{\prime}$ is the vector that should approximate $K\nabla f(v_{3})$, $K\nabla f(v_{3})$ and $K\nabla f(v_{3})$. Thus, in order to find the discrete version $W^{\prime}$ of the anisotropic flux vector $K\nabla f$ on the triangle, we need to solve the system (12)

The system (12) has a unique solution. Indeed, since

[TABLE]

then

[TABLE]

i.e.

[TABLE]

Since $w_{1}$ and $w_{2}$ are linearly independent

[TABLE]

i.e.

[TABLE]

Thus, making the appropriate substitutions, we see that the third equation in (12) is dependent on the first two independent equations, and there is a unique vector $W^{\prime}$ that solves the system.

For the sake of completeness, the values of the parameters are:

[TABLE]

where

[TABLE]

is the $90^{\circ}$ counter-clockwise rotation, and

[TABLE]

which is, in fact, the image under $K$ of the discrete isotropic flux and shows the consistency of our local reasoning. Also observe that this formula is the same as that of the FEML anisotropic flux.

3.3 Anisotropy on primal 1-forms

The system (12) can be rewritten in matrix form as follows

[TABLE]

Recalling the Remark at the end of Subsection 3.1, the matrix identity (39) states that the primal 1-form dual to the anisotropic flux vector $W^{\prime}$ is given by the local DEC discretization of the anisotropy tensor

[TABLE]

acting on the primal 1-form $D_{0}[f]$.

Remark. The matrix $K^{DEC}$ is the local DEC discretization on primal 1-forms of the pullback operator $K^{*}$ on 1-forms. In this case, the discretization of $K^{*}df:=df\circ K$.

3.3.1 Geometric interpretation of the entries of $K^{DEC}$

Let us examine $\lambda_{1}$ in the anisotropic case. Consider the following figure

We have

[TABLE]

where $A^{\prime}$ is the area of the red triangle and we have used a well known formula for the area of a triangle in terms of an inner angle. Thus $\lambda_{1}$ is the negative of the quotient of the area $A^{\prime}$ of the red triangle and the area $A$ of the original triangle. The calculations for the other entries are similar.

3.3.2 Isotropic case

Now, let us assume $K=k\,{\rm Id}_{2\times 2}$ on the triangle. The previous calculations show that

[TABLE]

Note that, in this case,

[TABLE]

4 2D anisotropic Poisson equation

In this section, we describe the local DEC discretization of the 2D anisotropic Poisson equation and compare it to that of FEML.

4.1 Local DEC discretization of the 2D anisotropic Poisson equation

The anisotropic Poisson equation reads as follows

[TABLE]

where $f$ and $q$ are two functions on a certain domain in $\mathbb{R}^{2}$. In terms of the exterior derivative $d$ and the Hodge star operator $\star$ it reads as follows

[TABLE]

where $K^{*}df:=df\circ K$ and $K=K^{T}$. Following the discretization of the discretized divergence operator [5], the corresponding local DEC discretization of the anisotropic Poisson equation is

[TABLE]

or equivalently

[TABLE]

In order to simplify the notation, consider the lengths and areas defined in the Figure 6.

Now, the discretized equation (42) looks as follows:

[TABLE]

The diffusive term matrix

[TABLE]

is actually symmetric (see Subsection 4.3.1).

4.2 Local FEML-Discretized 2D anisotropic Poisson equation

The diffusive elemental matrix in FEM (frequently called “stiffness matrix”) on an element $e$ is given by

[TABLE]

where $D$ is the matrix representing the anisotropic diffusion tensor $K$ in this paper, and the matrix $B$ is given explicitly by

[TABLE]

Since the matrix $B$ is constant on an element of the mesh, the integral is easy to compute. Thus, the difussive matrix $K_{e}$ for a linear triangular element (FEML) is given by

[TABLE]

Now, let us consider the first diagonal entry of the local FEML anisotropic Poisson diffusive matrix $K_{e}$,

[TABLE]

where

[TABLE]

is the $90^{\circ}$ counter-clockwise rotation. In this notation, the diffusive term in local FEML is given as follows

[TABLE]

4.3 Comparison between local DEC and FEML discretizations

For the sake of brevity, we are only going to compare the entries of the first row and first column of each formulation. Consider the various lengths, areas and angles given in the triangle of Figure 7.

We have the following:

[TABLE]

4.3.1 The diffusive term

We claim that

[TABLE]

Indeed,

[TABLE]

Thus

[TABLE]

All we have to do is show that

[TABLE]

Note that, since $J(v_{3}-v_{2})$ is orthogonal to $v_{3}-v_{2}$, $J(v_{3}-v_{2})$ must be parallel to $c-{v_{2}+v_{3}\over 2}$. Thus,

[TABLE]

Now we are going to express $c$ in terms of $v_{1},v_{2},v_{3}$. Let us consider

[TABLE]

where $a,b$ are coefficients to be determined. Taking inner products with $(v_{2}-v_{1})$ and $(v_{3}-v_{1})$ we get the two equations

[TABLE]

Solving for $a$ and $b$

[TABLE]

Substituting all the relevant quantities in (49) we have, for instance, that the coefficient of $(v_{2}-v_{1})$ is

[TABLE]

and similarly for the coefficient of $(v_{1}-v_{3})$. The calculations for the remaining entries are similar.

Thus, the local DEC and FEML diffusive terms of the 2D anisotropic Poisson equation coincide.

4.3.2 The source term

As already observed in [5], the right hand sides of the local DEC and FEML systems are different

[TABLE]

While FEML uses a barycentric subdivision to calculate the areas associated to each node/vertex, DEC uses a circumcentric subdivision. Eventually, this leads the DEC discretization to a better approximation of the solution (on coarse meshes).

5 Some remarks about DEC quantities

5.1 The discrete Hodge star quantities revisited

The numbers appearing in the local DEC matrices can be expressed both in terms of determinants and in terms of trigonometric functions. More precisely,

[TABLE]

These expressions are valid regardless of the location of the circumcenter and can, indeed, take negative values. he angles that are measured in the scheme can be negative as in the obtuse triangle of Figure 8

and some quantities can be zero or negative. For instance, if

[TABLE]

hen

[TABLE]

5.2 Area weights assigned to vertices

In order to understand how local DEC assigns area weights to vertices differently from FEML, let us consider the obtuse triangle shown in Figure 8 Let $p_{1},p_{2},p_{3}$ be the middle points of the segments $[v_{2},v_{3}],[v_{2},v_{3}],[v_{3},v_{1}]$ respectively. As shown in Figure 9, the triangle $[v_{1},p_{3},c]$ lies completely outside of the triangle $[v_{1},v_{2},v_{3}]$. Geometrically, this implies that its area must be assigned a negative sign, which is confirmed by the determinant formulas of Subsection 5.1. On the other hand, the triangle $[v_{1},p_{1},c]$ will have positive area. Thus, their sum gives us the area $A_{1}$ in Figure 9.

The area $A_{3}$ is computed similarly, where the triangle $[p_{3},v_{3},c]$ is assigned negative area (see Figure 10).

Note that for $A_{2}$, the two triangles $[p_{1},v_{3},c]$ and $[v_{2},p_{2},c]$ both have positive areas (see Figure 11).

6 Numerical Examples

In this section, we present three examples in order to illustrate the performance of DEC resulting from the local formulation and its implementation. In all cases, we solve the anisotropic Poisson equation. The FEML methodology that we have used in the comparison can be consulted [10, 11, 1].

6.1 First example: Heterogeneity

This example is intended to highlight how Local DEC deals effectively with heterogeneous materials. Consider the region in the plane given in Figure 12.

• •

The difussion constant for the region labelled mat1 is $k=12$ and its source term is $q=20$.

• •

The difussion constant for the region labelled mat2 is $k=6$ and its source term is $q=5$.

The meshes used in this example are shown in Figure 13 and vary from coarse to very fine.

The numerical results for the maximum temperature value are exemplified in Table 1.

[TABLE]

The temperature and flux-magnitude distribution fields are shown in Figure 14.

Figure 15 shows the graphs of the temperature and the flux-magnitude along a horizontal line crossing the inner circle for the first two meshes.

6.2 Second example: Anisotropy

Let us solve the Poisson equation in a circle of radius one centered at the origin $(0,0)$ under the following conditions (see Figure 16):

• •

heat anisotropic diffusion constants $K_{x}=1.5,K_{y}=1.0$;

• •

material angle $30^{\circ}$;

• •

source term $q=1$;

• •

Dirichlet boundary condition $u=10$.

The meshes used in this example are shown in Figure 17 and vary from very coarse to very fine.

The numerical results for the maximum temperature value ($u(0,0)=10.2$) are exemplified in Table 2 where a comparison with the Finite Element Method with linear interpolation functions (FEML) is also shown.

[TABLE]

The temperature distribution and Flux magnitude fields for the finest mesh are shown in Figure 18.

Figures 19(a), 19(b) and 19(c) show the graphs of the temperature and flux magnitude values along a diameter of the circle for the different meshes of Figures 17(a), 17(b) and 17(c) respectively.

6.3 Third example: Heterogeneity and anisotropy

Let us solve the Poisson equation in a circle of radius on the following domain (see Figure 20) with various material properties. The geometry of the domain is defined by segments of ellipses passing through the given points which have centers at the origin $(0,0)$.

[TABLE]

• •

The Dirichlet boundary condition is $u=10$ and material properties (anisotropic heat diffusion constants, material angles and source terms) are given according to Figure 21 and the table below.

[TABLE]

The meshes used in this example are shown in Figure 22.

The numerical results for the maximum temperature value ($u(0,0)=10.2$) are exemplified in Table 3 where a comparison with the Finite Element Method with linear interpolation functions (FEML) is also shown.

[TABLE]

The temperature distribution and Flux magnitude fields for the finest mesh are shown in Figure 23.

Figure 24 shows the graphs of the temperature and flux magnitude values along a diameter of the circle for different meshes of Figure 22

Remark. As can be seen from the previous examples, DEC behaves well on coarse meshes. As expected, the results of DEC and FEML are similar for fine meshes. We would also like to point out the the computational costs of DEC and FEML are very similar.

7 Conclusions

DEC is a relatively recent discretization scheme for PDE’s which takes into account the geometric and analytic features of the operators and the domains involved. The main contributions of this paper are the following:

We have made explicit the local formulation of DEC, i.e. on each triangle of the mesh. As is customary, the local pieces can be assembled, which facilitates the implementation of DEC by the interested reader. Furthermore, the profiles of the assembled DEC matrices are equal to those of assembled FEML matrices. 2. 2.

Guided by the local formulation, we have deduced a natural way to approximate the flux/gradient vector of a discretized function as well as the anisotropic flux vector. We have shown that the formulas defining the flux in DEC and FEML coincide. 3. 3.

We have deduced how the anisotropy tensor acts on primal 1-forms. 4. 4.

We have deduced the local DEC formulation of the 2D anisotropic Poisson equation, and have proved that the DEC and FEML diffusion terms are identical, while the source terms are not – due to the different area-weight allocation for the nodes. 5. 5.

Local DEC allows a simple treatment of heterogeneous material properties assigned to subdomains (element by element), which eliminates the need of dealing with it through ad hoc modifications of the global discrete Hodge star operator matrix.

On the other hand we would like to point the following features:

• •

The area weights assigned to the nodes of the mesh when solving the 2D anisotropic Poisson equation can even be negative (when a triangle has an inner angle greater that $120^{\circ}$), in stark contrast to the FEML formulation.

• •

The computational cost of DEC is similar to that of FEML. While the numerical results of DEC and FEML on fine meshes are virtually identical, the DEC solutions are better than those of FEML on coarse meshes. Furthermore, DEC solutions display numerical convergence.

Our future work will include the DEC discretization of convective terms and DEC on 2-dimensional simplicial surfaces in 3D. Preliminary results on both problems are promising and competitive with FEML.