## 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.

## 1. Introduction

The theory of Discrete Exterior Calculus (DEC) is a relatively recent discretization [8] of the classical theory of Exterior Differential Calculus developed by E. Cartan [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 [9], the Navier-Stokes equations [10], the simulation of elasticity, plasticity and failure of isotropic materials [5], some comparisons with the finite differences and finite volume methods on regular flat meshes [7], as well as applications in digital geometry processing [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 [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 [12] and Bossavit [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, 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, we deduce the discretization of the pullback operator on primal 1-forms. In Section 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, 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, we re-examine the geometry of some of the local DEC quantities. In Section 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, we summarize the contributions of this paper.

## 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 ${\textstyle v_{1},v_{2},v_{3}}$ (Figure 1).

 Figure 1. Triangle ${\textstyle v_{1},v_{2},v_{3}}$

Such a mesh has

• one oriented 2-dimensional face
 ${\displaystyle [v_{1},v_{2},v_{3}];}$
• three oriented 1-dimensional edges
 ${\displaystyle [v_{1},v_{2}],\quad [v_{1},v_{3}],\quad [v_{2},v_{3}];}$
• and three 0-dimensional vertices
 ${\displaystyle [v_{1}],\quad [v_{2}],\quad [v_{3}].}$

### 2.1 Boundary operator

There is a well known boundary operator

 ${\displaystyle \partial _{2,1}[v_{1},v_{2},v_{3}]=[v_{2},v_{3}]-[v_{1},v_{3}]+[v_{1},v_{2}],}$
(1)

which describes the boundary of the triangle as an alternated sum of its ordered oriented edges ${\textstyle [v_{1},v_{2}]}$, ${\textstyle [v_{1},v_{3}]}$ and ${\textstyle [v_{2},v_{3}]}$.

Similarly, one can compute the boundary of each edge

 ${\displaystyle \partial _{1,0}[v_{1},v_{2}]=[v_{2}]-[v_{1}],}$ ${\displaystyle \partial _{1,0}[v_{1},v_{3}]=[v_{3}]-[v_{1}],}$ (2) ${\displaystyle \partial _{1,0}[v_{2},v_{3}]=[v_{3}]-[v_{2}].}$

If we consider

• the symbol ${\textstyle [v_{1},v_{2},v_{3}]}$ as a basis vector of a 1-dimensional vector space,
• the symbols ${\textstyle [v_{1},v_{2}]}$, ${\textstyle [v_{1},v_{3}]}$, ${\textstyle [v_{2},v_{3}]}$ as an ordered basis of a 3-dimensional vector space,
• the symbols ${\textstyle [v_{1}]}$, ${\textstyle [v_{2}]}$, ${\textstyle [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

 ${\displaystyle [\partial _{2,1}]=\left({\begin{array}{r}1\\-1\\1\end{array}}\right),}$

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

 ${\displaystyle [\partial _{1,0}]=\left({\begin{array}{rrr}-1&-1&0\\1&0&-1\\0&1&1\end{array}}\right).}$

### 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 [6,8]. More precisely, suppose we have a function discretized by its values at the vertices

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

Its discrete derivative, according to DEC, is

 ${\displaystyle \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)}$ ${\displaystyle =\left({\begin{array}{l}f_{2}-f_{1}\\f_{3}-f_{1}\\f_{3}-f_{2}\end{array}}\right).}$

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

 ${\displaystyle f_{2}-f_{1}\approx df_{v_{1}}(v_{2}-v_{1}).}$

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

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

Let

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

### 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 ${\textstyle [v_{1},v_{2},v_{3}]}$ will correspond the 0-dimensional point given by the circumcenter ${\textstyle c}$ of the triangle (Figure 2).
 Figure 2. Circumcenter ${\textstyle c}$ of the triangle ${\textstyle [v_{1},v_{2},v_{3}]}$
• To the 1-dimensional edge ${\textstyle [v_{1},v_{2}]}$ will correspond the 1-dimensional straight line segment ${\textstyle [p_{1},c]}$ joining the midpoint ${\textstyle p_{1}}$ of the edge ${\textstyle [v_{1},v_{2}]}$ to the circumcenter ${\textstyle c}$ (Figure 3). Similarly for the other edges.
 Figure 3. Dual segment ${\textstyle [p_{1},c]}$ of the edge ${\textstyle [v_{1},v_{2}]}$
• To the 0-dimensional vertex/node ${\textstyle [v_{1}]}$ will correspond the oriented ${\textstyle 2}$-dimensional quadrilateral ${\textstyle [v_{1},p_{1},c,p_{3}]}$ (Figure 4).
 Figure 4. Dual quadrilateral ${\textstyle [v_{1},p_{1},c,p_{3}]}$ of the vertex ${\textstyle [v_{1}]}$

### 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 ${\textstyle M_{1}}$ applied to the discrete differential of a discretized function ${\textstyle f\sim (f_{1},f_{2},f_{3})}$ is given as follows:
• the value ${\textstyle f_{2}-f_{1}}$ assigned to the edge ${\textstyle [v_{1},v_{2}]}$ is changed to the new value
 ${\displaystyle {{length}[p_{1},c] \over {length}[v_{1},v_{2}]}(f_{2}-f_{1})}$

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

• the value ${\textstyle f_{3}-f_{1}}$ assigned to the edge ${\textstyle [v_{1},v_{3}]}$ is changed to the new value
 ${\displaystyle {{length}[p_{3},c] \over {length}[v_{1},v_{3}]}(f_{3}-f_{1})}$

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

• the value ${\textstyle f_{3}-f_{2}}$ assigned to the edge ${\textstyle [v_{2},v_{3}]}$ is changed to the new value
 ${\displaystyle {{length}[p_{2},c] \over {length}[v_{2},v_{3}]}(f_{3}-f_{2})}$

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

In other words,

 ${\displaystyle 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).}$
• Similarly, the discrete Hodge star map ${\textstyle M_{0}}$ on values on vertices is given as follows:
• the value ${\textstyle f_{1}}$ assigned to the vertex ${\textstyle [v_{1}]}$ is changed to the new value
 ${\displaystyle {Area}[v_{1},p_{1},c,p_{3}]f_{1}}$

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

• the value ${\textstyle f_{2}}$ assigned to the vertex ${\textstyle [v_{2}]}$ is changed to the new value
 ${\displaystyle {Area}[v_{2},p_{2},c,p_{1}]f_{2}}$

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

• the value ${\textstyle f_{3}}$ assigned to the vertex ${\textstyle [v_{3}]}$ is changed to the new value
 ${\displaystyle {Area}[v_{3},p_{3},c,p_{2}]f_{2}}$

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

In other words,

 ${\displaystyle 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).}$

### 2.5 Continuous 1-forms

A differential 1-form ${\textstyle \alpha }$ on ${\textstyle \mathbb {R} ^{2}}$ can be thought of as a function whose arguments are a point in ${\textstyle \mathbb {R} ^{2}}$ and a vector in ${\textstyle \mathbb {R} ^{2}}$, i.e.

 ${\displaystyle \alpha :\mathbb {R} ^{2}\times \mathbb {R} ^{2}\longrightarrow \mathbb {R} ,}$

and which is linear on the vector argument.

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

 ${\displaystyle (p,Y)\mapsto \alpha (p,Y):=X(p)\cdot Y,}$

where ${\textstyle Y\in \mathbb {R} ^{2}}$. On the other hand, from a continuos 1-form we can build a vector field as follows

 ${\displaystyle X(p):=\alpha (p,e_{1})e_{1}+\alpha (p,e_{2})e_{2}=(\alpha (p,e_{1}),\alpha (p,e_{2})),}$

where

 ${\displaystyle 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).}$

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

Given a tensor ${\textstyle S}$ (a linear transformation) and the vector field ${\textstyle X}$ associated to a 1-form ${\textstyle \alpha }$ as above, consider the vector field ${\textstyle S(X)}$. For any vector ${\textstyle Y}$, the value of the corresponding continuous 1-form is given by the product

 ${\displaystyle S(X(p))\cdot Y=X(p)\cdot S^{T}Y=\alpha (p,S^{T}Y)=(S^{T})^{*}\alpha (p,Y),}$

which is the value of the 1-form called the pullback of ${\displaystyle \alpha }$ by the tensor ${\displaystyle S^{T}}$.

#### 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 ${\textstyle v_{i}}$ to ${\textstyle v_{j}}$ we have

 ${\displaystyle \alpha _{ij}:=\int _{0}^{1}\alpha (v_{i}+t(v_{j}-v_{i}),v_{j}-v_{i})dt.}$

This is equivalent to calculating the line integral (circulation) of the corresponding vector field ${\textstyle X}$ along the path ${\textstyle \gamma _{ij}:=v_{i}+t(v_{j}-v_{i})}$, ${\textstyle 0\leq t\leq 1}$, parametrizing the edge ${\textstyle [v_{i},v_{j}]}$

 ${\displaystyle \alpha _{ij}:=\int _{\gamma _{ij}}X.}$

## 3. Pullback operator on primal 1-forms

In Section 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 [12], and Bossavit [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)

 ${\displaystyle N_{1}={1 \over 2A}[(y_{2}-y_{3})x+(x_{3}-x_{2})y+x_{2}y_{3}-x_{3}y_{2}],}$ ${\displaystyle N_{2}={1 \over 2A}[(y_{3}-y_{1})x+(x_{1}-x_{3})y+x_{3}y_{1}-x_{1}y_{3}],}$ (3) ${\displaystyle N_{3}={1 \over 2A}[(y_{1}-y_{2})x+(x_{2}-x_{1})y+x_{1}y_{2}-x_{2}y_{1}],}$

and their differentials

 ${\displaystyle dN_{1}={1 \over 2A}[(y_{2}-y_{3})dx+(x_{3}-x_{2})dy],}$ ${\displaystyle dN_{2}={1 \over 2A}[(y_{3}-y_{1})dx+(x_{1}-x_{3})dy],}$ (4) ${\displaystyle dN_{3}={1 \over 2A}[(y_{1}-y_{2})dx+(x_{2}-x_{1})dy].}$

The Whitney interpolation 1-forms are

 ${\displaystyle \phi _{12}=N_{1}dN_{2}-N_{2}dN_{1},}$ ${\displaystyle \phi _{13}=N_{1}dN_{3}-N_{3}dN_{1},}$ (5) ${\displaystyle \phi _{23}=N_{2}dN_{3}-N_{3}dN_{2}.}$

These differential 1-forms are such that, along the edge ${\textstyle [v_{k},v_{l}]}$

 ${\displaystyle \int _{v_{k}}^{v_{l}}\phi _{ij}=\delta _{ik}\delta _{jl}-\delta _{il}\delta _{jk},}$

i.e. ${\textstyle 0}$ or ${\textstyle \pm 1}$.

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

 ${\displaystyle \beta _{12},\quad \beta _{13},\quad \beta _{23}.}$

We form a continuous 1-form as follows

 ${\displaystyle \beta :=\beta _{12}\phi _{12}+\beta _{13}\phi _{13}+\beta _{23}\phi _{23}}$

which is such that, for ${\textstyle 1\leq i

 ${\displaystyle \int _{v_{i}}^{v_{j}}\beta =\beta _{ij}.}$

Now we will calculate

 ${\displaystyle \int _{v_{i}}^{v_{j}}S^{*}\beta }$

for an arbitrary linear transformation

 ${\displaystyle S:=\left({\begin{array}{cc}a&b\\c&d\end{array}}\right)}$

Let

 ${\displaystyle w_{1}=v_{2}-v_{1},}$ ${\displaystyle w_{2}=v_{3}-v_{1},}$ ${\displaystyle w_{3}=v_{3}-v_{2},}$

and

 ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|.}$

As it turns out,

 ${\displaystyle (S^{*}\beta )_{12}=\int _{v_{1}}^{v_{2}}S^{*}\beta }$ ${\displaystyle ={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],}$ ${\displaystyle (S^{*}\beta )_{13}=\int _{v_{1}}^{v_{3}}S^{*}\beta }$ ${\displaystyle ={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],}$ ${\displaystyle (S^{*}\beta )_{23}=\int _{v_{2}}^{v_{3}}S^{*}\beta }$ ${\displaystyle ={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],}$

where

 ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|,}$ ${\displaystyle \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|.}$

Anyhow, these quantities are areas of the parallelograms formed by one fixed edge of the original triangle and one edge transformed by ${\textstyle S}$.

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

 ${\displaystyle 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)}$
(6)

When the material is isotropic, i.e. ${\textstyle S}$ is a constant multiple of the ${\textstyle 2\times 2}$ identity matrix, we get a constant multiple of the ${\textstyle 3\times 3}$ 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 ${\textstyle f:\mathbb {R} ^{2}\longrightarrow \mathbb {R} }$ at a point ${\textstyle p\in \mathbb {R} ^{2}}$ in the direction of ${\textstyle w\in \mathbb {R} ^{2}}$ is defined by

 ${\displaystyle df_{p}(w):=\lim _{t\rightarrow 0}{f(p+tw)-f(p) \over t}=\nabla f(p)\cdot w.}$

Thus, we have three Vector Calculus identities

 ${\displaystyle df_{v_{1}}(v_{2}-v_{1})=\nabla f(v_{1})\cdot (v_{2}-v_{1}),}$ ${\displaystyle df_{v_{2}}(v_{3}-v_{2})=\nabla f(v_{2})\cdot (v_{3}-v_{2}),}$ ${\displaystyle df_{v_{3}}(v_{1}-v_{3})=\nabla f(v_{3})\cdot (v_{1}-v_{3}).}$

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

 ${\displaystyle df_{v_{1}}(v_{2}-v_{1})\approx f_{2}-f_{1},}$ ${\displaystyle df_{v_{2}}(v_{3}-v_{2})\approx f_{3}-f_{2},}$ ${\displaystyle df_{v_{3}}(v_{1}-v_{3})\approx f_{1}-f_{3}.}$

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

 ${\displaystyle W_{1}\cdot (v_{2}-v_{1})=f_{2}-f_{1}}$ (7) ${\displaystyle W_{1}\cdot (v_{3}-v_{1})=f_{3}-f_{1}.}$ (8)

If

 ${\displaystyle v_{1}=(x_{1},y_{1}),}$ ${\displaystyle v_{2}=(x_{2},y_{2}),}$ ${\displaystyle v_{3}=(x_{3},y_{3}),}$

then

 ${\displaystyle W_{1}=\left({f_{1}y_{2}-f_{1}y_{3}-f_{2}y_{1}+f_{2}y_{3}+f_{3}y_{1}-f_{3}y_{2} \over x_{1}y_{2}-x_{1}y_{3}-x_{2}y_{1}+x_{2}y_{3}+x_{3}y_{1}-x_{3}y_{2}},-{f_{1}x_{2}-f_{1}x_{3}-f_{2}x_{1}+f_{2}x_{3}+f_{3}x_{1}-f_{3}x_{2} \over x_{1}y_{2}-x_{1}y_{3}-x_{2}y_{1}+x_{2}y_{3}+x_{3}y_{1}-x_{3}y_{2}}\right)^{T}}$

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

 ${\displaystyle W_{2}\cdot (v_{1}-v_{2})=f_{1}-f_{2}}$ (9) ${\displaystyle W_{2}\cdot (v_{3}-v_{2})=f_{3}-f_{2}.}$

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

 ${\displaystyle f_{3}-f_{1}=W_{1}\cdot (v_{3}-v_{1})=W_{1}\cdot (v_{3}-v_{2}+v_{2}-v_{1})}$ ${\displaystyle =W_{1}\cdot (v_{3}-v_{2})+W_{1}\cdot (v_{2}-v_{1})=W_{1}\cdot (v_{3}-v_{2})+f_{2}-f_{1},}$

so that

 ${\displaystyle W_{1}\cdot (v_{3}-v_{2})=f_{3}-f_{2}.}$
(10)

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

 ${\displaystyle (W_{1}-W_{2})\cdot (v_{2}-v_{1})=0.}$
(11)

Subtracting (8) from (10) we get

 ${\displaystyle (W_{1}-W_{2})\cdot (v_{3}-v_{2})=0.}$
(12)

Since ${\textstyle v_{2}-v_{1}}$ and ${\textstyle v_{3}-v_{2}}$ are linearly independent and the two inner products in (11) and (12) vanish,

 ${\displaystyle W_{1}-W_{2}=0.}$

Analogously, the corresponding gradient vector ${\textstyle W_{3}}$ of ${\textstyle f}$ at the vertex ${\textstyle v_{3}}$ is equal to ${\textstyle W_{1}}$. This means that the three approximate gradient vectors at the three vertices coincide, i.e. we have a unique discrete flux vector ${\textstyle W=W_{1}=W_{2}=W_{3}}$ on the given element. Note that the discrete flux ${\textstyle W}$ satisfies

 ${\displaystyle W\cdot (v_{2}-v_{1})=f_{2}-f_{1},}$ ${\displaystyle W\cdot (v_{3}-v_{1})=f_{3}-f_{1},}$ ${\displaystyle W\cdot (v_{3}-v_{2})=f_{3}-f_{2}.}$

This means that the primal 1-form discretizing ${\textstyle df}$ can be obtained by performing the dot products of the discrete flux vector ${\textstyle W}$ 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 ${\textstyle f}$ in FEML is given by

 ${\displaystyle \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),}$

where ${\textstyle N_{1},N_{2}}$ and ${\textstyle N_{3}}$ are the basis functions given in (3). Explicitly

 ${\displaystyle \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)}$

where

 ${\displaystyle A={1 \over 2}[(x_{2}y_{3}-x_{3}y_{2})-(x_{1}y_{3}-x_{3}y_{1})+(x_{1}y_{2}-x_{2}y_{1})]}$

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

 ${\displaystyle \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},}$

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 ${\textstyle K}$ denote the symmetric anisotropy tensor

 ${\displaystyle K=\left({\begin{array}{cc}k_{11}&k_{12}\\k_{12}&k_{22}\end{array}}\right)}$

and recall the anisotropic Poisson equation

 ${\displaystyle -\nabla \cdot (K\,\nabla f)=q.}$

First recall that, in Exterior Differential Calculus, for any ${\textstyle w\in \mathbb {R} ^{2}}$,

 ${\displaystyle (K\nabla f(p))\cdot w=\nabla f(p)\cdot (K^{T}w)=\nabla f(p)\cdot (Kw)}$ ${\displaystyle =df_{p}(Kw)=(df_{p}\circ K)(w)=:(K^{*}df_{p})(w),}$

where ${\textstyle K^{*}df_{p}}$ is the pullback of ${\textstyle df_{p}}$ by ${\textstyle K}$.

Since we already have a discrete candidate ${\textstyle W}$ for ${\textstyle \nabla f}$ on the given element, the product of the matrix ${\textstyle K}$ and ${\textstyle W}$ gives us a candidate ${\textstyle W':=KW}$ for the anisptropic flux vector on such an element

 ${\displaystyle W':=KW=\left({\begin{array}{cc}k_{11}&k_{12}\\k_{12}&k_{22}\end{array}}\right)\left({\begin{array}{c}{[(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}\end{array}}\right)}$ ${\displaystyle =\left({\begin{array}{c}{k_{11}[f_{1}(y_{2}-y_{3})+f_{2}(y_{3}-y_{1})+f_{3}(y_{1}-y_{2})]+k_{12}[f_{1}(x_{3}-x_{2})+f_{2}(x_{1}-x_{3})+f_{3}(x_{2}-x_{1})] \over x_{1}y_{2}-x_{1}y_{3}-x_{2}y_{1}+x_{2}y_{3}+x_{3}y_{1}-x_{3}y_{2}}\\{k_{12}[f_{1}(y_{2}-y_{3})+f_{2}(y_{3}-y_{1})+f_{3}(y_{1}-y_{2})]+k_{22}[f_{1}(x_{3}-x_{2})+f_{2}(x_{1}-x_{3})+f_{3}(x_{2}-x_{1})] \over x_{1}y_{2}-x_{1}y_{3}-x_{2}y_{1}+x_{2}y_{3}+x_{3}y_{1}-x_{3}y_{2}}\end{array}}\right).}$

In order to produce a primal 1-form from ${\textstyle W'}$, all we have to do is assume that ${\textstyle K}$ is constant on the given triangle and take dot products with the edges of the triangle, i.e.

 ${\displaystyle \left({\begin{array}{c}W'\cdot w_{1}\\W'\cdot w_{2}\\W'\cdot w_{3}\end{array}}\right).}$

Nevertheless, in order to better understand the structure of the local system (as in FEML), it is desirable to have a factorization of the result as a matrix product

 ${\displaystyle \left({\begin{array}{c}W'\cdot w_{1}\\W'\cdot w_{2}\\W'\cdot w_{3}\end{array}}\right)=K^{DEC}D_{0}\left({\begin{array}{l}f_{1}\\f_{2}\\f_{3}\end{array}}\right)}$

where the matrix ${\textstyle K^{DEC}}$ is a real ${\textstyle 3\times 3}$ matrix depending on the geometry of the triangle and the matrix ${\textstyle K}$, and such that it is a multiple of the identity matrix when the domain is isotropic. In order to do this, we will use the discretization of the pullback opertator of Exterior Differential Calculus for arbitrary 1-forms from Section 3. Thus, using (6), ${\textstyle K^{DEC}}$ is the following ${\textstyle 3\times 3}$ matrix where we have set ${\textstyle S=K}$, i.e.

 ${\displaystyle K^{DEC}={1 \over 4A}\left({\begin{array}{ccc}\eta _{21}^{K}-\eta _{11}^{K}&\eta _{31}^{K}&\eta _{31}^{K}\\\eta _{22}^{K}&\eta _{32}^{K}-\eta _{12}^{K}&-\eta _{22}^{K}\\-\eta _{13}^{K}&-\eta _{13}^{K}&\eta _{33}^{K}-\eta _{23}^{K}\end{array}}\right)}$

where

 ${\displaystyle \eta _{kl}^{K}=\eta _{k}(K(w_{l})).}$

#### 4.2.1 Geometric interpretation of the entries of ${\displaystyle K^{DEC}}$4.2.1 Geometric interpretation of the entries of K D E C {\displaystyle K^{DEC}}

Consider the Figure 5, where ${\textstyle J}$ denotes the ${\textstyle 90^{\circ }}$ anti-clockwise rotation

 ${\displaystyle J=\left({\begin{array}{cc}0&-1\\1&0\end{array}}\right).}$
 Figure 5. Geometric interpretation of the entries of the anisotropy tensor discretization ${\displaystyle K^{DEC}}$

We have

 ${\displaystyle \eta _{21}^{K}=-{J(w_{2})\cdot K(w_{1}) \over 2A}=-{1 \over 2A}|J(w_{2})||K(w_{1})|\cos(\beta )}$ ${\displaystyle =-{1 \over 2A}|w_{2}||K(w_{1})|\cos(\alpha {+\pi }/2)={1 \over 2A}|w_{2}||K(w_{1})|\sin(\alpha )={A' \over A},}$

where ${\textstyle A'}$ is the area of the red triangle. Thus, the numbers ${\textstyle \eta _{kl}^{K}}$ are quotients of areas of transformed triangles divided by the area of the original triangle.

## 5. Anisotropic Poisson equation in 2D

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

### 5.1 Local DEC discretization of the 2D anisotropic Poisson equation

The anisotropic Poisson equation reads as follows

 ${\displaystyle -\nabla \cdot (K\,\nabla f)=q,}$

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

 ${\displaystyle -\star d\star (K^{*}df)=q}$

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

 ${\displaystyle -M_{0}^{-1}\,\left(-D_{0}^{T}\right)\,M_{1}\,K^{DEC}\,D_{0}\,[f]=[q],}$

or equivalently

 ${\displaystyle D_{0}^{T}\,M_{1}\,K^{DEC}\,D_{0}\,[f]=M_{0}\,[q].}$
(13)

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

 Figure 6. Triangle

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

 ${\displaystyle {1 \over 4A}\left({\begin{array}{ccc}1&-1&0\\-1&0&-1\\0&1&1\end{array}}\right)\left({\begin{array}{ccc}{l_{1} \over L_{1}}&0&0\\0&{l_{3} \over L_{3}}&0\\0&0&{l_{2} \over L_{2}}\end{array}}\right)\left({\begin{array}{ccc}\eta _{21}^{K}-\eta _{11}^{K}&\eta _{31}^{K}&\eta _{31}^{K}\\\eta _{22}^{K}&\eta _{32}^{K}-\eta _{12}^{K}&-\eta _{22}^{K}\\-\eta _{13}^{K}&-\eta _{13}^{K}&\eta _{33}^{K}-\eta _{23}^{K}\end{array}}\right)\left({\begin{array}{rrr}-1&1&0\\-1&0&1\\0&-1&1\end{array}}\right)\left({\begin{array}{c}f_{1}\\f_{2}\\f_{3}\end{array}}\right)=\left({\begin{array}{c}A_{1}q_{1}\\A_{2}q_{2}\\A_{3}q_{3}\end{array}}\right).}$

The diffusive term matrix is actually symmetric (see Subsection 5.3.1).

### 5.2 Local FEML-Discretized 2D anisotropic Poisson equation

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

 ${\displaystyle K_{e}=\int B^{t}DBdA,}$

where ${\textstyle D=K}$ is the matrix representing the anisotropic diffusion tensor, and the matrix ${\textstyle B}$ is given explicitly by

 ${\displaystyle B=\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).}$

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

 ${\displaystyle K_{e}=\int B^{T}DBdA=B^{T}DBA_{e}}$ ${\displaystyle ={1 \over 4A_{e}}\left({\begin{array}{ccc}y_{2}-y_{3}&x_{3}-x_{2}\\y_{3}-y_{1}&x_{1}-x_{3}\\y_{1}-y_{2}&x_{2}-x_{1}\end{array}}\right)\left({\begin{array}{cc}k_{11}&k_{12}\\k_{12}&k_{22}\end{array}}\right)\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)}$

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

 ${\displaystyle (K_{e})_{11}={1 \over 4A}(k_{11}(y_{2}-y_{3})^{2}+(k_{12}+k_{12})(y_{2}-y_{3})(x_{3}-x_{2})+k_{22}(x_{3}-x_{2})^{2})}$ ${\displaystyle ={1 \over 4A}({\begin{array}{cc}-(y_{3}-y_{2}),&x_{3}-x_{2}\end{array}})\left({\begin{array}{cc}k_{11}&k_{12}\\k_{12}&k_{22}\end{array}}\right)\left({\begin{array}{c}-(y_{3}-y_{2})\\x_{3}-x_{2}\end{array}}\right)}$ ${\displaystyle ={1 \over 4A}({\begin{array}{cc}x_{3}-x_{2},&y_{3}-y_{2}\end{array}})\left({\begin{array}{cc}0&1\\-1&0\end{array}}\right)\left({\begin{array}{cc}k_{11}&k_{12}\\k_{12}&k_{22}\end{array}}\right)\left({\begin{array}{cc}0&-1\\1&0\end{array}}\right)\left({\begin{array}{c}x_{3}-x_{2}\\y_{3}-y_{2}\end{array}}\right)}$ ${\displaystyle ={1 \over 4A}(J(v_{3}-v_{2}))^{T}KJ(v_{3}-v_{2})={1 \over 4A}J(v_{3}-v_{2})\cdot K(J(v_{3}-v_{2})),}$

where ${\textstyle J}$ is the ${\textstyle 90^{\circ }}$ anticlockwise rotation,

 ${\displaystyle J=\left({\begin{array}{cc}0&-1\\1&0\end{array}}\right).}$

In this notation, the diffusive term in local FEML is given as follows

 ${\displaystyle {1 \over 4A}\left({\begin{array}{ccc}J(v_{3}-v_{2})\cdot K(J(v_{3}-v_{2}))&J(v_{3}-v_{2})\cdot K(J(v_{1}-v_{3}))&J(v_{3}-v_{2})\cdot K(J(v_{2}-v_{1}))\\J(v_{1}-v_{3})\cdot K(J(v_{3}-v_{2}))&J(v_{1}-v_{3})\cdot K(J(v_{1}-v_{3}))&J(v_{1}-v_{3})\cdot K(J(v_{2}-v_{1}))\\J(v_{2}-v_{1})\cdot K(J(v_{3}-v_{2}))&J(v_{2}-v_{1})\cdot K(J(v_{1}-v_{3}))&J(v_{2}-v_{1})\cdot K(J(v_{2}-v_{1}))\end{array}}\right).}$

### 5.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 labeled in Figures 6 and 7.

 Figure 7. Circumscribed triangle.

We have the following identities:

 ${\displaystyle \pi =2(\alpha _{1}+\alpha _{2}+\alpha _{3}),}$ ${\displaystyle {2l_{i} \over L_{i}}=\tan(\alpha _{i}),}$ ${\displaystyle {l_{i} \over R}=\sin(\alpha _{i}),}$ ${\displaystyle {L_{i} \over 2R}=\cos(\alpha _{i}),}$ ${\displaystyle A_{1}={L_{1}l_{1} \over 4}+{L_{3}l_{3} \over 4},}$ ${\displaystyle A_{2}={L_{1}l_{1} \over 4}+{L_{2}l_{2} \over 4},}$ ${\displaystyle A_{3}={L_{2}l_{2} \over 4}+{L_{3}l_{3} \over 4}.}$

#### 5.3.1 The diffusive term

For instance, we claim that

 ${\displaystyle {J(v_{3}-v_{2})\cdot K(J(v_{3}-v_{2})) \over 4A}={1 \over 4A}\left({l_{1}(\eta _{3,1}^{K}+\eta _{2,1}^{K}-\eta _{1,1}^{K}) \over L_{1}}+{l_{3}(\eta _{3,2}^{K}+\eta _{2,2}^{K}-\eta _{1,2}^{K}) \over L_{3}}\right)}$

Indeed,

 ${\displaystyle \eta _{3,1}^{K}+\eta _{2,1}^{K}-\eta _{1,1}^{K}=\left|{\begin{array}{cc}x_{2}-x_{1}&y_{2}-y_{1}\\dx(K(w_{1}))&dy(K(w_{1}))\end{array}}\right|+\left|{\begin{array}{cc}x_{1}-x_{3}&y_{1}-y_{3}\\dx(K(w_{1}))&dy(K(w_{1}))\end{array}}\right|-\left|{\begin{array}{cc}x_{3}-x_{2}&y_{3}-y_{2}\\dx(K(w_{1}))&dy(K(w_{1}))\end{array}}\right|}$ ${\displaystyle =-2\left|{\begin{array}{cc}x_{3}-x_{2}&y_{3}-y_{2}\\dx(K(w_{1}))&dy(K(w_{1}))\end{array}}\right|=-2J(v_{3}-v_{2})\cdot K(v_{2}-v_{1})}$ ${\displaystyle \eta _{3,2}^{K}+\eta _{2,2}^{K}-\eta _{1,2}^{K}=\left|{\begin{array}{cc}x_{2}-x_{1}&y_{2}-y_{1}\\dx(K(w_{2}))&dy(K(w_{2}))\end{array}}\right|+\left|{\begin{array}{cc}x_{1}-x_{3}&y_{1}-y_{3}\\dx(K(w_{2}))&dy(K(w_{2}))\end{array}}\right|-\left|{\begin{array}{cc}x_{3}-x_{2}&y_{3}-y_{2}\\dx(K(w_{2}))&dy(K(w_{2}))\end{array}}\right|}$ ${\displaystyle =-2\left|{\begin{array}{cc}x_{3}-x_{2}&y_{3}-y_{2}\\dx(K(w_{2}))&dy(K(w_{2}))\end{array}}\right|=-2J(v_{3}-v_{2})\cdot K(v_{3}-v_{1})}$

Thus,

 $\begin{array}{ll}& \frac{1}{4A}\left(\frac{{l}_{1}\left({\eta }_{3,1}^{K}+{\eta }_{2,1}^{K}-{\eta }_{1,1}^{K}\right)}{{L}_{1}}+\frac{{l}_{3}\left({\eta }_{3,2}^{K}+{\eta }_{2,2}^{K}-{\eta }_{1,2}^{K}\right)}{{L}_{3}}\right)=\\ & \phantom{\rule{2em}{0ex}}=-\frac{1}{2A}\left(\frac{{l}_{1}}{{L}_{1}}J\left({v}_{3}-{v}_{2}\right)\cdot K\left({v}_{2}-{v}_{1}\right)+\frac{{l}_{1}}{{L}_{1}}J\left({v}_{3}-{v}_{2}\right)\cdot K\left({v}_{3}-{v}_{1}\right)\right)\\ & \phantom{\rule{2em}{0ex}}=\frac{J\left({v}_{3}-{v}_{2}\right)\cdot K\left({v}_{1}-{v}_{3}\right)}{2A}\frac{\mathrm{tan}\left({\alpha }_{3}\right)}{2}-\frac{}{}\end{array}$