# A Simple and Complete Discrete Exterior Calculus on General Polygonal Meshes

Lenka Ptáčková111Corresponding author, email<EMAIL_ADDRESS>Luiz
Velho

###### Abstract

Discrete exterior calculus (DEC) offers a coordinate–free discretization of
exterior calculus especially suited for computations on curved spaces. In this
work, we present an extended version of DEC on surface meshes formed by
general polygons that bypasses the need for combinatorial subdivision and does
not involve any dual mesh. At its core, our approach introduces a new
polygonal wedge product that is compatible with the discrete exterior
derivative in the sense that it satisfies the Leibniz product rule. Based on
the discrete wedge product, we then derive a novel primal–to–primal Hodge star
operator. Combining these three ‘basic operators’ we then define new discrete
versions of the contraction operator and Lie derivative, codifferential and
Laplace operator. We discuss the numerical convergence of each one of these
proposed operators and compare them to existing DEC methods. Finally, we show
simple applications of our operators on Helmholtz–Hodge decomposition,
Laplacian surface fairing, and Lie advection of functions and vector fields on
meshes formed by general polygons.

###### keywords:

Discrete exterior calculus , Geometry processing with polygonal meshes ,
Polygonal wedge product , Polygonal Hodge star operator , Discrete Laplace
operator , Discrete Lie advection

††journal: CAGD 2021

Figure 1: Comparison of implicit mean curvature flows on a general polygonal
mesh (29k vertices) after 10 iterations with time step $t=10^{-4}$. On the far
left is the original mesh. Our method (center left) and the algorithm of Alexa
and Wardetzky [2011] with their combinatorially enhanced Laplacian (center
right) produce visually well–smoothed meshes. However, their method with
purely geometric Laplacian (far right) exhibits some undesirable artifacts on
the ears, neck, and tail of the kitten.

##  1 Introduction

The discretization of differential operators on surfaces is fundamental for
geometry processing tasks, ranging from remeshing to vector fields
manipulation. Discrete exterior calculus (DEC) is arguably one of the
prevalent numerical frameworks to derive such discrete differential operators.
However, the vast majority of work on DEC is restricted to simplicial meshes,
and far less attention has been given to meshes formed by arbitrary polygons,
possibly non–planar and non–convex.

In this work, we propose a new discretization for several operators commonly
associated to DEC that operate directly on polygons without involving any
subdivision. Our approach offers three main practical benefits. First, by
working directly with polygonal meshes, we overcome the ambiguities of
subdividing a discrete surface into a triangle mesh. Second, our construction
operates solely on primal elements, thus removing any dependency on dual
meshes. Finally, our method includes the discretization of new differential
operators such as the contraction operator and Lie derivatives.

We concisely expose our framework, describe each of our operators and compare
them to existing DEC methods. We examine the accuracy of our numerical scheme
by a series of convergence tests on flat and curved surface meshes. We also
demonstrate the applicability of our method for Helmholtz–Hodge decomposition
of vector fields, surface fairing, and Lie advection of vector fields and
functions.

##  2 Related Work and Preliminaries

There is a vast literature on DEC on triangle meshes, e.g., [Hirani, 2003,
Desbrun et al., 2005, 2006, Crane et al., 2013] – all these publications have
in common that they deal with purely simplicial meshes and use a dual mesh to
define operators, we will refer to their approach as to the classical DEC.

As announced, unlike the classical DEC, our method works with general
polygonal meshes and does not involve any dual meshes. However, our operators
differ also in other aspects, e.g., support (see Figures 3–5). Next we briefly
introduce several basic DEC notions and point out the key differences between
our approach and existing schemes, principally the classical DEC.

###  2.1 Discrete differential forms and the exterior derivative

We strictly stick to the convention, common to previous DEC literature, that a
discrete $q$–form is located on $q$–dimensional cells of the given mesh.

Discrete differential forms are usually denoted by small Greek letters and
sometimes we add a number superscript to emphasize the degree of the form,
i.e., a $q$–form $\alpha$ can be denoted as $\alpha^{q}$.

A polygonal mesh $S$ is made of a set of vertices (0–dimensional cells),
oriented edges (1–dimensional cells), and oriented faces (2–dimensional
cells). A real discrete differential $q$–form $\alpha^{q}$ on $S$ is a
$q$–cochain, i.e., a real number assigned to each $q$–dimensional cell $c^{q}$
of $S$. E.g., if $(e_{0},e_{1},\dots,e_{n})$ is the vector of all edges of
$S$, then an 1–form $\alpha^{1}$ is a vector of real values

$\alpha^{1}=(\alpha(e_{0}),\dots,\alpha(e_{n})).$

The discrete exterior derivative $d$ is the coboundary operator and it holds:

$(d\alpha)(c^{q+1})=\alpha(\partial c^{q+1})=\sum_{c^{q}\in
S}[c^{q+1}:c^{q}]\alpha(c^{q}),$

where $\partial$ is the boundary operator and $[c^{q+1}:c^{q}]$ denotes the
incidence relation between cells $c^{q+1}$ and $c^{q}$, as depicted in the
example bellow.

$\displaystyle[f_{0}:e_{0}]$
$\displaystyle=1,[f_{0}:e_{1}]=1,\dots,[f_{0}:e_{4}]=-1$
$\displaystyle\partial f_{0}$ $\displaystyle=\sum\limits_{e_{i}\in
S}[f_{0}:e_{i}]e_{i}\;=\;e_{0}+e_{1}+e_{2}+e_{3}-e_{4}$
$\displaystyle\alpha^{1}$ $\displaystyle=(\alpha(e_{0}),\dots,\alpha(e_{n}))$
$\displaystyle d\alpha(f_{0})$ $\displaystyle=\sum\limits_{e_{i}\in
S}[f_{0}:e_{i}]\alpha(e_{i})\;=\;\alpha(e_{0})+\alpha(e_{1})+\alpha(e_{2})+\alpha(e_{3})-\alpha(e_{4})$

The boundary of a face $f_{0}$ is a sum of incident oriented edges, where we
take in account the orientation of the boundary edges with respect to the
given face. The discrete exterior derivative of a 1–form $\alpha^{1}$ (stored
on edges) is a 2–form $d\alpha$ located on faces and it is the “oriented sum”
of the values of $\alpha$ on boundary edges of $f_{0}$.

###  2.2 The cup product and the wedge product

We consider the wedge product on polygons to be the main building block of our
theory. On smooth manifolds, the wedge product allows for building higher
degree forms from lower degree ones. Similarly in algebraic topology of
pseudomanifolds, a cup product is a product of two cochains of arbitrary
degree $p$ and $q$ that returns a cochain of degree $p+q$ located on
$(p+q)$–dimensional cells. Thus we consider the cup product to be the
appropriate discrete version of the wedge product.

The cup product was introduced by J. W. Alexander, E. Čech, and H. Whitney in
1930’s and it became a well–studied notion in algebraic topology, mainly in
the simplicial setting. Later, the cup product was extended also to $n$–cubes
[Massey, 1991, Arnold, 2012].

In graphics, Gu and Yau [2003] presented a wedge product of two discrete
1–forms on triangulations, which turns to be equivalent to the cup product of
two 1–cochains on triangle complexes well studied in algebraic topology, e.g.,
[Whitney, 1957].

On triangles our discrete wedge product is equivalent to the cup product of
Whitney [1957], see also the Figure 2. On quadrilaterals our discrete wedge
product is equivalent to the cup product of Massey [1991].

Figure 2: The wedge product of two 1–forms on a triangle: the product of two
1–forms is a 2–form located on faces (far left).

In common to previous approaches, the discrete wedge product is
metric–independent and satisfies core properties such as the Leibniz product
rule, skew–commutativity, and associativity on closed forms.

###  2.3 The Hodge star operator

The most common discretization of the Hodge star operator on triangle meshes
is the so called diagonal approximation, see, e.g., [Desbrun et al., 2006],
which is computed based on the ratios between the volumes of primal simplices
and their dual cells. In contrast, we propose a Hodge star operator that does
not use a dual mesh. Since our dual forms are again located on primal
elements, we can compute the wedge product of primal and dual forms and hence
define further operators. However, this primal–primal definition brings some
drawbacks as well, we discuss them in Section 3.2.

###  2.4 The Hodge inner product

On smooth manifolds the Hodge star operator together with the wedge product
define the Hodge inner product, our definition is derived in the same fashion.

On the contrary, in classical DEC, e.g, [Desbrun et al., 2006], the Hodge star
is actually derived from a previously given inner product. Alexa and Wardetzky
[2011] in Lemma 3 also present a discrete version of inner product matrices,
that become the building blocks of their theory.

###  2.5 The codifferential

On a Riemannian $n$–manifold, the Hodge star operator is employed to define
the codifferential operator $\delta(\alpha^{k})=(-1)^{n(k-1)+1}\star
d\star\alpha.$ It is a linear operator that maps $k$–forms to $(k-1)$– forms.
On 1–forms it is also called the divergence operator.

In classical DEC the discrete codifferential operator on triangle meshes is
defined using the diagonal approximation of the Hodge star operator. Alexa and
Wardetzky [2011] in Section 3 hint at a codifferential of 1–forms on general
polygonal meshes. The main difference between these and our codifferentials is
in the support, see Figures 3 and 4.

Figure 3: Comparison of the support of the codifferential of 1–forms between
the classical DEC (L) and our method (R). The codifferential of a 1–form
$\alpha$ is a 0–form located on vertices. The value of $\delta\alpha$ on the
red vertex $v$ is a linear combination of values of $\alpha$ on edges colored
green. The edge thickness reflects the weight of the corresponding edge values
$\alpha$ on $\delta\alpha(v)$.  Figure 4: Comparison of the support of the
codifferential of 2–forms between the classical DEC (L) and our approach (R).
The codifferential of a 2–form $\beta$ is a 1–form $\delta\beta$ located on
edges. The value of $\delta\beta$ on the red edge $e$ is a linear combination
of the values of $\beta$ on faces colored green. The color intensity of faces
reflects their weight of influence on $\delta\beta(e)$.

###  2.6 The Laplace operator

In exterior calculus, the Laplace operator is given by $\Delta:=\delta
d+d\delta,$ where $\delta$ is the codifferential and $d$ the exterior
derivative. The Laplacian is defined in this way also in the classical DEC and
we follow this convention. The classical Laplacian on 0–forms is also called
the cotan Laplace operator. Even though many different approaches lead to the
cotan–formula, MacNeal [1949] was the first to derive it.

Discrete Laplacians of 0–forms on general polygonal meshes were introduced by
Alexa and Wardetzky [2011]. They derive Laplacians with a geometric and a
combinatorial part that improves their mesh processing methods, see also
Figure 1. By Theorem 2 therein, on triangle meshes their Laplacians reduce to
the cotan–formula. In Figure 5 we compare the support of their purely
geometric Laplacians to ours.

Figure 5: Comparison of support of the Laplacian of 0–forms between the
classical DEC (far left) and ours on triangle meshes (center left), Laplacian
of Alexa and Wardetzky [2011] for their $\lambda=0$ (center right) and ours on
polygonal meshes (far right). The Laplacian of a 0–form $\alpha$ is a 0–form
$\Delta\alpha$ located again on vertices. The value of $\Delta\alpha$ on the
red vertex $v$ is a linear combination of values of $\alpha$ on vertices
colored green. The support of our Laplacian is always larger, the point size
reflects the weight of respective $\alpha$s on $\Delta\alpha(v)$. We also
color yellow the faces whose vertices carry the $\alpha$s that enter as
variables for $\Delta\alpha(v)$.

###  2.7 The contraction operator and the Lie derivative

The Lie derivative can be thought of as an extension of a directional
derivative of a function to derivative of tensor fields (such as vector fields
or differential forms) along a vector field. It is invariant under coordinate
transformations, which makes it an appropriate version of a directional
derivative on curved manifolds. It evaluates the change of a tensor field
along the flow of a vector field and is widely used in mechanics.

Our discretization of Lie derivative of functions (0–forms) corresponds to the
functional map framework of Azencot et al. [2013], but now generalized to
polygonal meshes. Our discrete Lie derivatives are thus linear operators on
functions that produce new functions, with the property that the derivative of
a constant function is 0. Furthermore, both theirs and our Lie derivative of a
vector field produces a vector field. However, whereas their framework is
built for triangle meshes, we work with general polygonal meshes.

While maintaining the discrete exterior calculus framework, our work can also
be interpreted as an extension of the Lie derivative of 1–forms presented by
Mullen et al. [2011] from planar regular grids to surface polygonal meshes in
space.

##  3 Primal-to-Primal Operators

This section contains actual results of our research – we present the theory
and numerically evaluate the quality of our approximations by setting our
results against analytical solutions. We also compare our methods to other DEC
schemes.

Not using dual meshes simplifies the definition of several operators on
polygonal meshes, which may be a difficult task otherwise. Moreover, it helps
to maintain the compatibility of our operators since both the initial and the
mapped discrete forms are located on primal elements. However, this approach
also brings some drawbacks, we discuss them in this section.

###  3.1 The Discrete Wedge Product

Just like the wedge product of differential forms, our discrete wedge product
is a product of two discrete forms of arbitrary degree $k$ and $l$ that
returns a form of degree $k+l$ located on primal $(k+l)$–dimensional cells
(see Figure 6).

Figure 6: The wedge product on a quadrilateral: the product of two 0–forms is
a 0–form located on vertices (far left). The product of a 0–form with a 1–form
is a 1–form located on edges (center left). The product of a 0–form with a
2–form is a 2–form located on faces (center right), and the product of two
1–forms is a 2–form located on faces (far right).

On triangle meshes our discrete wedge product is identical to the cup product
given by Whitney [1957] and on quadrilaterals it is equivalent to the cubical
cup product of Arnold [2012]. Further, the wedge product of differential forms
satisfies the Leibniz product rule with exterior derivative and is
skew–commutative. The discrete wedge product must satisfy these properties as
well, we thus appropriately extend the discrete wedge product from triangles
and quads to general polygons and derive the following formulas:

###### Definition 3.1

Let $S$ be a surface mesh (pseudomanifold) whose faces (2–cells) are polygons.
The polygonal wedge product of two discrete forms $\alpha^{k},\beta^{l}$ is a
$(k+l)$–form $\alpha^{k}\wedge\beta^{l}$ defined on each $(k+l)$–cell
$c^{k+l}\in S$. Let $v$ be a vertex, $e=(v_{i},v_{j})$ an edge, and
$f=(v_{0},\dots,v_{p-1})$ a $p$–polygonal face with boundary edges
$e_{0},\dots,e_{p-1}$ having the same orientation as $f$. The polygonal wedge
product is given for each degree by:

$\displaystyle(\alpha^{0}\wedge\beta^{0})(v)$ $\displaystyle=$
$\displaystyle\alpha(v)\beta(v),$
$\displaystyle(\alpha^{0}\wedge\beta^{1})(e)$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\alpha(v_{i})+\alpha(v_{j}))\beta(e),$
$\displaystyle(\alpha^{0}\wedge\beta^{2})(f)$ $\displaystyle=$
$\displaystyle\frac{1}{p}\bigg{(}\sum_{i=0}^{p-1}\alpha(v_{i})\bigg{)}\beta(f),$
$\displaystyle(\alpha^{1}\wedge\beta^{1})(f)$ $\displaystyle=$
$\displaystyle\sum_{a=1}^{\lfloor\frac{p-1}{2}\rfloor}\bigg{(}\frac{1}{2}-\frac{a}{p}\bigg{)}\sum_{i=0}^{p-1}\alpha(e_{i})(\beta(e_{i+a})-\beta(e_{i-a})),\text{where
all indices are modulo }p.$

The polygonal wedge product is illustrated in Figures 2 and 6. It is a
skew–commutative bilinear operation:
$\alpha^{k}\wedge\beta^{l}=(-1)^{kl}\beta^{l}\wedge\alpha^{k},$ matching its
continuous analog, and it satisfies the Leibniz product rule with discrete
exterior derivative:
$d(\alpha^{k}\wedge\beta^{l})=d\alpha\wedge\beta+(-1)^{k}\alpha\wedge d\beta.$

The wedge product of three 0–forms is trivially associative (it reduces to
multiplication of three scalars). Unfortunately for higher degree forms it is
not associative in general, only if one of the 0–forms involved is constant.
This is a common drawback of discrete wedge products, see e.g. [Hirani, 2003,
Remark 7.1.4.].

In matrix form, our polygonal wedge product reads:

$\displaystyle\alpha^{0}\wedge\epsilon^{0}$
$\displaystyle=\alpha^{0}\odot\epsilon^{0},$
$\displaystyle\alpha^{0}\wedge\beta^{1}$
$\displaystyle=(\operatorname{\mathsf{B}}\alpha^{0})\odot\beta^{1},$
$\displaystyle\alpha^{0}\wedge\omega^{2}$
$\displaystyle=(\operatorname{\mathsf{fv}}\alpha^{0})\odot\omega^{2},$
$\displaystyle(\beta^{1}\wedge\gamma^{1})|_{f}$
$\displaystyle=(\beta^{1}|_{f})^{\top}\operatorname{\mathsf{R}}(\gamma^{1}|_{f}),$

where $\odot$ is the Hadamard (element–wise) product, $\beta|_{f}$ denotes the
restriction of $\beta$ to a $p$–polygonal face $f$, and the matrices
$\operatorname{\mathsf{B}}\in\mathbb{R}^{|E|\times|V|}$,
$\operatorname{\mathsf{fv}}\in\mathbb{R}^{|F|\times|V|}$, and
$\operatorname{\mathsf{R}}\in\mathbb{R}^{p\times p}$ per $f$ read:

$\displaystyle\operatorname{\mathsf{B}}[i,j]=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{1}{2}&\mbox{if }v_{j}\prec
e_{i},\\\ 0&\mbox{otherwise.}\end{array}\right.$ (3)
$\displaystyle\operatorname{\mathsf{fv}}[i,j]=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{1}{p_{i}}&\mbox{if }v_{j}\prec
f_{i},f_{i}\text{ is a }p_{i}\text{--gon},\\\
0&\mbox{otherwise.}\end{array}\right.$ (6)
$\displaystyle\operatorname{\mathsf{R}}=$
$\displaystyle\sum_{a=1}^{\lfloor\frac{p-1}{2}\rfloor}\bigg{(}\frac{1}{2}-\frac{a}{p}\bigg{)}\operatorname{\mathsf{R}}_{a},\quad\operatorname{\mathsf{R}}_{a}[k,j]=\left\\{\begin{array}[]{ll}1&\mbox{if
}e_{j}\text{ is $(k+a)$--th halfedge of }f,[f:e_{j}]=1,\\\ -1&\mbox{if
}e_{j}\text{ is $(k-a)$--th halfedge of }f,[f:e_{j}]=1,\\\ 0&\mbox{otherwise.
}\end{array}\right.$ (10)

####  3.1.1 Numerical evaluation

We perform the numerical evaluation of our polygonal wedge product as an
approximation to the continuous wedge product on a given mesh $S$ over a
smooth surface in the following fashion:

##### 1

We integrate each differential $l$–form over all $l$–dimensional cells of the
mesh $S$ and thus define discrete forms $\alpha^{0}$, $\beta^{1}$,
$\gamma^{1}$, and $\omega^{2}$:

$\alpha^{0}(v)=A(v),\;\;\beta^{1}(e)=\int\limits_{e}B,\quad\gamma^{1}(e)=\int\limits_{e}\Gamma,\quad\omega^{2}(f)=\int\limits_{f}\Omega,$

where Greek capital letters denote the respective continuous differential
forms. In practice, we integrate the continuous differential 2–form $\Omega$
over a set of triangles $(C,v_{i},v_{i+1})$ that approximate the possibly
non–planar face $f=(v_{0},\dots,v_{p-1})$, where $C$ is the centroid of $f$,
i.e.,

$\omega^{2}(f)=\sum_{i=0}^{p-1}\int\limits_{(C,v_{i},v_{i+1})}\Omega,\;\text{
where }C=\frac{1}{p}\sum_{i=0}^{p-1}v_{i}.$

##### 2

Next we compute the polygonal wedge products $(\alpha^{0}\wedge\beta^{1})(e)$,
$(\alpha^{0}\wedge\omega^{2})(f)$, $(\beta^{1}\wedge\gamma^{1})(f)$ for each
edge and face of the mesh using our formulas.

##### 3

We also calculate analytical solutions of the (continuous) wedge products and
discretize (integrate) these solutions.

##### 4

We then compute the $L^{\infty}$ and $L^{2}$ errors of our approximation. So
let $\xi^{k}$ denote our solution (a discrete $k$–form) and $\Xi^{k}$ the
respective discretized analytical solution, we compute:

$L^{2}\text{ error
}=\Big{(}\xi^{k}-\Xi^{k}\Big{)}^{\top}M_{k}\Big{(}\xi^{k}-\Xi^{k}),\quad\quad
L^{\infty}\text{ error
}=||\xi^{k}-\Xi^{k}||_{\infty}=\max_{c^{k}}(|\xi^{k}(c^{k})-\Xi^{k}(c^{k})|),$

where $M_{k}$ are inner product matrices. Concretely,
$M_{2}\in\mathbb{R}^{|F|^{2}}$, $M_{0}\in\mathbb{R}^{|V|^{2}}$ are diagonal
matrices given by

$M_{2}[i,i]=\frac{1}{|f_{i}|},\quad M_{0}[i,i]=\sum_{f_{j}\succ
v_{i}}\frac{|f_{j}|}{p_{j}},$ (11)

and $M_{1}$ is the inner product of two 1–forms of Alexa and Wardetzky [2011],
i.e., for two 1–forms $\epsilon$ and $\lambda$, $M_{1}$ is defined in the
sense that

$\epsilon^{\top}M_{1}\lambda=\sum_{f}\epsilon|_{f}^{\top}M_{f}\lambda|_{f},\quad
M_{f}:=\frac{1}{|f|}B_{f}B_{f}^{\top},$ (12)

where $\epsilon|_{f}$ again denotes the restriction of $\epsilon$ to a
$p$–polygonal face $f$ and $B_{f}$ denotes a $p\times 3$ matrix with edge
midpoint positions as rows (we take the centroid of each face as the center of
coordinates per face).

##### 5

To evaluate the numerical convergence behavior, we refine the mesh over the
given smooth underlying surface. The smooth surfaces used for tests are: unit
sphere, torus azimuthally symmetric about the z-axis, and planar square. To
create unstructured meshes, we randomly eliminate a given percentage of edges
of an initially regular mesh.

We also use jittering to evaluate the influence of irregularity of a mesh on
the experimental convergence. When jittering, we start with a regular mesh and
displace each vertex in a random tangent direction to distance $r\cdot|e|$,
where $|e|$ is the shortest edge length, and then project all thus displaced
vertices on a given underlying smooth surface. That is why we use simple
surfaces such as spheres and tori for our tests.

If not stated otherwise, all graphs use $\log_{10}$ scales on both the
horizontal and vertical axes.

Figure 7: Convergence of the wedge products on a set of unstructured polygonal
meshes on a planar square to analytical solutions in $L^{2}$ norm (far left)
and $L^{\infty}$ norm (center left). Both axes are in $\log_{10}$ scale. The
differential forms tested are trigonometric forms
$\alpha^{0}=\sin(x)\cos(y)+1$,
$\beta^{1}=(\sin^{2}(x)-1)dx+(3\cos(x+2)+\sin(y))dy$,
$\gamma^{1}=(\cos(x)\sin(y)+3)dx+\cos(y)dy$,
$\omega^{2}=(\sin(xy)+\cos(1))dx\wedge dy$. On the right are samples of tested
meshes, both over a planar $[-1,1]^{2}$ square.

We have tested quadratic and trigonometric differential forms on flat and
curved surface meshes (with non–planar faces) and our polygonal wedge products
exhibit at least linear convergence to the respective analytical solutions,
both in $L^{2}$ and $L^{\infty}$ norm. In Figure 7 we give an example.

###  3.2 The Hodge Star Operator

We define a discrete Hodge star operator as a homomorphism (linear operator)
from the group of $k$–forms to $(2-k)$–forms. But since we do not employ any
dual mesh and there is no isomorphism between the groups of $k$– and
$(2-k)$–dimensional cells, our Hodge star is not an isomorphism (invertible
operator), unlike its continuous counterpart and diagonal approximations.

On the other hand, thanks to the dual forms being located on elements of our
primal mesh, we can compute discrete wedge products of primal and dual forms
and thus define a discrete inner product and discrete contraction operator
later on.

Moreover, thanks to the Hodge star operating on primal meshes, we circumvent
the ambiguity of defining dual meshes of unstructured general polygonal
meshes. The idea of defining a Hodge star operator without using a dual mesh
was borrowed from Arnold [2012], where the author suggests metric–independent
Hodge star operators on simplicial and cubical complexes.

Our formulas for discrete Hodge star operators $\star$ are further motivated
by the condition that the Hodge dual of constant discrete forms on planar
surfaces is exact, hence $\star\mu=1$ and $\star 1=\mu$, where $\mu$ is the
volume form on a given Riemannian manifold, for details see [Abraham et al.,
1988, Section 6.5].

Figure 8: On the left, the Hodge dual of a 2–form $\omega$ is a 0–form
$\star\omega$, which value on a vertex $v$ (colored red) is a linear
combination of values of $\omega$ on adjacent faces (colored green). On the
right, the Hodge dual of a 0–form $\alpha$ is a 2–form $\star\alpha$, the
value of $\star\alpha$ on a face $f$ (colored red) is a linear combination of
values of $\alpha$ on vertices (green) of that face.

The Hodge star operator on 2–forms takes in account the degree $p_{i}$ of
$p_{i}$–polygonal faces $f_{i}$ and their vector areas $|f_{i}|$. If
$\omega^{2}$ is a 2–form, then the 0–form $\star\omega$ on a vertex $v$ is
given by

$(\star_{2}\omega)(v)=\frac{1}{\displaystyle\sum_{f_{i}\succ
v}\frac{|f_{i}|}{p_{i}}}\cdot\sum_{f_{i}\succ v}\frac{\omega(f_{i})}{p_{i}},$
(13)

i.e., it is a linear combination of values of $\omega$ on faces adjacent to
$v$, see Figure 8 left.

The Hodge star on an 1–form $\beta^{1}$ is first defined per halfedges of a
$p$–polygonal face $f$ as:

$\star_{1}\beta=\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top}\beta,$
(14)

where $\operatorname{\mathsf{R}}$ is the matrix defined in (10) and
$\operatorname{\mathsf{W}}_{1}$ is a symmetric $p\times p$ matrix given by:

$\operatorname{\mathsf{W}}_{1}[i,j]=\frac{\langle e_{i},e_{j}\rangle}{|f|},$

for $e_{k}$ the halfedges incident to and having the same orientation as the
face $f$, where $\langle.,.\rangle$ denotes the Euclidean dot product.

If an edge $e$ is not on boundary, it has two adjacent faces and two
halfedges, thus we compute the values of $\star\beta$ on corresponding
halfedges, sum their values with appropriate orientation sign and divide the
result by 2, see an example in Figure 9.

Figure 9: Let $\beta\in C^{1}$ and $e=(v_{0},v_{1})$ be the edge with
$e_{0},e_{4}$ as the corresponding halfedges, then
$\star\beta(e)=\frac{\star\beta(e_{0})-\star\beta(e_{4})}{2}$, where
$\star\beta(e_{4})$ is a linear combination of values of $\beta$ on dashed
orange halfedges and $\star\beta(e_{0})$ is a linear combination of values of
$\beta$ on dashed blue halfedges, concretely
$\star\beta(e_{0})=\frac{1}{4|f_{0}|}\Big{(}(\langle
e_{0},e_{1}\rangle-\langle
e_{0},e_{3}\rangle)(\beta(e_{0})-\beta(e_{2}))+(\langle
e_{0},e_{0}\rangle-\langle
e_{0},e_{2}\rangle)(\beta(e_{3})-\beta(e_{1})\Big{)}.$

The Hodge dual of a 0–form $\alpha$ is a 2–form $\star\alpha$ defined per a
$p$–polygonal face $f$ by:

$(\star_{0}\alpha)(f)=\frac{|f|}{p}\sum_{v_{i}\succ f}\alpha(v_{i}),$ (15)

and it is simply the arithmetic mean of the values of $\alpha$ on vertices of
the given face $f$ multiplied by the vector area $|f|$.

In matrix form, the discrete Hodge star operators read

$\displaystyle\star_{0}$
$\displaystyle=\operatorname{\mathsf{W}}_{F}\operatorname{\mathsf{fv}},$
$\displaystyle\star_{1}$
$\displaystyle=\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top},$
$\displaystyle\star_{2}$
$\displaystyle=\operatorname{\mathsf{W}}_{V}\operatorname{\mathsf{fv}}^{\top},$

where $\operatorname{\mathsf{fv}}$ and $\operatorname{\mathsf{R}}$ are defined
in equations (6) and (10), resp., and
$\operatorname{\mathsf{W}}_{F}\in\mathbb{R}^{|F|\times|F|}$,
$\operatorname{\mathsf{W}}_{V}\in\mathbb{R}^{|V|\times|V|}$,
$\operatorname{\mathsf{A}}\in\mathbb{R}^{|E|\times|E|}$ are given by

$\operatorname{\mathsf{W}}_{F}[i,i]=|f_{i}|,\quad\operatorname{\mathsf{W}}_{V}[i,i]=\frac{1}{\sum\limits_{f_{k}\succ
v_{i}}\frac{|f_{k}|}{p_{k}}},\quad\operatorname{\mathsf{A}}[i,j]=\left\\{\begin{array}[]{cl}1&\mbox{if
}i=j,\;e_{i}\text{ is on boundary},\\\ \frac{1}{2}&\mbox{if }i=j,\;e_{i}\text{
is not on boundary},\\\ -\frac{1}{2}&\mbox{if }e_{i}=-e_{j},\\\
0&\mbox{otherwise.}\end{array}\right.$

Although our Hodge star matrices are not diagonal, they are highly sparse and
thus computationally efficient. We have performed numerical tests on linear,
quadratic, and trigonometric forms on planar and curved meshes and they
exhibit the same at least linear convergence rate. We give an example in
Figure 10.

Figure 10: The approximation errors of the discrete Hodge star on a set of
irregular polygonal meshes on the torus. We have chosen
$\alpha^{0}=x^{2}+y^{2}$, $\beta^{1}=\mathbb{X}^{\flat}$, where
$\mathbb{X}=(-y,x,0)$ is a tangent vector field, and $\omega^{2}=\mu$ is the
area element on the torus. Thus $\star\mu=1$, $\star\alpha=(x^{2}+y^{2})\mu$,
and $\star\beta=\mathbb{Y}^{\flat}$, where
$\mathbb{Y}=2(-xz,-yz,x^{2}+y^{2}-\sqrt{x^{2}+y^{2}})$ is a tangent vector
field orthogonal to $\mathbb{X}$. On the right are two examples of meshes on
the torus with 5k vertices and 20k vertices.

###  3.3 The Hodge Inner Product

The $L^{2}$–Hodge inner product of differential forms $\Gamma^{k},\Omega^{k}$
on a Riemannian manifold $M$ is defined as:

$(\Gamma^{k},\Omega^{k}):=\int_{M}\Gamma\wedge\star\Omega.$

We define a discrete $L^{2}$–Hodge inner product of two discrete forms
$\alpha^{k},\beta^{k}$ on a mesh $S$ by:

$(\alpha^{k},\beta^{k}):=\sum_{f\in
S}\big{(}\alpha\wedge\star\beta\big{)}(f)=\alpha^{\top}\operatorname{\mathsf{M}}_{k}\beta,\;k=0,1,2,$

where $\operatorname{\mathsf{M}}_{k}$ are the discrete Hodge inner product
matrices that read:

$\displaystyle\operatorname{\mathsf{M}}_{0}=$
$\displaystyle\operatorname{\mathsf{fv}}^{\top}\operatorname{\mathsf{W}}_{F}\operatorname{\mathsf{fv}},$
$\displaystyle\operatorname{\mathsf{M}}_{1}=$
$\displaystyle\operatorname{\mathsf{R}}\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top},$
$\displaystyle\operatorname{\mathsf{M}}_{2}=$
$\displaystyle\operatorname{\mathsf{fv}}\operatorname{\mathsf{W}}_{V}\operatorname{\mathsf{fv}}^{\top}.$

It can be shown that our inner product of 1–forms restricted to a single face
$f$ is identical to the one of [Alexa and Wardetzky, 2011, Lemma 3]:
$\operatorname{\mathsf{R}}\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top}|_{f}=M_{f},$
where $M_{f}$ is defined as in equation (12). However, if a given mesh $S$ is
not just a single face, they differ in general, i.e., for 1–forms $\beta^{1}$,
$\gamma^{1}$:

$\beta^{\top}\operatorname{\mathsf{M}}_{1}\gamma=\beta^{\top}\operatorname{\mathsf{R}}\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top}\gamma\;\neq\;\beta^{\top}\operatorname{\mathsf{R}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top}\gamma=\sum_{f\in
S}\beta^{\top}M_{f}\gamma.$

To numerically evaluate our inner products, we calculate our discrete
$L^{2}$–Hodge norms of forms $\alpha^{0},\beta^{1},\omega^{2}$ over a mesh $S$
and compare them to their respective analytical $L^{2}$ norms. That is, if
$\Gamma^{k}$ is a differential $k$–form and $\gamma^{k}$ the corresponding
discrete $k$–form, we compute the error of approximation as:

$\int\displaylimits_{S}\Gamma\wedge\star\Gamma-\sum_{f\in
S}\gamma\wedge\star\gamma=\int\displaylimits_{S}\Gamma\wedge\star\Gamma-\gamma^{\top}\operatorname{\mathsf{M}}_{k}\gamma.$

An example of numerical evaluation of our $L^{2}$–Hodge inner products and
numerical evaluation of inner products $M_{0}$ and $M_{1}$ of Alexa and
Wardetzky [2011], see also the equations (11 – 12), is given in Figure 11. The
experimental convergence rate of our discrete $L^{2}$ norms is at least linear
on all tested forms on compact manifolds with or without boundary.

Figure 11: The influence of irregularity on experimental convergence of
discrete Hodge inner products to analytically computed solutions. Concretely,
on a set of jittered meshes with $r=0.2$ (far left) and $r=0.4$ (center left).
We show a sample mesh with 5k vertices for $r=0.2$ (center right) and $r=0.4$
(far right). We set $\alpha^{0}=x^{2}+y^{2}$, $\beta^{1}=-xzdx-
yzdy+(x^{2}+y^{2})dz$, $\omega^{2}=xdy\wedge dz+ydz\wedge dx+zdx\wedge dy$.
Here $\alpha\wedge\star\alpha$ denotes our Hodge inner product on $\alpha$,
and similarly for $\beta$ and $\omega$. $\alpha\mathsf{AW}$ denotes the inner
product of 0–forms and $\beta\mathsf{AW}$ the product of 1–forms of Alexa and
Wardetzky [2011].

###  3.4 The Contraction Operator

The contraction operator $\operatorname{\mathbf{i}}_{X}$, also called the
interior product, is the map that sends a $k$–form $\omega$ to a $(k-1)$–form
$\operatorname{\mathbf{i}}_{X}\omega$ such that
$(\operatorname{\mathbf{i}}_{X}\omega)(X_{1},\dots,X_{k-1})=\omega(X,X_{1},\dots,X_{k-1})$
for any vector fields $X_{1},\dots,X_{k-1}$. The following property holds
[Hirani, 2003, Lemma 8.2.1]:

###### Lemma 3.1

Let $M$ be a Riemannian $n$–manifold, $X\in\mathfrak{X}(M)$ a vector field,
then for the contraction of a differential $k$–form $\alpha$ with a vector
field $X$ holds:

$i_{X}\alpha=(-1)^{k(n-k)}\star(\star\alpha\wedge X^{\flat}),$

where $\flat:\mathfrak{X}(M)\rightarrow\Omega(M)$ is the flat operator.

Since we already have discrete wedge and Hodge star operators that are
compatible with each other, we can employ the lemma to define our discrete
contraction operator $\operatorname{\mathbf{i}}_{X}:C^{k}(S)\rightarrow
C^{k-1}(S)$ on a polygonal mesh $S$ by:

$\operatorname{\mathbf{i}}_{X}\alpha=(-1)^{k(2-k)}\star(\star\alpha\wedge
X^{\flat}),\;\;\alpha\in C^{k}(S),\;k=1,2,$ (16)

where the discrete flat operator on a vector field $X$ is given by
discretizing its value over all edges of $S$. Let $e=(v_{0},v_{1})$ be an edge
of $S$, then $e=e(t)=v_{0}+(v_{1}-v_{0})t,t\in[0,1]$, and we set:

$X^{\flat}(e)=\int\limits_{e}\langle
e^{\prime},X\rangle=\int\limits_{0}^{1}\langle e^{\prime}(t),X(e(t))\rangle
dt.$ (17)

Thus the discrete contraction operator is a linear operator that maps
$k$–forms located on $k$–dimensional primal cells to $(k-1)$–forms located on
$(k-1)$–dimensional primal cells.

Our discrete contraction of differential 2–forms with respect to different
vector fields exhibit linear convergence to the analytically computed
solutions, both in $L^{\infty}$ and $L^{2}$ norms. On 1–forms, the errors of
approximation decrease linearly in $L^{2}$ and with slope $0.5$ in
$L^{\infty}$ norm, see two examples in Figure 12.

Figure 12: The contraction operator on a unit sphere (L) and a torus (R). For
the sphere we have used the same set of jittered meshes with $r=0.4$ as in
Figure 11, and contracted the same forms as therein with respect to vector
field $X=(-y,x,0)$. For the torus (R) we have contracted the differential
forms of Figure 10 with respect to vector field
$X=2(-xz,-yz,x^{2}+y^{2}-\sqrt{x^{2}+y^{2}})$ on the same set of meshes as
therein. $\operatorname{\mathbf{i}}_{X}\beta\;\mathsf{L2}$ denotes the $L^{2}$
error approximation of the contraction operator on the 1–form $\beta$, whereas
$\operatorname{\mathbf{i}}_{X}\beta\;\mathsf{Inf}$ denotes the $L^{\infty}$
error approximation on $\beta$, and similarly for the 2–form $\omega$.

###  3.5 The Lie Derivative

We define the discrete Lie derivative
$\operatorname{\mathsf{L}}_{X}:C^{k}(S)\rightarrow C^{k}(S)$ using Cartan’s
magic formula:

$\operatorname{\mathsf{L}}_{X}\alpha=\operatorname{\mathbf{i}}_{X}d\alpha+d\operatorname{\mathbf{i}}_{X}\alpha,\;\alpha\in
C^{k}(S),\;k=0,1,2.$ (18)

Unfortunately, the Leibniz product rule of our contraction operator and Lie
derivative with discrete exterior derivative is not satisfied in general.
Concretely

$\displaystyle\operatorname{\mathbf{i}}_{X}(\alpha^{k}\wedge\beta^{l})$
$\displaystyle=$
$\displaystyle(\operatorname{\mathbf{i}}_{X}\alpha^{k})\wedge\beta^{l}+(-1)^{k}\alpha^{k}\wedge(\operatorname{\mathbf{i}}_{X}\beta^{l}),$
$\displaystyle\operatorname{\mathsf{L}}_{X}(\alpha^{k}\wedge\beta^{l})$
$\displaystyle=$
$\displaystyle(\operatorname{\mathsf{L}}_{X}\alpha^{k})\wedge\beta^{l}+\alpha^{k}\wedge(\operatorname{\mathsf{L}}_{X}\beta^{l}),$

holds only if $\alpha$ or $\beta$ is a closed 0–form. Already Hirani [2003]
noticed that the Leibniz rule for Lie derivative might not hold due to the
discrete wedge product not being associative in general. We confirm the
observation of Desbrun et al. [2005] that the Leibniz rule may be satisfied
only for closed forms.

Figure 13: The influence of jittering on experimental convergence of discrete
Lie derivatives on regular polygonal meshes on a sphere (far and center left)
and jittered meshes with vertex displacement by $0.4\times$ shortest edge
length (center and far right). We use the same set of jittered meshes as in
Figure 11. We employ the same forms and a vector field as in Figure 12.

The Lie derivatives exhibit converging behavior on all tested forms on regular
meshes, planar and non–planar. However, the $L^{2}$ error of approximation of
Lie derivatives of 1– and 2–forms on irregular meshes stays rather constant,
see an example on a set of regular versus jittered meshes on a unit sphere in
Figure 13. In this figure we can see that the $L^{2}$ error of the Lie
derivative of a 1–form $\beta$ and a 2–form $\omega$ on regular meshes
decreases with slope $-0.5$, whereas on very irregular meshes it stays
constant.

Although our discrete Lie derivative of 1–forms on irregular meshes does not
converge, in general, to analytically computed solutions, it can still be
employed for Lie advection (Section 4.3) of vector fields on irregular meshes
and produce visually satisfying results, as we demonstrate in Figure 20.

Figure 14: Lie derivative used for advection of a color function (far left)
encoded as a 0–form $\beta$ and advected using equation (23). The advected
function after 2000 (center left), 4000 (center right), and 5000 iterations
(right) with time step $10^{-2}$.

###  3.6 The Codifferential Operator

Just like on Riemannian $n$–manifolds, we define our discrete codifferential
operator $\delta$ as

$\delta_{k}\beta^{k}=(-1)^{n(k-1)+1}\star d\star\beta,\;\text{ $\beta$ a
discrete $k$-form}.$

Thus in matrix form our codifferential operators read:

$\displaystyle\delta_{1}=$
$\displaystyle-\operatorname{\mathsf{W}}_{V}\operatorname{\mathsf{fv}}^{\top}\,d_{1}\,\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top},$
$\displaystyle\delta_{2}=$
$\displaystyle-\operatorname{\mathsf{A}}\operatorname{\mathsf{W}}_{1}\operatorname{\mathsf{R}}^{\top}\,d_{0}\,\operatorname{\mathsf{W}}_{V}\operatorname{\mathsf{fv}}^{\top}.$

If $M$ is a compact manifold without boundary or if $\alpha$ or $\star\beta$
has zero boundary values, then the codifferential is the adjoint operator of
the exterior derivative with respect to the $L^{2}$–Hodge inner product:
$(d\alpha,\beta)=(\alpha,\delta\beta)\quad\forall\alpha\in\Omega^{k-1}(M),\beta\in\Omega^{k}(M).$
Alexa and Wardetzky [2011] use this equation to derive their discrete
codifferential operator on 1–forms, that reads

$\delta_{1}=M_{0}^{-1}d_{0}^{\top}M_{1},$

where $M_{0}$ and $M_{1}$ are as in equations (11–12). This codifferential
reduces to the classical codifferential, e.g., [Desbrun et al., 2006], in the
case of a pure triangle mesh.

Figure 15: The influence of jittering on experimental convergence of
codifferentials of
$\beta^{1}=(\sin(2x)+\cos(\frac{y}{2}))dx+(3\sin(x)-\cos(y))dy$ and
$\kappa^{2}=(\sin(\frac{x+1}{4})+\cos(1-\frac{y}{3}))dx\wedge dy$ on a set of
planar quadrilateral jittered meshes with $r=0.01$ (far left) and $r=0.2$
(center left). $\delta\beta$ denotes the approximation error of our
codifferential of $\beta$, and similarly for $\delta\kappa$.
$\beta\mathsf{AW}$ stands for the $L^{2}$ error of the codifferential of Alexa
and Wardetzky [2011]. On the center and far right are samples of such jittered
meshes.

In Figure 15 we test numerically our discretization and compare it to the
codifferential of Alexa and Wardetzky [2011]. We observe that the $L^{2}$
errors become constant. Concretely, for the jittered meshes with $r=0.4$, we
get circa $5.69\cdot 10^{-2}$ for $\delta\beta$ and $2.82\cdot 10^{-1}$ for
$\beta\mathsf{AW}$, i.e., our approximation error is roughly $5\times$
smaller. We have seen this difference on more or less irregular planar and
non–planar meshes for trigonometric, linear, and quadratic forms.

###  3.7 The Laplace Operator

The Laplace–deRham operator $\Delta$ takes differential $k$–forms to $k$–forms
and is defined as $\Delta=d\delta+\delta d$, where $\delta$ is the
codifferential and $d$ is the exterior derivative. On 0–forms (functions), it
simplifies to

$\Delta=\delta d.$

We define our discrete Laplacian in the same manner, using our codifferential.

Our Laplace operator is linearly precise, i.e., it is zero on linear 0–forms
in plane. In Figure 16 we depict the numerical behavior of our Laplacian on a
trigonometric 0–form and compare it to the combinatorially enhanced (so called
$\lambda$–simple choices with $\lambda=1$, $\lambda=2$) and purely geometric
($\lambda=0$) Laplacians of Alexa and Wardetzky [2011]. We note that all the
$L^{2}$ errors become constant. Our and the purely geometric Laplacian give a
better approximation to the analytical Laplacian than the combinatorially
enhanced Laplacians. We have observed this pattern also on different quadratic
and trigonometric 0–forms on more or less irregular polygonal meshes.

Figure 16: The discrete Laplacian of a trigonometric 0–form
$\alpha^{0}=sin(x-1)-cos(2y)$ on two sets of planar quadrilateral jittered
meshes with $r=0.01$ (far and center left) and $r=0.2$ (center and far right).
We use the same set of meshes as in Figure 15. $\Delta\alpha$ denotes our
Laplacian, $\alpha\,\mathsf{AW}\,\lambda=0$ the purely geometric,
$\alpha\,\mathsf{AW}\,\lambda=1$ and $\alpha\,\mathsf{AW}\,\lambda=2$ the
combinatorially enhanced Laplacians of Alexa and Wardetzky [2011]. The graphs
in the center left and on the far right are in arithmetic scale.

##  4 Applications

In this section we show some basic applications of our operators on general
polygonal meshes.

###  4.1 Implicit Mean Curvature Flow

One of the widely used methods for smoothing a surface is the implicit mean
curvature flow. If $f$ is a discrete 0–form representing vertex positions,
then $\Delta f$ give us the direction and magnitude in which we should move
each point in order to smooth the given mesh, see, e.g., [Desbrun et al.,
1999].

Let $f_{0}$ denote the initial state and $f_{t}$ the configuration after a
mean curvature flow of some duration $t>0$. We employ the backward Euler
scheme to calculate $f_{t}$ by solving the linear system:

$(I-t\Delta)f_{t}=f_{0},$

where $I$ is the identity matrix. To solve this system, we use the mldivide
algorithm of MATLAB.

In Figure 1 we show smoothing of general polygonal meshes and compare our
method to the one of Alexa and Wardetzky [2011] with purely geometric
Laplacians ($\lambda=0$) and combinatorially enhanced Laplacians
($\lambda=1$). After testing also other meshes and several other parameters
$\lambda$, time steps, and number of iterations, we conclude that our results
are visually comparable to theirs if $\lambda\in[1,2]$, and that our scheme
does not create as many undesirable artifacts as theirs for $\lambda=0$.

###  4.2 Helmholtz–Hodge Decomposition

By the Hodge Decomposition Theorem [Abraham et al., 1988, Theorem 7.5.3], if
$M$ is a compact oriented Riemannian manifold without boundary and
$\omega^{k}\in\Omega^{k}(M)$, then there exist uniquely determined forms
$\alpha^{k-1}$, $\beta^{k+1}$, $\gamma^{k}$ ($\gamma$ harmonic, i.e.,
$\Delta\gamma=0$) such that

$\omega=d\alpha+\delta\beta+\gamma.$ (19)

If instead of forms, we think about a sufficiently smooth vector field
$X=(\omega^{1})^{\sharp}$, where $\sharp$ is the sharp operator, then an
analogous Helmholtz theorem states that any vector field $X$ can be decomposed
into an irrotational vector field (corresponding to $d\alpha$), a
divergence–free component (analogous to $\delta\beta$), and a both
irrotational and divergence–free vector field (corresponding to $\gamma$).
Thus the equation (19) is also referred to as to Helmholtz–Hodge decomposition
(HHD).

If $X$ is a divergence–free vector field (also known as solenoidal), we can
find its two–component HHD, i.e., decompose $X$ into a rotational and
irrotational part. In terms of differential forms, for $\omega^{1}=X^{\flat}$
we get

$\omega=\delta\beta+\gamma,$ (20)

where $\gamma$ is a harmonic 1–form and thus $\gamma^{\sharp}$ is an
irrotational vector field, and $(\delta\beta)^{\sharp}$ is a rotational vector
field. The two–component HHD is used for decomposition of vector fields of
incompressible flows.

We use our codifferential operator to find our discrete two–component
Helmholtz–Hodge decomposition as in equation (20) by performing these steps:

  1. 1.

Discretize a given vector field $X$ with discrete flat operator (17) and
define discrete 1–form $\omega^{1}=X^{\flat}$.

  2. 2.

Find the 2–form $\beta$ by solving the equation $d\delta\beta=d\omega$.

  3. 3.

Set $\gamma=\omega-\delta\beta$.

Figure 17: The discrete sharp operator on a vertex $v$ restricted to a face
$f$. Let $X^{\flat}$ be a 1–from computed by applying a discrete flat operator
on a constant vector field $X$, then the orthogonal projection of $X$ on the
unit direction vector of the edge $e_{1}$ equals $\frac{\langle
X,e_{1}\rangle}{|e_{1}|}=\frac{X^{\flat}(e_{1})}{|e_{1}|}$. Similarly the
orthogonal projection of $X$ on the unit direction vector of $e_{2}$ is
$\frac{\langle X,e_{2}\rangle}{|e_{2}|}=\frac{X^{\flat}(e_{2})}{|e_{2}|}$.
Reconstructing the vector field $X$ from $X^{\flat}$, i.e., applying the sharp
operator on $X^{\flat}$ as in equation (21), yields vector
$(X^{\flat})^{\sharp}|_{f}=\frac{\langle
X,e_{2}\rangle}{|e_{2}|}\frac{n_{f}\times e_{1}}{|e_{1}|}-\frac{\langle
X,e_{1}\rangle}{|e_{1}|}\frac{n_{f}\times e_{2}}{|e_{2}|}$ that has the same
direction as $X$ and approximates its magnitude.

We can then map the discrete 1–forms $\delta\beta$ and $\gamma$ to discrete
vector fields by applying discrete sharp operator $\sharp$ defined on an
1–form $\epsilon$ and per a vertex $v$ by:

$\epsilon^{\sharp}(v)=\frac{1}{\rho(v)}\sum_{f\succ
v}\bigg{(}\frac{\epsilon(e_{2})}{|e_{2}|}\frac{n_{f}\times
e_{1}}{|e_{1}|}-\frac{\epsilon(e_{1})}{|e_{1}|}\frac{n_{f}\times
e_{2}}{|e_{2}|}\bigg{)},$ (21)

where $\rho(v)$ is the number of faces adjacent to $v$. Further
$e_{1},e_{2}\prec f$, $e_{1}$ is the edge which endpoint is $v$, $e_{2}$ is
the edge with $v$ as the starting point, see Figure 17, and $n_{f}$ is a unit
normal vector of the face $f=(v_{0},\dots,v_{n-1})$ computed as:

$n_{f}=\frac{\hat{n}_{f}}{|\hat{n}_{f}|},\quad\hat{n}_{f}=\frac{1}{2}\sum_{j=0}^{n-1}(v_{j}\times
v_{j+1}),\;\text{indices modulo }n.$

In Figure 18 we give an example of our HHD of an incompressible vector field
on a general polygonal mesh of a torus. In Figure 19 we then employ the HHD to
remove vortices of an arbitrary vector field.

Figure 18: HHD of a solenoidal vector field $X$ on a torus centered at the
origin (a general polygonal mesh with 20k vertices). The decomposed vector
field is $X=X_{H}+X_{R}$, where $X_{H}=(-y,x,0)$ is a harmonic field on the
torus, and $X_{R}$ is a rotational vector field given by
$X_{R}=\nabla(\exp(-(x-x_{1})^{2}-(y-y_{1})^{2}-(z-z_{1})^{2})-\exp(-(x-x_{2})^{2}-(y-y_{2})^{2}-(z-z_{2})^{2}))\times
n,$ where $n$ is the unit normal vector of the torus. We have chosen the
center of CCW rotation $(x_{1},y_{1},z_{1})=(\frac{3}{2},0,0)$, where the
vector potential $\beta^{\sharp}$ reaches its maximum, and the center of CW
rotation at
$(x_{2},y_{2},z_{2})=(\frac{-\sqrt{2}}{2},\frac{\sqrt{2}}{2},\frac{1}{2})$,
where the vector potential $\beta^{\sharp}$ has its minimum (is negative). The
vector field $X$ is shown on the far left. Our discrete decomposition gives
approximate expected results: the harmonic part $\gamma^{\sharp}$ calculated
by our method is in the center left and the rotational part
$(\delta\beta)^{\sharp}$ in the the center right. On the far right we
visualize in pseudocolors the vector potential $\beta^{\sharp}$ computed by
our algorithm.

Figure 19: HHD applied to remove the vortices of a vector field $\omega$ on a
mesh of Spot (model created by Keenan Crane). On the left is the original
vector field $\omega=\delta\beta+\phi$, on the right is its curl–free part
$\phi$.

###  4.3 Lie Advection

The Lie derivative finds its application in dynamical systems. In computer
graphics the Lie advection of differential forms (including scalar and vector
fields) is used for tasks ranging from fluid flow simulation [McKenzie, 2007]
to authalic parametrization of surfaces [Zou et al., 2011].

In Figure 20 we employ our Lie derivative to perform a simple discrete Lie
advection of a tangent vector field $Y$ along a tangent vector field
$X=(-y,x,0)$ on a torus azimuthally symmetric about the $z$ axis. $Y$ is a
vorticial vector field given as

$Y=-\nabla\exp\bigg{(}-\Big{(}x+\frac{\sqrt{2}}{2}\Big{)}^{2}-\Big{(}y-\frac{\sqrt{2}}{2}\Big{)}^{2}-\Big{(}z-\frac{1}{2}\Big{)}^{2}\bigg{)}\times
n,$ (22)

where $n$ is the unit normal vector of the torus. To advect $Y$, we discretize
it as a 1–form $\beta=Y^{\flat}$ and denote this initial state as $\beta_{0}$.
We then iterate over our discrete solutions using a simple forward Euler
method:

$\beta_{k+1}=\beta_{k}-t\operatorname{\mathsf{L}}_{X}\beta_{k},\;k=0,\dots,$
(23)

where $t$ is the time step, $k$ is the number of iterations, and each
$\operatorname{\mathsf{L}}_{X}\beta_{k}$ is computed using our discrete Lie
derivative (18). Note that the vector field $X$ is also discretized as a
discrete 1–form.

Figure 20: The Lie advection on a polygonal mesh on a torus (8k vertices), the
mesh is shown in the top left corner. We advect an 1–form $\beta=Y^{\flat}$
along the flow of a tangent vector field $X=(-y,x,0)$, where $Y$ is given in
equation (22). The second picture from the left depicts $\beta_{0}^{\sharp}$.
We apply time steps of length $10^{-3}$. From left to right and top to bottom,
we plot $Y=\beta^{\sharp}_{k}$ after 1000, 2000,…, 5000 iterations. At the
bottom right we see $\beta_{k}^{\sharp}$ for $k=6283$. Because the domain is
periodic, $\beta$ should be advected back to its original state $\beta_{0}$
after $2\pi\cdot 10^{3}\approx 6283$ iterations. We can see that
$\beta^{\sharp}_{6283}$ gets close to $\beta^{\sharp}_{0}$, but some
undesirable artifacts appear – especially at the small irregularity around the
spot with coordinate position
$(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},\frac{1}{2})$.

The Lie derivative can be employed also for advection of a function by a
vector field. In Figure 14 we advect a color function on a mesh of a vase.

##  5 Discussion

Geometry processing with general polygonal meshes is a new developing area. We
propose various discrete operators that act directly on meshes made of
arbitrary polygons, possibly non–planar and non–convex, and thus open the
possibility to perform many geometry processing tasks directly on these
meshes.

Tangent vector fields on surfaces are used in many applications in computer
graphics and other areas. We propose to represent vector fields as 1–forms and
we provide methods for their design such as the Helmholtz–Hodge decomposition
or the Lie advection. However, further applications are now ready to be
explored, e.g., finding pairs of vector fields with zero Lie derivative for
surface parameterization or employing Le advection of volume forms for area
preserving parametrization.

Furthermore, we present a novel discrete Laplace operator that is numerically
comparable to the purely geometric Laplacian of Alexa and Wardetzky [2011],
but results in a better mesh smoothing. On the other hand, our Laplacian gives
a better numerical approximation to the analytically computed solutions than
their combinatorially enhanced Laplacians, yet performs as well as theirs in
smoothing of tested general polygonal meshes.

## References

  * Abraham et al. [1988] Abraham, R., Marsden, J.E., Ratiu, T., 1988.  Manifolds, Tensor Analysis, and Applications: 2Nd Edition.  Springer-Verlag New York, Inc., New York, NY, USA.
  * Alexa and Wardetzky [2011] Alexa, M., Wardetzky, M., 2011\.  Discrete Laplacians on general polygonal meshes, in: ACM SIGGRAPH 2011 Papers, ACM, New York, NY, USA. pp. 102:1–102:10.  URL: http://doi.acm.org/10.1145/1964921.1964997, doi:10.1145/1964921.1964997.
  * Arnold [2012] Arnold, R.F., 2012.  The Discrete Hodge Star Operator and Poincaré Duality.  Ph.D. thesis. Virginia Polytechnic Institute and State University.
  * Azencot et al. [2013] Azencot, O., Ben-Chen, M., Chazal, F., Ovsjanikov, M., 2013\.  An operator approach to tangent vector field processing, in: Proceedings of the Eleventh Eurographics/ACMSIGGRAPH Symposium on Geometry Processing, Eurographics Association, Aire-la-Ville, Switzerland, Switzerland. pp. 73–82.  URL: http://dx.doi.org/10.1111/cgf.12174, doi:10.1111/cgf.12174.
  * Crane et al. [2013] Crane, K., de Goes, F., Desbrun, M., Schröder, P., 2013\.  Digital geometry processing with discrete exterior calculus, in: ACM SIGGRAPH 2013 Courses, ACM, New York, NY, USA. pp. 7:1–7:126.  URL: http://doi.acm.org/10.1145/2504435.2504442, doi:10.1145/2504435.2504442.
  * Desbrun et al. [2005] Desbrun, M., Hirani, A.N., Leok, M., Marsden, J.E., 2005\.  Discrete exterior calculus.  URL: http://arXiv.org/math.DG/0508341. preprint.
  * Desbrun et al. [2006] Desbrun, M., Kanso, E., Tong, Y., 2006.  Discrete differential forms for computational modeling, in: ACM SIGGRAPH 2006 Courses, ACM, New York, NY, USA. pp. 39–54.  URL: http://doi.acm.org/10.1145/1185657.1185665, doi:10.1145/1185657.1185665.
  * Desbrun et al. [1999] Desbrun, M., Meyer, M., Schröder, P., Barr, A.H., 1999\.  Implicit fairing of irregular meshes using diffusion and curvature flow, in: Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques, ACM Press/Addison-Wesley Publishing Co., USA. p. 317–324.  URL: https://doi.org/10.1145/311535.311576, doi:10.1145/311535.311576.
  * Gu and Yau [2003] Gu, X., Yau, S.T., 2003.  Global conformal surface parameterization, in: Proceedings of the 2003 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing, Eurographics Association, Aire-la-Ville, Switzerland, Switzerland. pp. 127–137.  URL: http://dl.acm.org/citation.cfm?id=882370.882388.
  * Hirani [2003] Hirani, A.N., 2003.  Discrete Exterior Calculus.  Ph.D. thesis. California Institute of Technology. Pasadena, CA, USA.  AAI3086864.
  * MacNeal [1949] MacNeal, R.H., 1949.  The solution of partial differential equations by means of electrical networks.  Ph.D. thesis. California Institute of Technology. Pasadena, CA, USA.
  * Massey [1991] Massey, W.S., 1991.  A Basic Course in Algebraic Topology.  Graduate Texts in Mathematics, Springer New York.
  * McKenzie [2007] McKenzie, A., 2007.  HOLA: A High-Order Lie Advection of Discrete Differential Forms, with Applications in Fluid Dynamics.  Master’s thesis. California Institute of Technology. Division of Engineering and Applied Science.
  * Mullen et al. [2011] Mullen, P., McKenzie, A., Pavlov, D., Durant, L., Tong, Y., Kanso, E., Marsden, J.E., Desbrun, M., 2011\.  Discrete Lie advection of differential forms.  Found. Comput. Math. 11, 131–149.  URL: http://dx.doi.org/10.1007/s10208-010-9076-y, doi:10.1007/s10208-010-9076-y.
  * Whitney [1957] Whitney, H., 1957.  Geometric Integration Theory.  Princeton University Press.
  * Zou et al. [2011] Zou, G., Hu, J., Gu, X., Hua, J., 2011.  Authalic parameterization of general surfaces using Lie advection.  IEEE Transactions on Visualization and Computer Graphics 17, 2005–2014.  URL: http://dx.doi.org/10.1109/TVCG.2011.171, doi:10.1109/TVCG.2011.171.


16 | PR100013_TG001703_V0200 | tra | 2021-11-11T22:59:21 | 9.1 | 298 | 54 | 2224 | 1.19
17 | PR100016_TG015001_V0200 | occ | 2021-12-05T09:14:25 | 11.5 | 384 | 55 | 2245 | 1.00
18 | PR100013_TG001704_V0200 | tra | 2021-12-05T20:55:21 | 10.5 | 361 | 57 | 2246 | 0.93
19 | PR100013_TG001705_V0200 | tra | 2021-12-07T00:38:23 | 10.6 | 362 | 56 | 2247 | 0.93
20 | PR100016_TG015002_V0200 | occ | 2021-12-07T11:28:21 | 11.1 | 373 | 55 | 2247 | 1.08
21 | PR100013_TG001706_V0200 | tra | 2021-12-15T16:39:21 | 12.2 | 424 | 58 | 2255 | 1.12
22 | PR100016_TG015003_V0200 | occ | 2021-12-16T05:36:20 | 12.4 | 441 | 59 | 2255 | 1.33
23 | PR100013_TG001707_V0200 | tra | 2021-12-16T18:51:21 | 12.2 | 425 | 58 | 2256 | 1.27
24 | PR100016_TG015004_V0200 | occ | 2021-12-17T07:12:21 | 12.6 | 407 | 53 | 2256 | 1.04
25 | PR100016_TG015005_V0200 | occ | 2021-12-24T22:55:20 | 10.7 | 379 | 58 | 2263 | 1.21
26 | PR100013_TG001708_V0200 | tra | 2021-12-26T14:35:21 | 9.8 | 327 | 55 | 2265 | 1.32
27 | PR100013_TG001709_V0200 | tra | 2021-12-29T21:36:21 | 11.7 | 394 | 55 | 2268 | 1.64
28 | PR100016_TG015006_V0200 | occ | 2021-12-30T09:32:22 | 13.6 | 479 | 58 | 2268 | 0.98
29 | PR100013_TG001710_V0200 | tra | 2021-12-30T23:58:21 | 10.8 | 394 | 60 | 2269 | 1.82
30 | PR100016_TG015007_V0200 | occ | 2021-12-31T11:43:21 | 14.9 | 522 | 58 | 2269 | 1.49
31 | PR100013_TG001711_V0200 | tra | 2022-01-04T08:15:21 | 9.8 | 331 | 56 | 2273 | 1.68
32 | PR100013_TG001712_V0200 | tra | 2022-01-06T12:37:21 | 12.1 | 421 | 58 | 2275 | 1.34
33 | PR100016_TG015008_V0200 | occ | 2022-01-07T00:53:22 | 11.7 | 404 | 57 | 2275 | 1.17
34 | PR100013_TG001713_V0200 | tra | 2022-01-09T19:41:22 | 12.1 | 429 | 58 | 2278 | 1.10
35 | PR100013_TG001714_V0200 | tra | 2022-01-19T15:49:21 | 9.8 | 309 | 52 | 2287 | 1.44
36 | PR100013_TG001715_V0200 | tra | 2022-01-26T04:51:21 | 9.8 | 331 | 56 | 2293 | 1.28
37 | PR100013_TG001716_V0200 | tra | 2022-01-29T11:12:21 | 9.0 | 314 | 58 | 2296 | 1.35
38 | PR100013_TG001717_V0200 | tra | 2022-02-03T23:05:22 | 10.8 | 375 | 58 | 2301 | 1.09
39 | PR100013_TG001718_V0200 | tra | 2022-02-04T23:49:21 | 9.1 | 325 | 59 | 2302 | 0.90
40 | PR100013_TG001719_V0200 | tra | 2022-02-08T06:57:21 | 10.0 | 328 | 54 | 2305 | 1.08
41 | PR100013_TG001720_V0200 | tra | 2022-02-18T03:33:21 | 10.8 | 337 | 51 | 2314 | 1.14
42 | PR100016_TG015601_V0200 | occ | 2022-02-23T00:04:21 | 11.4 | 365 | 53 | 2318 | 1.00
43 | PR100016_TG015602_V0200 | occ | 2022-02-24T02:07:21 | 10.5 | 352 | 56 | 2319 | 1.22
44 | PR100016_TG015603_V0200 | occ | 2022-11-22T19:53:22 | 10.7 | 374 | 58 | 2568 | 1.68
45 | PR100016_TG015604_V0200 | PC | 2022-11-23T22:12:21 | 24.0 | 787 | 54 | 2569 | 1.21
46 | PR330093_TG000201_V0200 | occ | 2022-12-21T04:47:22 | 8.1 | 270 | 55 | 2594 | 1.00
47 | PR100016_TG015605_V0200 | occ | 2022-12-24T11:46:21 | 10.6 | 390 | 61 | 2597 | 1.53

Table 5: TESS observation log of WASP-12  Sector |  | Start date
---
[UTC]
| Duration
---
[days]
| Data
---
points [#]
| Exp.
---
time [s]
$\beta_{w}\beta_{r}$
20 | 2019-12-24 | 27 | 16552 | 120 | 0.58
43 | 2021-09-16 | 27 | 15577 | 120 | 0.88
44 | 2021-11-12 | 27 | 15777 | 120 | 0.97
45 | 2021-11-06 | 27 | 16085 | 120 | 1.23
Table 6: Spitzer transit observation log of WASP-12  | pass
---
band
| Start date
---
[UTC]
| Duration
---
[hrs]
| Data
---
points [#]
| Exp.
---
time [s]
$\beta_{w}\beta_{r}$
3.6$\mu$m | 2010-11-17 | 9.6 | 264 | 128 | 1.01
4.5$\mu$m | 2010-12-11 | 9.6 | 266 | 128 | 1.03
3.6$\mu$m | 2013-12-12 | 9.6 | 264 | 128 | 1.01
4.5$\mu$m | 2013-12-15 | 9.6 | 266 | 128 | 0.95
Table 7: Derived power-2 LDC priors in the different passbands and the posterior from the ellipsoidal and spherical planet model fits. Parameter | Prior | Posterior
---|---|---
|  |  | Spherical planet
---
| Ellipsoidal planet
---
| $c_{{}_{\mathrm{CHEOPS}}}$
---
$\alpha_{{}_{\mathrm{CHEOPS}}}$
|  $\mathcal{N}$(0.714, 0.013)
---
$\mathcal{N}$(0.631, 0.016)
| $0.706\pm{0.010}$
---
$0.624\pm{0.010}$
| $0.714\pm{0.010}$
---
$0.630\pm{0.011}$
| $c_{{}_{\mathrm{TESS}}}$
---
$\alpha_{{}_{\mathrm{TESS}}}$
|  $\mathcal{N}$(0.634,0.018)
---
$\mathcal{N}$(0.554,0.019)
| $0.627\pm 0.013$
---
$0.542\pm{0.014}$
| $0.636\pm 0.013$
---
$0.548\pm{0.014}$
| $c_{{}_{\mathrm{Spitzer3.6}}}$
---
$\alpha_{{}_{\mathrm{Spitzer3.6}}}$
|  $\mathcal{N}$(0.313, 0.011)
---
$\mathcal{N}$(0.351, 0.015)
| $0.316\pm 0.008$
---
$0.348\pm 0.011$
| $0.317\pm 0.010$
---
$0.346\pm 0.011$
| $c_{{}_{\mathrm{Spitzer4.5}}}$
---
$\alpha_{{}_{\mathrm{Spitzer4.5}}}$
|  $\mathcal{N}$(0.245, 0.011)
---
$\mathcal{N}$(0.398, 0.022)
| $0.241\pm 0.010$
---
$0.406\pm 0.020$
| $0.240\pm 0.096$
---
$0.406\pm 0.019$

Table 8: Derived mid-transit times for the individual CHEOPS transit observations of WASP-12 b. | Visit
---
[#]
Type | Epoch |  | Timing
---
[BJDTBD – 2459000]
$\sigma_{T}$ [d]
13 | tra | 1946 | 227.183907 | 0.000225
14 | tra | 2214 | 519.683809 | 0.000351
15 | tra | 2223 | 529.506215 | 0.000294
16 | tra | 2224 | 530.598008 | 0.000305
18 | tra | 2246 | 554.608948 | 0.000286
19 | tra | 2247 | 555.700197 | 0.000309
21 | tra | 2255 | 564.431242 | 0.000353
23 | tra | 2256 | 565.523282 | 0.000336
26 | tra | 2265 | 575.345658 | 0.000318
27 | tra | 2268 | 578.620503 | 0.000449
29 | tra | 2269 | 579.711388 | 0.000458
31 | tra | 2273 | 584.076978 | 0.000373
32 | tra | 2275 | 586.259887 | 0.000212
34 | tra | 2278 | 589.534364 | 0.000313
35 | tra | 2287 | 599.356962 | 0.000307
36 | tra | 2293 | 605.905245 | 0.00024
37 | tra | 2296 | 609.180016 | 0.000348
38 | tra | 2301 | 614.636613 | 0.000316
39 | tra | 2302 | 615.728374 | 0.000211
40 | tra | 2305 | 619.002623 | 0.000373
41 | tra | 2314 | 628.824968 | 0.000283
45b | tra | 2570 | 908.228003 | 0.000327
S1 | occ | 1900 | 177.523576 | 0.002121
S2 | occ | 2272 | 583.530597 | 0.002471
S3 | occ | 2582 | 921.869782 | 0.002904

##  Appendix B Figures

Figure 7: CHEOPS transit light curves of WASP-12b labeled according to the
visit number. Left: The best-fit transit and systematics model is overplotted
on the data. Middle: Systematics detrended flux with transit model
overplotted. The obtained transit time uncertainty for each visit is shown in
seconds. Right:  Residuals after subtraction of best-fit transit and
systematics model. The 30-min bins and the root-mean-square (rms) of each
visit are also shown. Figure 8: CHEOPS occultation light curves of WASP-12b
labeled according to their visit numbers. Left: The best-fit occultation and
systematics model is overplotted on the data. Middle: Systematics detrended
flux with occultation model overplotted. Right: Residuals after subtraction of
best-fit occultation and systematics model. The 30-min bins and the root-mean-
square (rms) of each visit are also shown. Figure 9:  Different approaches
adopted to model the stellar limb darkening using model stellar intensity
profiles from the PHOENIX (red) and ATLAS (blue) libraries. The plots
illustrate the approaches for the Spitzer 4.5$\mu$m passband observations.
Panel (a) illustrates our ”fiducial” approach where LDCs are generated to be
used as priors in the transit model fitting. The black curve shows the fit of
the power-2 law to the combined model intensity profiles. The parameter space
(gray) allowed by the 1$\sigma$ uncertainty of the obtained LDCs encompasses
both intensity profiles and associated 1$\sigma$ uncertainties. Panel (b)
shows the ”alternative-1” approach which merges the PHOENIX and ATLAS model
profiles to create a new joint intensity profile (cyan) whose 1$\sigma$
uncertainties at each $\mu$ encompasses the 1$\sigma$ uncertainty of the
individual profiles. This is illustrated in the inset for $\mu=0.4$. The new
joint profile is fitted with an LD law alongside the transit observation at
each passband. Panel (c) shows the ”alternative-2” approach which similarly
combines both intensity profiles but creates a uniform bound (cyan) spanning
the median of both profiles such that LD profile points within the bound have
equal likelihood but decrease as the profile points deviate from the bound.
Figure 10: Posterior distribution of the retrieval using Pyrat Bay with the
median and $1\sigma$.

Figure 11:  Corner plot showing the posterior probability distributions of the
fitted parameters for the spherical planet model (red) and ellipsoidal planet
model. Vertical colored lines indicate the median of the distributions while
shaded regions define the 68%credible intervals. The different instruments are
labeled: CHEOPS (ch), TESS (tess), Spitzer 3.6$\mu$m (sp1), and Spitzer
4.5$\mu$m (sp2).