A note on the rate of convergence of integration schemes for closed surfaces

In this paper, we issue an error analysis for integration over discrete surfaces using the surface parametrization presented in [PS22] as well as prove why even-degree polynomials exhibit a higher convergence rate than odd-degree polynomials. Additionally, we provide some numerical examples that illustrate our findings and propose a potential approach that overcomes the problems associated with the original one.


Introduction
Many applications, including mathematical physics and mathematical biology [YQ09], require accurate approximation of integrals on curved surfaces.The concept of surface integration is a fundamental procedure in a wide range of numerical methods, including the boundary integral method, the finite element method, the surface finite element method, and the finite volume method.As a result of piecewise linear approximations of the surface and integrand, integration over discrete surfaces (such as surface triangulation) is typically only of first-or second-order accuracy [Geo98].To improve the order of accuracy, using the approach presented in [PS22] a polynomial approximation of the geometry of a smooth closed embedded hypersurface M is considered, where a local interpolation polynomial of the closest point projection map π for each element T ∈ T 1 h will be constructed, with T 1 h denoting a conforming triangulation of M.
The contribution of this paper is twofold.First, we deliver a parametrization of smooth, closed surfaces, enabling to numerically approximate surface integrals based on discrete triangulations.That is approximately providing a decomposition M = n i=1 V i into non-overlapping regions V i , given as images of regular maps, see Fig. (1) where σ is a reference simplex defined as σ := {(s, t) , 0 ≤ s ≤ 1, 0 ≤ t ≤ 1 − s} in R 2 .Consequently, the surface integration problem becomes: where g i (s, t) = det(J T J) with J = [∂ s ψ i (s, t), ∂ t ψ i (s, t)] and ∂ s ψ i (s, t), ∂ t ψ i (s, t) denoting the partial derivatives of ψ i with respect to s and t.We propose to approximate the function f and each parametrization ψ i by polynomials of degree k ∈ N. Consequently, the integral is approximated by numerically computing where gi (s, t) = det( JT J) with J = [∂ s Q ψ i ,k (s, t), ∂ t Q ψ i ,k (s, t)], and Q ψ i ,k denoting a k−th order polynomial approximating the map ψ i , whereas Q f,k is a k−th order polynomial approximating the integrand f .The second part of our contribution is motivated by recent results of Ray et al. [RWJG12], where a stabilized least squares approximation, a blending procedure based on linear shape functions, and high-degree quadrature rules are combined into a method for integration over discrete surfaces.This method has an accuracy order of O h k + h 5 with k denoting the polynomial degree and h > 0 the mesh size ([RWJG12], Theorem 1).Based on their numerical experiments, they observe that evendegree polynomials exhibit a higher convergence rate than predicted by their theory.The aim of this paper is to publicize this phenomenon and expand our understanding of it with the aid of new numerical experiments (section 4) and a new Theorem 3.3, which also provides a detailed analysis of errors for integration over discrete triangulated surfaces.
The article is structured as follows.Section 2 gives a mathematical description of the surface parametrization introduced in [PS22].In Section 3, we analyze the error using this approach and justify why even-degree polynomials exhibit higher convergence rates.In section 4 we conclude with some numerical examples to illustrate our findings and propose an alternative approach that overcomes some of the disadvantages of the original method.Codes for reproducing the results of this manuscript are summarized and available in the repository https://github.com/zavala92/code_paper.

Polynomial Approximation for Closed Surfaces
Let M be a smooth connected, orientable closed hypersurface, with smooth hypersurface we mean a C ∞ topological manifold that is second-countable, Hausdorff, and locally Euclidean of dimension 2, which is embedded in an ambient space of dimension 3.According to the Jordan-Brouwer decomposition theorem [Lim88], M divides R 3 into an interior and an exterior domain and we denote by d the signed distance function to M oriented in such a way that d > 0 in the exterior, d < 0 in the interior of M. Let us denote the outward unit normal of M with n(x) = ∇d(x), where ∇ is the standard gradient in R 3 (see [Dem09] for more details).Given M ⊂ R 3 for each point x ∈ M the tangent space T x M ⊆ R 3 is a linear subspace.The disjoint union of tangent spaces to M is the tangent bundle of M: The normal space of M in R 3 at x is the orthogonal complement of the tangent space T x M, namely Let δ : M −→ R + be a positive, continuous function.Consider the following open tubular neighborhood of the normal bundle: where y x is normal vector attached at x.The map F : N M −→ R 3 (x, y) → x + y is smooth and there exists a δ such that the restriction F| N δ becomes a diffeomorphism onto its image [Lee13].Consequently, open, smooth, embedded submanifold of R 3 that forms a tubular neighborhood of M.
We will assume that we have a polyhedral surface M h in Euclidean three-space, which is defined to be a compact subset M h ⊆ R 3 homeomorphic to M. This discrete surface is composed of finitely many triangles whose vertices are located in M, ensuring that each edge is contained in a certain (affine) line and each face is contained in a certain (affine) plane: where each triangle T i is parameterized over a reference and T 1 h is a collection of flat triangles with mesh size h = max h is called a conforming triangulation if for any T i , T j ∈ T 1 h with T i ̸ = T j the intersection T i ∩ T j is either empty or a proper k-sub-simplex of T i (k < 2).Let assume that M h is contained in the tubular neighborhood N δ .Under these conditions, we define a unique nonlinear closest point projection map: which assigns to every x ∈ M h the closest point on M, so that (x − π(x)) ⊥ T π(x) M, ∀x ∈ M h .We assume that N δ has the additional property that for each point x ∈ M h ⊆ N δ there is a unique point π(x) on M that minimizes the distance from M to x.In other words, the computation of the closest point projection (2.1) is a local minimizer problem [BV21] in the sense: where the nonlinear projection maps a point x ∈ M h to the point on M that minimizes the distance to x.The accuracy of standard surface integration methods is limited to only first or second order due to the use of piecewise linear approximations of the surface geometry and the integrand.To obtain high order accuracy, we must construct a high order approximation, both, to the geometry of the surface and to the integrand.
By relying on [PS22] we give now the construction of a surface M k h , which is locally parametrised over the reference simplex σ by polynomials of degree k, interpolating the smooth surface M. As aforementioned, we first consider a piecewise flat triangulation M h of the smooth surface with vertices lying on M. To simplify the notation, the subscripts from the transformation map (1.1) and the triangulation of M h are dropped.For any triangle T ∈ T 1 h with vertices q 1 , q 2 , q 3 ∈ R 3 we define the following map to be the affine linear parametrization which maps each vertex pi of σ to the vertex q i of T. We denote the images of the non-vertex nodes of σ on the simplex T by qi = ξ ( pi ) , where N (2, k) is the dimension of the vector space P 2,k (σ) of bivariate polynomials of degree k.
be the local Lagrange basis functions of degree k on σ corresponding to the nodal points p1 , . . .pN(2,k) .Set and define Q ψ,k to be a k−th order polynomial interpolation of the mapping ψ:

Now,
defines a polynomial mapping Q ψ,k by interpolating the points p i ∈ M with the Lagrange-polynomials . Thus, for every simplex T ∈ T 1 h we compute the projection π(q i ) and define an isoparametric simplex T k by applying Lagrange interpolation of order k to the coordinates of the projected equidistant nodes (see Fig ( 2)).Furthermore, if the base-triangulation T 1 h is fine enough, then the map Q ψ,k is a diffeomorphism.By differentiating the interpolation polynomial of the map ψ, we obtain: (2.4) The Jacobian of the transformation is calculated using the equation (2.4).Given that Q ψ,k is a diffeomorphism ∀ T ∈ T 1 h , then by union of non-overlapping mapped elements: (2.5) Therefore, we have successfully obtained a k-th order discrete approximation M k h of the continuous surface M. Fig. (3) depicts an illustration of this procedure applied to a torus and a sphere.In the specific scenario where k = 1, we denote the discrete approximation as

Accuracy of Integration
In this section, we will prove our main theorem about the accuracy of integration when replacing M by its k−th order polygonal approximation M k h and show that symmetric triangles of the mesh combined with even-degree polynomials improve global errors.This motivates the following definition.4) are called symmetric with respect to their common vertex x 1 , if they satisfy the following property

Definition 3.1 (symmetric triangles). Given a collection of triangles
(3.1) x 5 x 4 x 1 Traditionally, the specific technique employed for triangulation refinement of a triangulation mesh T 1 h is often regarded as unimportant, provided that the grid size h −→ 0 as n −→ ∞.However, a more intricate picture emerges when the influence of integration process is taken into account.The choice of the appropriate type of triangulation can potentially result in a serendipitous cancellation of errors.Now, when examining a conforming triangulation T 1 h , we refine each triangle T ∈ T 1 h into four smaller triangles using bisection refinement, which involves connecting the midpoints of its edges with straight lines.The new elements are all congruent, and they are similar to T .Repeating the bisection refinement m− times, the total number of triangles at that level will be approximately O( 2 2m ).An advantage of this form of refinement is that each set of mesh points contains those mesh points at the preceding level.During this refinement procedure, we obtain different triangular elements in size.For example, in Fig. 4, you can observe three distinct triangular elements of varying sizes.
Consider C m , m ∈ N as the set of triangles with equal size within T ∈ T 1 h at the m-th refinement level, and let n m represent the number of triangles in this set, then we can write For example, from the illustration in Fig. 4, we have Observe that during each bisection refinement step, at least two triangles fail to meet the Definition 3.1.As an illustration, these triangles can be discerned based on their blue shading in Fig. 4. As a result, the density of symmetric pairs of triangles within T ∈ T 1 h after the m-th refinement level is O (ñ m ) while the remaining triangles exhibit a density of O √ ñm .Here, ñm denotes the total number of triangles within T ∈ T 1 h after the m-th refinement.Remark 3.2 (symmetric triangulation).Any (smoothed closed) surface triangulation mesh that has undergone the refinement process previously described will be referred to as a symmetric triangulation.
In light of these facts, we state the main theorem of this article, showing that even degree approximation is superior to odd degree approximation for symmetric triangulations.
Theorem 3.3.Let M be a smooth closed embedded hypersurface and f ∈ C k+2 (M, R).Consider a piecewise linear triangulation M h with mesh size h of the smooth surface having vertices lie on M and let M k h be the k−th order approximation of the smooth surface constructed using local fittings of π ∈ C k+3 M h , R 3 .Consider a symmetric triangulation, consisting of O (n) = O h −2 symmetric triangles, whereas the number of non-symmetric triangles is bounded by where As we approach the proof of Theorem 3.3, we need three lemmata.
Lemma 3.4.Let T be a triangle in T 1 h and assume that ψ ∈ C k+1 (T ) and k ≥ 0, then we have with h = diam(T ).The constant c depends on k, but it is independent of both ψ and T .
Proof.In order to prove this lemma, we will use Taylor's formula in several variables.Let's denote with α = (α 1 , α 2 ), where As a first step, let's consider the k− order interpolation of the map ψ Using Taylor's formula for ψ around the point (0, 0) we obtain: dµ, by interpolating each term in the Eq.(3.4) using a polynomial of order k and using the fact that the interpolation of the first term is exact because it is a polynomial of degree ≤ k, namely Subtracting Eq.(3.6) from Eq.(3.4), we have The right hand side of Eq.(3.7) can be written where E (µ; s, t) has the following form (3.12) Combining Eq. (3.7) with Eq. (3.12), yields Eq. (3.3).
Next, we recall the following Lemma.

Lemma 3.5 ([Chi93]
).Let k be an even integer.Let ψ (p) ∈ P k+1 be a polynomial of degree k + 1 on σ, and let Q ψ,k (p) be its interpolant of degree k.Then, for each p ∈ σ Lemma 3.6.Let M be a smooth closed hypersurface and consider a piecewise linear triangulation M h with mesh size h of the smooth surface M and let M k h be the k−th order approximation of the smooth surface constructed using local fittings of π ∈ C k+3 M h , R 3 .Consider a symmetric triangulation, consisting of O (n) = O h −2 symmetric triangles, whereas the number of non-symmetric triangles is bounded by Proof.First, let us write every term in Eq. (3.14) over a reference simplex In the same manner, as in Lemma 3.4, we obtain where Analogously, an expansion can be given for the errors in the partial derivatives of ∂ α 1 ψ (p) , ∂ α 2 ψ (p), amounting to apply the derivative on (3.16).Therefore, we have Using Taylor's theorem and expanding about (s = 0, t = 0), we obtain where E h k+2 ; (x k − x j ) , (x ℓ − x j ) and E h k+3 ; (x k − x j ) , (x ℓ − x j ) , describes the collection of terms with order k+2, k+3 in h, whose coefficients are combinations of the vertices of triangles with appropriate indices j, ℓ, and k.For example for a triangle T x 1 ,x 2 ,x 3 the coefficients are combination of (x 3 − x 1 ) and (x 2 − x 1 ) both for E (k + 2) and E (k + 3).From above computation we see that |g i − gi | is at least of order O h k+2 , confirming also the result obtained in ([RWJG12], Theorem 7).However, this result can be further improved assuming that the triangulation mesh comprises symmetric triangles.Thus, integrating Eq. (3.19) over a reference simplex, we have then, the left hand side of the Eq.(3.20) is at least of order O h k+2 for every T ∈ T 1 h .Writing Eq. (3.20) with respect to a pair of symmetric triangles T x 1 ,x 2 ,x 3 and T x 1 ,x 4 ,x 5 , we obtain From Eq. (3.21) and Eq.(3.22), we have the following This is due to the fact that the first right hand side integrand of the Eq.(3.21) and Eq.(3.22), are the collection of terms which are of order k +2 in h, where the coefficients are the combination of the vertices of triangles, thus each integrand is an odd function.In general for a triangle T x 1 ,x 2 ,x 3 we may write: In the same manner for Using the property (3.1),Eq. (3.24) and Eq.(3.25), we obtain The first part of inequality(3.14) is obtained by leveraging the fact that the number of symmetric triangles in the mesh M h is of order , in combination with Eq. (3.20) and Eq.(3.23).At this point, we notice that if k is even, then σ E (k + 2) dsdt = 0.This is due to the Lemma 3.5 because E (k + 2) represents errors when integrating a polynomial of order k + 1.Thus, we have (3.26) Again, writing Eq. (3.26) with respect to a pair of symmetric triangles, we obtain T x 1 ,x 2 ,x 3 and Now, if k is even and triangles are pairwise symmetric then the error contributed is This is attributed to the fact that for k even, the following holds true.
As in the odd case, combining Eq. (3.26) and Eq.(3.29), it yields the second part of inequality(3.14).
We can now prove Theorem 3.3 Poof of theorem 3.3.As in Lemma 3.6, expanding the function f using Taylor's formula around the point (0, 0), yields where R k+1 has the form Integrate both side of the Eq.(3.30), we have Using |g i |= O h 2 , the integral term on the right hand side is at least of order O h k+3 .If k is even due to the presence of symmetric triangles the first term on the right-hand side integral is 0, so we can compute the left-hand integral with an accuracy of order O h k+4 , while for k odd, we have the following T k i and the smooth closed surface, parameterized using (1.1) Let us rewrite Eq. (3.33) over a reference simplex where the quadrature rules are defined By making use of Lemma 3.6, we have For the last inequality, we have used Eq.(3.31) and Eq.(3.32).

Results and Discussion
In order to validate our findings, we developed tests employing the Gauss-Bonnet theorem [PPCT01], [Spi99] for a series of classic smooth closed surfaces given by the following equations: The analytic expressions for the Gaussian Curvature are 2) Torus K Gauss = cos v r(R+r cos v) , where we used toric coordinates: (x, y, z) = (R + r cos θ) cos φ, (R + r cos θ) sin φ, r sin θ , φ, θ ∈ [0, 2π).
3) Sphere In all experiments we utilize the algorithm of Persson and Strang [PS04] to generate Delaunay triangulations1 serving as the initial M 1 h surface approximation and a Gaussian quadrature rule of deg = 12 on each triangle, employing 32 quadrature points [Dun85].In Fig. 5 and Fig 6 we show the relative errors under mesh refinement.In order to calculate the relative error, we integrate the Gaussian curvature on the manifold and compare it with the predictions of the Gauss-Bonnet Theorem.The convergence rates shown on the plots coincide, confirming our theoretical results.
However, as can be seen from Fig. 7, the approach presented in [PS22] has some limitations, it runs into Runge's phenomenon.This is due to the fact that the Lebesgue constant Λ for the set of equidistant points grows exponentially [MS92].To mitigate this effect, we propose an alternative approach.It is   well known that under a proper map, the operations (e.g., interpolation, numerical differentiation, and quadrature) on a triangular element can be performed on the reference square.Below we will state and prove a theorem that we hope will pave the way for developing another powerful numerical integration method for closed surfaces.Prior to presenting the theorem, we first recall the modulus of continuity from [TIM63] (Chapter 3): Let δ 1 ≥ 0 and δ 2 ≥ 0 be given.The modulus of continuity for a continuous function f is represented by ω (f ; δ 1 , δ 2 ), which is defined as: The function ω (f ; δ 1 , δ 2 ) is semi-additive, i.e, To address Runge's phenomenon, we propose the use of Chebyshev-Lobatto nodes as the interpolation points.Specifically, when working within a one-dimensional interpolation domain of [−1, 1], we define  In light of the multivariate extension of Jackson's theorem [BBL02], we have the following result.

Figure 1 :
Figure 1: Representation of a smooth surface parametrization, where every region Vi forms a curved triangle.

Figure 2 :
Figure 2: Construction of the second order approximation of the smooth surface M 2 h (blue line).A simplex of a 'base' triangulation M h (green line) is shown.The interpolation nodes, here the center qi of an edge, are projected (grey line) onto the smooth surface M (red line) via the projection π.The projected nodal points π ( qi) and the vertices of M h are then interpolated, giving the second order approximation of the smooth surface M 2 h .

Figure 3 :
Figure 3: Lagrange parametrization for a torus and a sphere using equidistant nodes on vertices and edges.

Figure 4 :
Figure 4: Illustration of bisection refinement and a duo of symmetric triangles obtained from the triangulation T 1 h .

Figure 5 :
Figure 5: Relative errors by integrating the Gaussian curvature over the torus with radii R = 2, r = 1 with the ideal convergence lines h n .

Figure 6 :
Figure 6: Relative errors by integrating the Gaussian curvature over the unit sphere with the ideal convergence lines h n .

Figure 7 :
Figure 7: Relative errors by integrating the Gaussian curvature over the torus and the ellipsoid using N∆ = 2528, N∆ = 6152 respectively.