High-Order Approximation of Gaussian Curvature with Regge Finite Elements

TL;DR

This paper introduces a high-order method for approximating Gaussian curvature on triangulated surfaces using Regge finite elements, providing error estimates and a generalization of classical angle defect techniques.

Contribution

It develops a novel high-order approximation framework for Gaussian curvature based on Regge finite elements, extending classical methods and providing rigorous error analysis.

Findings

Error estimates in Sobolev norms for high-order approximations

Relation between angle defect linearization and div-div operator discretization

Generalization of angle defect to higher order methods

Abstract

A widely used approximation of the Gaussian curvature on a triangulated surface is the angle defect, which measures the deviation between $2\pi$ and the sum of the angles between neighboring edges emanating from a common vertex. We show that the linearization of the angle defect about an arbitrary piecewise constant Regge metric is related to the classical Hellan-Herrmann-Johnson finite element discretization of the div-div operator. Integrating this relation leads to an integral formula for the angle defect which is well-suited for analysis and generalizes naturally to higher order. We prove error estimates for these high-order approximations of the Gaussian curvature in $H^k$-Sobolev norms of integer order $k \ge -1$.

Tables

Table 1. Table 1: Errors ‖ κ h ​ ( g h ) − κ ​ ( g ) ‖ L 2 ​ ( Ω ) subscript norm subscript 𝜅 ℎ subscript 𝑔 ℎ 𝜅 𝑔 superscript 𝐿 2 Ω \|\kappa_{h}(g_{h})-\kappa(g)\|_{L^{2}(\Omega)} for the metric ( 1 ) on the square ( − 1 , 1 ) × ( − 1 , 1 ) 1 1 1 1 (-1,1)\times(-1,1) .

		$r=0$				$r=1$				$r=2$
$h$		Error		Order		Error		Order		Error		Order
$1$	$.654\cdot {10}^{-1}$	$1$	$.535\cdot {10}^{-1}$			$4$	$.479\cdot {10}^{-2}$			$2$	$.887\cdot {10}^{-3}$
$8$	$.606\cdot {10}^{-2}$	$1$	$.812\cdot {10}^{-1}$	$-0$	$.25$	$5$	$.236\cdot {10}^{-2}$	$-0$	$.24$	$9$	$.119\cdot {10}^{-4}$	$1$	$.76$
$4$	$.442\cdot {10}^{-2}$	$1$	$.787\cdot {10}^{-1}$	$0$	$.02$	$5$	$.152\cdot {10}^{-2}$	$0$	$.02$	$3$	$.523\cdot {10}^{-4}$	$1$	$.44$
$2$	$.260\cdot {10}^{-2}$	$1$	$.780\cdot {10}^{-1}$	$0$	$.01$	$5$	$.165\cdot {10}^{-2}$	$0$	$.00$	$1$	$.592\cdot {10}^{-4}$	$1$	$.17$

Equations

$$g(x,y)=\left(\begin{array}{cc}1+{\left(\frac{\partial f}{\partial x}\right)}^2& \frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\\ \frac{\partial f}{\partial x}\frac{\partial f}{\partial y}& 1+{\left(\frac{\partial f}{\partial y}\right)}^2\end{array}\right),$$

$$\kappa (g)(x,y)=\frac{81(1-x^2)(1-y^2)}{{\left(9+x^2(x^2-3{)}^2+y^2(y^2-3{)}^2\right)}^2}.$$

$$({\mathrm{Hess}}_{g}v{)}_{ij}=\frac{{\partial }^{2}v}{\partial x^i\partial x^j}-{\Gamma }_{ij}^k\frac{\partial v}{\partial x^k}.$$

$${\Gamma }_{ij}^k=\frac{1}{2}{g}^{k\ell }\left(\frac{\partial g_{\ell i}}{\partial x^j}+\frac{\partial g_{\ell j}}{\partial x^i}-\frac{\partial g_{ij}}{\partial x^{\ell }}\right),$$

$${\Delta }_{g}v={\mathrm{div}}_g{\nabla }_{g}v.$$

$$({\mathrm{div}}_g\sigma {)}^i=\frac{\partial {\sigma }^{ij}}{\partial x^j}+{\Gamma }_{jk}^i{\sigma }^{kj}+{\Gamma }_{jk}^j{\sigma }^{ik},$$

$$\kappa (g)=\frac{1}{2}R(g)=\frac{1}{2}{g}^{ij}\left(\frac{\partial {\Gamma }_{ij}^k}{\partial x^k}-\frac{\partial {\Gamma }_{ik}^k}{\partial x^j}+{\Gamma }_{ij}^{\ell }{\Gamma }_{k\ell }^k-{\Gamma }_{ik}^{\ell }{\Gamma }_{j\ell }^k\right).$$

$$\mathrm{Ric}(g)=\kappa (g)g.$$

$$D(\kappa \mu )(g)\cdot \sigma =\frac{1}{2}({\mathrm{div}}_g{\mathrm{div}}_g{S}_g\sigma )\mu (g),$$

$$S_g\sigma =\sigma -g\mathrm{Tr}(g^{-1}\sigma ).$$

$$D\kappa (g)\cdot \sigma$$

$$D\mu (g)\cdot \sigma$$

$$⟨u,v{⟩}_g={\int }_{\Omega }uv\mu (g)$$

$${\frac{d}{dt}|}_{t=0}⟨\kappa (g+t\sigma ),v{⟩}_{g+t\sigma }=\frac{1}{2}⟨{\mathrm{div}}_g{\mathrm{div}}_g{S}_g\sigma ,v{⟩}_g.$$

$$⟨\kappa (g),v{⟩}_g=\frac{1}{2}{\int }_0^{1}b((1-t)\delta +tg;g-\delta ,v)dt,$$

$$b(g;\sigma ,v)=⟨{\mathrm{div}}_g{\mathrm{div}}_g{S}_g\sigma ,v{⟩}_g.$$

$$\underset{K\in {\mathcal{T}}_h}{\max }\frac{h_K}{{\rho }_K}$$

$$\underset{K\in {\mathcal{T}}_h}{\max }\frac{h}{h_K}$$

$${\tau }_g=\frac{1}{\sqrt{{\tau }^{T}g\tau }}\tau ,n_g=\frac{Jg\tau }{\sqrt{{\tau }^{T}g\tau }\sqrt{\det g}}.$$

$$n_g^T(S_g\sigma )n_g=-{\tau }_g^T\sigma {\tau }_g,$$

$$\left\llbracket v\right\rrbracket=\left.v\right|_{\partial K_{1}}+\left.v\right|_{\partial K_{2}}.$$

$$\left\llbracket v\right\rrbracket=\left.v\right|_{\partial\Omega}.$$

$$V=\{v\in H_0^1(\Omega )\mid {v|}_K\in H^2(K),\forall K\in {\mathcal{T}}_h\}$$

$$\begin{array}{r}\Sigma =\{\sigma \in L^2(\Omega )\otimes \mathbb{S}\mid {\sigma |}_K\in H^1(K)\otimes \mathbb{S},\forall K\in {\mathcal{T}}_h,\text{ and }\\ {\tau }^T\sigma \tau \text{ is continuous across every }e\in {\mathring{\mathcal{E}}}_h\}.\end{array}$$

$$b_h(g;\sigma ,v)$$

$$⟨u,v{⟩}_{g,e}$$

$$⟨\sigma ,\rho {⟩}_{g,K}$$

$$\frac{\partial v}{\partial n_g}=n_g^{T}g{\nabla }_{g}v.$$

$$V_h=\{v\in V\mid {v|}_K\in {\mathcal{P}}_q(K),\forall K\in {\mathcal{T}}_h\}$$

$${\Sigma }_h=\{\sigma \in \Sigma \mid {\sigma |}_K\in {\mathcal{P}}_r(K)\otimes \mathbb{S},\forall K\in {\mathcal{T}}_h\},$$

$${\Sigma }_{h+}=\{\sigma \in {\Sigma }_h\mid \sigma (x)\succ 0,\forall x\in \mathring{K},\forall K\in {\mathcal{T}}_h\}.$$

$$⟨{\kappa }_h(g_h),v_h{⟩}_{g_h}=\frac{1}{2}{\int }_0^1{b}_h((1-t)\delta +t{g}_h;g_h-\delta ,v_h)dt,\forall v_h\in V_h.$$

Full text

{AMS}

65N30, 65N15, 65D18, 83C27

1 Introduction

One of the most widely used approximations of the Gaussian curvature on a triangulated surface is the angle defect: $2\pi$ minus the sum of the angles between neighboring edges emanating from a common vertex. This approximation (and its generalization to higher dimensions) is used in several applications, including discrete analogues of Ricci flow [10, 21], discrete theories of relativity [25, 1], discrete differential geometry [19, 28], and computer graphics algorithms [22, 13, 27]. Despite its widespread use, the angle defect leaves much to be desired if one is interested in accurately approximating the curvature of a smooth surface (or smooth Riemannian manifold) with a discretization thereof. It is manifestly a low-order approximation of the curvature, relying in essence on piecewise constant approximations of the underlying smooth metric tensor.

In this paper, we introduce and analyze a family of high-order approximations of the Gaussian curvature using piecewise polynomial approximations of the metric tensor. The cornerstone of our construction is an integral formula for the angle defect that mimics a certain integral formula for the Gaussian curvature which is valid in the smooth setting. In the discrete setting, the integral formula follows from the observation that the linearization of the angle defect about an arbitrary piecewise constant metric (more precisely, a piecewise constant Regge metric) is related to the classical Hellan-Herrmann-Johnson finite element discretization of the $\operatorname{div}\operatorname{div}$ operator. This observation generalizes one made by Christiansen [11], who derived the linearization of the angle defect about the Euclidean metric and related it to the jumps in the tangential-normal components of the metric perturbation (see Remark 3.5 for more details).

To generalize the angle defect to higher order, we rely on the Regge finite element spaces recently developed by Li [24], which have their origins in the work of Christiansen [11]. These finite element spaces consist of piecewise polynomial $(0,2)$-tensor fields with continuous tangential-tangential components across element interfaces. In the lowest order setting, a positive definite Regge finite element realizes a piecewise flat triangulation whose (squared) edge lengths correspond to the degrees of freedom for the finite element space.

The advantages of high-order approximation of the Gaussian curvature are easily illustrated with an example. Consider the square $\Omega=(-1,1)\times(-1,1)$ equipped with the Riemannian metric

[TABLE]

where $f(x,y)=\frac{1}{2}x^{2}-\frac{1}{12}x^{4}+\frac{1}{2}y^{2}-\frac{1}{12}y^{4}$. This is nothing more than the induced metric for the surface $z=f(x,y)$ in $\mathbb{R}^{3}$. The exact Gaussian curvature of $g$ is

[TABLE]

Figure 1 plots the approximate Gaussian curvature of $g$, as computed using Definition 3.1 below with piecewise polynomial approximations of $g$ of degree $r=0,1,2$ on the triangulation of $\Omega$ depicted in Figure 1a. As we will show later in this paper, the approximate curvature produced by Definition 3.1 in the case $r=0$ is precisely the angle defect, normalized by the consistent mass matrix for piecewise linear finite elements. (In the interest of fairness, we normalized by a lumped mass matrix to produce Figure 1b; the appearance of Figure 1b worsens if the consistent mass matrix is used.) Notice that the results for $r=0$ and $r=1$ are not particularly satisfactory. This should come as no surprise; one expects the second derivatives of a degree-$r$ polynomial approximation of $g$ to converge in $W^{r-2,p}$-Sobolev norms under mesh refinement, but not in stronger norms. A chief goal of this paper is to verify this intuition with an error estimate.

Our error analysis complements a number of related results in the literature on scalar curvature approximation. Cheeger, Müller, and Schrader [8] prove that the angle defect converges in the sense of measures to the (densitized) scalar curvature if a smooth Riemannian manifold (not necessarily of dimension 2) is approximated with a suitable sequence of triangulations. Christiansen [12] proves a dual result: if a piecewise constant metric is approximated with a sequence of mollifications thereof, then the exact (densitized) scalar curvature of the mollified metric converges in the sense of measures to the angle defect. Other analyses of the angle defect appear in [5, 29] and the references therein; many of these analyses are guided by Taylor expansions on parametrized surfaces and impose special conditions on the triangulation. Our error analysis will focus not on the angle defect but instead on its higher order generalizations.

There appear to be relatively few studies on Gaussian curvature approximation that offer quantitative error bounds in Sobolev norms like the bounds in Theorem 4.1 below. Fritz [17] has proven bounds of this type for a curvature approximation which uses isoparametric approximations of surfaces embedded in $\mathbb{R}^{3}$. His results apply more generally to Ricci tensor approximation on hypersurfaces, and have been used to discretize Ricci flow [18]. However, they are inapplicable if the manifold does not admit a codimension-one embedding into Euclidean space. Note that studies related to mean curvature approximation are more widespread; see, for instance, [20, 4, 14, 23].

This paper is organized as follows. In Section 2, we introduce our notation and point out an integral formula for the Gaussian curvature. In Section 3, we discretize this integral formula with the aid of Regge finite elements, and we show that it reduces to the angle defect in the lowest-order setting. We state and prove error estimates for the aforementioned discretization in Section 4. We conclude with numerical examples in Section 5.

For simplicity, we perform much of the forthcoming analysis on a triangulated polygonal domain in $\mathbb{R}^{2}$ equipped with a non-Euclidean metric. Note, however, that the curvature approximations we introduce are coordinate-free and can be readily applied to two-dimensional orientable simplicial complexes with more general topology. We refer the reader to [24] for the appropriate definitions of the Regge finite elements in this general setting. Importantly, our curvature approximations do not rely on (nor assume the existence of) an embedding of the manifold under consideration in $\mathbb{R}^{3}$.

2 Curvature in the Smooth Setting

2.1 Preliminaries

Let $\Omega\subset\mathbb{R}^{2}$ be a polygonal domain. We use standard notation for the Sobolev spaces $W^{s,p}(\Omega)$ of differentiability index $s\in\mathbb{R}$ and integrability index $p\in[1,\infty]$. We denote $H^{s}(\Omega)=W^{s,2}(\Omega)$ and $H^{1}_{0}(\Omega)=\{v\in H^{1}(\Omega)\mid\left.v\right|_{\partial\Omega}=0\}$.

Let $g$ be a smooth Riemannian metric on $\Omega$. Fix coordinates $(x^{1},x^{2})$ on $\mathbb{R}^{2}$ so that $g$ may be regarded as a map from $\Omega$ to $\mathbb{S}=\{\sigma\in\mathbb{R}^{2\times 2}\mid\sigma=\sigma^{T}\}$. If $v$ is a scalar field on $\Omega$, then we denote $\nabla v=\left(\frac{\partial v}{\partial x^{1}},\frac{\partial v}{\partial x^{2}}\right)^{T}$ and $\nabla_{g}v=g^{-1}\nabla v$. If $w$ is a vector field on $\Omega$, then we regard it as a column vector and denote $\operatorname{div}w=\frac{\partial w^{1}}{\partial x^{1}}+\frac{\partial w^{2}}{\partial x^{2}}$ and $\operatorname{div}_{g}w=\frac{1}{\sqrt{\det g}}\operatorname{div}(\sqrt{\det g}\,w)$.

The Riemannian Hessian of a scalar field $v$ is denoted $\operatorname{Hess}_{g}v$. In coordinates,

[TABLE]

Here, the Einstein summation convention is adopted, and $\Gamma^{k}_{ij}$ are the Christoffel symbols of the second kind. That is,

[TABLE]

where $g_{ij}$ denotes the $(i,j)$-component of $g$, and $g^{ij}$ denotes the $(i,j)$-component of $g^{-1}$. The Laplacian of $v$ is

[TABLE]

If $\sigma$ is a symmetric $(0,2)$-tensor field on $\Omega$, then we regard it as a map from $\Omega$ to $\mathbb{S}$ and denote its components by $\sigma_{ij}$. In a slight abuse of notation, we denote by $\operatorname{div}_{g}\sigma$ the vector field with components

[TABLE]

where $\sigma^{ij}=g^{ik}\sigma_{k\ell}g^{\ell j}$. We omit the subscript $g$ if $g$ is the Euclidean metric $\delta\equiv\left(\begin{smallmatrix}1&0\\ 0&1\end{smallmatrix}\right)$. Thus, $\operatorname{div}\sigma=\operatorname{div}_{\delta}\sigma$, $\Delta v=\Delta_{\delta}v$, etc.

The Gaussian curvature of $g$ is denoted $\kappa(g)$; it is half the scalar curvature $R(g)$:

[TABLE]

Because $\Omega$ is two-dimensional, the Ricci tensor $\operatorname{Ric}(g)$ is proportional to the metric:

[TABLE]

2.2 Linearization of the Curvature

Let $\mu(g)=\sqrt{\det g}\,dx$ be the volume form on $\Omega$ induced by $g$. The linearization of $(\kappa\mu)(g)=\kappa(g)\mu(g)$ about $g$ will be of fundamental importance in this paper. For any smooth $\sigma:\Omega\rightarrow\mathbb{S}$, we have

[TABLE]

where

[TABLE]

This can be seen by combining the well-known relations ([16, Lemma 2], [9, Equation 2.11])

[TABLE]

with (2), noting that $\operatorname{div}_{g}\operatorname{div}_{g}(gv)=\Delta_{g}v$ for any scalar field $v$.

Another way of writing (3) is as follows. Let

[TABLE]

be the $L^{2}$-inner product on $\Omega$ induced by $g$. Then, for any scalar function $v$,

[TABLE]

Since the Euclidean metric $\delta$ has zero curvature, integrating the above relation leads to the integral formula

[TABLE]

where

[TABLE]

Our strategy for discretizing $\kappa(g)$ will be to discretize the integral formula (6).

3 Discretization

Let $\{\mathcal{T}_{h}\}_{h>0}$ be a shape-regular, quasi-uniform family of triangulations of $\Omega$ parametrized by $h=\max_{K\in\mathcal{T}_{h}}h_{K}$, where $h_{K}=\operatorname{diam}(K)$ denotes the diameter of a triangle $K$. In other words, there are constants $C_{1}$ and $C_{2}$ such that for every $h$,

[TABLE]

where $\rho_{K}$ denotes the inradius of $K$.

Let $e\subset\partial K$ be an edge of a triangle $K\in\mathcal{T}_{h}$. The outward unit normal vector to $K$ along $e$ relative to the Euclidean metric is denoted $n$, and the unit tangent vector relative to the Euclidean metric is $\tau=-Jn$, where $J=\left(\begin{smallmatrix}0&1\\ -1&0\end{smallmatrix}\right)$. Relative to $g$, the unit tangent and normal vectors are

[TABLE]

One checks that $\tau_{g}^{T}g\tau_{g}=n_{g}^{T}gn_{g}=1$ and $\tau_{g}^{T}gn_{g}=0$ since $gJg=(\det g)J$. We also note that if $\sigma:\Omega\rightarrow\mathbb{S}$, then

[TABLE]

owing to the definition (4) of $S_{g}$ and the identity $\tau_{g}\tau_{g}^{T}+n_{g}n_{g}^{T}=g^{-1}$.

Let $\mathcal{E}_{h}$ denote the set of edges of triangles in $\mathcal{T}_{h}$, and let $\mathring{\mathcal{E}}_{h}$ denote the set of interior edges; these are the edges $e\in\mathcal{E}_{h}$ with $e\not\subset\partial\Omega$. Let $v$ be a scalar field. Along any edge $e=K_{1}\cap K_{2}\in\mathring{\mathcal{E}}_{h}$, we denote

[TABLE]

If $e\in\mathcal{E}_{h}$ is on the boundary of $\Omega$, we denote

[TABLE]

Let

[TABLE]

and

[TABLE]

Note that $V$ and $\Sigma$ depend on $h$, but we have omitted a subscript $h$ to emphasize that they are infinite-dimensional spaces. We define a metric-dependent bilinear form $b_{h}(g;\cdot,\cdot):\Sigma\times V\rightarrow\mathbb{R}$ by

[TABLE]

where

[TABLE]

and

[TABLE]

The bilinear form $b_{h}$ is well-studied. In view of (9), it is (up to the appearance of $S_{g}$) a non-Euclidean generalization of the bilinear form used to discretize the $\operatorname{div}\operatorname{div}$ operator in the classical Hellan-Herrmann-Johnson mixed discretization of the biharmonic equation [3, 2, 7]. See [24, Section 4.2] for more insight into this connection in the Euclidean setting.

Now let $q\in\mathbb{N}$ and $r\in\mathbb{N}_{0}$, and define finite-dimensional subspaces

[TABLE]

and

[TABLE]

where $\mathcal{P}_{r}(K)$ denotes the space of polynomials of degree $\leq r$ on $K$. The space $\Sigma_{h}$ is the space of Regge finite elements of degree $r$ [24, 11].

Of particular importance to us will be the space of positive definite Regge finite elements,

[TABLE]

We define the discrete Gaussian curvature of a metric $g_{h}\in\Sigma_{h+}$ as follows.

Definition 3.1**.**

Let $q\in\mathbb{N}$, $r\in\mathbb{N}_{0}$, and $g_{h}\in\Sigma_{h+}$. The discrete Gaussian curvature of $g_{h}$ is the unique function $\kappa_{h}(g_{h})\in V_{h}$ satisfying

[TABLE]

Note that for each $v_{h}\in V_{h}$, the value of $\langle\kappa_{h}(g_{h}),v_{h}\rangle_{g_{h}}$ depends only on the values of $g_{h}$ in $\operatorname{supp}(v_{h})$. Thus, Definition 3.1 extends readily to orientable triangulations with more general topology: for each function $v_{h}$ in the canonical basis for $V_{h}$, one computes the spatial integral in (10) over $\operatorname{supp}(v_{h})$, which (generically) is a patch of triangles admitting a Euclidean metric.

One of our main reasons for favoring this definition is that when $(r,q)=(0,1)$, $\kappa_{h}(g_{h})$ reduces to the popular angle defect. In detail, let $\{y^{(i)}\}_{i=1}^{N}\subset\Omega$ be the vertices of $\mathcal{T}_{h}$, enumerated in such a way that $y^{(i)}\notin\partial\Omega$ if and only if $1\leq i\leq N_{0}<N$. Let $\{\phi_{i}\}_{i=1}^{N_{0}}$ be the basis for $V_{h}$ defined by

[TABLE]

Theorem 3.2**.**

Let $r=0$, $q=1$, and $g_{h}\in\Sigma_{h+}$. For every $i=1,2,\dots,N_{0}$, we have

[TABLE]

*where $\omega_{i}$ is the set of triangles in $\mathcal{T}_{h}$ sharing vertex $i$, and $\theta_{iK}$ is the interior angle of $K$ at vertex $i$ as measured by $g_{h}$. *

The preceding theorem is a consequence of the following identity that mimics (5).

Lemma 3.3**.**

Let $r=0$, $q=1$, $g_{h}\in\Sigma_{h+}$, $\sigma_{h}\in\Sigma_{h}$, and $i\in\{1,2,\dots,N_{0}\}$. For each $K\in\omega_{i}$ and $t$ sufficiently small, let $\theta_{iK}(t)$ be the interior angle of $K$ at vertex $i$ as measured by $g_{h}+t\sigma_{h}$. Then

[TABLE]

Proof 3.4**.**

Let $K\in\omega_{i}$ be a triangle with edges $e_{a}$, $e_{b}$, and $e_{c}$ of length $a(t)$, $b(t)$, and $c(t)$ relative to $g_{h}+t\sigma_{h}$. Assume that $e_{c}$ is opposite the angle $\theta_{iK}(t)$. Differentiating the law of cosines

[TABLE]

with respect to $t$ at $t=0$ and solving for $\dot{\theta}_{iK}=\dot{\theta}_{iK}(0)$ gives

[TABLE]

where $A=\frac{1}{2}ab\sin\theta_{iK}$ is the area of $K$ relative to $g_{h}$ and $a=a(0)$, $\dot{a}=\dot{a}(0)$, etc. On the other hand, we have

[TABLE]

and

[TABLE]

Indeed, the relations in (12) follow from the fact that relative to $g_{h}$, $\nabla_{g_{h}}\phi_{i}$ is a vector of length $\frac{c}{2A}$ that is perpendicular to the edge $e_{c}$. The relations in (13) follow from differentiating the identity

[TABLE]

and its counterpart for the edges $e_{b}$ and $e_{c}$. Summing over all $K\in\omega_{i}$ and noting that $\phi_{i}=0$ on each $K\notin\omega_{i}$, we conclude that

[TABLE]

*Since $\operatorname{Hess}_{g_{h}}\phi_{i}$ vanishes on each $K\in\mathcal{T}_{h}$, the relation (11) follows. *

Remark 3.5**.**

*Let us clarify the distinction between the preceding lemma and the results of Christiansen [11]. Christiansen works in three dimensions and computes the first- and second-order variation of the Regge action $\sum_{e}\ell_{e}\kappa_{e}$ ($\kappa_{e}$ being the angle defect at an edge $e$, and $\ell_{e}$ being its length) around the Euclidean metric. He notes that the first variation vanishes (a fact proven by Regge [25]), while the second variation is related to the distributional $\operatorname{curl}T\operatorname{curl}$ operator (which coincides with the operator $-\operatorname{div}\operatorname{div}S_{\delta}$ in two dimensions). Along the way, he relates the linearization of $\kappa_{e}$ around the Euclidean metric to a summation of jumps of $\tau^{T}\sigma n$ (see [11, Proposition 2]). With some manipulation, this relation can be restated in the form of Theorem 3.3 with $g_{h}=\delta$. Note that the case $g_{h}\neq\delta$ is not addressed in [11]; this missing ingredient plays a crucial role in our work. *

4 Error Estimates

In this section, we prove error estimates for $\kappa_{h}(g_{h})-\kappa(g)$ in $H^{k}$-Sobolev norms of integer order $k\geq-1$.

Denote

[TABLE]

and $\|v\|_{H^{1}(\Omega,g)}^{2}=\|v\|_{L^{2}(\Omega,g)}^{2}+|v|_{H^{1}(\Omega,g)}^{2}$. Also let

[TABLE]

Note that if the eigenvalues of $g$ are bounded above and below by positive constants on $\overline{\Omega}$, then the norms $\|v\|_{L^{2}(\Omega,g)}$ and $\|v\|_{H^{1}(\Omega,g)}$ are equivalent to the norms $\|v\|_{L^{2}(\Omega)}$ and $\|v\|_{H^{1}(\Omega)}$, respectively. These facts imply that (15) is equivalent to the norm obtained by replacing $\|u\|_{H^{1}(\Omega)}$ with $\|u\|_{H^{1}(\Omega,g)}$ in the denominator of  (15).

Our analysis will make use of the broken Sobolev norms

[TABLE]

with the obvious modifications for $p=\infty$. We denote $\|\cdot\|_{H^{s}_{h}(\Omega)}=\|\cdot\|_{W^{s,2}_{h}(\Omega)}$, and we use analogous notation $|\cdot|_{W^{s,p}_{h}(\Omega)}$ and $|\cdot|_{H^{s}_{h}(\Omega)}$ for the corresponding broken Sobolev semi-norms.

For $p\in[1,\infty]$ and $s>2/p$, we denote

[TABLE]

Note that the condition $g(x)\succ 0$ is meaningful for $s>2/p$, since the Sobolev embedding theorem implies that elements of $W^{s,p}(\Omega)$ with $s>2/p$ are continuous on $\overline{\Omega}$. It is known that for $p\in(1,\infty)$ and $s>2/p+1$, the curvature operator $\kappa$ maps $\mathcal{M}^{s,p}(\Omega)$ into $W^{s-2,p}(\Omega)$ [16, Lemma 1].

Theorem 4.1**.**

Let $g\in\mathcal{M}^{2,\infty}(\Omega)$. Suppose that $\{g_{h}\in\Sigma_{h}\}_{h>0}$ is a sequence satisfying $\lim_{h\rightarrow 0}\|g_{h}-g\|_{L^{\infty}(\Omega)}=0$, $\lim_{h\rightarrow 0}h^{-1}\log h^{-1}\|g_{h}-g\|_{L^{2}(\Omega)}=0$, and $C_{0}:=\sup_{h>0}\|g_{h}\|_{W^{1,\infty}_{h}(\Omega)}<\infty$. Then there exists a constant $C$ depending on $\Omega$, $\|g\|_{W^{1,\infty}(\Omega)}$, $\|g^{-1}\|_{L^{\infty}(\Omega)}$, $\|\kappa(g)\|_{L^{2}(\Omega)}$, $r$, $q$, $C_{0}$, and the shape regularity and quasi-uniformity constants $C_{1}$ and $C_{2}$ such that for every $h$ sufficiently small,

[TABLE]

Furthermore, if $\kappa(g)\in H^{m}(\Omega)\cap H^{1}_{0}(\Omega)$ with $m\in\mathbb{N}$, then for each $k=0,1,\dots,q$ and each $\ell=k,k+1,\dots,\min\{q+1,m\}$,

[TABLE]

*for every $h$ sufficiently small. *

Note that the theorem is vacuous when $r=0$, since we generally cannot expect $\|g_{h}-g\|_{L^{2}(\Omega)}$ to decay faster than $O(h^{r+1})$. It should also be noted that the assumptions in the theorem statement guarantee that $g_{h}\in\Sigma_{h+}$ for $h$ sufficiently small, so that $\kappa_{h}(g_{h})$ is well-defined. See Section 4.2 for more details.

Our proof will be structured as follows. First, we verify in Lemma 4.6 that the exact Gaussian curvature satisfies

[TABLE]

for every $v\in V$. Next, we consider an arbitrary $v\in H^{1}_{0}(\Omega)$ and write

[TABLE]

with $v_{h}\in V_{h}$. The three terms above are estimated in Propositions 4.14, 4.16, and 4.18. Combining them will lead to Theorem 4.1.

4.1 Consistency of $b_{h}$

We begin by recalling two integration-by-parts formulas.

Lemma 4.2**.**

Let $K\in\mathcal{T}_{h}$. For every $v\in H^{1}(K)$ and every $w\in H^{1}(K)\otimes\mathbb{R}^{2}$, we have

[TABLE]

In addition, for every $v\in H^{2}(K)$ and every $\sigma\in H^{1}(K)\otimes\mathbb{S}$, we have

[TABLE]

Proof 4.3**.**

The (Euclidean) divergence theorem gives

[TABLE]

where we have used the fact that

[TABLE]

To prove (20), note first that

[TABLE]

which can be verified by expanding both sides in coordinates. The formula (20) now follows from the observation that

[TABLE]

*where we have used (21) again. *

Lemma 4.4**.**

For every $\sigma\in H^{2}(\Omega)\otimes\mathbb{S}$ and every $v\in V$,

[TABLE]

Proof 4.5**.**

Let $\widetilde{\sigma}=S_{g}\sigma$. Recalling that $\tau_{g}^{T}\sigma\tau_{g}=-n_{g}^{T}\widetilde{\sigma}n_{g}$, we have

[TABLE]

Using the fact that $\tau_{g}\tau_{g}^{T}+n_{g}n_{g}^{T}=g^{-1}$, the first and third terms can be combined to give

[TABLE]

Integrating the second term by parts, we obtain

[TABLE]

Since $v\in V$, both $v$ and its tangential derivative are continuous across element interfaces, and they vanish on $\partial\Omega$. In particular, $\llbracket\tau_{g}^{T}g\nabla_{g}v\rrbracket=\llbracket\tau^{T}\nabla v/\sqrt{\tau^{T}g\tau}\rrbracket=0$ for every $e\in\mathcal{E}_{h}$. On the other hand, since $\widetilde{\sigma}\in H^{2}(\Omega)\otimes\mathbb{S}$, $n_{g}^{T}\widetilde{\sigma}\tau_{g}$ is continuous across element interfaces, and $\llbracket(\operatorname{div}_{g}\widetilde{\sigma})^{T}gn_{g}\rrbracket=0$ for every $e\in\mathring{\mathcal{E}}_{h}$. It follows that the summations over $e\in\mathcal{E}_{h}$ in (23) vanish. Thus,

[TABLE]

Lemma 4.6**.**

For every $v\in V$,

[TABLE]

Proof 4.7**.**

Let $G(t)=(1-t)\delta+tg$ and $\sigma=g-\delta$. From (5), we obtain

[TABLE]

so

[TABLE]

4.2 Basic Estimates

We now collect a few basic estimates in preparation for our estimation of $\langle\kappa_{h}(g_{h})-\kappa(g),v\rangle_{g}$.

Throughout what follows, we make use of the fact that the eigenvalues of $g$, being positive and continuous on $\overline{\Omega}$, are bounded above and below by positive constants $C_{3}$ and $C_{4}$ that depend on $\|g\|_{L^{\infty}(\Omega)}$ and $\|g^{-1}\|_{L^{\infty}(\Omega)}$. Hence, for any $w\in\mathbb{R}^{2}$, we have

[TABLE]

It follows also that $\|\sqrt{\det g}\|_{L^{\infty}(\Omega)}$ and $\|1/\sqrt{\det g}\|_{L^{\infty}(\Omega)}$ are each bounded above by constants depending on $\|g\|_{L^{\infty}(\Omega)}$ and $\|g^{-1}\|_{L^{\infty}(\Omega)}$. Differentiation then reveals that $\|\sqrt{\det g}\|_{W^{1,\infty}(\Omega)}$ and $\|1/\sqrt{\det g}\|_{W^{1,\infty}(\Omega)}$ are each bounded above by constants depending on $\|g\|_{W^{1,\infty}(\Omega)}$ and $\|g^{-1}\|_{L^{\infty}(\Omega)}$.

We will use this information to establish analogous estimates for $g_{h}$, under the assumption that $C_{0}:=\sup_{h>0}\|g_{h}\|_{W^{1,\infty}_{h}}<\infty$ and $\lim_{h\rightarrow 0}\|g_{h}-g\|_{L^{\infty}(\Omega)}=0$. From this point forward, we use the letter $C$ to denote a constant which is not necessarily the same at each occurrence and may depend on $\Omega$, $\|g\|_{W^{1,\infty}(\Omega)}$, $\|g^{-1}\|_{L^{\infty}(\Omega)}$, $\|\kappa(g)\|_{L^{2}(\Omega)}$, $r$, $q$, $C_{0}$, and the shape regularity and quasi-uniformity constants $C_{1}$ and $C_{2}$, but not on $h$.

Lemma 4.8**.**

Let $p\in[1,\infty]$. For every $h$ sufficiently small,

[TABLE]

Proof 4.9**.**

The identity

[TABLE]

shows that for $h$ sufficiently small,

[TABLE]

Thus, for every $K\in\mathcal{T}_{h}$ and every $x\in\mathring{K}$, the operator norm of $g_{h}(x)^{-1}-g(x)^{-1}$ satisfies

[TABLE]

It follows that

[TABLE]

*for $h$ sufficiently small. Taking $h$ small enough so that (say) $\|g-g_{h}\|_{L^{\infty}(\Omega)}\leq\frac{1}{2\|g^{-1}\|_{L^{\infty}(\Omega)}}$ completes the proof. *

From the above lemma and the assumption that $\sup_{h>0}\|g_{h}\|_{W^{1,\infty}_{h}}<\infty$, we conclude that for $h$ sufficiently small,

[TABLE]

In particular, the eigenvalues of $g_{h}$ are bounded above and below by positive constants independent of $h$, and estimates of the form

[TABLE]

and

[TABLE]

hold for $h$ sufficiently small. In the remainder of this section, we will tacitly assume that $h$ is small enough that the preceding estimates hold.

Lemma 4.10**.**

For every $k\in\{0,1\}$ and $p\in[1,\infty]$,

[TABLE]

Proof 4.11**.**

Let $g=\left(\begin{smallmatrix}a&b\\ c&d\end{smallmatrix}\right)$ and $g_{h}=\left(\begin{smallmatrix}a_{h}&b_{h}\\ c_{h}&d_{h}\end{smallmatrix}\right)$. The identity

[TABLE]

shows that

[TABLE]

and the identity

[TABLE]

shows that

[TABLE]

To obtain the $W^{1,p}_{h}(\Omega)$ estimate, we compute

[TABLE]

which gives

[TABLE]

From this, (25), and Lemma 4.8, we conclude that

[TABLE]

4.2.1 The $L^{2}(\Omega,g)$-Orthogonal Projector

Let $P_{h}:L^{2}(\Omega)\rightarrow V_{h}$ denote the $L^{2}(\Omega,g)$-orthogonal projector onto $V_{h}$, so that

[TABLE]

Lemma 4.12**.**

We have

[TABLE]

Proof 4.13**.**

For any $u\in L^{2}(\Omega)$ and any $u_{h}\in V_{h}$, we have

[TABLE]

The equivalence of the norms $\|\cdot\|_{L^{2}(\Omega,g)}$ and $\|\cdot\|_{L^{2}(\Omega)}$ implies

[TABLE]

so $\|P_{h}u-u\|_{L^{2}(\Omega)}\leq C\|u_{h}-u\|_{L^{2}(\Omega)}$. Since $u_{h}$ was arbitrary, (27) follows. It follows also that $\|P_{h}u-u\|_{L^{2}(\Omega)}\leq Ch|u|_{H^{1}(\Omega)}$ if $u\in H^{1}_{0}(\Omega)$. Now let $Q_{h}:H^{1}(\Omega)\rightarrow V_{h}$ denote the Scott-Zhang interpolation operator [26]. Using an inverse estimate, interpolation estimates, and the stability of $Q_{h}$ in $H^{1}(\Omega)$, we obtain

[TABLE]

Finally, for any $v\in H^{1}_{0}(\Omega)$,

[TABLE]

so

[TABLE]

4.3 Proof of Theorem 4.1

We are now ready to estimate the three terms in (18).

Proposition 4.14**.**

If $v_{h}=P_{h}v$, then

[TABLE]

Proof 4.15**.**

For any $u_{h}\in V_{h}$, we have

[TABLE]

so

[TABLE]

*Taking $u_{h}=P_{h}\kappa(g)$ and invoking the $H^{-1}(\Omega,g)$-estimate (26) completes the proof. *

Proposition 4.16**.**

For $h$ sufficiently small, we have

[TABLE]

Proof 4.17**.**

By Lemma 4.10 and the inverse estimate $\|v_{h}\|_{L^{\infty}(\Omega)}\leq C\log h^{-1}\|v_{h}\|_{H^{1}(\Omega)}$ [6, Lemma 4.9.2], we have

[TABLE]

Now note that

[TABLE]

The second term above is bounded by a constant (a multiple of $\|\kappa(g)\|_{L^{2}(\Omega)}$) owing to (27). For the first term, an inverse estimate gives

[TABLE]

Thus, using the $H^{-1}(\Omega,g)$-estimate (26), we obtain

[TABLE]

for $h$ sufficiently small, since $\lim_{h\rightarrow 0}\inf_{u_{h}\in V_{h}}\|\kappa(g)-u_{h}\|_{L^{2}(\Omega)}=0$ [15, Equation 1.99]. It follows that

[TABLE]

*Combining this with (29) completes the proof. *

We will now estimate the first term in (18). We state the result below, and prove it with a series of lemmas.

Proposition 4.18**.**

We have

[TABLE]

To prove Proposition 4.18, let $G(t)=(1-t)\delta+tg$ and $G_{h}(t)=(1-t)\delta+tg_{h}$, so that

[TABLE]

Note that $G_{h}(t)$, being a convex combination of $\delta$ and $g_{h}$, has eigenvalues bounded above and below by positive constants independent of $h$ (and $t$), and all of the estimates we established for the norms of $g_{h}-g$, $g_{h}^{-1}-g^{-1}$, $\sqrt{\det g_{h}}-\sqrt{\det g}$, etc. carry over easily to $G_{h}-G$, $G_{h}^{-1}-G^{-1}$, $\sqrt{\det G_{h}}-\sqrt{\det G}$, etc.

Lemma 4.19**.**

For every $t\in[0,1]$,

[TABLE]

Proof 4.20**.**

Let $\sigma=g_{h}-\delta$. We have

[TABLE]

where

[TABLE]

Along any edge $e\in\mathcal{E}_{h}$, we may write

[TABLE]

so, with $\widetilde{G}=G(\det G)^{-1/2}(\tau^{T}G\tau)^{-1}$ and $\widetilde{G}_{h}=G_{h}(\det G_{h})^{-1/2}(\tau^{T}G_{h}\tau)^{-1}$, we have

[TABLE]

Writing

[TABLE]

shows that for each $K_{j}\in\mathcal{T}_{h}$ with $e\subset\partial K_{j}$,

[TABLE]

The trace inequality then gives, if $e\in\mathring{\mathcal{E}}_{h}$,

[TABLE]

For edges $e$ on $\partial\Omega$, the same estimate holds with terms involving $K_{2}$ replaced by zero.

Turning now to $B_{K}$, observe that

[TABLE]

Observe also that

[TABLE]

where

[TABLE]

and $\Gamma_{ij}^{m}$ are the Christoffel symbols associated with $G$. If $a_{h}^{\ell k}$ is defined similarly with $G_{h}$ in place of $G$, then it is not hard to see that

[TABLE]

and

[TABLE]

so

[TABLE]

Using a similar argument for the second term in (32), we get

[TABLE]

Hence,

[TABLE]

*The proof is completed by summing (30) over all edges $e$, summing (33) over all triangles $K$, and noting that $\|\sigma\|_{L^{\infty}(\Omega)}=\|g_{h}-\delta\|_{L^{\infty}(\Omega)}\leq C$ and $G_{h}-G=t(g_{h}-g)$. *

Lemma 4.21**.**

For every $t\in[0,1]$,

[TABLE]

Proof 4.22**.**

Let $\sigma=g_{h}-g$, so that

[TABLE]

As was shown in the proof of Lemma 4.19, we have

[TABLE]

where $\widetilde{G}=G(\det G)^{-1/2}(\tau^{T}G\tau)^{-1}$. Thus, if $e=K_{1}\cap K_{2}\in\mathring{\mathcal{E}}_{h}$, then the trace inequality gives

[TABLE]

and similarly for edges $e$ on $\partial\Omega$. On the other hand,

[TABLE]

*The conclusion follows from summing over all edges $e$ and all triangles $K$. *

We will now finish the proof of Theorem 4.1. Taking $v_{h}=P_{h}v$ and combining Propositions 4.14, 4.16, and 4.18 with the inverse estimate $|v_{h}|_{H^{2}_{h}(\Omega)}\leq Ch^{-1}\|v_{h}\|_{H^{1}(\Omega)}$, the stability estimate $\|v_{h}\|_{H^{1}(\Omega)}\leq C\|v\|_{H^{1}(\Omega)}$, and the upper bound $\log h^{-1}\leq h^{-1}$, we obtain

[TABLE]

Taking the supremum over $v\in H^{1}_{0}(\Omega)$ and rearranging gives

[TABLE]

Since we have assumed that $\lim_{h\rightarrow 0}h^{-1}\log h^{-1}\|g_{h}-g\|_{L^{2}(\Omega)}=0$, we may divide by $\left(1-Ch^{-1}\log h^{-1}\|g_{h}-g\|_{L^{2}(\Omega)}\right)$ for $h$ sufficiently small to arrive at the $H^{-1}(\Omega,g)$-estimate in (16).

To obtain the $L^{2}(\Omega)$-estimate in (16), observe that for any $v\in L^{2}(\Omega)$, we have

[TABLE]

Now since

[TABLE]

we deduce that

[TABLE]

Invoking the $H^{-1}(\Omega,g)$-estimate in (16) gives

[TABLE]

Next, assume $\kappa(g)\in H^{m}(\Omega)\cap H^{1}_{0}(\Omega)$, let $0\leq k\leq\min\{q+1,m\}$, and let $v_{h}\in V_{h}$ be the Scott-Zhang interpolant of $\kappa(g)$. Then, using standard inverse estimates and interpolation error estimates, we obtain

[TABLE]

for each $\ell=k,k+1,\dots,\min\{q+1,m\}$. Combining this with the $L^{2}(\Omega)$-estimate (34) gives (17). This completes the proof of Theorem 4.1.

5 Numerical Examples

In this section, we present numerical experiments to illustrate the convergence rates predicted by Theorem 4.1. We focus on the $L^{2}(\Omega)$-error, which, upon taking $k=0$ and $\ell=q+1$ in (17), satisfies

[TABLE]

assuming that $g\in H^{r+1}(\Omega)\otimes\mathbb{S}$, $\kappa(g)\in H^{q+1}(\Omega)$, and $g_{h}$ is taken equal to a suitable interpolant of $g$ [24, Theorem 2.6].

In our numerical implementation, we took $g_{h}$ equal to a subdivision-based interpolant of $g$ defined as follows. For each $K\in\mathcal{T}_{h}$, let $K_{1},K_{2},\dots,K_{(r+1)^{2}}$ be the triangles formed by partitioning $K$ along the lines given in barycentric coordinates $(\lambda^{0},\lambda^{1},\lambda^{2})$ by $\lambda^{i}=j/(r+1)$, $i=0,1,2$, $j=1,2,\dots,r$, as depicted in Figure 2. Let $\mathcal{E}_{h}^{r}(K)$ be the union of the edges of $K_{1},K_{2},\dots,K_{(r+1)^{2}}$. Let $\mathcal{E}_{h}^{r}=\cup_{K\in\mathcal{T}_{h}}\mathcal{E}_{h}^{r}(K)$. For each $e\in\mathcal{E}_{h}^{r}$, let $z^{(e)}\in\Omega$ denote the midpoint of $e$. It is known that the linear functionals

[TABLE]

form a basis for the dual of $\Sigma_{h}$, where $\tau$ is a unit vector (relative to the Euclidean metric) tangent to $e$ [24, pp. 38-39]. Let $\{\psi_{e}\}_{e\in\mathcal{E}_{h}^{r}}\subset\Sigma_{h}$ denote the basis for $\Sigma_{h}$ dual to these functionals:

[TABLE]

We took

[TABLE]

Table 1 shows the errors $\|\kappa_{h}(g_{h})-\kappa(g)\|_{L^{2}(\Omega)}$ for the metric (1) on $\Omega=(-1,1)\times(-1,1)$. We computed the discrete Gaussian curvature using Definition 3.1 with $q=1$ and $r=0,1,2$ on triangulations with maximum element diameter $h\in[0.02,0.17]$. We constructed each triangulation by randomly perturbing the interior vertices of a uniform triangulation of $\Omega$, as depicted in Figure 1a. Observe that the $L^{2}(\Omega)$-error converges linearly when $r=2$ and does not converge when $r<2$, in agreement with the estimate (35).

Acknowledgements

I am grateful to Melvin Leok for many helpful discussions.