Orientational Order on Surfaces: The Coupling of Topology, Geometry, and Dynamics

Abstract

We consider the numerical investigation of surface bound orientational order using unit tangential vector fields by means of a gradient flow equation of a weak surface Frank–Oseen energy. The energy is composed of intrinsic and extrinsic contributions, as well as a penalization term to enforce the unity of the vector field. Four different numerical discretizations, namely a discrete exterior calculus approach, a method based on vector spherical harmonics, a surface finite element method, and an approach utilizing an implicit surface description, the diffuse interface method, are described and compared with each other for surfaces with Euler characteristic 2. We demonstrate the influence of geometric properties on realizations of the Poincaré–Hopf theorem and show examples where the energy is decreased by introducing additional orientational defects.

Appendices

Appendix A: Thin-Film Limit of Penalized Frank–Oseen Energy

Considering a thin shell $\Omega _{\delta } = \mathcal {S}\times [-\delta /2,\delta /2]$ around the surface $\mathcal {S}$ with thickness $\delta$, the local coordinates $\theta$ and $\varphi$ of the surface immersion $\mathbf {x}$ and an additional coordinate $\xi$, which acts along the surface normal $\varvec{\nu }$, lead to a thin shell parametrization $\widetilde{\mathbf {x}}:U_{\delta } \rightarrow \mathbb {R}^{3}$ for the parameter domain $U_\delta := U\times [-\delta /2,\delta /2]$, with $\widetilde{\mathbf {x}}$ defined by

$$\begin{aligned} \widetilde{\mathbf {x}}(\theta ,\varphi ,\xi ) = \mathbf {x}(\theta ,\varphi ) + \xi \varvec{\nu }(\theta ,\varphi ) . \end{aligned}$$

(81)

The thickness $\delta$ is sufficiently small to guarantee the injectivity of the pushforward, see Napoli and Vergori (2012b).

For a better readability, we denote indices which mark all three components $\{\theta ,\varphi ,\xi \}$ by capital letters. The indices for the surface components $\{\theta ,\varphi \}$ are denoted by small letters. The metric tensor $\widetilde{\mathbf {g}}$ of the thin shell is given by its components $\widetilde{g}_{IJ} = \partial _{I}\widetilde{\mathbf {x}}\cdot \partial _{J}\widetilde{\mathbf {x}}$, i.e.,

$$\begin{aligned} \begin{aligned} \widetilde{g}_{ij} = g_{ij} - 2\xi \mathcal {B}_{ij} + \mathcal {O}(\xi ^{2})_{ij} = g_{ij} + \mathcal {O}(\xi )_{ij} , \quad \widetilde{g}_{\xi \xi } = 1 \text { and } \widetilde{g}_{i\xi } = \widetilde{g}_{\xi i} = 0 . \end{aligned} \end{aligned}$$

(82)

The pure formal indices on $\mathcal {O}$ extend the asymptotic polynomial behavior to tensor context and preserve summation conventions. Hence, for the Christoffel symbols $\widetilde{\varGamma }_{I J}^{K} = \frac{1}{2}\widetilde{g}^{KL}\left( \partial _{I}\widetilde{g}_{JL} + \partial _{J}\widetilde{g}_{IL} - \partial _{L}\widetilde{g}_{IJ} \right)$, we obtain

$$\begin{aligned} \begin{aligned} \widetilde{\varGamma }_{i j}^{k}&= \varGamma _{i j}^{k} + \mathcal {O}(\xi )_{ij}^{k}, \quad \widetilde{\varGamma }_{i j}^{\xi } = \mathcal {B}_{ij} + \mathcal {O}(\xi )_{ij}, \quad \widetilde{\varGamma }_{i \xi }^{k} = \widetilde{\varGamma }_{\xi i}^{k} = -{{\mathcal {B}}_{i}}^{k} + \mathcal {O}(\xi )_{i}^{k}, \\ \widetilde{\varGamma }_{\xi \xi }^{K}&= \widetilde{\varGamma }_{I \xi }^{\xi } = \widetilde{\varGamma }_{\xi I}^{\xi } = 0 . \end{aligned} \end{aligned}$$

(83)

We can approximate the square root of the determinant $|\widetilde{\mathbf {g}}|$ on $\mathcal {S}$ by $\sqrt{|\widetilde{\mathbf {g}}|} = \sqrt{\widetilde{g}_{\xi \xi }|\mathbf {g}|} + \mathcal {O}(\xi )= \left( 1 + \mathcal {O}(\xi )\right) \sqrt{|\mathbf {g}|}$. Therefore, the volume element becomes

$$\begin{aligned} \,\text {d}{V}= \sqrt{|\widetilde{\mathbf {g}}|}d\xi \wedge d\theta \wedge d\varphi = \left( 1 + \mathcal {O}(\xi )\right) d\xi \wedge \,\text {d}{\mathcal {S}}. \end{aligned}$$

(84)

The 3-tensor, with the same qualities as the volume element, is the Levi–Civita tensor

$$\begin{aligned} \widetilde{E}_{IJK} = \,\text {d}{V}\left( \partial _{I}\widetilde{\mathbf {x}}, \partial _{J}\widetilde{\mathbf {x}}, \partial _{K}\widetilde{\mathbf {x}}\right) = \sqrt{|\widetilde{\mathbf {g}}|}\varepsilon _{IJK} = \sqrt{|\mathbf {g}|}\varepsilon _{IJK} + \mathcal {O}(\xi )_{IJK}, \end{aligned}$$

(85)

with the common Levi–Civita symbols $\varepsilon _{IJK}\in \{-1,0,1\}$. With the Levi–Civita tensor $\mathbf {E}$ on the surface, defined by $E_{ij} = \,\text {d}{\mathcal {S}}\left( \partial _{i}\mathbf {x}, \partial _{j}\mathbf {x}\right) = \sqrt{|\mathbf {g}|}\varepsilon _{ij}$, and the fact that all non-vanishing components of the Levi–Civita tensor $\mathbf {\widetilde{E}}$ in the thin shell have exactly one $\xi$-index, we obtain

$$\begin{aligned} \widetilde{E}_{\xi ij} = - \widetilde{E}_{i \xi j} = \widetilde{E}_{ij \xi } = E_{ij} + \mathcal {O}(\xi )_{ij} . \end{aligned}$$

(86)

For a better distinction, we use a semicolon in the thin shell and a straight line on the surface to mark the components of the covariant derivative, i.e., for the vector fields $\widetilde{\mathbf {p}}\in C^{1}\left( \Omega _{\delta }, \mathsf {T}\Omega _{\delta } \right)$ and $\mathbf {p}\in C^{1}\left( \mathcal {S}, \mathsf {T}\mathcal {S}\right)$, we write

$$\begin{aligned} {{\widetilde{\text {p}}}^{I}}_{;J}&= \partial _{J}\widetilde{\text {p}}^{I} + \widetilde{\varGamma }_{J K}^{I}\widetilde{\text {p}}^{K} \text { and} \end{aligned}$$

(87)

$$\begin{aligned} {{p}^{i}}_{|j}&= \partial _{j}p^{i} + \varGamma _{j k}^{i}p^{k}. \end{aligned}$$

(88)

The contravariant derivatives are given by $\widetilde{\text {p}}^{I;J} = \widetilde{g}^{JK}{{\widetilde{\text {p}}}^{I}}_{;K}$ and $p^{i|j} = g^{jk}{{p}^{i}}_{|k}$. Henceforward, we assume that $\widetilde{\mathbf {p}}\in \mathsf {T}\Omega _{\delta }$ is an extension of $\mathbf {p}$, i.e., $\widetilde{\mathbf {p}}\big |_{\mathcal {S}} = \mathbf {p}\in \mathsf {T}\mathcal {S}$, and $\widetilde{\mathbf {p}}$ is parallel and length preserving in direction of $\varvec{\nu }$, i.e., ${{\widetilde{\text {p}}}^{I}}_{;\xi } = 0$ as a consequence. Therefore, the Taylor approximation on the surface of the contravariant tangential components becomes

$$\begin{aligned} \begin{aligned} \widetilde{\text {p}}^{i}&= p^{i} + \xi \partial _{\xi }\widetilde{\text {p}}^{i}\big |_{\mathcal {S}} + \mathcal {O}(\xi ^{2})^{i} = p^{i} +\xi \left( {{\widetilde{\text {p}}}^{i}}_{;\xi } - \widetilde{\varGamma }_{\xi K}^{i}\widetilde{\text {p}}^{K}\big |_{\mathcal {S}} \right) + \mathcal {O}(\xi ^{2})^{i} \\&= p^{i} +\xi {{\mathcal {B}}_{k}}^{i}p^{k} + \mathcal {O}(\xi ^{2})^{i} . \end{aligned} \end{aligned}$$

(89)

It holds $\widetilde{\text {p}}^{\xi } = 0$, because $\widetilde{\text {p}}^{\xi }\big |_{\mathcal {S}} = 0$ and $\partial _{\xi }\widetilde{\text {p}}^{\xi } = {{\widetilde{\text {p}}}^{\xi }}_{;\xi } - \widetilde{\varGamma }_{\xi K}^{\xi }\widetilde{\text {p}}^{K} = 0$, but nonetheless, we get non-vanishing covariant tangential derivatives

$$\begin{aligned} {{\widetilde{\text {p}}}^{\xi }}_{;j} = \widetilde{\varGamma }_{j K}^{\xi }\widetilde{\text {p}}^{K} = \mathcal {B}_{jk}p^{k} + \mathcal {O}(\xi )_{j} . \end{aligned}$$

(90)

All remaining covariant derivatives can be approximated by

$$\begin{aligned} {{\widetilde{\text {p}}}^{i}}_{;j}&= \partial _{j}\widetilde{\text {p}}^{i} + \widetilde{\varGamma }_{j K}^{i}\widetilde{\text {p}}^{K} = \partial _{j}\text {p}^{i} + \varGamma _{j k}^{i}\text {p}^{k} + \mathcal {O}(\xi )^{i}_{j} = {{p}^{i}}_{|j} + \mathcal {O}(\xi )^{i}_{j} . \end{aligned}$$

(91)

The divergence of a vector field is the trace of its covariant derivative reads

$$\begin{aligned} \nabla \cdot \widetilde{\mathbf {p}}&= {{\widetilde{\text {p}}}^{I}}_{;I} = {{\widetilde{\text {p}}}^{i}}_{;i} = {{p}^{i}}_{|i} + \mathcal {O}(\xi )= \mathrm{div}\,\mathbf {p}+ \mathcal {O}(\xi ). \end{aligned}$$

(92)

The covariant curl of a vector field can be obtained by a double contraction of the Levi–Civita tensor and the contravariant derivative, i.e.,

$$\begin{aligned} \left[ \nabla \times \widetilde{\mathbf {p}} \right] _{I} = -\widetilde{E}_{IJK}\widetilde{\text {p}}^{J;K} . \end{aligned}$$

(93)

With (86), the $\xi$-component of the curl can be approximated by

$$\begin{aligned} \left[ \nabla \times \widetilde{\mathbf {p}} \right] _{\xi }&= -E_{jk}\widetilde{g}^{kL} {{\widetilde{\text {p}}}^{j}}_{;L} + \mathcal {O}(\xi )= -E_{jk}g^{kl} {{p}^{j}}_{|l} + \mathcal {O}(\xi )= \mathrm{rot}\,\mathbf {p}+ \mathcal {O}(\xi )\end{aligned}$$

(94)

and the covariant tangential components by

$$\begin{aligned} \begin{aligned} \left[ \nabla \times \widetilde{\mathbf {p}} \right] _{i}&= -\left( \widetilde{E}_{ij\xi }\widetilde{\text {p}}^{j;\xi } + \widetilde{E}_{i\xi j}\widetilde{\text {p}}^{\xi ; j}\right) = E_{ij}\widetilde{g}^{jK}{{\widetilde{\text {p}}}^{\xi }}_{;K} + \mathcal {O}(\xi )_{i} \\&= E_{ij}{{\mathcal {B}}^{j}}_{l}p^{l} + \mathcal {O}(\xi )_{i} = -[*(\mathcal {B}\mathbf {p})^{\flat }]_{i} + \mathcal {O}(\xi )_{i} , \end{aligned} \end{aligned}$$

(95)

where we use, that for a every $\mathbf {q}\in \mathsf {T}\mathcal {S}$

$$\begin{aligned} *\mathbf {q}^{\flat } = \mathbf {i}_{\mathbf {q}}(\,\text {d}{\mathcal {S}}) = \sqrt{|\mathbf {g}|}\left( -q^{\theta }d\varphi + q^{\varphi }d\theta \right) = -\mathbf {E}\mathbf {q}\end{aligned}$$

(96)

is valid on $\mathcal {S}$, see Abraham et al. (1988). The Hodge star operator is length preserving and the metric $\mathbf {\widetilde{g}}$ induces the common norm in the thin shell; therefore, it holds

$$\begin{aligned} \left\| \nabla \times \widetilde{\mathbf {p}} \right\| ^{2}_{\Omega _{\delta }}&= \left\| -*(\mathcal {B}\mathbf {p})^{\flat } \right\| ^{2}_{\mathcal {S}} + \widetilde{g}^{\xi \xi }\left( \mathrm{rot}\,\mathbf {p}\right) ^{2} + \mathcal {O}(\xi )= \left\| \mathcal {B}\mathbf {p}\right\| ^{2}_{\mathcal {S}} + \left( \mathrm{rot}\,\mathbf {p}\right) ^{2} + \mathcal {O}(\xi ). \end{aligned}$$

Finally, with $\left\| \widetilde{\mathbf {p}} \right\| ^{2}_{\Omega _{\delta }} = \left\| \mathbf {p}\right\| ^{2}_{\mathcal {S}} + \mathcal {O}(\xi )$, (84), (92), (94), and (95), we can approximate the penalized Frank–Oseen energy (3) in the thin shell $\Omega _{\delta }$ by

$$\begin{aligned} {F}_\mathrm {{\omega _\mathrm{n}}}\left[ \widetilde{\mathbf {p}}, \Omega _{\delta } \right]= & {} \int _{\mathcal {S}} \int _{-\delta /2}^{\delta /2} \frac{{K}}{2} \left( \left( \mathrm{div}\,\mathbf {p}\right) ^{2} + \left( \mathrm{rot}\,\mathbf {p}\right) ^{2} \right. \\&\left. + \,\left\| \mathcal {B}\mathbf {p}\right\| ^{2}_{\mathcal {S}} \right) + \frac{{\omega _\mathrm{n}}}{4}\left( \left\| \mathbf {p}\right\| ^{2}_{\mathcal {S}} - 1 \right) ^{2} + \mathcal {O}(\xi )\; d\xi \wedge \,\text {d}{\mathcal {S}}\\= & {} \,\delta \left( {F}_\mathrm {\omega _\mathrm{n}}^\mathcal {S}[\mathbf {p}] + \mathcal {O}(\delta ) \right) \end{aligned}$$

for $\widetilde{\mathbf {p}}\in H^{\text {DR}}( \Omega _{\delta };\, \mathsf {T}\Omega _{\delta })$ and $\mathbf {p}\in H^{\text {DR}}( \mathcal {S};\, \mathsf {T}\mathcal {S})$.

Appendix B: Integral Theorems

The exterior derivative $\mathbf {d}$ is the $L^{2}$-adjoint of $(-* \mathbf {d}*)$. This allows to obtain some frequently used integral identities for the tangential vector field $\mathbf {p}=\varvec{\alpha }^{\sharp }:\mathcal {S}\rightarrow \mathsf {T}\mathcal {S}$ on a closed surface $\mathcal {S}$ and also for its $\mathbb {R}^{3}$ extension $\widehat{\mathbf {p}}:\mathcal {S}\rightarrow \mathbb {R}^3$, with $\mathbf {p}= \pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {p}}$. We get

$$\begin{aligned} -\int _{\mathcal {S}} \left\langle {\mathrm{grad}\,f , \widehat{\mathbf {p}}}\right\rangle \,\text {d}{\mathcal {S}}&= -\int _{\mathcal {S}} \left\langle {\mathrm{grad}\,f , \mathbf {p}}\right\rangle \,\text {d}{\mathcal {S}}= -\int _{\mathcal {S}} \left\langle {\mathbf {d}f , \varvec{\alpha }}\right\rangle \,\text {d}{\mathcal {S}}\\&= \int _{\mathcal {S}} f * \mathbf {d}* \varvec{\alpha }\,\text {d}{\mathcal {S}}= \int _{\mathcal {S}} f \mathrm{div}\,\mathbf {p}\,\text {d}{\mathcal {S}}\nonumber \\&= \int _{\mathcal {S}} f \mathrm{div}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {p}}) \,\text {d}{\mathcal {S}}= \int _{\mathcal {S}} f \mathrm{div}\,\widehat{\mathbf {p}}- \mathcal {H}\left( \widehat{\mathbf {p}}\cdot \varvec{\nu }\right) \,\text {d}{\mathcal {S}}\nonumber \end{aligned}$$

and

$$\begin{aligned} -\int _{\mathcal {S}} \left\langle {\mathrm{rot}\,f , \widehat{\mathbf {p}}}\right\rangle \,\text {d}{\mathcal {S}}&= -\int _{\mathcal {S}} \left\langle {\mathrm{rot}\,f , \mathbf {p}}\right\rangle \,\text {d}{\mathcal {S}}= -\int _{\mathcal {S}} \left\langle {*\mathbf {d}f , \varvec{\alpha }}\right\rangle \,\text {d}{\mathcal {S}}\\&= \int _{\mathcal {S}} \left\langle {\mathbf {d}f , *\varvec{\alpha }}\right\rangle \,\text {d}{\mathcal {S}}= -\int _{\mathcal {S}} f * \mathbf {d}**\varvec{\alpha }\,\text {d}{\mathcal {S}}= \int _{\mathcal {S}} f \mathrm{rot}\,\mathbf {p}\,\text {d}{\mathcal {S}}\nonumber \\&= \int _{\mathcal {S}} f \mathrm{rot}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {p}}) \,\text {d}{\mathcal {S}}= \int _{\mathcal {S}} f \mathrm{rot}\,\widehat{\mathbf {p}}\,\text {d}{\mathcal {S}}\nonumber . \end{aligned}$$

Note that $**\alpha = -\alpha$ and the inner product is invariant with respect to $*$, $\flat$, and $\sharp$, applied to both arguments of the product simultaneously, see Abraham et al. (1988). Hence, we obtain for the Laplace–deRham operator

$$\begin{aligned} \int _{\mathcal {S}} \left\langle {\varvec{\Delta }^{\text {dR}}\mathbf {p}, \widehat{\mathbf {q}}}\right\rangle \,\text {d}{\mathcal {S}}&= \int _{\mathcal {S}} \left\langle {\varvec{\Delta }^{\text {dR}}\mathbf {p}, \mathbf {q}}\right\rangle \,\text {d}{\mathcal {S}}= -\int _{\mathcal {S}} \left\langle {\mathrm{grad}\,\mathrm{div}\,\mathbf {p}, \mathbf {q}}\right\rangle + \left\langle {\mathrm{rot}\,\mathrm{rot}\,\mathbf {p}, \mathbf {q}}\right\rangle \,\text {d}{\mathcal {S}}\\&= \int _{\mathcal {S}} (\mathrm{div}\,\mathbf {p})(\mathrm{div}\,\mathbf {q}) + (\mathrm{rot}\,\mathbf {p})(\mathrm{rot}\,\mathbf {q}) \,\text {d}{\mathcal {S}}\nonumber \\&= \int _{\mathcal {S}} \mathrm{div}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {p}})\mathrm{div}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {q}}) + \mathrm{rot}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {p}})\mathrm{rot}\,(\pi _{\mathsf {T}\mathcal {S}}\widehat{\mathbf {q}}) \nonumber \\&= \int _{\mathcal {S}} (\mathrm{div}\,\widehat{\mathbf {p}}- \mathcal {H}\left( \widehat{\mathbf {p}}\cdot \varvec{\nu }\right) )(\mathrm{div}\,\widehat{\mathbf {q}}- \mathcal {H}\left( \widehat{\mathbf {q}}\cdot \varvec{\nu }\right) ) + (\mathrm{rot}\,\widehat{\mathbf {p}})(\mathrm{rot}\,\widehat{\mathbf {q}}) \nonumber . \end{aligned}$$

Appendix C: Convergence Study of the Laplace–deRham Approximation

To justify the approximation $\varvec{\Delta }^{\text {dR}}\mathbf {p}\approx \widehat{\varvec{\Delta }}^{\text {dR}}\widehat{\mathbf {p}}+ {\omega _\mathrm{t}}\varvec{\nu }\left( \varvec{\nu }\cdot \widehat{\mathbf {p}}\right)$, we set up a test case consisting of a vector-valued Helmholtz equation on an ellipsoidal surface $\mathcal {S}^{E}$ (major axis: 1.0, 0.5, and 1.5)

$$\begin{aligned} -\varvec{\Delta }^{\text {dR}}\mathbf {p}+ \mathbf {p}= -\varvec{\Delta }^{\text {dR}}\mathbf {p}_s + \mathbf {p}_s =: \mathbf {f} \quad \text{ on } \mathcal {S}^{E}\end{aligned}$$

(97)

with given analytical solution $\mathbf {p}_{s} = \left[ -2y, 0.5 x, 0 \right] ^\mathrm{T} \in C( \mathcal {S}^{E}, \mathsf {T}\mathcal {S}^{E})$. We solve

$$\begin{aligned} -\widehat{\varvec{\Delta }}^{\text {dR}}\widehat{\mathbf {p}}+ \widehat{\mathbf {p}}+ {\omega _\mathrm{t}}\varvec{\nu }\left( \varvec{\nu }\cdot \widehat{\mathbf {p}}\right) = \mathbf {f} \quad \text{ on } \mathcal {S}^{E}\end{aligned}$$

(98)

using sFEM on a conforming triangulation $\mathcal {S}^{E}_h$ of $\mathcal {S}^{E}$ with piecewise linear Lagrange elements $\mathbb {V}_h(\mathcal {S}^{E}_h) = \{ v_h \in C^0(\mathcal {S}^{E}_h) \; : \; v_h|_T \in \mathbb {P}^1, \, \forall \, T \in \mathcal {T}\}$ as trial and test space for all components $\hat{\text {p}}_{i}$. This leads to a sequence of linear discrete equations

$$\begin{aligned}&\int _{\mathcal {S}^{E}_h} \nabla _{\mathcal {S}}\cdot \widehat{\mathbf {p}}D_i \psi + \nabla _{\mathcal {S}}\cdot \left( \widehat{\mathbf {p}}\times \varvec{\nu }\right) \nabla _{\mathcal {S}}\cdot \left( \mathbf {e}_i \psi \times \varvec{\nu }\right) \,\text {d}{\mathcal {S}}\nonumber \\&\quad +\int _{\mathcal {S}^{E}_h}\hat{\text {p}}_{i} \psi \,\text {d}{\mathcal {S}}+ {\omega _\mathrm{t}}\int _{\mathcal {S}^{E}_h} \nu _{i}\left( \varvec{\nu }\cdot \widehat{\mathbf {p}}\right) \psi \,\text {d}{\mathcal {S}}= \int _{\mathcal {S}^{E}_h} f_i \psi \,\text {d}{\mathcal {S}}. \end{aligned}$$

(99)

To assemble and solve the resulting system, we use the FEM toolbox AMDiS (Vey and Voigt 2007; Witkowski et al. 2015).

Figure 9 shows the $L^2$-error $\epsilon _{L^2}(\mathbf {p}) = \left( \int _{\mathcal {S}^{E}} \sum _{i=1} (\widehat{\mathbf {p}}_i - \text {p}_{s,i})^2 \,\text {d}{\mathcal {S}}\right) ^{1/2}$ vs ${\omega _\mathrm{t}}$ and linear convergence, which is only limited by the mesh quality.

Fig. 9

$L^2$-error for $\widehat{\varvec{\Delta }}^{\text {dR}}$ approximation (solid lines) for two well-centered triangulations of $\mathcal {S}^{E}$ with 25k and 100k vertices. The black dashed line indicates linear rate of convergence. The dash doted line shows the result for a componentwise approximation of $\varvec{\Delta }^{\text {dR}}$ in (100) (Color figure online)

As a complementary result and to emphasize the delicate nature of the coupling between curvature and spatial derivatives, we also show in Fig. 9 the $L^2$-error of a componentwise approximation of $\varvec{\Delta }^{\text {dR}}$

$$\begin{aligned} \varvec{\Delta }^{\text {dR}}\mathbf {p}\approx \sum _{i=i}^3 \nabla _{\mathcal {S}}\cdot \nabla _{\mathcal {S}}\hat{\text {p}}_{i} \mathbf {e}_i + {\omega _\mathrm{t}}\varvec{\nu }\left( \varvec{\nu }\cdot \widehat{\mathbf {p}}\right) . \end{aligned}$$

(100)

As clearly visible in Fig. 9, this approximation fails for any values of ${\omega _\mathrm{t}}$ to reproduce the $\varvec{\Delta }^{\text {dR}}$ behavior on $\mathcal {S}^{E}$.

Appendix D: DEC—Notations and Details

1.1 D.1: Notations

We often use the strict order relation $\succ$ and $\prec$ on simplices, where $\succ$ is proverbial the “contains” relation, i.e., $e \succ v$ means: the edge e contains the vertex v. Correspondingly is $\prec$ the “part of” relation, i.e., $v \prec T$ means: the vertex v is part of the face $T$. Hence, we can use this notation also for sums, like $\sum _{T\succ e}$, i.e., the sum over all faces $T$ containing edge e, or $\sum _{v\prec e}$, i.e., the sum over all vertices v being part of edge e. Sometimes we need to determine this relation for edges more precisely with respect to the orientation. Therefore, a sign function is introduced,

$$\begin{aligned} s_{T,e}&:= {\left\{ \begin{array}{ll} +1 &{}\quad \text {if } e\prec T\text { and } T\text { is on left side of } e \\ -1 &{}\quad \text {if } e\prec T\text { and } T\text { is on right side of } e \\ 0 &{}\quad e \nprec f, \end{array}\right. }\end{aligned}$$

(101)

$$\begin{aligned} s_{v,e}&:= {\left\{ \begin{array}{ll} +1 &{}\quad \text {if } v\prec e \text { and } e \text { points to } v \\ -1 &{}\quad \text {if } v\prec e \text { and } e \text { points away from } v \\ 0 &{}\quad v \nprec e, \end{array}\right. } \end{aligned}$$

(102)

to describe such relations between faces and edges, or vertices and edges, respectively. Figure 10 gives a schematic picture.

Fig. 10

Left this simple example mesh leads to $s_{T_{1},e} = +1$, $s_{T_{2},e} = -1$, $s_{v_{1},e} = -1$ and $s_{v_{2},e} = +1$. Right the vertex v (green) and its Voronoi cell $\star v$ (semitransparent green); the edge e (blue) and its Voronoi edge $\star e$ (blue); the face $T$ (semitransparent red) and its Voronoi vertex (red) (Color figure online)

The property of primal mesh to be well centered ensures the existence of a Voronoi mesh (dual mesh), which is also an orientable manifold-like simplicial complex, but not well centered.

The basis of the Voronoi mesh is not simplices, but chains of them. To identify these basic chains, we apply the (geometrical) star operator $\star$ on the primal simplices, i.e., $\star v$ is the Voronoi cell corresponding to the vertex v and inherits its orientation from the orientation of the polytope $\left| \mathcal {K}\right|$. $\star v$ is, from a geometric point of view, the convex hull of circumcenters $c(T)$ of all triangles $T\succ v$. The Voronoi edge $\star e$ of an edge e is a connection of the right face $T_{2}\succ e$ with the left face $T_{1}\succ e$ over the midpoint c(e) . The Voronoi vertex $\star T$ of a face $T$ is simply its circumcenter $c(T)$ (see Fig. 10). For greater details and a more mathematical discussion, see, e.g., Hirani (2003), VanderZee et al. (2010).

The boundary operator $\partial$ maps simplices (or chains of them) to the chain of simplices that describes its boundary, with respect to its orientation (see Hirani 2003), e.g., $\partial (\star v)=-\sum _{e\succ v} s_{v,e}(\star e)$ (formal sum for chains) and $\partial e = \sum _{v\prec e} s_{v,e} v$.

The expression $\left| \cdot \right|$ measures the volume of a simplex, i.e., $\left| T\right|$ the area of the face $T$, $\left| e \right|$ the length of the edge e, and the zero-dimensional volume $\left| v \right|$ are set to be 1. Therefore, the volume is also defined for chains and the dual mesh, since the integral is a linear functional.

1.2 D.2: Laplace Operators

With the Stokes theorem and the discrete Hodge operator defined in Hirani (2003), we can develop a DEC discretized Rot-Rot-Laplace for a discrete 1-form $\alpha \in \varLambda _{h}^{1}(\mathcal {K})$ by

$$\begin{aligned} \varvec{\Delta }^{\text {RR}}_{h}\alpha (e):= & {} \left( *\mathbf {d}*\mathbf {d}\alpha \right) (e) = -\frac{\left| e \right| }{\left| \star e \right| } \left( \mathbf {d}* \mathbf {d}\alpha \right) (\star e) \nonumber \\= & {} -\frac{\left| e \right| }{\left| \star e \right| } \left( * \mathbf {d}\alpha \right) (\partial \star e) = -\frac{\left| e \right| }{\left| \star e \right| } \sum _{T\succ e} s_{T,e} \left( * \mathbf {d}\alpha \right) (\star T)\\= & {} -\frac{\left| e \right| }{\left| \star e \right| } \sum _{T\succ e} \frac{s_{T,e}}{\left| T\right| } \left( \mathbf {d}\alpha \right) (T) = -\frac{\left| e \right| }{\left| \star e \right| } \sum _{T\succ e} \frac{s_{T,e}}{\left| T\right| }\alpha (\partial T)\nonumber \\= & {} -\frac{\left| e \right| }{\left| \star e \right| } \sum _{T\succ e} \frac{s_{T,e}}{\left| T\right| } \sum _{\tilde{e}\prec T} s_{T,\tilde{e}} \alpha (\tilde{e})\nonumber \end{aligned}$$

(103)

and a DEC discretized Grad-Div-Laplace by

$$\begin{aligned} \begin{aligned} \varvec{\Delta }^{\text {GD}}_{h}\alpha (e)&:= \left( \mathbf {d}* \mathbf {d}* \alpha \right) (e) = \left( * \mathbf {d}*\alpha \right) (\partial e) \\&= \sum _{v\prec e} s_{v,e} \left( * \mathbf {d}* \alpha \right) (v) = \sum _{v\prec e} \frac{s_{v,e}}{\left| \star v \right| } \left( \mathbf {d}* \alpha \right) (\star v)\\&= \sum _{v\prec e} \frac{s_{v,e}}{\left| \star v \right| } \left( * \alpha \right) (\partial \star v) = -\sum _{v\prec e} \frac{s_{v,e}}{\left| \star v \right| } \sum _{\tilde{e}\succ v} s_{v,\tilde{e}} \left( * \alpha \right) (\star \tilde{e}) \\&= -\sum _{v\prec e} \frac{s_{v,e}}{\left| \star v \right| } \sum _{\tilde{e}\succ v} s_{v,\tilde{e}} \frac{\left| \star \tilde{e} \right| }{\left| \tilde{e} \right| }\alpha (\tilde{e}) . \end{aligned} \end{aligned}$$

(104)

Hence, we obtain the DEC discretized Laplace–deRham operator by

$$\begin{aligned} \varvec{\Delta }^{\text {dR}}_h\alpha (e) = -\varvec{\Delta }^{\text {RR}}_h\alpha (e) - \varvec{\Delta }^{\text {GD}}_h\alpha (e). \end{aligned}$$

1.3 D.3: Conflate Linear Operators and Its Hodge Dual to a PD-(1, 1)-Tensor

For a linear operator $\mathbf {M}:\mathsf {T}^{*}\mathcal {S}\rightarrow \mathsf {T}^{*}\mathcal {S}$ pointwise defined as a mixed co- and contravariant (1,1)-tensor with components ${{M}_{i}}^{j}$, we discretize the 1-form $\mathbf {M}\varvec{\alpha }$ on an edge $e\in \mathcal {E}$ by definition (27) and approximate the operator on the projected midpoint of the edge, i.e.,

$$\begin{aligned} \left( \mathbf {M}\varvec{\alpha }\right) _{h}(e)&= \int _{\pi (e)}{{M}_{i}}^{j}\alpha _{j}\ \hbox {d}x^{i} \approx \left[ \mathbf {M}(e)\right] _{ik}g^{kj}\int _{\pi (e)}\alpha _{j} \hbox {d} x^{i}, \end{aligned}$$

(105)

with $\mathbf {M}(e) := \mathbf {M}|_{\pi (c(e))}$. With respect to an orthogonal basis $\{\partial _{i}\mathbf {x},\partial _{j}\mathbf {x}\}$ with metric tensor $\mathbf {g}=g_{i}(dx^{i})^{2}$, we obtain for the 1-form $\varvec{\alpha }=\alpha _{i}dx^{i}$ the Hodge dual

$$\begin{aligned} *\varvec{\alpha }&= [*\alpha ]_{1}dx^{1} + [*\alpha ]_{2}dx^{2} = -\sqrt{\frac{g_{1}}{g_{2}}}\alpha _{2}dx^{1} + \sqrt{\frac{g_{2}}{g_{1}}}\alpha _{1}dx^{2} . \end{aligned}$$

(106)

Hence, we can replace the 1-forms beneath the integrals by

$$\begin{aligned} \begin{bmatrix} \alpha _{1}dx^{1}&\quad \alpha _{2}dx^{1} \\ \alpha _{1}dx^{2}&\quad \alpha _{2}dx^{2} \end{bmatrix}&=\begin{bmatrix} \alpha _{1}dx^{1}&\quad -\sqrt{\frac{g_{2}}{g_{1}}}[*\alpha ]_{1}dx^{1} \\ \sqrt{\frac{g_{1}}{g_{2}}}[*\alpha ]_{2}dx^{2}&\quad \alpha _{2}dx^{2} \end{bmatrix} . \end{aligned}$$

(107)

Now, we use the basis $\{\mathbf {e},\mathbf {e}_{\star }\}$ defined in Sect. 4.2 on the polytope $|\mathcal {K}|$ and the resulting metric (34), i.e., $g_{1}=|e|^{2}$ and $g_{2} = |\star e|^{2}$. This leads to an approximation of $\left( \mathbf {M}\varvec{\alpha }\right) _{h} \in \varLambda ^{1}_{h}(\mathcal {K})$ as a linear combination of $\alpha _{h}, (*\alpha )_{h} \in \varLambda ^{1}_{h}(\mathcal {K})$, or rather, evaluated on an edge $e\in \mathcal {E}$

$$\begin{aligned} \left( \mathbf {M}\varvec{\alpha }\right) _{h}(e)&\approx \frac{1}{|e|^{2}}M_{\mathbf {e},\mathbf {e}}(e)\alpha _{h}(e) - \frac{1}{|e||\star e|}M_{\mathbf {e},\mathbf {e}_{\star }}(*\alpha )_{h}(e) \end{aligned}$$

(108)

and, in general, for $\mathbf {v},\mathbf {w}\in \mathrm{Span}\,\left\{ \mathbf {e},\mathbf {e}_{\star } \right\}$ is $M_{\mathbf {v},\mathbf {w}}(e) = \mathbf {v}\cdot \mathbf {M}(e)\cdot \mathbf {w} = v^{i}\left[ \mathbf {M}(e)\right] _{ij}w^{j}$ the evaluation of the complete covariant tensor $\mathbf {M}(e)$ in direction $\mathbf {v}$ and $\mathbf {w}$. Note, if $\mathbf {M}\in \mathsf {T}\mathcal {S}\times \mathsf {T}\mathcal {S}$ is formulated in Euclidean $\mathbb {R}^{3}$ coordinates, so that $\mathbf {M}(e)\in \mathbb {R}^{3 \times 3}$, there is no distinction between co- and contravariant components of $\mathbf {M}(e)$. Furthermore, if we use the approximation $\left( *\mathbf {M}\varvec{\alpha }\right) _{h}(e) \approx -\frac{|e|}{|\star e|}\left( \mathbf {M}\varvec{\alpha }\right) _{h}(\star e)$, we get with respect to (105) and (107)

$$\begin{aligned} \left( *\mathbf {M}\varvec{\alpha }\right) _{h}(e)&\approx - \frac{1}{|e||\star e|}M_{\mathbf {e}_{\star },\mathbf {e}}\alpha _{h}(e) + \frac{1}{|\star e|^{2}}M_{\mathbf {e}_{\star },\mathbf {e}_{\star }}(*\alpha )_{h}(e) . \end{aligned}$$

(109)

Finally, we can summarize (108) and (109) with the PD-1-form $\underline{\varvec{\alpha }}\in \varLambda ^{1}_{h}(\mathcal {K};\mathfrak {T}^{*}\mathcal {E})$ on every edge $e\in \mathcal {E}$ to

$$\begin{aligned} \underline{\underline{\mathbf {M}}}\cdot \underline{\varvec{\alpha }}&:= \begin{bmatrix} \frac{1}{|e|^{2}}M_{\mathbf {e},\mathbf {e}}&\quad - \frac{1}{|e||\star e|}M_{\mathbf {e},\mathbf {e}_{\star }} \\ - \frac{1}{|e||\star e|}M_{\mathbf {e}_{\star },\mathbf {e}}&\quad \frac{1}{|\star e|^{2}}M_{\mathbf {e}_{\star },\mathbf {e}_{\star }} \end{bmatrix} \cdot \underline{\varvec{\alpha }}\approx \begin{bmatrix} \left( \mathbf {M}\varvec{\alpha }\right) _{h} \\ \left( *\mathbf {M}\varvec{\alpha }\right) _{h} \end{bmatrix} , \end{aligned}$$

(110)

where the evaluation argument e is omitted for a better readability.