Analysis of curvature approximations via covariant curl and incompatibility for Regge metrics

The metric tensor of a Riemannian manifold can be approximated using Regge finite elements and such approximations can be used to compute approximations to the Gauss curvature and the Levi-Civita connection of the manifold. It is shown that certain Regge approximations yield curvature and connection approximations that converge at a higher rate than previously known. The analysis is based on covariant (distributional) curl and incompatibility operators which can be applied to piecewise smooth matrix fields whose tangential-tangential component is continuous across element interfaces. Using the properties of the canonical interpolant of the Regge space, we obtain superconvergence of approximations of these covariant operators. Numerical experiments further illustrate the results from the error analysis.


Introduction
This paper is concerned with the finite element approximation of the Gauss curvature K of a two-dimensional Riemannian manifold.As shown by Gauss's Theorema Egregium, K is an intrinsic quantity of the manifold.It can be computed solely using the metric tensor of the manifold.Therefore, when a finite element approximation of the metric tensor is given, it is natural to ask if an approximation to K can be computed.The answer was given in the affirmative by the recent work of [27], assuming that the metric is approximated using Regge finite elements, and further improved by [10].The convergence theorems of this paper are heavily based on these works.We prove that the resulting curvature and connection approximations converge at a higher rate than previously known for the approximation given by the canonical Regge interpolant.Our method of analysis is different and new.In particular, we show how covariant curl and incompatibility can be approximated using appropriate finite element spaces, given a nonsmooth Regge metric.These operators arise in a myriad of other applications, so our intermediate results regarding them are of independent interest.Notions of curvature while gluing together piecewise smooth metrics has been a preoccupation in varied fields far away from computing, as early as [29] to recent years [45], so we note at the outset that we approach the topic with numerical computation in mind.
The Regge finite element takes its name from Regge calculus, originally developed for solving Einstein field equations in general relativity.It discretizes the metric tensor through edge-length specifications, allowing the curvature to be approximated by means of angle deficits [40].Regge calculus was established in theoretical and numerical physics and routinely finds applications in relativity and quantum mechanics.In [52,41,9] a comprehensive overview of the development of Regge calculus over the last fifty years can be found.Just as Whitney forms [51] can be interpreted as finite elements, it was observed that Regge's approach of prescribing quantities on edges is equivalent to defining a piecewise constant metric tensor whose tangential-tangential components are continuous across element interfaces [44].The first rigorous proof of convergence of Regge's angle deficits to the scalar curvature, for a sequence of appropriate triangulations in the sense of measures, was accomplished in [15].Later, it was also shown [19] that for a given metric in the lowest order Regge finite element space, a sequence of mollified metrics converges to the angle deficit in the sense of measures.Methods based on angle deficits for approximating the Gauss curvature on triangulations consisting of piecewise flat triangles are well-established in discrete differential geometry and computer graphics.On specific triangulations satisfying certain conditions, convergence in the L 8 -norm up to quadratic order was proven, but for a general irregular grid there is no reason to expect convergence [12,53,54].In [35], Regge's concept of angle deficits has been extended to quadrilateral meshes.Notable among the results applicable for higher dimensional manifolds is the proof of convergence for approximated Ricci curvatures of isometrically embedded hypersurfaces Ă R n`1 presented in [24], and used later for Ricci flows [25].
Another natural perspective to place the modern developments on the Regge finite element is within the emergence of finite element exterior calculus (FEEC) [6,5].Finite element structures for Regge calculus were developed in [17,18] and the resulting elements became popular in FEEC under the name Regge finite elements [34].Regge elements approximating metric and strain tensors were extended to arbitrary polynomial order on triangles, tetrahedra, and higher dimensional simplices in [34], and for quadrilaterals, hexahedra, and prisms in [36].The utility of Regge elements when discretizing parts of the Kröner complex, or the elasticity complex, was considered in [7,18,28].Properties of Regge elements were exploited to construct a method avoiding membrane locking for general triangular shell elements [37].
In this backdrop, the recent work of [27] provides an interesting application of Regge elements by developing a high-order Gauss curvature approximation formula based on higher degree Regge elements.(It was applied to Ricci and Ricci-DeTurck flow [26].)The key is an integral formulation of the angle deficit, extendable to higher orders.Using it, rigorous proofs of convergence at specific rates were proved in [27].Even more recently, in [10], this approach has been reformulated in terms of a nonlinear distributional curvature and connection 1-form (Levi-Civita connection), using the element-wise Gauss curvature, jump of geodesic curvature at edges, and angle deficits at vertices as sources of curvature [47,46].The authors show that L 2 -convergence of the approximated curvature can be obtained if Regge elements using polynomials of degree at least two are used to approximate the metric.This is in line with the rule of thumb that a second order differential operator approximated using polynomials of degree k leads to convergence rates of order k ´2.Nonetheless, convergence rates better than this rule of thumb have often been observed in compatible discretizations in FEEC.One of our goals in this paper is to establish such an improved rate for the curvature and connection approximations, as well as for the intermediate covariant operators arising in our analysis, such as the curl and incompatibility.
In a later section, we extend the ideas in [27,10] by exploiting certain orthogonality properties for the error in the canonical interpolation by Regge finite elements to obtain one extra order of convergence for the curvature approximation.This extra order is comparable to super-convergence properties of mixed methods [11,21] and has been observed for the Hellan-Herrmann-Johnson method for the biharmonic plate and shell equation [8,49].
Another difference in our analysis, in comparison to [27,10], is the use of the intrinsic (or covariant) incompatibility operator (which we define using the covariant curl on the manifold).It is now well known that linearizing the curvature operator around the Euclidean metric gives a first order term involving the incompatibility operator [18] and we exploit this relationship in the analysis of the curvature approximation.On Euclidean manifolds, the incompatibility operator is well known to be the natural differential operator for Regge elements in any dimension.By showing that the curvature approximation can be analyzed via the incompatibility operator, we hope to generate new ideas for computing and analyzing approximations of the intrinsic curvature tensor of higher dimensional manifolds.The incompatibility operator also arises in modeling elastic materials with dislocations [1,2], another potential area of application.The key insight on which we base our definition of these covariant operators for Regge metrics is revealed by the essential role played by a glued smooth structure (described in §4.2).Since coordinates in this glued smooth structure are generally inaccessible for computations, we detail how to compute these operators in the coordinates in which the Regge metric is given as input.This paper can be read linearly, but we have structured it so a numerical analyst can also directly start with the error analysis in Section 6-where the main convergence theorems appear in §6.1-referring back to the previous sections as needed.(Only coordinate-based formulas are used in §6; their derivations from intrinsic geometry are in the previous sections.)The next section ( §2) establishes notation and introduces geometric preliminaries and finite element spaces.Section 3 defines the curvature approximation formula and details coordinate formulas we use for numerical computations.In §4, covariant differential operators on piecewise smooth metric tensors are derived, concentrating on the covariant curl and incompatibility operator, and how they arise from linearization of curvature.Section 5 is devoted to the approximation of the connection 1-form.Section 6 is devoted to the numerical analysis of the errors in the method.The analysis is performed by first proving optimal convergence rates for the distributional covariant curl and inc, and then for the approximations of the Gauss curvature and connection 1-form.Numerical examples illustrating the theoretical results are presented in §7.

Notation and preliminaries
This section provides definitions that we use throughout.We give intrinsic definitions of quantities on a manifold, but in view of our computational goals, we also make extensive use of coordinate expressions.We use the Einstein summation convention, by which a term where the same integer index appears twice, as both an upper and a lower subscript, is tacitly assumed to be summed over the values of that index in t1, 2u.Summation convention does not apply when a repeated index is not an integer (such as when a subscript or a superscript represents a vector field or other non-integer quantities).
2.1.Spaces on the manifold.Let M denote a two-dimensional oriented manifold with or without boundary.Endowed with a smooth metric ḡ, let pM, ḡq be a Riemannian manifold.Let the unique Levi-Civita connection generated by ḡ be denoted by ∇.Let XpM q, Ź k pM q, and T k l pM q denote, respectively, the sets of smooth vector fields on M , k-form fields on M , and pk, lq-tensor fields on M .The value of a tensor ρ P T k l pM q acting on k vectors X i P XpM q and l covectors µ j P Ź 1 pM q is denoted by ρpµ 1 , . . ., µ l , X 1 , . . ., X k q.Note that Ź 1 pM q " T 1 0 pM q and XpM q " T 0 1 pM q.Note also that it is standard to extend the Levi-Civita connection ∇ from vector fields to tensor fields (see e.g., [32,Lemma 4.6]) so that Leibniz rule holds.
For coordinate computations, we use a chart to move locally to a Euclidean domain with coordinates x 1 , x 2 .Let the accompanying coordinate frame and coframe be denoted by B i and dx i .We assume these coordinates preserve orientation, so the orientation of M is given by the ordering pB 1 , B 2 q.Let SpM q " tσ P T 2 0 pM q : σpX, Y q " σpY, Xq for X, Y P XpM qu and S `pM q " tσ P SpM q : σpX, Xq ą 0 for 0 ‰ X P XpM qu.They represent the subspace of symmetric tensors in T 2 0 pM q, whose elements σ can be expressed in coordinates as σ " σ ij dx i b dx j with smoothly varying coefficients σ ij satisfying σ ij " σ ji and are additionally positive definite, respectively.
We will use standard operations on 2-manifold spaces such as the Hodge star ‹ : Ź k pM q Ñ Ź 2´k pM q, the exterior derivative d k : Ź k pM q Ñ Ź k`1 pM q, the tangent to cotangent isomorphism 5 : XpM q Ñ Ź 1 pM q, and the reverse operation 7 : Ź 1 pM q Ñ XpM q.Their definitions can be found in standard texts [32,38,48].

2.2.
Curvature.The exact metric ḡ is an element of S `pM q.We define the Riemann curvature tensor R P T 4 0 pM q of the manifold following [32], If X and Y are linearly independent, the Gauss curvature of M can be expressed by whose value is well known to be independent of the choice of the basis (see, e.g., [32, p. 144] or [14,Ch. 4,Proposition 3.1]).We will also need the geodesic curvature along a curve Γ in the manifold pM, ḡq.To recall its standard definition (see [48, p. 140] or [32]), we let 0 ă s ă a be the ḡ-arclength parameter so that Γ is described by µpsq for some smooth µ : r0, as Ñ M and its ḡ-unit tangent vector is tpsq " dµ{ds.Let npsq be such that ptpsq, npsqq is a ḡ-orthonormal set of two vectors in the tangent space whose orientation is the same as that of M , i.e., dx 1 ^dx 2 ptpsq, npsqq ą 0. Then κpḡq " ḡp ∇tpsq tpsq, npsqq (2.3) gives the geodesic curvature at the point µpsq of Γ .
2.3.Approximate metric.We are interested in approximating Kpḡq when the metric is given only approximately.We assume that M has been subdivided into a geometrically conforming triangulation T .The edges of T may be curved, but do not necessarily consist of geodesics.
On each element T P T , we are given an approximation g| T P SpT q of ḡ| T .When the approximation is sufficiently good, g will also be positive definite since ḡ is.Then each T P T can be considered to be a Riemannian manifold pT, g| T q with g| T as its metric.Since g| T is smooth within each element T (not across BT ), we use the unique Levi-Civita connection ∇ generated by g| T to compute covariant derivatives within T .(Constraints on g across element boundaries are clarified below in (2.12)).We drop the accent ¯in any previous definition to indicate that it pertains to the manifold pT, g| T q instead of pM, ḡq, e.g., R refers to the Riemann curvature tensor computed using g and ∇ in place of ḡ and ∇ in (2.1).
A point p P T can be viewed either as a point in the manifold M or as a point in the manifold T .Irrespective of the two viewpoints, the meanings of coordinate frame B i , coframe dx i , and the tangent space T p M at p are unchanged.In coordinates, g ij " gpB i , B j q, g ij " g ´1pdx i , dx j q (2.4) may be viewed as entries of symmetric positive definite matrices.Christoffel symbols of the first kind (Γ ijk ) and the second kind They can alternately be expressed, using (2.4), as Later, we will also use Γ ijl pσq to denote 1 2 pB i σ jl `Bj σ li ´Bl σ ij q for other tensors σ in T 2 0 pM q. 2.4.Tangents and normals on element boundaries.Throughout this paper, we use τ to denote a tangent vector (not g-normalized; cf.(3.15) later) along an element boundary BT for any T P T .The orientation of τ is aligned with the boundary orientation of BT (inherited from the orientation of T , which is the same as the orientation of M ).For any p P BT, define ν P T p M by gpν, Xq " pdx 1 ^dx 2 qpτ, Xq, for all X P T p M. (2.6) It is easy to see from (2.6) that the ordered basis pτ, νq has the same orientation as pB 1 , B 2 q since pdx 1 ^dx 2 qpτ, νq ą 0, and moreover, gpν, τ q " 0, and gpν, νq detpgq " gpτ, τ q. (2.7) In particular, defining we obtain a g-orthonormal basis pτ , νq of normal and tangent vectors along every element boundary BT , whose orientation matches the manifold's orientation.E.g., if M is the unit disc in R 2 with the Euclidean metric g " δ and the standard orientation, then τ is oriented counterclockwise and ν points inward.In (2.8) and throughout, we use σ uv to denote σpu, vq for any vectors u, v P T p M and σ P SpM q.Note that g uv is not to be confused with the g ij introduced in (2.4) where the indices are integers (which trigger the summation convention) rather than vectors.
To write ν in coordinates, it is useful to introduce the alternating symbol ε ij whose value is 1, ´1, or 0 according to whether pi, jq is an even permutation, odd permutation, or not a permutation of p1, 2q, respectively.The value of the symbols ε ij , ε j i , and and furthermore (e.g., using (2.9), (2.7), and the Hodge star in coordinates, ‹dx i " ?det gg ij ε jk dx k ) that p‹ωqpνq " ωpτ q, p‹ωqpτ q " ´ωpνq, for any ω P Ź 1 pM q.
2.5.Finite element spaces.Let C 8 pT q denote the space of piecewise smooth functions on T , by which we mean functions that are infinitely smooth within each mesh element and continuous up to (including) the boundary of each mesh element T .A notation in §2.1 with T in place of M indicates the piecewise smooth analogue, e.g., XpT q " tX " X i B i : X i P C 8 pT qu, SpT q " tσ ij dx i b dx j : σ ij " σ ji P C 8 pT qu, etc.Note that a σ P SpT q need not be continuous across the element interfaces.Let E " BT `XBT ´denote an interior mesh edge (possibly curved) shared between elements T ˘P T .Let T p E denote the one-dimensional tangent space of the curve E at any one of its points p. (Note that the tangent space at p from either element T ˘coincides with T p M and T p E Ă T p M .)We say that a σ P SpT q has "tangential-tangential continuity" or that σ is tt-continuous if at all p P E, and for every interior mesh edge E. Let RpT q " tσ P SpT q : σ is tt-continuousu, (2.12a) R `pT q " tσ P RpT q : σpX, Xq ą 0 for all X P T p M u. (2.12b) The approximate metric g is assumed to be in R `pT q.For the numerical analysis later, we will additionally assume that it is in R k h defined below.In finite element computations, we use a reference element T, the unit triangle, and the space P k p T q of polynomials of degree at most k on T .Let T denote a Euclidean triangle with possibly curved edges that is diffeomorphic to T via Φ : T Ñ T .For finite element computations on manifolds, we need charts so that each whole element T P T of the manifold is covered by a single chart giving the coordinates x i on T .The chart identifies the parameter domain of T as the (possibly curved) Euclidean triangle T diffeomorphic to T .Let Φ : T Ñ T denote the diffeomorphism.Then Φ T " Φ´1 ˝Φ : T Ñ T maps diffeomorphically to the reference element where P k p T q is defined.We use its pullback Φ T below.
Define the Regge finite element space of degree k on the manifold M by R k h " tσ P RpT q : for all T P T , σ| T " σ ij dx i b dx j with σ ij " Φ T q for some q P P k p T qu. (2.13) The subscript h indicates a mesh size parameter, e.g., on meshes whose elements are close to straight-edged triangles, one may set h " max T PT diampT q.Let VpT q " tu P Ź 0 pT q : u is continuous on M u, VΓ pT q " tu P VpT q : u| Γ " 0u, where Γ denote a subset of the boundary BM of positive boundary measure.
The Lagrange finite element space on M and its subspaces with essential boundary conditions are defined by V k h " tu P VpT q : for all T P T , u| T " Φ T û for some û P P k p T qu, Vk h,Γ " tu P V k h : u| Γ " 0u and Vk h " Vk h,BM . ( The previous definitions in this subsection were independent of the metric.We will now introduce a metric-dependent space of normal-continuous vector fields.First, we introduce the following notation surrounding an interior mesh edge E shared by two adjacent elements in T , In this context, the g-orthonormal tangent and normal vectors introduced above along BT ˘are denoted by τ ˘and ν˘, respectively.For a collection of scalar functions, tf BT pνq : T P T u, each depending on the normal ν at an element boundary, we define the jump on E by f pνq " f BT `p ν`q `fBT ´p ν´q . (2.15b) Thus the jump function f pνq is well defined (and single-valued) on the union of all interior mesh edges, excluding the mesh vertices.The jump of an element boundary function dependent on τ is defined similarly.We say that a piecewise smooth vector field W P XpT q has "g-normal continuity" across element interfaces if gpW, νq " 0. Define W g pT q " tW P XpT q : gpW, νq " 0u, Wg,Γ pT q " tW P W g pT q : gpW | Γ , νq " 0u, Wg pT q " Wg,BM pT q.
(2.16) Also define their polynomial subspaces W k g,h " tW P W g pT q : for all T P T , W | T " Φ T Ŵ for some Ŵ P P k p T, R 2 qu, Wk g,h,Γ " tW P W k g,h : gpW | Γ , νq " 0u, Wk g,h " Wk g,h,BΩ .
2.6.Integrals over the manifold's triangulation.On every element T P T , in order to integrate a scalar function f P Ź 0 pT q, adopting the notation of [33], we tacitly use the unique Riemannian volume form vol T,g to convert it to a 2-form and then pullback to integrate over the Euclidean parameter domain T , i.e., where detpDΦq denotes the Jacobian determinant of Φ, we have used the standard extension of pullback to forms, and we have appended an area measure notation "da" to emphasize that the right most integral is a standard Lebesgue integral over the Euclidean domain T .For v, w with the understanding that the right hand sides above must be evaluated using (2.17).
In order to integrate along the boundary curve BT , we use the one-dimensional analogue of the formula in (2.17) to compute on the Euclidean domain B T " ΦpBT q, namely ż pBT,gq f " where t is a tangent vector along B T of unit Euclidean length-and to emphasize that the last integral is a standard Euclidean integral, we have appended the length measure "dl".We use w to simplify notation for sum of integrals over element boundaries.

Curvature approximation
In this section we give the curvature approximation formula and discuss a few nontrivial computational details on curved elements.In order to approximate the Gauss curvature Kpḡq, one may consider computing Kpg| T q on each element T P T using the given approximation g of the exact metric ḡ.However, this alone cannot generally be a good approximation to Kpḡq because discontinuities of g across elements generate additional sources of curvature on the edges and vertices of the mesh.Below we provide a curvature approximation incorporating these extra sources.Since it coincides with the formula given in a recent work [10] for a specific case, we opt for a brief description, expanding only on aspects complementary to that work.

3.1.
A finite element curvature approximation.Given a metric g P R `pT q approximating ḡ we identify three sources of curvature, modeled after similar terms in the Gauss-Bonnet formula, and define them as the following linear functionals acting on ϕ P VpT q: xK T g , ϕy VpT q " ż pT,gq Kpgq ϕ, xK T E,g , ϕy VpT q " ż pE,gq κpgq ϕ, where T V denotes the set of all elements of T which have V as a vertex.Here Kpgq and κpgq are defined by (2.2) and (2.3) after replacing ḡ by g, and T V p¨q denotes the interior angle at a vertex V of T determined using the metric in its argument (computable using (3.5) below).Throughout, we use xf, ϕy H to denote duality pairing on a vector space H that gives the action of a linear functional f P H 1 acting on a ϕ P H. Also, V and E denote the set of mesh vertices and edges, respectively (so in (3.1), V P V and E P E ).Define K g P VpT q 1 by Here E T denotes the set of three edges of BT .In addition to g, suppose that we are also given boundary curvature data in essential (Dirichlet) or natural (Neumann) forms for manifolds with boundary.The former type of data arises when we know that M is a submanifold of a larger manifold whose Gauss curvature is known outside of M .To accommodate such information only on a part of the boundary of M , we split BM into two non-overlapping parts Γ D and Γ N .One of these can be empty.In case none of them is empty, both must have positive length.On Γ D , we assume that we are given K D " Kpḡq| Γ D and that K D is in the trace of the Lagrange finite element space V k h .E.g., when a manifold is flat around Γ D (i.e., Kpḡq vanishes in a neighborhood of Γ D ), we may set homogeneous Dirichlet data K D " 0 on Γ D .The other type of boundary data, in the form of a natural (or Neumann) boundary condition, is motivated by the Gauss-Bonnet theorem, and provides geodesic curvature data at the boundary.These natural Neumann-type boundary data is given in the form of a data functional for ϕ P VpT q, where V N Ă V X Γ N is the subset of the manifold's vertices contained in the interior of Γ N and ˜ N V denotes the exterior angle measured by the edges of Γ N at V .(If V is part of a smooth boundary, such an angle amounts to π, whereas at kinks of the boundary the angle has to be provided as input data.)The action of functional (3.3) on the finite-dimensional subspace V k h Ă VpT q is computable if we are given the exact metric ḡ on and near Γ N .For manifolds without boundary, there is no need to provide any boundary data.Definition 3.1.Let g P R `pT q and k ě 0 be an integer.The finite element curvature approximation K h pgq of degree k `1 is the unique function in V k`1 h determined by requiring that K h pgq| Γ D " K D on Γ D and for all u h P Vk`1 h,Γ D , ż T K h pgq u h " xK g , u h y VpT q ´xκ N , u h y VpT q .
(3.4) 3.2.Implementation issues.We now discuss how to numerically compute the quantities in (3.4) in the given computational coordinates x i .Recall (from §2.4) that the tangent vector τ along the boundary BT is aligned with the boundary orientation of BT .Let V T denotes the set of three vertices of an element T P T .At any vertex V P V T , the tangent τ undergoes a change in direction, between an incoming and an outgoing tangent vector, which we denote by τ ´and τ `, respectively.The angle at V with respect to the metric g is then computed by This is what we use to calculate the angle deficit functional K T V,g (3.1).Next, consider the interior source term K T g , defined using Kpgq, and related to the Riemann curvature by (2.2).By (2.1), where Γ k ij and Γ ijk are as in (2.5).For two-dimensional manifolds the Gauss curvature can be expressed [14]  It thus remains to discuss the computation of the edge sources K T E,g using the definition of the κpgq in (2.3).In finite element computations, we usually do not have ready access to the g-arclength parameter s used there.But κpgq can be computed using the readily accessible τ and ν of §2.4, as shown below.Let γptq be an orientation-preserving parametrization that gives an oriented mesh edge E Ă BT as E " tγptq : t 0 ď t ď t 1 u.Parametrizing scalar functions a on E by t, we abbreviate da{dt to 9 a.Note that the components of τ " τ k B k are given by τ k " 9 γ k and their derivatives by where u " u i B i , v " v i B i , w P XpT q and recall that ν was defined in (2.6).
Returning to the edge source term, using Lemma 3.2 and (2.18), Thus, through (3.5), (3.7), and (3.14), we have shown it is possible to easily compute all the terms in the curvature approximation (3.4) using standard finite element tools.

3.3.
A model case for further analysis.Having shown how curvature of general manifolds can be computed, we now focus our analysis on the following model case for the remainder of the paper.We assume that the manifold M , as a set, equals Ω Ă R 2 , a bounded open connected domain, that Γ N " H, and that ḡ is compactly supported in Ω (so Dirichlet boundary data is zero).The set Ω forms the full parameter domain of M for the single trivial chart Φ : M Ñ Ω with Φ equaling the identity map.The triangulation T of M is now a conforming finite element mesh in the planar domain Ω and its elements are (possibly curved) triangles.
In this setting, an element T P T can be considered as either the Euclidean manifold pT, δq, equipped with the identity metric δ, or the Riemannian manifold pT, gq.Let us reconsider the tangent vector τ on an element boundary BT with the orientation described in §2.4,previously left un-normalized.Henceforth, we assume that 1 " δpτ, τ q. (3.15) We computed the normal vector ν from τ by (2.6) for the manifold pT, gq.For the manifold pT, δq, a generally different Euclidean normal vector arises at any p in BT and we denote it by ν.It can be computed by simply replacing g with δ in (2.6), i.e., ν P T p M satisfies δpν, Xq " dx 1 ^dx 2 pτ, Xq for all X P T p M. (3.16) Analogous to (2.9), we now have the accompanying coordinate expression, Note that (3.15) implies that δpν, νq " 1.These identities, together with (2.8) guide us move between the δ-orthonormal tangent and normal vectors (τ and ν) and gorthonormal tangent and normal vectors (τ and ν), while preserving the orientation.Jumps of functions of ν and τ on the Euclidean element boundaries pBT, δq, T P T , are defined in analogy to (2.15b).
Lemma 3.3.The geodesic curvature along any mesh edge E is given by for any X P T p M. Using (3.20) in Lemma 3.2, the result follows.
In this model case of a single parameter domain, since the identity metric δ is well defined throughout, the angle deficit in (3.1) can be reformulated using δ as a sum of element contributions at each vertex.Define When summing T V δ over all V P V T , we clearly obtain 2π.Hence ÿ Proposition 3.4.In this model case, equation (3.4) implies that for all If g P S `pM q, then all terms vanish except the last.
Proof.Equation (3.22) follows from (3.21), (3.7), (3.14), and Lemma 3.3, once Φ is set to the identity.To prove the last statement, divide the set of mesh vertices V into the set of boundary and interior vertices V bnd " V X BM and V int " V zV bnd .At every interior vertex V P V int , a rearrangement gives ÿ which vanishes because the smoothness of g implies ř The element boundary integrals can be rewritten using the set of interior mesh edges E int , as since u h " 0 on Γ D and that the trace of ?det g{g τ τ is well defined (single valued) on E due to the given smoothness of g P S `pM q.It is easy to see that 9 τ has the same value from adjacent elements of E while τ and ν changes sign, so 9 τ ν " 0 and Γ ν τ τ " 0. Hence (3.23) vanishes.

Covariant derivatives using the nonsmooth metric
The objective of this section is to formulate a covariant incompatibility operator that can be applied to our situation with piecewise smooth metrics.To this end, we first define several covariant derivatives in the smooth case, restricting ourselves to a single element, i.e., the smooth Riemannian manifold pT, g| T q.Then we proceed to consider the changes needed due to the jumps of the metric across element interfaces.4.1.Covariant curl and incompatibility for smooth metric.For a µ P Ź 1 pT q, the exterior derivative d 1 µ P Ź 2 pT q is characterized in terms of the connection [38] by for any X, Y P XpT q.Next, for a σ P T 2 0 pT q writing σ Z pY q " σpZ, Y q for any X, Y, Z P XpT q, we define an operation analogous to (4.1) on σ by Since the expressions in (4.1) and (4.2)-holding Z fixed-are skew-symmetric in X and Y , they may be viewed as elements of Ź 2 pT q.Then we may use the Hodge star ‹ operation to convert them to 0-forms, since T is two dimensional.Doing so, define curl g µ " ‹pd 1 µq, µ P T 1 0 pT q " Ź 1 pT q, (4.3a) pcurl g σqpZq " ‹pd 1 σ Z q, σ P T 2 0 pT q, Z P XpT q. (4.3b) The latter, due to linearity in Z, is in Ź 1 pT q, while the former curl g µ is in Ź 0 pT q.Combining these operations in succession, we define the covariant incompatibility, on two-dimensional manifolds.Clearly inc g σ is in Ź 0 pT q.We will now quickly write down coordinate expressions of these tensors.Expanding the right hand side of (4.2) using the Leibniz rule p∇ X σqpZ, Y q " XσpZ, Y q σp∇ X Z, Y q ´σpZ, ∇ X Y q twice, substituting Z " B i , X " B j , Y " B k , and simplifying using (2.5a), for σ " σ jk dx j b dx k P T 2 0 pT q.Note that equation (4.4) can be rewritten as det g for any scalar field f P Ź 0 pT q, we arrive at Similarly, one obtains coordinate expressions for the remaining operators in (4.3), namely In the derivation of (4.5c), we have employed the useful formula It is useful to contrast the expressions in (4.5) with the standard Euclidean curl and inc.To this end we use matrix and vector proxies, rσs P R 2ˆ2 and rµs P R 2 , of µ P Ź 1 pT q and σ P T 2 0 pT q, respectively [3].These proxies are made up of coefficients in the coordinate basis expansion, specifically rσs is the matrix whose pi, jqth entry is σ ij " σpB i , B j q, and rµs is the Euclidean vector whose ith component, denoted by rµs i , equals µ i " µpB i q.Then the standard two-dimensional curl operator applied to the vector rµs gives curlrµs " ε ij B i µ j .When this operator is repeated row-wise on a matrix, we get the standard row-wise matrix curl, namely rcurlrσss i " ε jk B j σ ik .The standard Euclidean incompatibility operator [4,1,16] in two dimensions is incrσs " ε qi ε jk B jq σ ik .Using these, we can rewrite the formulas in (4.5) as rcurl g µs " 1 ?det g curlrµs, (4.7a) det g prcurlrσss i ´εjk Γ m ji σ mk q, (4.7b) Note how the expressions for covariant curl and covariant inc contain, but differ from their Euclidean analogues.
Other useful covariant operators include rot g f " ´p‹d 0 f q 7 , f P Ź 0 pT q, (4.8a) pT q, X P XpT q. (4.8b) Clearly, rot g f is in XpT q, while rot g X P T 0 2 pT q.Their coordinate expressions are where, in the latter expressions, the proxy notation has been extended to XpT q and T 0 2 pT q in an obvious fashion to use the Euclidean matrix and vector rot .
It is easy to see that the following integration by parts formula ż pT,gq pcurl g σqpZq " ż pT,gq σprot g Zq `żpBT,gq σpZ, τ q (4.10) holds for all σ P T 2 0 pT q and Z P XpT q. (This can be seen, e.g., using the coordinate expressions (4.9b) and (4.7b) and standard integration by parts on the Euclidean parameter domain.)Here τ is the unit tangent defined in (2.8), the integrals are computed as indicated in (2.17), and σprot g Zq denotes the result obtained by acting the T 2 0 pT q-tensor σ on the T 0 2 pT q-tensor rot g Z. Equation (4.10) shows that rot g can be interpreted as the adjoint of curl g .A similar integration by parts formula for φ P Ź 0 pT q and µ P Ź 1 pT q ż pT,gq connects the other curl g and rot g defined in (4.3a) and (4.8a).

4.2.
Covariant curl in the Regge metric.We proceed to extend the definitions of the covariant operators to the case when the metric g is only tt-continuous (see (2.11)) across element interfaces.Let M denote the open set obtained from M by removing all the mesh vertices (after its triangulation by T ).The topological manifold M can be endowed with a natural glued smooth structure based on the tt-continuity of g, as alluded to in the literature [15,20,30,31], [34,Theorem 3.3] and [50].This glued smooth structure is different from that given by the coordinates x i (see §3.3) in which we plan to conduct all computations.A striking difference is that while g is only tt-continuous in x i , it is fully continuous in the natural glued smooth structure.The glued smooth structure can be understood using the following coordinate construction around any interior mesh edge E. Let z P E and let U z denote an open neighborhood of z not intersecting any other mesh edge or mesh vertex.Let d g p¨, ¨q denote the distance function generated by g on the manifold M .For any p P U z , let πppq " arg min qPE d g pq, pq.We use T ˘, ν˘, τ ˘introduced in (2.15).Let E p denote the submanifold of E connecting z to πppq oriented in the τ `direction.Then for any p P U z , define new coordinates x1 ppq " ˘dg pπppq, pq if p P T ˘, x2 ppq " Denote the coordinate frame of xi by Bi .It can be shown [30] that B1 " ν`a nd B2 " τ àt points in U z XE, and that gp Bi , Bj q is continuous across U z XE for all i, j.Augmenting the set M with the maximal atlas giving such coordinates, we obtain a manifold with the glued smooth structure, which we continue to denote by M .Moreover, p M , gq is a Riemannian manifold with piecewise smooth and globally continuous metric g.
For the next result, we need the subspace of smooth vector fields Because the transformations B i " pBx j {Bx i q Bj and dx i " pBx i {Bx j qdx j are smooth within mesh elements, previously defined piecewise smooth spaces like RpT q carry over to the glued smooth structure.
Boundary mesh edges do not contribute to the last integral since ϕ is compactly supported.Across an interior mesh edge, since g is continuous in the glued smooth structure of M , and since σ is tt-continuous, the contributions to the last integral from adjacent elements cancel each other.Hence (4.13) follows.
We use p M , gq to extend the definition of covariant curl.Recall that the adjoint of curl g is rot g , as shown by (4.10).Hence, as in the theory of Schwartz distributions, a natural extension would be to consider curl g σ, for σ P RpT q, as a linear functional on X c p M q defined by for all ϕ P X c p M q, where we have used Proposition 4.1 in the second equality.The next key observation is that we may extend the above functional to act on a piecewise smooth vector field W P XpT q instead of the smooth ϕ, provided W is g-normal continuous (an intrinsically verifiable property on the manifold).
Definition 4.2.For any σ P RpT q, define curl g σ as a linear functional on Wg pT q, the space defined in (2.16), by for all W P Wg pT q, where the first term on the right hand side is evaluated using the smooth case in (4.3b).The rationale for this definition is that there are functions ϕ in X c p M q approaching W P Wg pT q in such a way that the right hand side of (4.14) converges to that of (4.15): see Proposition A.1 in Appendix A. Furthermore, because of the g-normal continuity of W , the last term in (4.15) vanishes when σ is globally smooth, so (4.15) indeed extends curl g on smooth functions.

Their finite element subspaces of interest are
W k h " tw P WpΩq : for all T P T , w| T " Φ T ŵ for some ŵ where the pull-back Φ T ŵ is the Piola transformation.For k ą 0, (4.16) coincides with the Brezzi-Douglas-Marini finite element space (BDM ) [13] on the parameter domain.
In practice, it is more convenient to work with the BDM space than (2.16).For any w " rw 1 , w 2 s P WpΩq, let .17) Proposition 4.3 (Extended covariant curl in computational coordinates).
(1) A vector field w on Ω is in WpΩq if and only if Q g w is in Wg pT q.
(2) For any σ P RpT q and w P WpΩq, Hence the continuity of g τ τ across element interfaces implies that gpQ g w, νq " 0 if and only if w ν " 0.
Similarly w ν | BΩ " 0 if and only if gpQ g w, νq " 0 vanishes on the boundary.Hence Q g w P Wg pT q if and only if w P WpΩq.
Integrating this over each element boundary using the measure ?g τ τ dl and summing over elements, the right hand sides of (4.15) and (4.18a) are seen to be the same.
To prove the second identity (4.18b), consider any W P Wg pT q.We start by applying the integration by parts formula (4.10) to the first term on the right hand side of the definition (4.15): `żBT σpW, τ q ´żBT gpW, τ q σpν, τ q " ż T σprot g W q `żBT gpW, τ q σ τ τ after simplifying using g-orthogonal decomposition W " gpW, τ qτ `gpW, νqν on element boundaries.Substituting W " Q g w, applying the quotient rule to compute rot g pQ g wq using (4.6), and expressing the result in x i coordinates, (4.18b) follows.
In analogy with the finite element curvature approximation, we may now also lift the functional curl g σ to a finite element space to get a computable representative of the covariant curl on the parameter domain Ω.Using the BDM space in (4.16), we define curl g,h σ, for any σ P RpT q, as the unique element in Wk h satisfying ż Ω δpcurl g,h σ, w h q da " xcurl g σ, Q g w h y WgpT q for all w h P Wk h , where the right hand side can be evaluated using either of the formulas in (4.18).

4.4.
Covariant incompatibility in the Regge metric.To extend the smooth covariant incompatibility defined in (4.3c), we use the space VpT q of piecewise smooth and globally continuous functions on M .For any u P VpT q, the vector field rot g u, by definition (4.8a), satisfies gprot g u, νq " gp´p‹d 0 uq 7 , νq " ´p‹d 0 uqpνq.Hence (2.10) implies gprot g u, νq " ´pd 0 uqpτ q, (4.22) so in particular, gprot g u, νq " 0 due to the continuity of u.Also note that in the formula (4.3c) for the smooth case, inc g σ " curl g curl g σ, the outer curl g 's adjoint is the rot g appearing in (4.11).These facts motivate us to use Definition 4.2 to extend inc g to tt-continuous σ as follows.
Definition 4.4.For any σ P RpT q, extend inc g σ as a linear functional on VpT q by xinc g σ, uy VpT q " xcurl g σ, rot g uy WgpT q (4.23) for all u P VpT q.Note that rot g u is an allowable argument in the right hand side pairing since it is g-normal continuous (and hence in Wg pT q) by (4.22).The next result shows that (4.23) indeed extends the smooth case.
Illustration of the vertex jump defined in (4.24).
The "vertex jump" of σ at a vertex V P V T of an element T P T (cf.[18,28] and see Figure 1), denoted by σ ν τ T V , represents the jump in the value of σpν, τ q across the vertex V when traversing BT in the τ -direction.Alternately, enumerating the vertices of V T as V 0 , V 1 , V 2 so that the indices increase while moving in the τ -direction, naming the edge opposite to V i as E i , and calculating the indices mod 3, we put (4.24) Proposition 4.5.For any σ P RpT q and u P VpT q, When σ is globally smooth, all terms on the right hand side vanish except the first.
When summing over the three edges E i Ă BT , the above vertex values of σ ν τ yield vertex jumps.Hence the first statement follows by summing over all T in T .
To prove the second statement, note that the edge integrals from adjacent triangles cancel each other.To show that the last term with vertex contributions also vanish, let E V denote the set of vertices connected to V P V by an edge and let τ V denote the gunit tangent vector along an E P E V pointing away from V .Then the jump σpν, τ V q on any edge E P E V is defined as before, using (2.15).Its limit as we approach a vertex V along any edge E P E V , is denoted by σ V ν τ .Using it, the last sum can be rearranged to ÿ where V int Ă V is the subset of mesh vertices in the interior of the domain.Each summand on the right hand side vanishes when σ is smooth.
As before, one may now lift the incompatibility functional into a finite element space to get a computable representative.Namely, for any σ P RpT q and g P R `pT q, let inc g,h σ be the unique function in Vk`1 where the right hand side can be computed using the formula in Proposition 4.5.
4.5.Linearization of curvature.Linearization of curvature was discussed in various forms by many previous authors.For example, the linearization of our vertex curvature sources can be guessed from the three-dimensional case presented in [18], while that of the edge and interior curvature sources were derived in [27] (making use of [23]) in terms of the covariant divergence operator.Here we revisit the topic to derive the linearization of edge and interior curvature sources directly in terms of the covariant incompatibility and curl operators.While doing so, we also present different and elementary proofs.
The variational derivative of a scalar function f : SpT q Ñ R in the direction of a σ P SpT q is given by D σ pf pρqq " lim tÑ0 f pρ `tσq ´f pρq t when the limit exists.We use this exclusively for scalar functions of the metric g (i.e., with ρ " g above).Note that the changes in ρppq and σppq as p varies in M are immaterial in the above definition.Hence, we may use Riemann normal coordinates [32] to prove the pointwise identities of tensorial quantities in the next result.Let xi denote the Riemann normal coordinates given by a chart covering a point p P T, Bi denote the corresponding coordinate frame, σij " σp Bi , Bj q for any σ P T 2 0 pT q, Rijkl " Rp Bi , Bj , Bk , Bl q, and let Γijk , Γk ij be defined by (2.5) with B i replaced by Bi .Then, by [32,Proposition 5.11 which greatly simplify calculations.As an example, consider the expression for covariant incompatibility of a σ P SpT q given by (4.7c) in coordinate proxies.It simplifies in normal coordinates, by virtue of (4.26), to rinc g σs " incrσs´ε ij ε kl σml Bi Γm jk .Expanding out the last term in terms of gij , using (2.5b), and simplifying, rinc g σs " incrσs ´1 2 trrσs incrgs.(4.27) Lemma 4.6 (Variations of curvature terms).Consider an element manifold pT, gq for T P T .Let p be an arbitrary point in T and let X, Y P T p M .Let q be any point in one of the edges E of BT and let τ P T q E. Then D σ pKpgqvol T,g pX, Y qq " ´1 2 vol T,g pX, Y q inc g σ, at the point p P T, (4.28a) D σ pκpgqvol E,g pτ qq " 1 2 pcurl g σ `d0 σ τ ν qpτ q, at the point q P E, (4.28b)To prove (4.28b), we start with its right hand side.By (4.7b), At any point on the edge E, without loss of generality, we may choose a Riemann normal coordinate system so that B1 " τ " τ and B2 " ν " ν.Then, using (4.26), the above expression becomes pcurl g σ `d0 σ τ ν qpτ q " rcurlrσss i τ i `Bτ σ τ ν " rcurlrσss i τ i `pB τ rσsq ν τ `pσ ν ν ´στ τ q 9 τν , ( where we used (3.19) to get the last equality.Now we work on the left hand side of (4.28b).Differentiating the expression for geodesic curvature from Lemma 3.3, we get ˙pΓ ν τ τ `9 τ ν q `τ i τ j ν l δ kl pg km Γ ijm pσq ´gka σ ab g bm Γ ijm q  in general coordinates.Specializing to the previously used Riemann normal coordinates, applying (4.26) and simplifying Γ ijm pσq, D σ pκ g pgqvol E,g pτ qq " ´1 2 trrσs ´rσs τ τ ¯9 τ ν `pB τ rσsq ν τ ´1 2 p Bν rσsq τ τ .
This coincides with the expression in (4.29), so (4.28b) is proved.

Connection approximation
In this section we approximate the Levi-Civita connection when only an approximation to the true metric is given, namely g P R `pT q.To do so, we assume we are given a g-orthonormal frame pe 1 , e 2 q in each T P T .Then, the connection is fully determined by a single connection form pg; ¨q " g P Ź 1 pT q, within each element T , given by g pXq " gpe 1 , ∇ X e 2 q " ´gp∇ X e 1 , e 2 q (5.1) for any X P XpM q.This section is largely based on [10] (so we will be brief), but we note that while they approximate the Hodge star of g , we approximate g directly (and also note that the orientation in their work is opposite to ours).
To extend the connection to accommodate the possible discontinuities of the frame pe 1 , e 2 q across element interfaces, let g pa, bq denote the counterclockwise angle from b to a measured in the g-inner product, for any two vectors a, b P T p M .This angle is well defined even for points p on a mesh edge E (excluding the vertices) since we use the glued smooth structure (see (4.12)) in which g is continuous across the edge.On each interior mesh edge E, let T ˘, ν˘, τ ˘be as in (2.15), orient the edge E by τ E " τ `, and put νE " ν`, e ˘,i " e i | T ˘.Let Θ E " g pe `,1 , e ´,1 q; see Figure 2. (It is possible to compute this angle without resorting to the coordinates in (4.12), as we will explain later in Appendix B.) This is the angle by which a vector must be rotated while parallel transporting it across the edge E in the νE direction.To account for this rotation, we extend g as follows: Definition 5.1.Given g P R `pT q and g-orthonormal piecewise smooth frame e 1 , e 2 P XpT q, define g P Wg pT q 1 by for all W P Wg pT q.
Angle between frames on different elements.
Within each element, the well-known identity d 1 g " Kpgqvol T,g holds.Equivalently, using the curl in (4.3a), curl g p g | T q " Kpg| T q for each T P T .To speak of curl g g for the functional g in (5.2), we must extend curl g .Motivated by (4.11), we define xcurl g µ, uy VpT q " xµ, rot g uy WgpT q , for all u P VpT q, µ P Wg pT q 1 . (5. 3) The right hand side is well defined for µ P Wg pT q 1 since rot g u P Wg pT q by (4.22).Next, for each V P V and E P E V , let s EV equal `1 if τ E points towards V and ´1 otherwise.Following [10], we assume that at each interior vertex V , the "consistency" condition ÿ holds.It can be seen from the proof of [10, Proposition 5.4] that the left hand side above always equals 2πm for some integer m.The condition (5.4) requires that e i be chosen so as to achieve m " 1. (We'll give a recipe for doing this shortly: see (5.7) below.) Proposition 5.2.Let K g be as in (3.2), g be as (5.2) for a g-orthonormal frame e i satisfying (5.4), and curl g g be as given by (5.3).Then curl g g " K g .
Proof.This was proved in [10], so we will merely indicate how to apply their result to xcurl g g , uy VpT q " x g , rot g uy WgpT q " ż T g prot g uq `ÿ EPE int ż pE,gq Θ E gprot g u, νE q.
A computable representative of the connection form is obtained by lifting the g into the BDM finite element space (defined in (4.16)) on the parameter domain, as follows.

Definition 5.3 (Connection 1-form approximation). Define h pgq as the unique function in Wk
h such that for all v h P Wk where Q g is as Proposition 4.3 and the right hand side is evaluated using (5.2).
In the remainder, we assume that h pgq is computed using a specific g-orthonormal frame pe 1 , e 2 q satisfying (5.4), that we describe now.We start with a globally smooth δ-orthonormal Euclidean basis pE 1 , E 2 q on the parameter domain (e.g., the standard unit basis on R 2 ).Then, this basis is continuously transformed to a g-orthonormal frame as in the next lemma.Let Gptq " δ `tpg ´δq.We consider the ordinary differential equation (ODE) Gptq ´1pg ´δq, up0q " I. (5.6) Lemma 5.4.The solution of ODE (5.6) is given by uptq " Gptq ´1{2 .The frame puptqE 1 , uptqE 2 q is Gptq-orthonormal, so at t " 1, it is g-orthonormal.
Proof.Let g " V ΛV ´1 be a diagonalization of g with eigenvalues λ i in Λ " diagpλ i q and eigenvectors in the orthogonal matrix V .Then Gptq " V pp1 ´tqδ `tΛqV ´1 and Gptq ´1{2 " V diag `1 `tpλ i ´1q ˘´1{2 V ´1.Using these expressions, the statements of the lemma can be easily verified.
We use the g-orthonormal frame e i " up1qE i " g ´1{2 E i (5.7) for computations.As stated in [10], since the frame E i obviously satisfies (5.4) with g " δ, the chosen e i obtained by the continuous deformation of the metric and frame, satisfies (5.4).In Appendix B we present a stable algorithm for computing the right hand side of (5.5) using the discontinuous metric g in the computational coordinates.

Error analysis
In this section, we prove a priori estimates for the error in the previously defined approximations.We restrict ourselves to the model case (see §3.3) and work on the parameter domain Ω, where we shall use the Euclidean dot product u ¨v " δpu, vq and the standard Frobenius inner product A : B between matrices A, B. We assume that the triangulation T consists of affine-equivalent elements, is shape-regular, and is quasi-uniform of meshsize h :" max T PT diampT q.
6.1.Convergence results.All results here concern the canonical interpolant into the Regge space R k h (defined in (2.13)) of degree k ě 0. This well-known interpolant [34], k σq τ τ q dl " ż E σ τ τ q dl for all q P P k pEq and edges E of BT , (6.1a) Note that when ρ is a skew-symmetric matrix, both sides of (6.1b) vanish, so (6.1b) is nontrivial only for symmetric ρ.Throughout, we use standard Sobolev spaces W s,p pΩq and their norms and seminorms for any s ě 0 and p P r1, 8s.When the domain is Ω, we omit it from the norm notation if there is no chance of confusion.We also use the element-wise norms }u} p W s,p h " ř T PT }u} p W s,p pT q , with the usual adaption for p " 8.When p " 2, we put Most of our results assume that k is an integer, k ě 0, g " I R k ḡ, ḡ P W 1,8 pΩ, S `q, ḡ ´1 P L 8 pΩ, S `q.(6.2) We use a b to indicate that there is an h-independent generic constant C, depending on Ω and the shape-regularity of the mesh T , such that a ď Cb.The C may additionally depend either on t}ḡ} W 1,8 , }ḡ ´1} L 8 , pḡq, Kpḡqu, or on t}g} W 1,8 h , }g ´1} L 8 u, depending on whether we assume (6.2) or not, respectively.We use p¨, ¨qD to denote the integral of the appropriate inner product (scalar, dot product, Frobenius product, etc.) of its arguments over a Euclidean measurable set D, e.g., pσ τ τ , qq E and pσ, ρq T equal the right hand sides of (6.1a) and (6.1b), respectively.Theorem 6.1 (Approximation of covariant curl).Suppose g P R `pT q, σ h P R k h , v h P Wk h (the BDM space in (4.16)), σ P H 1 pΩ, Sq X C 0 pΩ, Sq, and let curl g,h be as defined in (4.21).Then, there exists an h 0 ą 0 such that for all h ă h 0 , and if (6.2) holds, Proofs of this and other theorems in this subsection are presented in later subsections.For now, let us note that on Euclidean manifolds with ḡ " δ, the expressions of our distributional covariant curl (either (4.15) or (4.18a)) reduce to It is easy to see from (6.1) (and integrating the right hand side above by parts) that pcurl δ,h pσ ´IR k σq, v h q Ω " 0 (6.6) for all v h P Wk h .This equality has the flavor of typical FEEC identities (also known as commuting diagram properties).On general manifolds however, it appears that we must trade this equality for the inequality (6.3).The remaining inequalities (6.4)-(6.5)bound the nonlinear changes in the covariant operator arising due to the perturbations in the metric.Theorem 6.1 directly implies error bounds in L 2 norm, while error bounds in stronger norms follow from it: Corollary 6.2.Under the assumptions of Theorem 6.1, for all 1 ď l ď k, Similar results can be proved for the incompatibility operator.
Theorem 6.3 (Approximation of covariant incompatibility operator).Suppose g P R `pT q, σ h P R k h , u h P Vk`1 h (the Lagrange space in (2.14)), σ P H 1 pΩ, Sq X C 0 pΩ, Sq, and let inc g,h be as defined in (4.25).Then, there exists an h 0 ą 0 such that for all h ă h 0 , and if (6.2) holds, Here again, as in the case of covariant curl, comparison with the Euclidean case is illuminating.On Euclidean manifolds, instead of (6.8), the stronger result pinc δ,h pσ ´IR k σq, u h q Ω " 0 (6.11) holds for the distributional incompatibility (which has element, edge, and vertex contributions: see Proposition 4.5).Indeed, (6.11) follows immediately from (6.6) and (4.23).The theorem also implies error bounds in stronger norms.
Corollary 6.4.Under the assumptions of Theorem 6.3, for all 0 ď l ď k, Our remaining results are for the approximations of connection and curvature.Let I denote the identity operator (on some space that will be obvious from context) and let Π W k and Π V k`1 denote the L 2 -orthogonal projection into Wk h and Vk`1 h , respectively.Theorem 6.5 (Approximation of Gauss curvature).Suppose (6.2) holds, ḡ P W k`1,8 pΩq, Kpḡq P H k pΩq, and K h pgq P Vk`1 h be as in (3.4).Then, there exists an h 0 ą 0 such that for all h ă h 0 , Corollary 6.6.Under the assumptions of Theorem 6.5, for all 0 ď l ď k, Theorem 6.7 (Approximation of Levi-Civita connection).Suppose (6.2) holds, ḡ P H k`1 pΩq, pḡq P H k`1 pΩq, and let h pgq P Wk h be as in (5.5).Then, there exists an h 0 ą 0 such that for all h ă h 0 , and when k ě 1, Corollary 6.8.Under the assumptions of Theorem 6.7, for all 1 ď l ď k, Since the curvature Kpḡq has second order derivatives of the metric ḡ, at first glance it may seem surprising that Theorem 6.5 gives H ´1-convergence of curvature approximations at the same rate as ~ḡ ´g~8.Even for the lowest order case k " 0 (while using piecewise constant metric approximations), where one might only expect convergence in the H ´2-norm, the theorem gives first order convergence of the curvature in the H ´1-norm.The convergence rates of Theorems 6.7 and 6.5 are both higher than those proved in [27,10].As we shall see, the reason behind these higher rates is a property of the Regge interpolant proved in the next subsection.6.2.Distributional Christoffel symbols of the first kind.In a neighborhood where the metric g is smooth, the Christoffel symbols of the first kind, Γ lmn pgq, are given by (2.5b).To see what further terms must be supplied to obtain their distributional version when the metric g is only tt-continuous across an element interface, consider a ψ lmn in the Schwartz test space DpΩq of smooth compactly supported functions.Since Γ lmn p¨q is a linear first order differential operator applied to a smooth metric argument, its distributional definition is standard: pψ νlm `ψlνm ´ψlmν q dl ˙.
We can split the integrand over BT into νν, τ ν, ντ , and τ τ components.When summing over BT for all T P T , the tt-continuity of g implies that the τ τ -terms cancel out.The remaining terms give the boundary contribution as ř ş BT pg νν ψ ννν `2g ντ ψ νντ q dl, so the result follows.Proposition 6.9 serves to motivate the introduction of Γpg, Σq :" for piecewise smooth g P R `pT q and Σ P C 8 pT , R 2ˆ2ˆ2 q.As we proceed to analyze distributional covariant operators, it is perhaps not a surprise that this quantity will reappear in our analysis with various arguments Σ, including those in P k pT , R 2ˆ2ˆ2 q " tΣ : Σ| T P P k pT, R 2ˆ2ˆ2 q for all T P T u.The next result gives a property of Γp¨, ¨q in connection with the Regge interpolation error.Lemma 6.10.If k, ḡ, g are as in (6.2), then for any Σ h P P k pT , R 2ˆ2ˆ2 q, Γpḡ ´g, Σ h q " 0. ( Proof.We start with (6.15) on T P T and integrate by parts The first and second inner products above vanish by (6.1b) and (6.1a), respectively.
6.3.Basic estimates.We need a number of preliminary estimates to proceed with the analysis.The approximation properties of the Regge elements are well understood.By the Bramble-Hilbert lemma, on any T P T , }pI ´IR 0 qḡ} L p pT q h|ḡ| W 1,p pT q (6.17a) for ḡ P W 1,p pT, Sq and p P r1, 8s, for the lowest order case (and certainly for the higher k ě 1).A similar estimate holds for the element-wise L 2 pT q-projection into the space of constants, which we denote by Π 0 : }pI ´Π0 qf } L p pT q ď hC|f | W 1,p pT q (6.17b) for f P W 1,p pT q.Since g " I R k ḡ approaches ḡ as h Ñ 0, we will tacitly assume throughout that h has become sufficiently small to guarantee that the approximated metric g is positive definite throughout.Let E Ă BT be an edge of T .We also need the following well-known estimates that follow from scaling arguments: for all u P H 1 pT q and for all u P P k pT q, }u} L 2 pEq h ´1{2 }u} L 2 pT q , }∇u} L 2 pT q h ´1}u} L 2 pT q .(6. 19) Staying in the setting of (6.2), the following estimates are a consequence of [27, Lemma 4.5 and Lemma 4.6]: for p P r1, 8s, }g ´1 ´ḡ ´1} L p }g ´ḡ} L p , (6.20) and for all x in the interior of any element T P T and for all u P R 2 , Better control of differences of some functions of the metric is possible through the next lemma.Let Lemma 6.11.In the setting of (6.2), for sufficiently small h, there exist smooth functions f 1 , f 2 P C 8 pS `, Rq such that at each point on BT, β 1 pḡq ´β1 pgq " pḡ ´gq τ τ f 1 pgq `η1 pḡ, gq ` 2 1 , β 2 pḡq ´β2 pgq " pḡ ´gq τ τ f 2 pgq `η2 pḡ, gq ` 2 2 , where maxp 1 , 2 q " Op}ḡ ´g}q for some Euclidean norm } ¨} at the point.
The proof of (6.26) is similar.6.4.Analysis of the covariant curl approximation.This subsection is devoted to proving Theorem 6.1.We will start with the first estimate of the theorem, which is easier to prove.The remaining inequalities will be proved using Lemma 6.10 afterward.Lemma 6.13.Suppose g P R `pT q, σ P H 1 pΩ, Sq X C 0 pΩ, Sq, σ h " I R k σ, and v h P Wk h .Then, pcurl g,h pσ ´σh q, v h q Ω ~σ ´σh ~2 }v h } L 2 .
Proof.Using (4.21) and (4.18b), ´`pσ ´σh q τ τ , v i h g iτ β 1 pgq{g τ τ ˘BT ı where β 1 is as in (6.24).The first term on the right is zero when k " 0. When k ě 1, we use (6.1b) to insert a projection to constants: where we also used (6.17b) and the inverse estimate (6.19).The second term in (6.27) is easily bounded by absorbing the maximum of g-dependent Γ k ij into a generic constant: Finally, for the last (boundary) term in (6.27), we use (6.25) of Lemma 6.12 for each edge E Ă BT , setting Ψ " β 1 pgqg iτ {g τ τ P W 1,8 pT q, after extending the constant tangent vector τ from E into the element T .Then `pσ ´σh q τ τ , v i h g iτ β 1 pgq{g τ τ ˘BT p}σ ´σh } L 2 pT q `h|σ ´σh | H 1 h pT q q}v h } L 2 pT q , where we have absorbed the norm }Ψ} W 1,8 pT q from the lemma into the generic gdependent constant in the inequality.Thus and the result follows by applying Cauchy-Schwarz and Young inequalities.
Proof of Theorem 6.1.Error estimate (6.3) directly follows from Lemma 6.13.Using the estimates of Lemma 6.15 to bound the Γ-terms in the inequalities of Lemma 6.14, the remaining error estimates of the theorem follow.
Proof of Corollary 6.2.The proof follows along the lines of [27, pp. 1818], where the error is compared to the weaker L 2 -norm using the Scott-Zhang interpolant, inverse estimates, and the triangle inequality.6.5.Analysis of the covariant incompatibility.The error analysis here can now be given by a simple argument using Theorem 6.1.
We conclude this short subsection with a few remarks on our analysis so far.To compare our analysis with [27], the easily spotted difference is that we use our operator inc g instead of the operator div g div g in [27].These two operators are closely related in two dimensions (see Appendix C).A more substantial difference is that while [27] separately estimates certain element terms and inter-element jumps (by applying a triangle inequality first), we do not.Instead, we kept such terms together, gathered terms of good convergence rates, and identified the remainder as a collection of terms that look like those arising from the formula for distributional Christoffel symbols of first kind.The next key insight was (6.16), which zeroed out the latter collection.Also notable from our analysis so far is the idea of splitting the error terms into a high-order ones and ones that might be sub-optimal in general, but vanishes in specific cases.Such an idea was also used in [49], where the convergence of a surface div div operator on an approximated triangulation is proven (in an extrinsic manner) to converge.6.6.Analysis of the curvature approximation.Now we turn to the proof of Theorem 6.5.The key extra ingredient here is a technique pioneered in [27] to represent the curvature approximation using an integral of its first variation, as stated in the next lemma.For any g P R `pT q and σ P RpT q, using inc g of Definition 4.4, let b h pg, σ, uq " ´xinc g σ, uy VpT q , Gptq " δ `tpg ´δq, Ḡptq " δ `tpḡ ´δq.Lemma 6.16.Let g P R `pT q, K h pgq P Vk`1 h be as in Definition 3.1 for some k ě 0. Also let ḡ be a smooth metric and let Kpḡq denote its smooth Gauss curvature.Then Proof.Since K h pδq " 0, Now, expand the inner integral using (3.4) and differentiate each term in the direction σ " G 1 ptq " g´δ, using each of the three identities of Lemma 4.6.Comparing the result, term by term, to the expression in Proposition 4.5, equation (6.35) is proved.The proof of (6.36) is similar, after noting that the global smoothness of ḡ ´δ implies that the jump terms in b h p Ḡptq, ḡ ´δ, uq vanish (by the last statement of Proposition 4.5).
To state a simple lemma before the error analysis, for any u, v, f P L where we have used the equivalence of L 2 pΩq-norm and } ¨}ḡ -norm.
Proof of Theorem 6.5.The general structure of the proof follows [27].Let u P H 1 0 pΩq and let u h P Vk`1 h .Then pK h pgq ´Kpḡq, uq ḡ " pK h pgq ´Kpḡq, u ´uh `uh q ḡ " s 1 `s2 `s3 where s 1 " pK h pgq, u h q g ´pKpḡq, u h q ḡ, s 2 " pK h pgq, u h q ḡ ´pK h pgq, u h q g , and s 3 " pK h pgq ´Kpḡq, u ´uh q ḡ.We proceed to estimate each s i .
6.7.Analysis of the connection approximation.The integral representation of the connection 1-form can be given analogously to the curvature case once we know the variation of the connection with respect to the metric.Recall that we compute the connection using the canonical g-orthonormal frame e i ptq " Gptq ´1{2 E i (see (5.7)) obtained using the flow (5.6) evaluated at t " 1.We only discuss the variation of g with respect to g along this flow (i.e., unlike Lemma 4.6, the following is valid only for a specific direction σ).
Let c h pg, σ, vq " ´xcurl g σ, Q g vy WgpT q for v P WpΩq, where Q g is as in (4.17) and the distributional covariant curl is as defined in (4.18).Lemma 6.19.Let g P R `pT q, k P N 0 .There holds for h pgq P Wk h from Definition 5.3, the exact metric ḡ and its connection 1-form pḡq Proof.The proof is similar to proof of Lemma 6.16: (6.41) follows from Lemma 6.18 and the fundamental theorem of calculus, and for (6.42), we note that the jump terms in c h p Ḡptq, ḡ ´δ, vq vanish due to the global smoothness of ḡ ´δ.
Proof of Theorem 6.7.For v P L 2 pΩ, R 2 q we add and subtract the L 2 -orthogonal projector v h " Π W k v P Wk h and write where s 1 " p h pgq´ pḡq, v´v h q Ω and s 2 " p h pgq´ pḡq, v h q Ω .Due to the properties of the L 2 -orthogonal projector we obtain For s 2 , we use Lemma 6.19 and (4.21) to get pcurl Ḡptq,h pg ´δq ´curl Gptq,h pg ´δq, v h q Ω ´1 2 pcurl Ḡptq,h pg ´ḡq, v h q Ω .
Proof of Corollary 6.8.This is analogous to proof of Corollary 6.2.Remark 6.20 (Case k " 0 for connection approximation).In the case k " 0, the space W0 h consists of piecewise constant vector fields with normal continuity.Thus every function in W0 h is exactly divergence-free.If the exact connection 1-form ḡ is not divergence-free, then it cannot generally be approximated by functions in W0 h .Thus, no convergence for the connection approximation in the lowest order case k " 0 should be expected.This is confirmed by numerical experiments in the next section.

Numerical examples
In this section we confirm, by numerical examples, that the theoretical convergence rates from Theorem 6.5 and Theorem 6.7 are sharp.All experiments were performed in the open source finite element software NGSolve1 [42,43], where the Regge elements are available.
To test also the case of non-homogeneous Dirichlet and Neumann boundary conditions we consider only one quarter Ω " p0, 1q ˆp0, 1q and define the right and bottom boundaries as Dirichlet and the remaining parts as Neumann boundary.To avoid possible super-convergence properties due to a structured grid, we perturb all internal points of the triangular mesh by a uniform distribution in the range r´h 4 , h 4 s, h denoting the maximal mesh-size of the originally generate mesh.The geodesic curvature on the left boundary is exactly zero, whereas on the top boundary we compute x 2 px 2 ´3q 2 `y2 py 2 ´3q 2 `9.
The vertex expressions K V at the vertices of the Neumann boundary can directly be computed by measuring the angle arccosp ḡpτ 1 ,τ 2 q }τ 1 }ḡ}τ 2 }ḡ q.To illustrate our theorems, we must use g " I R k ḡ.In implementing the Regge interpolant, the moments on the edges have to coincide exactly: see (6.1).To this end,  we use a high enough integration rule for interpolating ḡ for minimizing the numerical integration errors.
We compute and report the curvature error in the L 2 -norm, namely }Kpḡq´K h pgq} L 2 .We also report the H ´1-norm of the error.It can be computed by solving for w P H 1 0 such that ´∆w " Kpḡq ´Kh pgq and observing that }Kpḡq ´Kh pgq} H ´1 " }w} H 1 .Of course the right hand side can generally be computed only approximately.To avoid extraneous errors, we approximate w using Lagrange finite elements of two degrees more, i.e., w h P V k`3 h when K h pgq P V k`1 h .We start by approximating ḡ by the lowest order piecewise constant Regge elements g P R 0 h .As shown in Figure 5 (left), we do not obtain convergence in the L 2 -norm, but do obtain linear convergence in the weaker H ´1-norm, in agreement with Theorem 6.5.When increasing the approximation order of Regge elements to linear and quadratic polynomials we observe the appropriate increase of the convergence rates: see Figure 5 (middle and right), confirming that the results stated in Theorem 6.5 and Corollary 6.6 are practically sharp.In Figure 6 snap-shots of the approximated Gauss curvature are shown.
AppA`9q ?A`6Aq ‚, where A " y 6 ´6y 4 `9 `x6 ´6x 4 `9x 2 , E i P R 2 are the Euclidean basis vectors, e i " ḡ´1 2 E i in accordance with (5.7), and the covariant derivative p ∇X Y q i " pp ∇XqY q i Γi jk X j Y k with Γk ij denoting the Christoffel symbol of second kind with respect to ḡ.For the results shown in Figure 8 we use BDM and also Raviart-Thomas [39]RT elements.The optimal L 2 -convergence rates stated in Theorem 6.7 are confirmed for k ą 0 to be sharp when using BDM elements for k ą 0. If we increase the test-space, however, to RT elements, which additionally include specific polynomials of one order higher than BDM , one order of convergence is lost, compare also Figure 8. Note, that to construct the finite element space W 0 h " BDM 0 , we consider the lowest order Raviart-Thomas elements RT 0 and lock the linear part by enforcing that divpRT 0 q " 0. As depicted in Figure 8 (left) the discrete solution converges linearly at the beginning,  however, after some refinements the error stagnates.This is in accordance with the explanation provided in Remark 6.20.Solution snap-shots are displayed in Figure 9.
Appendix A. Rationale for the g-normal continuity This section briefly presents a justification for the usage of g-normal continuous vector fields in Definition 4.2.We show that there are smooth functions ϕ in X c p M q approaching a g-normal continuous W P Wg pT q in such a way that the right hand side of (4.14) converges to that of (4.15).For any mesh vertex V P V , let B ε pV q " tq P M : d g pq, V q ď εu.Then put D ε " Y V PV B ε pV q and M ε " M zD ε .Let U i , Φi : U i Ñ R 2 denote a chart of the glued smooth structure.In the parameter domain Φi pU i q, using the Euclidean divergence operator, define W p pdiv, Φi pU i qq " tw P L p p Φi pU i qq, divpwq P L p p Φi pU i qqu for any 1 ď p ă 8 with its natural Euclidean norm.This norm and the duality pairings defined in (4.14) and (4.15) feature in the next result.
Proposition A.1.Let σ P RpT q and W P Wg pT q.For any given ε 1 ą 0, there exists a p ą 2, an ε 2 ą 0, finitely many charts tpU i , Φi q : i P Iu covering M ε 2 in the glued element boundary terms, ż BT gpχ ε 2 W ´ϕ, νq σpν, τ q " ÿ iPI ż BT gpW i ´ϕi , νq σpν, τ q, (A.7) we focus, as before, on a neighborhood U i intersecting an edge E " BT ´XBT `.On BT `, by (A.4), gpW i ´ϕi , νq " W 1 i ´ϕ1 i yields the normal component of the pushforward p Φi q ˚pW i ´ϕi q on Y -axis in the parameter domain.Since the latter converges to zero in W p pdiv, Φi pU i qq by (A.5), its normal trace converges to zero in W ´1{p,p p Φi pE X U i qq.Choose p ą 2 and q such that 1{p `1{q " 1.When σpν, τ q is mapped to Φi pE X U i q and extended to Y -axis by zero, is in W 1´1{q,q p Y q since 1 ´1{q ă 1{2.Hence the contribution from U i X E to the right hand side of (A.7) vanishes.Repeating this argument on other charts, (A.6) is proved.

Appendix B. Angle computation for connection approximation
In this section we discuss and present a stable angle computation used in the connection 1-form approximation.Since g ij is in general discontinuous across an edge E (in the computational coordinates x i ), we cannot use it to directly compute the angle between frame vectors on two different triangles.Instead, we compute angles the frame makes with an intermediate vector element by element and then use it to compute Θ E as explained below.
Before going into details we make the following observation: if the metric g approximates a smooth metric ḡ, we expect that the frame e i fixed by (5.7) will be such that their restrictions to adjacent elements, pe `,1 , e `,2 q and pe ´,1 , e ´,2 q, will differ by a small angle, say less than π.This is the case in Figure 10 (left), where the angle difference Θ E is negative and |Θ E | ă π.
On each interior mesh edge E, let T ˘, ν˘, τ ˘be as in (2.15), orient the edge E by τ E " τ `, and put νE " ν`, e ˘,i " e i | T ˘.Let Θ E ˘" g pe ˘,1 , ˘ν ˘q.Clearly, Θ E ˘can be computed using g| T ˘, specifically using its components g ij in the computational coordinates x i on either triangle.Then Θ E " Θ E `´Θ E ´in most cases (and certainly in the case illustrated in left drawing of Figure 10).In some cases however, such as that in the middle drawing of Figure 10, although Θ E is negative and |Θ E | ă π, the number Θ E `´Θ E ´is positive and larger than π.Thus setting Θ E " Θ E `´Θ E ´would be incorrect and would lead to a bad numerical approximation of the connection 1-form.
To cure this problem we change the choice of the starting angle on the fly.First, we compute Θ E ˘" g pe ˘,1 , ˘ν ˘q on every edge as a pre-processing step.On each edge, set ´.In other words, we compute after reversing the sign of the artifically introduced g-normal vector ν˘o n both the adjacent elements of an edge if the modulus of the pre-computed angle is larger than π.This is illustrated in Figure 10 (right), where the sign change of the normal vector is depicted in red.The following computational formula is easy to prove.By symmetrization of (5.1) we have pg; e i q " 1 2 g jk pB i e j 2 `Γj li e l 2 e k 1 ´Bi e j 1 ´Γj li e l 1 e k 2 q and the claim follows.

Appendix C. Relation between distributional covariant inc and divdiv
In this section we show that the distributional covariant inc from Proposition 4.5 and the covariant distributional divdiv operator of a rotated sigma used in [10] coincide in the sense xdiv g div g S g σ, uy VpT q " ´xinc g σ, uy VpT q for all u P VpT q, (C.1) where the covariant divergence is defined below.This is in common with the Euclidean identity divdivSσ " divdivg ´∆ trpσq " ´inc σ, trp¨q denoting the trace of a matrix.The distributional covariant divdiv reads xdiv g div g S g σ, uy VpT q " ż T u div g div g S g σ `żBT u ppdiv g S g σq 5 pνq `pd 0 σ ν τ qpτ qq where S g σ " σ ´tr g pσqg with tr g pσq " σ ij g ij .(Note that the authors in [10] used pν, τ q as positively oriented frame, whereas we use pτ , νq such that the signs in the boundary and vertex terms differ.Also, a different orientation in the vertex jump is used.) The covariant divergence is defined as the L 2 -adjoint of the covariant gradient.For f P Ź 0 pΩq its covariant gradient is given by the equation det g B i p a det gv i q, v P XpΩq.
Lemma C.1.There holds for g P R `pT q and σ P RpT q (1) ‹pdiv g S g σq 5 " ´curl g σ and pdiv g S g σq 5 pνq " pcurl g σqpτ q on BT , (2) div g div g S g σ " div g div g σ ´∆g tr g pσq " ´inc g σ on T , where ∆ g f :" div g grad g f , f P Ź 0 pΩq denotes the Laplace-Beltrami operator.
Proof.In the first identity only first order derivatives of g are involved.Thus, in normal coordinates, xi (see §4.5), pdiv g S g σq 5 becomes the Euclidean version divpσ ´trpσqδq, which reads in components pdiv g S g σq 5 " ˆB 2 σ12 ´B 1 σ22 B1 σ12 ´B 2 σ11 ˙.
The identity div g div g S g σ " div g div g σ ´∆g tr g pσq is well known [27], so we focus on proving its relationship with inc, by means of normal coordinates.The Laplace-Beltrami operator becomes r∆ g tr g pσqs " rdiv g ∇ g tr g pσqs " rdiv g pg ij B j tr g pσqqs " B2 i trpσg ´1q " ∆ trpσq ´trpσ∆gq and the divdiv part becomes rdiv g div g σs " divdivrσs ´2B 2 ij gik σkj `B i Γlji σlj ´B i Γjlj σil .
Note that we abused notation and summed over repeated indices all of which are subscripts (forgivable while using normal coordinates).Furthermore, by inserting the

b
h pGptq, g ´δ, u h q dt, for all u h P Vk`1 h , Ḡptq, ḡ ´δ, uq dt, for all u P VpT q.(6.36)

Figure 3 .
Figure 3. Left: Embedded surface, color indicates to z-component.Right: Exact Gauss curvature as graph over the domain Ω.

Figure 4 .
Figure 4. Left: Domain with Dirichlet and Neumann boundaries.Middle and right: perturbed unstructured triangular mesh grids.

Figure 5 .
Figure 5. Convergence of Gauss curvature with respect to number of degrees of freedom (ndof) in different norms for Regge elements g P R k h