==Anisotropy in 2D Discrete Exterior Calculus==

​

''''''

-->

==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 (assigning material properties element by element). It also allows us to deduce, in a principled  manner, anisotropic fluxes and the DEC discretization of the pullback of 1-forms by the anisotropy tensor, i.e. we deduce the discrete action of the anisotropy tensor 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 interpolating 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.

<!--

<span id="fn-1"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) Centro de Investigación en Matemáticas, Calle Jalisco s/n,  Guanajuato, GTO 36023, México. Email: esqueda,rherrera,botello@cimat.mx </span>

​

<span id="fn-2"></span>

<span style="text-align: center; font-size: 75%;">([[#fnc-2|<sup>2</sup>]]) Departamento de  Matemáticas, Universidad de Guanajuato,  Guanajuato, GTO 36023, México. Email: valerocar@gmail.com</span>

-->

​

==1. Introduction==

​

The theory of Discrete Exterior Calculus (DEC) is a relatively recent discretization  <span id='citeF-8'></span>[[#cite-8|[8]]] of the classical theory of  Exterior Differential Calculus developed by  E. Cartan <span id='citeF-3'></span>[[#cite-3|[3]]], which is a fundamental tool in Differential Geometry and Topology.  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 <span id='citeF-9'></span>[[#cite-9|[9]]],  the Navier-Stokes equations <span id='citeF-10'></span>[[#cite-10|[10]]], the simulation of  elasticity, plasticity and failure of isotropic materials <span id='citeF-5'></span>[[#cite-5|[5]]], some comparisons with  the finite differences and finite volume methods on regular flat meshes <span id='citeF-7'></span>[[#cite-7|[7]]], as well as applications in digital geometry processing <span id='citeF-4'></span>[[#cite-4|[4]]].

​

In this paper, we describe a local formulation of DEC which is reminiscent of that of the Finite Element Method  (FEM). Indeed,  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 if one is dealing dealing with an entire mesh. Besides, we believe the local description to DEC will be accesible to a wide readership. We will, therefore, take a local approach when recalling all the objects required by 2D DEC <span id='citeF-6'></span>[[#cite-6|[6]]]. Our main results are the following:

​

* We present 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 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).

​

* In order to understand how to discretize the anisotropic Poisson equation, we develop the discretization of the pull-back operator of 1-forms under an arbitrary linear trasformation or tensor. The pull-back operator is one of the basic ingredients in Exterior Differential Calculus in the smooth setting and is essential in all of its applications  in Topology [bott-Tu].  We discretize the ''pullback operator on primal 1-forms induced by an arbitrary tensor'' using Whitney interpolation (the Whitney map). The Whitney interpolation forms were introduced by Hassler Whitney in 1957 (see <span id='citeF-12'></span>[[#cite-12|[12]]] p. 139) and  Bossavit <span id='citeF-1'></span>[[#cite-1|[1]]] explained their relevance in “mixed methods” of finite elements.  To our knowledge, this is the first time that the discrete version of this operator is presented in the DEC literature.   This has allowed us to discretize, in a principled manner, the anisotropic heat equation.

​

*  We carry out an analytic 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 DEC-solutions exhibit numerical convergence (see the error measurement tables) and a competitive performance, 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 Preliminaries on DEC from a local viewpoint|2]], we describe the local versions of the discrete derivative operator, the dual mesh, the discrete Hodge star operator and the meaning of a continuous 1-form on the plane.  In Section [[#3 Pullback operator on primal 1-forms|3]], we deduce the discretization of the pullback operator on primal 1-forms. In Section [[#4 Flux and anisotropy|4]], we deduce the natural way of computing flux vectors in DEC  (which turns out to be equivalent to the FEML result), as well as the anisotropic flux vectors.  In Section [[#5 Anisotropic Poisson equation in 2D|5]],  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 in both schemes, while the source terms are differentl due to a different area-weight assignment for the nodes.  In Section [[#6 Some remarks about DEC quantities|6]], we re-examine the geometry of some of the local DEC quantities. In Section [[#7 Numerical Examples|7]], 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 [[#8 Conclusions|8]], 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==

​

In this section, we will recall the basic operators of Dicrete Exterior Calculus  restricting ourselves to a mesh made up of one simplex/triangle. The local results we derive in the  paper can be assembled, just as in the Finite elemnt Method,  due to the additivity of both differentiation and integration.

​

Let us consider a primal mesh made up of a single (positively oriented) triangle with vertices <math display="inline">v_1,v_2,v_3</math>.

​

<div id='img-1'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"

|-

|style="padding:10px;"|[[Image:Esqueda_et_al_2020a-Triangle01.png|228px|Triangle ]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" style="padding-bottom:10px;"| '''Figure 1'''. Triangle 

|}

​

Such a mesh has

​

* one oriented 2-dimensional face

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>[v_1,v_2,v_3];</math>

|}

|}

* three oriented 1-dimensional edges

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>[v_1,v_2],\quad [v_1,v_3],\quad [v_2,v_3];</math>

|}

|}

* and three 0-dimensional vertices

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>[v_1],\quad [v_2],\quad [v_3].</math>

|}

|}

​

===2.1 Boundary operator===

​

There is a well known boundary operator

​

<span id="eq-1"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\partial _{2,1} [v_1,v_2,v_3]=[v_2,v_3]-[v_1,v_3]+[v_1,v_2], </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1)

|}

​

which describes the boundary of the triangle as an alternated sum of its ordered oriented edges <math display="inline">[v_1,v_2]</math>, <math display="inline">[v_1,v_3]</math> and <math display="inline">[v_2,v_3]</math>.

​

Similarly, one can compute the boundary of each edge

​

<span id="eq-2"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\partial _{1,0} [v_1,v_2] = [v_2]-[v_1], </math>

|-

| style="text-align: center;" | <math> \partial _{1,0} [v_1,v_3] = [v_3]-[v_1], </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|-

| style="text-align: center;" | <math> \partial _{1,0} [v_2,v_3] = [v_3]-[v_2].  </math>

|}

|}

​

If we consider

​

* the symbol <math display="inline">[v_1,v_2,v_3]</math> as a basis vector of a 1-dimensional vector space,

* the symbols <math display="inline">[v_1,v_2]</math>, <math display="inline">[v_1,v_3]</math>, <math display="inline">[v_2,v_3]</math> as an ordered basis of a 3-dimensional vector space,

* the symbols <math display="inline">[v_1]</math>, <math display="inline">[v_2]</math>, <math display="inline">[v_3]</math> as an ordered basis of a 3-dimensional vector space,

​

​

then the map ([[#eq-1|1]]), which sends the oriented triangle to a sum of its oriented edges, is represented by the matrix

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>[\partial _{2,1}]=\left(\begin{array}{r} 1 \\  -1 \\  1 \end{array}\right),</math>

|}

|}

​

while the map ([[#eq-2|2]]), which sends the oriented edges to sums of their oriented vertices, is represented by the matrix

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>[\partial _{1,0}]=\left(\begin{array}{rrr} -1 & -1 & 0 \\  1 & 0 & -1 \\  0 & 1 & 1 \end{array} \right).</math>

|}

|}

​

===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 <span id='citeF-8'></span><span id='citeF-6'></span>[[#cite-8|[8,6]]]). More precisely, suppose we have a function discretized by its values at the vertices

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>f\sim \left(\begin{array}{l} f_1\\ f_2\\ f_3               \end{array} \right).</math>

|}

|}

​

Its discrete derivative, according to DEC, is

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left(\begin{array}{rrr}-1 & -1 & 0 \\  1 & 0 & -1 \\  0 & 1 & 1 \end{array} \right)^T \left(\begin{array}{l}f_1\\ f_2\\ f_3               \end{array} \right) =\left(\begin{array}{rrr}-1 & 1 & 0 \\  -1 & 0 & 1 \\  0 & -1 & 1 \end{array} \right) \left(\begin{array}{l}f_1\\ f_2\\ f_3               \end{array} \right)</math>

|-

| style="text-align: center;" | <math> =\left(\begin{array}{l}f_2-f_1\\ f_3-f_1\\ f_3-f_2               \end{array} \right). </math>

|}

|}

​

Indeed, such differences are rough approximations of the directional derivatives of <math display="inline">f</math> along the oriented edges.  For instance, <math display="inline">f_2-f_1</math> is a rough approximation of the directional derivative of <math display="inline">f</math> at <math display="inline">v_1</math> in the direction of the vector <math display="inline">v_2-v_1</math>, i.e.

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>f_2-f_1 \approx df_{v_1}(v_2-v_1) .</math>

|}

|}

​

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

​

* the value <math display="inline">f_2-f_1</math> is assigned to the edge <math display="inline">[v_1,v_2]</math>,

* the value <math display="inline">f_3-f_1</math> is assigned to the edge <math display="inline">[v_1,v_3]</math>,

* and the value <math display="inline">f_3-f_2</math> is assigned to the edge <math display="inline">[v_2,v_3]</math>.

​

Let

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>D_0:=\left(\begin{array}{rrr} -1 & 1 & 0 \\  -1 & 0 & 1 \\  0 & -1 & 1 \end{array} \right). </math>

|}

|}

​

===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 <math display="inline">[v_1,v_2,v_3]</math> will correspond the 0-dimensional point given by the circumcenter <math display="inline">c</math> of the triangle.

​

<div id='img-2'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"

|-

|style="padding:10px;"|[[Image:Esqueda_et_al_2020a-Triangle02.png|240px|Circumcenter ]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" style="padding-bottom:10px;"| '''Figure 2'''. Circumcenter 

|}

​

* To the 1-dimensional edge <math display="inline">[v_1,v_2]</math> will correspond the 1-dimensional straight line segment <math display="inline">[p_1,c]</math> joining the midpoint <math display="inline">p_1</math> of the edge <math display="inline">[v_1,v_2]</math> to the circumcenter <math display="inline">c</math>.   Similarly for the other edges.

​

<div id='img-3'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;"

|-

|style="padding:10px;"|[[Image:Esqueda_et_al_2020a-Triangle03.png|240px|Dual segment ]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" style="padding-bottom:10px;"| '''Figure 3'''. Dual segment 

|}

​

* To the 0-dimensional vertex/node <math display="inline">[v_1]</math> will correspond the oriented <math display="inline">2</math>-dimensional quadrilateral <math display="inline">[v_1,p_1,c,p_3]</math>.

​

<div id='img-4'></div>

{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: bottom;"

|-

|style="padding:10px;"|[[Image:Esqueda_et_al_2020a-Triangle04.png|240px|Dual quadrilateral  ]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" style="padding-bottom:10px;"| '''Figure 4'''. Dual quadrilateral  

|}

​

===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 <math display="inline">M_1</math> applied to the discrete differential of a discretized function <math display="inline">f\sim (f_1,f_2,f_3)</math>  is given as follows:

​

* the value <math display="inline">f_2-f_1</math> assigned to the edge <math display="inline">[v_1,v_2]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{{length}[p_1,c]\over {length}[v_1,v_2]}(f_2-f_1)</math>

|}

|}

​

assigned to the segment <math display="inline">[p_1,c]</math>;

* the value <math display="inline">f_3-f_1</math> assigned to the edge <math display="inline">[v_1,v_3]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{{length}[p_3,c]\over {length}[v_1,v_3]}(f_3-f_1)</math>

|}

|}

​

assigned to the edge <math display="inline">[p_3,c]</math>.

* the value <math display="inline">f_3-f_2</math> assigned to the edge <math display="inline">[v_2,v_3]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{{length}[p_2,c]\over {length}[v_2,v_3]}(f_3-f_2)</math>

|}

|}

​

assigned to the segment <math display="inline">[p_2,c]</math>;

​

In other words,

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>M_1=\left(\begin{array}{ccc} \displaystyle{{length}[p_1,c]\over {length}[v_1,v_2]} & 0 & 0 \\  0 & \displaystyle{{length}[p_3,c]\over {length}[v_1,v_3]} & 0 \\  0 & 0 & \displaystyle{{length}[p_2,c]\over {length}[v_2,v_3]} \end{array} \right).</math>

|}

|}

​

* Similarly, the discrete Hodge star map <math display="inline">M_0</math> on values on vertices is given as follows:

​

* the value <math display="inline">f_1</math> assigned to the vertex <math display="inline">[v_1]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{Area}[v_1,p_1,c,p_3]f_1</math>

|}

|}

​

assigned to the quadrilateral <math display="inline">[v_1,p_1,c,p_3]</math>;

* the value <math display="inline">f_2</math> assigned to the vertex <math display="inline">[v_2]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{Area}[v_2,p_2,c,p_1]f_2</math>

|}

|}

​

assigned to the quadrilateral <math display="inline">[v_2,p_2,c,p_1]</math>;

* the value <math display="inline">f_3</math> assigned to the vertex <math display="inline">[v_3]</math> is changed to the new value

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{Area}[v_3,p_3,c,p_2]f_2</math>

|}

|}

​

assigned to the quadrilateral <math display="inline">[v_3,p_3,c,p_2]</math>.

​

In other words,

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>M_0=\left(\begin{array}{ccc} {Area}[v_1,p_1,c,p_3] & 0 & 0 \\  0 & {Area}[v_2,p_2,c,p_1] & 0 \\  0 & 0 & {Area}[v_3,p_3,c,p_2] \end{array} \right).</math>

|}

|}

​

===2.5 Continuous 1-forms===

​

A differential 1-form <math display="inline">\alpha </math> on <math display="inline">\mathbb{R}^2</math> can be thought of as a function whose arguments are  a point in <math display="inline">\mathbb{R} ^2</math> and a vector in <math display="inline">\mathbb{R} ^2</math>, i.e.

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\alpha : \mathbb{R}^2 \times \mathbb{R}^2 \longrightarrow \mathbb{R},</math>

|}

|}

​

and which is linear on the vector argument.

​

The dot product of <math display="inline">\mathbb{R}^2</math> allows us to see continuous 1-forms in a more familiar form as vector fields in the following way: If <math display="inline">X</math> is a vector field on <math display="inline">\mathbb{R}^2</math>, at each point <math display="inline">p\in \mathbb{R}^2</math> we can consider the  1-form

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>(p,Y)\mapsto   \alpha (p,Y):=X(p)\cdot Y, </math>

|}

|}

​

where <math display="inline">Y\in \mathbb{R}^2</math>. On the other hand, from a continuos 1-form we can build a vector field as follows

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>X(p) := \alpha (p,e_1) e_1 + \alpha (p,e_2) e_2 = (\alpha (p,e_1),\alpha (p,e_2)),</math>

|}

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> e_1=\left(\begin{array}{c} 1 \\  0 \end{array} \right) \quad  \mbox{and} \quad  e_2=\left(\begin{array}{c} 0 \\  1 \end{array} \right). </math>

|}

|}

​

====2.5.1 Pullback of a conitnuous 1-form====

​

Given a tensor <math display="inline">S</math> (a linear transformation) and the vector field <math display="inline">X</math> associated to a 1-form <math display="inline">\alpha </math> as above, consider the vector field <math display="inline">S(X)</math>. For any vector <math display="inline">Y</math>, the value of the corresponding continuous 1-form is given by the product

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>S(X(p))\cdot Y     =X(p)\cdot S^T Y </math>

|-

| style="text-align: center;" | <math>    =\alpha (p, S^T Y) </math>

|-

| style="text-align: center;" | <math>    =(S^T)^*\alpha (p, Y),  </math>

|}

|}

​

which is the value of the 1-form called ''the pullback of <math>\alpha </math> by the tensor <math>S^T</math>''.

​

====2.5.2 Discretization of continuous 1-forms====

​

In order to produce a primal 1-form from such a continuous 1-form we need to inegrate it along the edges of the mesh (the de Rham map). In our local description, this produces a collection of 3 numbers associated to the 3 oriented edges of the triangle, i.e. for the oriented edge joining <math display="inline">v_i</math> to <math display="inline">v_j</math> we have

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\alpha _{ij}:=\int _0^1 \alpha (v_i+t(v_j-v_i), v_j-v_i) dt.</math>

|}

|}

​

This is equivalent to calculating the line integral (circulation) of the corresponding vector field <math display="inline">X</math> along the path <math display="inline">\gamma _{ij}:=v_i+t(v_j-v_i)</math>, <math display="inline">0\leq t\leq 1</math>, parametrizing the edge <math display="inline">[v_i,v_j]</math>

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\alpha _{ij}:=\int _{\gamma _{ij}} X. </math>

|}

|}

​

==3. Pullback operator on primal 1-forms==

​

In Section [[#2 Preliminaries on DEC from a local viewpoint|2]], we dealt with the primal discretization of a continuous 1-form. Now we wish to understand the pullback of a primal 1-form. In order to do this, we will interpolate a primal 1-form using the Whitney interpolation forms to produce a  continuous 1-form (Whitney map) to which we can apply the pullback operator by a tensor (linear transformation), and then intergrate along the edges (de Rham map).

​

The Whitney interpolation forms were introduced by  Whitney in 1957 <span id='citeF-12'></span>[[#cite-12|[12]]], and Bossavit <span id='citeF-1'></span>[[#cite-1|[1]]] explained their relevance to “mixed methods” of finite elements.  For the sake of simplicity, we will only define the Whitney 0 and 1-forms.

​

First consider the Finite Element linear basis functions (or Whitney interpolation functions)

​

<span id="eq-3"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>N_1 = {1\over 2A}[(y_2-y_3)x+(x_3-x_2)y+x_2y_3-x_3y_2] ,</math>

|-

| style="text-align: center;" | <math> N_2 = {1\over 2A}[(y_3-y_1)x+(x_1-x_3)y+x_3y_1-x_1y_3] ,</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|-

| style="text-align: center;" | <math> N_3 = {1\over 2A}[(y_1-y_2)x+(x_2-x_1)y+x_1y_2-x_2y_1] , </math>

|}

|}

​

and their differentials

​

<span id="eq-4"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>dN_1 = {1\over 2A}[(y_2-y_3)dx+(x_3-x_2)dy] ,</math>

|-

| style="text-align: center;" | <math> dN_2 = {1\over 2A}[(y_3-y_1)dx+(x_1-x_3)dy] ,</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|-

| style="text-align: center;" | <math> dN_3 = {1\over 2A}[(y_1-y_2)dx+(x_2-x_1)dy] . </math>

|}

|}

​

The Whitney interpolation 1-forms are

​

<span id="eq-5"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\phi _{12}=N_1dN_2-N_2dN_1,</math>

|-

| style="text-align: center;" | <math> \phi _{13}=N_1dN_3-N_3dN_1 ,</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (5)

|-

| style="text-align: center;" | <math> \phi _{23}=N_2dN_3-N_3dN_2. </math>

|}

|}

​

These differential 1-forms are such that, along the edge <math display="inline">[v_k,v_l]</math>

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _{v_k}^{v_l}\phi _{ij} = \delta _{ik}\delta _{jl}-\delta _{il}\delta _{jk},</math>

|}

|}

​

i.e. <math display="inline">0</math> or <math display="inline">\pm 1</math>.

​

Now, suppose we have a primal 1-form given on our triangle by the collection of numbers

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\beta _{12},\quad \beta _{13},\quad \beta _{23}.</math>

|}

|}

​

We form a continuous 1-form as follows

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\beta :=\beta _{12}\phi _{12} + \beta _{13}\phi _{13} + \beta _{23}\phi _{23} </math>

|}

|}

​

which is such that, for <math display="inline">1\leq i <j \leq 3</math>

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _{v_i}^{v_j} \beta = \beta _{ij}. </math>

|}

|}

​

Now we will calculate

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _{v_i}^{v_j} S^*\beta </math>

|}

|}

​

for an arbitrary linear transformation

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>S:= \left( \begin{array}{cc} a & b \\  c & d \end{array}  \right)</math>

|}

|}

​

Let

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>w_1 = v_2-v_1,</math>

|-

| style="text-align: center;" | <math> w_2 = v_3-v_1,</math>

|-

| style="text-align: center;" | <math> w_3 = v_3-v_2, </math>

|}

|}

​

and

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\eta _1  =           (y_2 - y_3)dx +  (x_3 - x_2)dy \quad =\quad  2A\,\, dN_1  \quad      =\quad \left|\begin{array}{cc}x_3-x_2 & y_3-y_2 \\      dx & dy     \end{array} \right|, </math>

|-

| style="text-align: center;" | <math>  \eta _2 =           (y_3- y_1 )dx +  (x_1 -  x_3)dy \quad =\quad  2A\,\, dN_2  \quad      =\quad \left|\begin{array}{cc}x_1-x_3 & y_1-y_3 \\      dx & dy     \end{array} \right|, </math>

|-

| style="text-align: center;" | <math>  \eta _3 =           (y_1 -  y_2)dx +  (x_2-x_1 )dy \quad =\quad  2A\,\, dN_3  \quad      =\quad \left|\begin{array}{cc}x_2-x_1 & y_2-y_1 \\      dx & dy     \end{array} \right|. </math>

|}

|}

​

As it turns out,

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>(S^*\beta )_{12}     =\int _{v_1}^{v_2} S^*\beta </math>

|-

| style="text-align: center;" | <math>    = {1\over 4A} \left[(\eta _2(S(w_1))-\eta _1(S(w_1)))\beta _{12} + \eta _3(S(w_1))\beta _{13} + \eta _3(S(w_1))\beta _{23}\right],</math>

|-

| style="text-align: center;" | <math> (S^*\beta )_{13}     =\int _{v_1}^{v_3} S^*\beta </math>

|-

| style="text-align: center;" | <math>    = {1\over 4A} \left[\eta _2(S(w_2))\beta _{12} + (\eta _3(S(w_2))-\eta _1(S(w_2)))\beta _{13}  - \eta _2(S(w_2))\beta _{23}\right],</math>

|-

| style="text-align: center;" | <math> (S^*\beta )_{23}     =\int _{v_2}^{v_3} S^*\beta </math>

|-

| style="text-align: center;" | <math>    = {1\over 4A} \left[- \eta _1(S(w_3))\beta _{12} - \eta _1(S(w_3))\beta _{13} (\eta _3(S(w_3))-\eta _2(S(w_3)))\beta _{23} \right], </math>

|}

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\eta _1(S(w_1))    = \left|\begin{array}{cc}x_3-x_2 & y_3-y_2 \\      a(x_2-x_1)+b(y_2-y_1) & c(x_2-x_1)+d(y_2-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _1(S(w_2))    = \left|\begin{array}{cc}x_3-x_2 & y_3-y_2 \\      a(x_3-x_1)+b(y_3-y_1) & c(x_3-x_1)+d(y_3-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _1(S(w_3))    = \left|\begin{array}{cc}x_3-x_2 & y_3-y_2 \\      a(x_3-x_2)+b(y_3-y_2) & c(x_3-x_2)+d(y_3-y_2)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _2(S(w_1))    = \left|\begin{array}{cc}x_1-x_3 & y_1-y_3 \\      a(x_2-x_1)+b(y_2-y_1) & c(x_2-x_1)+d(y_2-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _2(S(w_2))    = \left|\begin{array}{cc}x_1-x_3 & y_1-y_3 \\      a(x_3-x_1)+b(y_3-y_1) & c(x_3-x_1)+d(y_3-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _2(S(w_3))    = \left|\begin{array}{cc}x_1-x_3 & y_1-y_3 \\      a(x_3-x_2)+b(y_3-y_2) & c(x_3-x_2)+d(y_3-y_2)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _3(S(w_1))    = \left|\begin{array}{cc}x_2-x_1 & y_2-y_1 \\      a(x_2-x_1)+b(y_2-y_1) & c(x_2-x_1)+d(y_2-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _3(S(w_2))    = \left|\begin{array}{cc}x_2-x_1 & y_2-y_1 \\      a(x_3-x_1)+b(y_3-y_1) & c(x_3-x_1)+d(y_3-y_1)     \end{array} \right|,</math>

|-

| style="text-align: center;" | <math> \eta _3(S(w_3))    = \left|\begin{array}{cc}x_2-x_1 & y_2-y_1 \\      a(x_3-x_2)+b(y_3-y_2) & c(x_3-x_2)+d(y_3-y_2)     \end{array} \right|. </math>

|}

|}

​

Anyhow, these quantities are areas of the parallelograms formed by one fixed edge of the original triangle and one edge transformed by <math display="inline">S</math>.

​

Thus, we have that the induced transformation on primal 1-forms is

​

<span id="eq-6"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> S^*\beta ={1\over 4A}\left(\begin{array}{ccc} \eta _2(S(w_1))-\eta _1(S(w_1)) & \eta _3(S(w_1)) & \eta _3(S(w_1)) \\  \eta _2(S(w_2)) & \eta _3(S(w_2))-\eta _1(S(w_2)) & -\eta _2(S(w_2)) \\  - \eta _1(S(w_3)) & - \eta _1(S(w_3)) & \eta _3(S(w_3))-\eta _2(S(w_3)) \end{array} \right) \left(\begin{array}{c} \beta _{12} \\  \beta _{13} \\  \beta _{23} \end{array} \right)   </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|}

​

When the material is isotropic, i.e. <math display="inline">S</math> is a constant multiple of the <math display="inline">2\times 2</math> identity matrix,  we get a constant multiple of the <math display="inline">3\times 3</math> identity matrix above as discrete pull-back operator.

​

==4. 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.

​

===4.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 <math display="inline">f:\mathbb{R}^2\longrightarrow \mathbb{R}</math> at a point <math display="inline">p\in \mathbb{R}^2</math> in the direction of <math display="inline">w\in \mathbb{R}^2</math> is defined by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>df_p(w):=\lim _{t\rightarrow 0}{f(p+tw)-f(p)\over t}=\nabla f(p)\cdot w.</math>

|-

<!--| style="text-align: center;" | <math> =\nabla f(p)\cdot w. </math>-->

|}

|}

​

Thus, we have three Vector Calculus identities

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>df_{v_1}(v_2-v_1) = \nabla f(v_1)\cdot (v_2-v_1)  ,</math>

|-

| style="text-align: center;" | <math>  df_{v_2}(v_3-v_2) = \nabla f(v_2)\cdot (v_3-v_2),</math>

|-

| style="text-align: center;" | <math>  df_{v_3}(v_1-v_3) = \nabla f(v_3)\cdot (v_1-v_3) . </math>

|}

|}

​

As in subsection [[#2.2 Discrete derivative|2.2]], the rough approximations to directional derivatives of a function <math display="inline">f</math> in the directions of the (oriented) edges are given as follows

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>df_{v_1}(v_2-v_1) \approx  f_2-f_1  ,</math>

|-

| style="text-align: center;" | <math>  df_{v_2}(v_3-v_2) \approx  f_3-f_2  ,</math>

|-

| style="text-align: center;" | <math>  df_{v_3}(v_1-v_3) \approx  f_1-f_3  . </math>

|}

|}

​

Thus, if we want to find a discrete gradient vector <math display="inline">W_1</math> of <math display="inline">f</math> at the point <math display="inline">v_1</math>,  we need to solve the equations

​

<span id="eq-7"></span>

<span id="eq-8"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_1\cdot (v_2-v_1) = f_2-f_1 </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|-

| style="text-align: center;" | <math> W_1\cdot (v_3-v_1) = f_3-f_1. </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (8)

|}

|}

​

If

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>v_1=(x_1,y_1),</math>

|-

| style="text-align: center;" | <math> v_2=(x_2,y_2),</math>

|-

| style="text-align: center;" | <math> v_3=(x_3,y_3), </math>

|}

|}

​

then

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_1 =\left({f_1y_2-f_1y_3-f_2y_1+f_2y_3+f_3y_1-f_3y_2\over x_1y_2-x_1y_3-x_2y_1+x_2y_3+x_3y_1-x_3y_2} ,  -{f_1x_2-f_1x_3-f_2x_1+f_2x_3+f_3x_1-f_3x_2\over x_1y_2-x_1y_3-x_2y_1+x_2y_3+x_3y_1-x_3y_2}\right)^T </math>

|}

|}

​

Now, if we were to find  a discrete gradient vector <math display="inline">W_2</math> of <math display="inline">f</math> at the point <math display="inline">v_2</math>,  we need to solve the equations

​

<span id="eq-9"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_2\cdot (v_1-v_2) = f_1-f_2</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|-

| style="text-align: center;" | <math> W_2\cdot (v_3-v_2) = f_3-f_2.</math>

|}

|}

​

The vectors <math display="inline">W_2</math> solving these equations is actually equal to <math display="inline">W_1</math>.  Indeed, consider

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>f_3-f_1 = W_1\cdot (v_3-v_1) </math>

|-

| style="text-align: center;" | <math> = W_1\cdot (v_3-v_2+v_2-v_1)</math>

|-

| style="text-align: center;" | <math> = W_1\cdot (v_3-v_2)+W_1\cdot (v_2-v_1)</math>

|-

| style="text-align: center;" | <math> = W_1\cdot (v_3-v_2)+f_2-f_1, </math>

|}

|}

​

so that

​

<span id="eq-10"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_1\cdot (v_3-v_2) = f_3-f_2. </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

Thus, adding up ([[#eq-7|7]]) and ([[#eq-9|9]]) we get

​

<span id="eq-11"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>(W_1-W_2)\cdot (v_2-v_1)=0. </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}

​

Subtracting  ([[#eq-8|8]]) from ([[#eq-10|10]]) we get

​

<span id="eq-12"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>(W_1-W_2)\cdot (v_3-v_2)=0. </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (12)

|}

​

Since <math display="inline">v_2-v_1</math> and <math display="inline">v_3-v_2</math> are linearly independent  and the two inner products in ([[#eq-11|11]]) and ([[#eq-12|12]]) vanish,

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_1-W_2=0.</math>

|}

|}

​

Analogously, the corresponding gradient vector <math display="inline">W_3</math> of <math display="inline">f</math> at the vertex <math display="inline">v_3</math> is equal to <math display="inline">W_1</math>. This means that the three approximate gradient vectors at the three vertices coincide, i.e. we have a unique discrete flux vector <math display="inline">W=W_1=W_2=W_3</math> on the given element.  Note that the discrete flux <math display="inline">W</math> satisfies

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W\cdot (v_2-v_1) = f_2-f_1,</math>

|-

| style="text-align: center;" | <math> W\cdot (v_3-v_1) = f_3-f_1,</math>

|-

| style="text-align: center;" | <math> W\cdot (v_3-v_2) = f_3-f_2  . </math>

|}

|}

​

This means that the primal 1-form discretizing <math display="inline">df</math> can be obtained by performing the dot products of the discrete flux vector <math display="inline">W</math>  with the vectors given by the triangle's oriented edges.

​

====4.1.1 Comparison of DEC and FEML local fluxes====

​

The local flux (gradient) of <math display="inline">f</math> in FEML is given by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> \left( \begin{array}{ccc} \displaystyle{\partial N_1\over \partial x} & \displaystyle\displaystyle{\partial N_2\over \partial x} & \displaystyle{\partial N_3\over \partial x} \\  \displaystyle{\partial N_1\over \partial y} & \displaystyle{\partial N_2\over \partial y} & \displaystyle{\partial N_3\over \partial y} \end{array}  \right) \left( \begin{array}{c} f_1 \\  f_2 \\  f_3 \end{array}  \right), </math>

|}

|}

​

where <math display="inline">N_1,N_2</math> and <math display="inline">N_3</math> are the basis functions given in ([[#eq-3|3]]). Explicitly

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left( \begin{array}{ccc}\displaystyle{\partial N_1\over \partial x} & \displaystyle{\partial N_2\over \partial x} & \displaystyle{\partial N_3\over \partial x} \\  \displaystyle{\partial N_1\over \partial y} & \displaystyle{\partial N_2\over \partial y} & \displaystyle{\partial N_3\over \partial y} \end{array}  \right) ={1\over 2A} \left( \begin{array}{ccc}y_2-y_3 & y_3-y_1 & y_1-y_2 \\  x_3-x_2 & x_1-x_3 & x_2-x_1 \end{array}  \right)  </math>

|}

|}

​

where

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>A = {1\over 2} [(x_2y_3-x_3y_2) -(x_1y_3-x_3y_1)+(x_1y_2-x_2y_1)] </math>

|}

|}

​

is the area of the triangle, so that the FEML flux is given by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> \left( {[(y_2-y_3)f_1+(y_3-y_1)f_2+(y_1-y_2)f_3]\over 2A} , {[(x_3-x_2)f_1+(x_1-x_3)f_2+(x_2-x_1)f_3]\over 2A}  \right)^T, </math>

|}

|}

​

and we can see that it coincides with the DEC flux.

​

===4.2 The anisotropic flux vector in local DEC===

​

We will now discuss how to discretize anisotropy in 2D DEC. Let <math display="inline">K</math> denote the symmetric anisotropy tensor

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>K=\left(\begin{array}{cc} k_{11} & k_{12} \\  k_{12} & k_{22} \end{array} \right)</math>

|}

|}

​

and recall the anisotropic Poisson equation

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>-\nabla \cdot (K\, \nabla f) = q.</math>

|}

|}

​

First recall that, in Exterior Differential Calculus, for any <math display="inline">w\in \mathbb{R}^2</math>,