High-Order Method with Moving Frames to Compute the Covariant Derivatives of Vectors on General 2D Curved Surfaces

The covariant derivative is a generalization of differentiating vectors. The Euclidean derivative is a special case of the covariant derivative in Euclidean space. The covariant derivative gathers broad attention, particularly when computing vector derivatives on curved surfaces and volumes in various applications. Covariant derivatives have been computed using the metric tensor from the analytically known curved axes. However, deriving the global axis for the domain has been mathematically and computationally challenging for an arbitrary two-dimensional (2D) surface. Consequently, computing the covariant derivative has been difficult or even impossible. A novel high-order numerical scheme is proposed for computing the covariant derivative on any 2D curved surface. A set of orthonormal vectors, known as moving frames, expand vectors to compute accurately covariant derivatives on 2D curved surfaces. The proposed scheme does not require the construction of curved axes for the metric tensor or the Christoffel symbols. The connectivity given by the Christoffel symbols is equivalently provided by the attitude matrix of orthonormal moving frames. Consequently, the proposed scheme can be extended to the general 2D curved surface. As an application, the Helmholtz‐Hodge decomposition is considered for a realistic atrium and a bunny.


Introduction
Modern scientific computing involves complex geometries with various materials.Stable and accurate numerical schemes are often required to simulate physical or biological phenomena in complex geometry.The covariant derivative emerges in a complex geometry such as a curved surface, which is distinguishably different from the widely-used Euclidean derivative in the classical numerical differentiation methods.

Covariant Formulation
Suppose that x i is the set of axes for a domain of two-dimensional (2D) surface.Let x i be the tangent vector of the axis.In the Euclidean domain, the magnitude and direction of the tangent vector x i can be made uniform across the surface.However, the tangent vector does not have the same magnitude and direction on a curved surface.Thus, comparing vectors on a curved domain should reflect the change of the axis.For this reason, a covariant derivative is needed on the curved surface where the axis changes accordingly.Therefore, the metric tensor and Christoffel symbols of the domain emerge because they represent the magnitude and rotation of the curved axes' tangent vectors, respectively.
Consider a vector on a 2D curved surface Ω .Suppose orthogonal curved axes x i exist on the surface of domain for i = 1, 2 .The vector is expressed as the expansion of the tan- gent vectors as follows: Differentiating Eq. ( 1) with respect to the x j -axis, we obtain The first component in Eq. (2) is the directional derivative.The second component in Eq. (2) is the sum of two components: a component in the tangent plane and another component along the surface normal vector, denoted by .Then, the second component in Eq. ( 2) is expressed as follows: where Γ k ij is known as the Christoffel symbol of the second kind to measure how the basis vector changes as it marches along each coordinate axis.By projecting Eq. ( 3) onto the surface, we obtain the following expression of the covariant derivative [30]: where we exchange indices i and k for the second component in Eq. (4) for convenience.We introduce the bold symbol j to represent the covariant derivative along the axis x j .u i x j can be viewed as the directional derivative of u i along the axis x j .Consider a vector expressed in the two axes of the domain as follows: The covariant derivative along the vector can be formulated using the linearity of the covariant derivative as follows: (1) = ũ1 x 1 + ũ2 x 2 . (2) x j x i . (3) Equation ( 5) is also referred to as the Levi-Civita connection which preserves the Riemann metric and is torsion-free [31].The primary difficulty when computing Eq. ( 4) for the general 2D surface is to find a continuous and differentiable curved axis x i .This task is computationally and analytically challenging, particularly in regions with various curvatures or nonuniform anisotropic conductivity.The coordinate system should be constructed individually for every domain, even for the same PDE.For example, the spherical coordinate system for PDEs on the sphere provides the analytically-exact directional derivatives and Christoffel symbols in Eq. ( 4).However, a deformed sphere requires the new parametrization of the coordinate system, which is challenging.More practical applications to complex 2D surfaces such as human faces or biological tissues may not have a continuous coordinate system, which could be beyond the scope of the general mathematical expression.
The second difficulty involves accurately computing the corresponding metric tensor g ij and the Christoffel symbol Γ k ij .The error of the geometric factors is often required to be sufficiently small compared to the discretization error.Otherwise, the error of the geometric factors significantly deteriorates the solution to the equations.The metric tensor is defined as the inner product of the tangent vectors as follows: where g ii represent the square of the magnitude of the tangent vector of the x i -axis.For example, g ij = 0 for i ≠ j for the orthogonal axis.Correspondingly, the Christoffel symbol can be computed by differentiating the metric tensor g ij as follows [2]: where the tensor g ij is the inverse of g ij .For example, g ii = 1∕g ii and g ij = 0 for i ≠ j for the orthogonal axis.Computing the metric tensor is also challenging because it requires computing the length of the curved axis.An invalid construction of the curved axis or inaccurate computation of the length of the axis yields the Christoffel symbols with non-negligible errors.Inaccurate metric tensors or Christoffel symbol functions appear as corrupted coefficients of the PDEs, which cause nonphysical dynamics.

Applications of Covariant Formulation
The covariant derivative has been a crucial tool in image processing and surface PDEs.Recent works on the computation and application of covariant derivatives include the computation of covariant derivatives by discrete connection on a triangulated 2-manifold [27], the application of covariant derivatives to image regularization [4], and the reformation of covariant derivatives in Cartesian coordinates in the context of finite element methods [32].
The previous numerical schemes, which solve the PDEs with covariant derivatives, can be summarized as follows.The first type of numerical schemes computes the geometric factors of g ij and Γ k ij directly [10,28,33].This method is similar to the classical Euclidean numerical scheme except for the geometric factors.For diffusion equations, it is equivalent to solving the Laplace-Beltrami operator.For shallow water equations, using covariant derivatives has been one of the most popular methods for the spherical coordinate system on the sphere in recent years [25,36,39].However, the stability and accuracy of the numerical scheme significantly depend on the accuracy and smoothness of the metric tensor and Christoffel symbols, particularly if the geometric factors cannot be derived analytically [42].The geometric singularity is another critical problem, particularly for closed 2D surfaces such as a sphere.The second type of numerical schemes for solving equations on curved surfaces is to use the surface's ambient space.The ambient space is the space surrounding the domain.Using the ambient space either increases the dimensionality of the space or expands the given domain [20,23,38,43].For this scheme, a set of fixed axes is constructed by increasing or expanding the surface's domain.Consequently, geometric factors are not required.The numerical scheme is similar to those in the Euclidean space, except for the extra operation that projects or extrapolates variables into the domain.However, the problem of dimensionality emerges for this type of numerical schemes.
Projecting vectors or the gradient onto the tangent plane has attracted considerable attention because of its simplicity.However, these projection methods seem to have intrinsic limitations in terms of the stability and accuracy depending on the surface type [5,6,22,24].Exterior calculus has been introduced to tackle this problem [26,35], especially with finite-element approaches, known as finite-element exterior calculus (FEEC) [3,21].The exterior calculus approach has been known to be efficient in analyzing the stability of a scheme and identifying covariant derivatives.Moreover, spin-weighted spherical harmonics (SWSH) has been used for PDEs on the sphere [1].SWSH, as the generalization of harmonics on the sphere, was implemented in the context of the spectral method to solve scalar and vector advection equations and 2 + 1 Maxwell's equations on a deformed sphere [7].

Advantages and Limitations
We propose a novel scheme to compute covariant derivatives in the context of spectral/hp and finite element methods.The proposed scheme uses orthonormal basis vectors constructed at each grid point, called the moving frame.By moving frame, it means a set of orthonormal vectors { 1 , 2 , 3 } at each grid point of the domain as follows.For 1 ⩽ i, j ⩽ 3 , the tangent vectors of each moving frame satisfy the following properties: where ij is the Kronecker delta.Moving frames are plural when representing orthonormal basis vectors in the domain, denoted by i by omitting the index i and the corresponding location.On a 2D surface, 1 and 2 lie on the tangent plane, whereas 3 is aligned along the surface normal vector .When one vector for a moving frame should be mentioned, i is called the ith tangent vector with a specific index i for 1 ⩽ i ⩽ 3 .We suppose that the domain is locally Euclidean such that any region of the domain can be approximated as a Euclidean domain.A moving frame can be built at every point in the domain.In addition, we suppose that i is differentiable in each curved element but may not be differentiable across the interfaces.Subsequently, moving frames do not need to be differentiable in the entire domain.
Moving frames have been successfully used as numerical solutions to various PDEs on 2D curved surfaces, such as conservation laws, diffusion equations, shallow water equations, and Maxwell's equations [13,14,16,17].In the previous works [13,14,16,17], the covariant derivative is computed in the weak Galerkin method, and the metric tensor or Christoffel symbols are not required.Integration by part either eliminates the geometric factors or replaces them with the compensation of Euclidean divergence or curl of moving frames.Therefore, the direct computation of covariant derivatives with moving frames should be naturally sought after without explicitly computing the metric tensor.This scheme can be similarly implemented in the context of the finite difference or finite volume methods for complex domains.
The proposed scheme has the following advantages.First, the covariant derivative can be computed without constructing any analytic curved axis.Thus, this proposed scheme can be effective for any complex domain where a curved axis cannot be analytically constructed.Complex geometry includes closed, irregular, or non-convex surfaces.Second, the new scheme employs a simple pre-processing to construct the discontinuous moving frames across the interfaces of curved or regular elements.Consequently, geometric singularities disappear to achieve almost-optimal exponential high-order convergence for any differential operator.
However, the computational efficiency of the proposed scheme is not significantly different from that of the covariant formulation.The most crucial factor in determining computational efficiency is the space dimension of the domain because it is critical to the scalability of the numerical scheme.For example, the covariant formulation and the proposed scheme with moving frames have the same basis dimension as 2D curved surfaces.Any scheme using the ambient space of the domain increases the basis dimension up to three.
This paper provides the convergence results of the sphere only for differential operators.Convergence results are compared with those by covariant formulation and spherical moving frames.The comparison for a more complex geometry is not included for two reasons.First, deriving the accurate metric tensor and Christoffel symbol for a complex domain is challenging, analytically and computationally, even for a relatively simple geometry other than the sphere and plane.Second, the proposed scheme is applied to a curved or regular element obtained from the tessellation of the domain.Therefore, the proposed scheme is not affected by the global shape of the domain.Instead, the accuracy of the scheme is only affected by the alignment of moving frames and the curvature in each element.From this perspective, the sphere can be not much different from a complex 2D surface in terms of the application and accuracy of the proposed scheme.In order to demonstrate this feature of the proposed scheme, the Helmholtz-Hodge decomposition (HHD) is considered for a complex 2D surface like an atrium or bunny.

Order of the Paper
This paper is organized as follows: Sect. 2 explains the mathematical scheme to compute covariant derivatives with moving frames.Section 3 provides examples of covariant derivatives on the sphere.Section 4 explains the algorithm used to construct local moving frames.In Sect.5, we perform several test cases on a sphere to demonstrate the exponential convergence of the proposed scheme.Section 6 demonstrates the application of the proposed scheme using the HHD scheme on a bunny and an atrium.In Sect.7, we conclude this paper with a discussion on the theory developed and the result of the proposed numerical scheme.

Covariant Derivative by Orthonormal Moving Frames
Consider a sufficiently smooth 2D curved surface, denoted by Ω .Let Ω e be the tessellation of the surface Ω , such that where the empty set ∅ means that the common area between different elements is zero, excepts lines or points in common.Suppose that each element is locally Euclidean.Moving frames i are constructed at each point such that each tangent vector i for 1 ⩽ i ⩽ 2 is differentiable within each element.Let xi be the axis corresponding to the tangent vector i .Consider a vector on an element Ω e .Expand the vector in the tangent vectors i orthogonal to the surface normal vector as follows: Remind that the moving frames are orthonormal.Thus, the magnitude of the curved axis' tangent vector is incorporated into the coefficients u i .Differentiating Eq. ( 7) with respect to the axis xj , we obtain The first component in Eq. ( 8) is the directional derivative, same as in the Euclidean space.The second component corresponds to the effects by the change of the coordinate system.The second component is similar to the covariant derivative formula with the Christoffel symbols, as shown in Eq. ( 4).However, a different measurement regarding changing axis is required.The orthonormal tangent vector i is not a holonomic basis for the domain because the magnitude of the tangent vector is incorporated into the coefficient u i , not the basis vectors i .
The change of moving frames can be measured by the same moving frames as follows [11,12,31].Consider a small change of moving frames d i .This change is a vector meaning that the change of each tangent vector has a direction and magnitude.Then, this change of moving frames can be expanded by the moving frame itself as follows: where d i is the vector-valued 1-form, which takes a vector to yield a vector.A scalarvalued 1-form takes a vector to yield a scalar.For example, d i ( j ) yields a scalar to meas- ure how i changes as it moves along j .Note that ij in Eqs. ( 9)-( 11) is also 1-form and is called the connection 1-form.The notation ij as used popularly in modern era [31,34] is equivalent to Cartan's classical notation of i j [12].The orthonormality condition of moving frames yields the following equalities for 1 ⩽ i, j ⩽ 3: Then, we have the matrix expression of Eqs. ( 9)-( 11) as follows: 1 3 or we have where the bracket represents that the corresponding quantity is a matrix.In the Cartesian coordinate system (x, y, z), moving frames are expanded as follows: In the matrix form, Eq. ( 16) is expressed as follows: or we have where ( , , ) is the tangent vector of the Cartesian coordinate (x, y, z), respectively.The notation T corresponds to the transpose of the vector or matrix.The matrix [A] consists of the coefficients of moving frames, known as the attitude matrix [31,34].Note that the coefficient (e x i , e y i , e z i ) represents only the angle of the tangent vector with respect to x, y, and z axes, respectively.Because the tangent vectors are orthonormal, The matrix [A]'s inverse, which exists all the times, is the same as its transpose, that is Then, Eq. ( 18) can be written as follows: In the following proposition, we propose that the covariant derivative can be equivalently computed without the metric tensor.

Proposition 1
The covariant derivative with the Christoffel symbol Γ i jk shown in Eq. ( 4) can be equivalently computed with the connection 1-form ij without the metric tensor g ij or Christoffel symbol Γ i jk as follows: The covariant derivative of along a vector is given as follows: where v i = ⋅ i .Correspondingly, the covariant derivative along an arbitrary vector com- puted by Eq. ( 4) is the same as that computed by Eq. (21).
Proof By differentiating both sides in Eq. ( 18), we obtain Using Eq. ( 20), Eq. ( 23) turns into the following equality: By comparing Eq. ( 15) with Eq. ( 24), we obtain the following equality: or, in the matrix form, we have In the component wise, Eq. ( 26) can be expressed as follows.For 1 ⩽ i, j ⩽ 3 , we have Note that if i = j , ii = for any vector.Note that (de x , de y , de z ) and ij are all 1-form, thus for a specific tangent vector k , we have the following equality: The component (de ) can be computed by the following equalities: Let us introduce the Jacobian matrix, defined as follows: Then, the above equations can be expressed as the following vector dot product: Then, the component ij for k is derived as follows: x j in Eq. ( 8) is obtained as follows: By substituting Eq. ( 35) into Eq.( 8), we obtain the expression of the covariant derivative along j in Eq. ( 21).
The magnitude of the axis' tangent vector is not required in Eq. ( 21).Instead, it is implicitly incorporated into the coefficient u i .The tangent vector is only required when initially computing the vector coefficient along each axis.The vector's coefficients can also be directly obtained from the actual data or analytic expression.Therefore, the Christoffel symbol as the derivative of the metric tensor is not generally required.Nevertheless, the connection 1-form ij as the derivative of moving frames emerges in Eq. ( 21).
The computation of the connection 1-form is more convenient than that of the Christoffel symbol Γ i jk .The connection 1-form is computed by the directional derivative along the unit tangent vector j , as shown in Eq. (34).On the other hand, in order to derive the Christoffel symbols, a continuous curved axis should be first constructed, followed by the differentiation of the metric tensor along the curved axis in Eq. (6).No metric tensor or Christoffel symbols enables the covariant derivative scheme by Eq. ( 21) to be applied to complex 2D surfaces such as irregular and non-convex surfaces.

Summary
A novel numerical scheme of a high-order direct covariant derivative is formulated on curved surfaces.Construct moving frames { i } on a smooth curved surface Ω such that 1 and 2 lie on the tangent plane when 3 is aligned along the surface normal vector.The velocity vectors and on the surface Ω are expanded by moving frames as follows: The covariant derivative of ∇ is transformed into the following expression using the lin- earity of the covariant derivative: The above equation can be written as follows: Or, it can be written as follows: ( Substituting Eq. ( 21) into the above equation (38), we obtain Each covariant derivative along the unit tangent vector i is provided as follows: For example, the covariant derivative of the unit tangent vector i is obtained by using Or, by solving the summation notation, we obtain Note that Eqs. ( 45) and (48), representing the magnitude changes of i along the same tangent vector, hold due to the constant magnitude of orthonormal moving frames.Suppose that the curved axis of a 2D surface has two orthonormal tangent vectors; the tangent vectors are in unit length and orthogonal.Then, the covariant derivative in Eq. ( 4) is the same as Eq. ( 21): For the general 2D curved surface, Eq. ( 49) is not generally true. (44) (45)

Mathematical Comparison on the Sphere
We should check whether the two different covariant derivatives in Eqs. ( 4) and ( 21) yield the same result.Let us choose a sphere for convenience.Consider the sphere of the unit radius with the following representation: where is the polar angle and is the azimuthal angle.The infinitesimal length is computed as ds 2 = d 2 + sin 2 d 2 .Let the tangent vectors of each axis be denoted by and .Let a velocity vector be expanded in the spherical coordinate system as follows: Note that ‖ ‖ = 1 and ‖ ‖ = sin .The metric tensor is obtained by the inner product of the tangent vector to obtain g = 1 and g = sin 2 .
For the construction of orthonormal moving frames { i } , let the first and second tangent vectors and be aligned along the and , respectively.Let them be normalized with the unit length.We call this set of moving frame as the spherical moving frames.Similar to Eq. (51), a vector can be expanded in spherical moving frames as follows: where the basis vectors and coefficients have the following relations: The covariant derivative along the -axis contains nontrivial Christoffel symbols and can be obtained as follows: where Γ on the sphere is computed as follows: The covariant derivative along the unit tangent vector can be also computed by Eq. ( 21): The corresponding connection form has only one nontrivial component ( ) , which can be computed as follows.Consider the unit tangent vector of the spherical moving frames as follows: (50) (x, y, z) = (sin cos , sin sin , cos ), (57) = (∇u ⋅ ) + ∇u ⋅ .
By using Eq. ( 34), the connection 1-form can be computed as follows: Swapping the indices of ij yields the following sign changes: Comparing Eq. ( 55) with Eq. (57) yields the following corollary.
Corollary 1 For the sphere of the unit radius, the covariant derivatives of along the -axis computed by using the Christoffel symbol in Eq. (55) are the same as those computed by using the connection 1-form in Eq. (57).
Proof Let us substitute Eqs. ( 53) and (54) into Eq.( 55).Then, Γ ũ cancels out to yield the following equation: Note that Eq. ( 57) is the same as Eq.(60), that is By similar calculus, we obtain the following covariant derivative along the -axis: where Γ is computed by Eq. (56).Similarly, Γ is derived as follows: The covariant derivative along is computed by Eq. ( 21) as follows: where ( ) and ( ) are computed by Eq. (58).Then, we obtain the following Corollary.
Corollary 2 For the sphere of the unit radius, the covariant derivatives of along the -axis computed by using the Christoffel symbol in Eq. (62) are the same as those computed by using the connection 1-form in Eq. (64) multiplied by sin .
sin cos sin sin cos cos cos cos sin − sin − sin cos 0 Proof By substituting Eqs. ( 53) and (54) into Eq.(62) and using the linearity of the covariant derivative, Eq. ( 62) is transformed into the following equation: Comparing Eq. (65) with Eq. (64), we obtain Using Corollaries 1 and 2, we have the following proposition.
Proposition 2 On the sphere of the unit radius, the covariant derivative of along a vector by using the Christoffel symbol in Eq. ( 4) is the same as that by using the connection 1-form in Eq. (21).
Proof Consider the following expansions of the vector : The covariant derivative of along the vector by the connection 1-form in Eq. ( 4) is provided as follows: Moreover, the covariant derivative of along the vector by the connection 1-form in Eq. ( 21) as follows: By substituting that ṽ = v  and ṽ = v  ∕ sin  , the difference between Eqs. (68) and ( 69) is given as follows: where we used Corollaries 1 and 2.
Proposition 2 confirms that the covariant derivative by using the Christoffel symbol is the same as that by using the connection 1-form on the sphere.The covariant derivative using the Christoffel symbol [2] is most popular in scientific computing.The Christoffel symbol requires a set of coordinate axis and their tangent vectors.Consequently, constructing a set of continuous coordinate axes is often required.However, this process is challenging for a complex surface.The connection 1-form requires only the surface normal vector, which is relatively easy to derive [29].No continuous coordinate axis is required.The connection 1-form is relatively cheap and can be conveniently derived for a complex surface without constructing its coordinate system.(66)

Local Moving Frames: Algorithm
The first advantage of the covariant derivative in moving frames is that the direction of the first unit tangent vector is arbitrary.Moving frames do not have to be aligned along the spherical coordinate axis as follows: where 1 and 2 are arbitrarily aligned.
A question arises on which direction of 1 and 2 would yield the most optimal accuracy.Varying curvature affects the overall error because the differentiation along the unit tangent vector on a curved element is an approximation to the differentiation along the unit tangent vector on a plane [13].However, moving frames for a constant or slowly-varying curvature would also be acceptable because the error caused by varying curvature is significantly smaller than the discretization error.Modern meshing software provides a sufficiently accurate approximation of an object using refinement such as smaller elements for a more highly curved area.Thus, varying curvature can be effectively dealt with in modern scientific computing.
The second critical problem is the geometric singularity, particularly for a closed 2D surface like a sphere.This singularity emerges due to a continuous curved axis built on the closed surface.For example, the north and the south poles are the location where differentiation along the -axis is impossible.To overcome this problem, a set of discontinuous and almost-Euclidean moving frames is constructed, called the LOC moving frames, short for local moving frames.
The spherical coordinate system around = π∕2 on the sphere is primarily Euclid- ean.Applying this coordinate system to each element means that an individual spherical coordinate system is constructed for each element using the different location of its pole.Consequently, this kind of construction eliminates the continuity and differentiability of moving frames across the interface of the elements.Let us call these frames as LOC moving frames.Compare LOC moving frames with the moving frames that are aligned along the global spherical coordinate axis with the unique pole on the sphere.Figure 1 illustrates the alignment of moving frames across the interfaces of the elements.The LOC moving frames are a discontinuous reference coordinate system.Thus, no geometric singularity exists even for a closed 2D surface like a sphere.Moreover, the overall covariant divergence of LOC moving frames is significantly less than that of spherical moving frames.
This paper proposes an efficient algorithm to construct LOC moving frames with low covariant divergence.Consequently, LOC moving frames produce a relatively smaller error by directional differentiation.The proposed algorithm constructs the LOC moving frames that are generally aligned at various oblique angles with the -axis of the spherical coordinate axis.The unit tangent vectors of LOC moving frames are not perfectly Euclidean.Thus, the divergence of the unit tangent vectors is not zero.However, the divergence of LOC moving frames is usually less than that of the spherical coordinate system.
Let us consider the mapping from the standard element to a curved triangular element as shown in Fig. 2. Suppose that = (r, s) represents the coordinate axis of grid points for an element.(r, s) is also a map because each grid point of Ω i e (r, s) , denoted by (r, s) , is mapped from each grid point of the standard element Ω st (r, s) .The dif- ferentiation of (r, s) with respect to the s-axis produces the tangent vector s .This vector is parallel to the line in Ω i e (r, s) , which is mapped from r = constant in Ω st (r, s) . (70 A similar argument can be applied to r .In general, s and r are not orthogonal.The inner product is not zero, that is Let E k represent each curved edge of the triangular element Ω i e for 1 ⩽ k ⩽ 3 .In Fig. 2, E 1 and E 2 are parametrized by r and s, respectively.For the construction of the edge E 3 , sup- pose that a new axis emanating from the vertex where (r, s) = (0, 1) .The axis is repre- sented such that = s + r .Then, the edge E 3 can be represented by = 1 .In summary, the unit tangent vectors along each edge are constructed as follows: For each kth edge, the covariant divergence of the unit tangent vector E k , or ∇ ⋅ E k , is com- puted for 1 ⩽ k ⩽ 3 .If the unit tangent vector's axis is less curved in Ω i e , then the divergence of E k is likely to be closer to zero.Thus, it is reasonable to choose E k , such that E k is the smallest divergence among the three in Eqs. ( 72)-(74).The selected E k becomes the first tangent vector 1 .Because the third unit tangent vector 3 is the same as the surface normal vector , the second unit tangent vector 2 is obtained by the cross product as 2 = × 1 .Repeat- ing this procedure for all the elements, we obtain a LOC moving frame at every grid point.Observe that the edge index k indicates that the edge to construct 1 could be different for each element.In summary, consider the following algorithm for LOC moving frames.Figure 3 presents the difference of ∇ ⋅ i between spherical and LOC moving frames for a tessellated spherical mesh with 498 elements and a 4.99 × 10 −8 mesh error.For spherical moving frames { , } , the maximum divergence magnitude of and in the L 2 norm is 2.5 and 0.25, respectively.This divergence can be expected from Eqs. ( 40)-(43) as follows: For LOC moving frames, the maximum divergence magnitude of 1 and 2 in the L 2 norm is 0.35 and 0.4, respectively.Table 1 displays the convergence of divergence as the element size (h) varies.As h decreases, the covariant divergence of decreases.The covariant divergence of remains significantly large.For LOC moving frames, the covariant divergence of 1 or 2 does not converge to zero.However, the covariant divergence of LOC moving frames remains significantly smaller than that of .Moreover, the difference between the covariant and Euclidean divergence provides a vital indicator of moving frames for computational accuracy.The covariant divergence is the same as the Euclidean divergence if the Christoffel symbol or the connection 1-form is zero.In Table 1, this difference between the divergences for LOC moving frames is approximately equal to or less than 0.1% of that for spherical moving frames, for both 1 and 2 .

Table 1 Comparison between the divergence of spherical and LOC moving frames
The divergence of is much larger than that of LOC moving frames.The divergence of converges to zero."diff" measures the difference between the covariant and Euclidean divergence.For the spherical MF, some elements around each pole are eliminated to remove geometric singularities

Test Cases on the Unit Sphere
This section displays several test cases to display that the LOC moving frames with the connection 1-form ( ij ) compute the covariant derivative with sufficient accuracy.The proposed scheme can be applied to any 2D curved surface, including irregular or non-convex surfaces.For the convenience convergence error evaluation, we use the unit sphere with the following metric in this section.Thus, the invariant square of an infinitesimal length, known as an interval, is obtained as follows: Two types of meshes of the sphere are used.The first mesh is a projected mesh.Place a cube of unit length inside the sphere such that each vertex lies on the sphere.The edges of the cube are equally dissected, respectively.Let each node in the kth edge be denoted by i for 1 ⩽ i ⩽ N k .The vertices of the nodes i are projected onto a sphere such that ‖ 2 i ‖ = 1 , i.e., the nodes of the mesh are exactly on the sphere.Let us call this mesh ProjMesh [18].However, the grid points other than the nodes may slightly lie outside the sphere.Because of this approximation error, the mesh error of ProjMesh does not converge as the order of the solution polynomial (p) increases.
The second mesh is a high-order mesh, called NekMesh [37,40].By the optimal distribution of all the grid points by energy minimization, the location of grid points is redistributed to minimize the mesh error, even p increases.However, the vertices are not exactly on the exact domain.The mesh error of each mesh is displayed in Fig. 4. Table 2 shows the mesh error and ( ) error versus p. NekMesh is generated by an open-source 3D high- order mesh generator NekMesh within the Nektar++ framework [9].For further references of the two meshes and downloading the files, refer to [18].
In the following test problems, we demonstrate that the proposed scheme for covariant derivatives converges on both meshes.In contrast to weak formulations, the differentiation based on the grid points strongly depends on the point-wise location of grid points and the relative location of grid points along with a directional differentiation.Thus, the proposed scheme may yield different accuracy depending on the direction of differentiation.For example, the spherical coordinate axis may yield better accuracy than the proposed scheme with LOC moving frames for some meshes.The grid points in this mesh are more (77) Fig. 4 Distribution of mesh error for ProjMesh and NekMesh accurately aligned along the spherical coordinate axis than the direction of LOC moving frames.Therefore, the comparison between the spherical coordinate system and LOC moving frames conveys little meaning-only the convergence of every mesh matters.
To minimize the effect of the mesh error, LOC moving frames are projected into the exact domain.This projection implies that the third unit tangent vector 3 of LOC moving frames is aligned along with the exact surface normal vector, not the surface normal vector of the domain with some degree of the mesh error.The first step is to replace the third tangent vector 3 with the unit surface normal vector of the exact domain.The second step requires the projection of 1 and 2 to the new tangent plane orthogonal to the new 3 .This type of moving frames is called LOCSPH, short for local spherical moving frames [19].In a high-order mesh with a negligible mesh error, LOC moving frames are almost the same as LOCSPH moving frames.

Covariant Derivative on the Sphere
Construct a flow modified from the steady zonal flow [41].The velocity vector of the steady zonal flow is multiplied by sin to remove the case that a component is divided by zero at = 0 .For a vector on the sphere, consider the following expansion of the vector in moving frames: The vector components in the spherical coordinate system are derived as follows: where is the angle between the velocity flow and the north pole.Let be arbitrary.The vector component in the Cartesian coordinate system is expressed as follows: The exact covariant derivative along each spherical coordinate axis is derived as follows: (79) u = −u 0 sin sin sin , (80) u = u 0 sin (cos cos + sin cos sin ),  where the directional derivatives are given as follows: Two s are considered: = 0 and = π∕2 .The velocity vector for each is displayed in Fig. 5. Table 3 presents the exponential error convergence in the L 2 norm in ProjMesh and NekMesh.For both meshes, the proposed scheme with LOCSPH moving frames shows almost exponential convergence.For ProjMesh, the mesh error does not converge as p increases as expected.However, the connection 1-form ( ) converges as p increases.For = 0 , the proposed scheme with moving frames shows a slightly better accuracy than that with spherical moving frames.For = π∕2 , it shows a different result; the scheme with spherical moving frames shows slightly better accuracy.The numerical accuracy of covariant derivative depends on the location accuracy of the grid points.This location accuracy varies from mesh to mesh, regardless of the mesh error.However, they all show an exponential convergence.For a sufficiently large p, the difference of errors for both is negligible.
For NekMesh, the mesh error converges as p increases, in contrast to ProjMesh.( ) also converges as p increases.The convergence with NekMesh is similar to that with Pro-jMesh.For = 0 , the scheme with LOCSPH shows significantly better accuracy than that with spherical moving frames.For = π∕2 , the scheme with spherical moving frames shows better accuracy than that with LOCSPH.In spite of significantly lesser mesh error, the grid points along LOCSPH moving frames are less accurately aligned than those along with the spherical coordinate system.For the remaining test problems, a more complex velocity field is used.This velocity field is similar to the Rossby-Haurwitz wave popularly used as a test problem for the shallow water equations: where the vector components are derived as follows: where = K = 7.848 × 10 −6 s −1 .Equations ( 82) and ( 83) are normalized such that the radius of the earth corresponds to the unit length.In the following tests, only NekMesh is considered to remove the effects of the mesh error.

Gradient
The gradients of a scalar variable on a 2D curved surface should be independent of the chosen axis.In the spherical coordinate axis of ( , ) , the gradient of a scalar variable v is expressed as follows:  where each component is computed such as where and K are the same as Eq. ( 81) and ( 82).The computation of the gradient in spherical or LOC moving frames should produce the same gradient, expressed as follows: Computationally, Eq. ( 85) should converge to Eq. ( 84), particularly as the space resolution increases.Table 4 confirms the the exponential convergence of the difference between the covariant differentiation and LOC moving frames for the gradient of v as p increases.

Divergence
In the spherical coordinate axis for the unit sphere, the divergence of the velocity vector in Eq. ( 81) can be computed by the following covariant formulation: On the other hand, the divergence in moving frames can be computed as follows: Note that the direction of the unit tangent vector 1 is arbitrary and possibly could be discontinuous across the interfaces of elements.Using the covariant derivatives in Eqs. ( 40)-(43), we obtain where the connection 1-form ij ( k ) is computed by Eq. (34).Note that the divergence of the velocity vector in Eq. ( 81) is analytically zero.Figure 6a displays the velocity vector on the sphere and the distribution of the divergence.Figure 6b demonstrates the exponential convergence by Eq. (86), Eq. (88) with spherical moving frames, and LOC moving frames.Figure 6b confirms that Eq. ( 88) with LOC moving frames has the optimal convergence order.The difference between LOC moving frames and the other methods becomes larger as p increases.
The geometric error's contribution to the overall discretization error seems to increase as p increases.However, there is no noticeable difference in computational efficiency among the three methods.Their difference only lies in the magnitude and direction of the unit tangent vectors.Table 5 presents the convergence order for the three methods of computing covariant divergence.Equation (88) with LOC moving frames is the most accurate with the optimal convergence order, closer to the theoretical expectation.The order should be theoretically the same as the order of p because it is the first derivative of a vector.However, the geometric error of the mesh undermines this order, ending up with 3.47, 5.3, 4.03 for covariant formulation and 3.27, 5.25, 4.04 for spherical moving frames.The order is increased to be nearly equivalent to the optimal order for LOC moving frames: 5.08, 4.79, 5.58.Consider that the divergence of the unit tangent vector becomes larger as it is closer to the pole.Thus, it is possible that the lesser accuracy of spherical moving frames is caused by the geometric singularities of the sphere.

Curl
Consider the projection of the curl of the vector onto the surface normal vector , denoted by ⋅ ∇ × .Note that the curl for the velocity vector of the Rossby-Haurwitz wave in Eq. ( 81) can be analytically obtained as follows: The covariant curl operator in the spherical coordinate axis is provided as follows: Similar to covariant divergence, the above projected curl can be computed with orthonormal moving frames as follows: Using the covariant derivatives in Eqs. ( 40)-( 43), we obtain Figure 7a displays the velocity vector on the unit sphere and the distribution of the projected curl ⋅ (∇ × ) . Figure 7b shows the similar exponential convergence as that of the divergence test by Eq. (89), Eq. ( 91) with spherical moving frames, and Eq. ( 91) with LOC moving frames.Similarly, Eq. ( 91) with LOC moving frames shows the optimal convergence.Table 6 presents the convergence order for the three methods.Equation (91) with LOC moving frames shows an improved convergence order of 5.077 5, 4.786 5, 5.584 0, closer to the optimal convergence of p. Similar to the divergence case, no geometric singularity for LOC moving frames seems to provide better accuracy.
6 Application: Helmholtz-Hodge Decomposition This section illustrates an example of the proposed scheme using a complex surface.As an application of the proposed scheme for covariant derivatives, the numerical scheme of the HHD was proposed for a complex curved surface.On a 2D curved surface Ω with either a Neumann boundary or no boundary, the HHD finds the three unique components of a vector field [8] ⋅ (∇ × ) = −2 cos + 30K sin 4 cos cos 4 .
(89)  (94) and (95) fails to locate the exact source of the flow that is represented as the irrotational and incompressible components.Equations ( 94) and (95) are solved by the Helmholtz solver, which was implemented at the open-source spectral/hp library and is referred to as Nektar++ [9], in the context of continuous or discontinuous Galerkin methods.For computational tests, a curl-less vector is derived from spherical moving frames and is multiplied by sin .A divergenceless-vector was also derived from the aforementioned Rossby-Haurwitz vector in Eq. (81).Figures 8 and 9 represent the HHD of a curl-less vector and a divergence-less vector.The error for each component in both HHD is shown in Table 7.
Two complex domains are introduced to demonstrate the HHD for complexly-curved 2D surfaces: the first is the surface model of the human atrium, and the second is the Stanford bunny.The initial vector was obtained by propagating a diffusion-reaction type wave, such as the two-variable Aliev-Panfilov model.The wave was initiated from a point and then moving frames are aligned along with the gradient of the action potential.The geometry has no holes, corresponding to zero Betti number [3].However, the wave propagation on the geometry is depicted by eliminating the initial and final stages of point-wise propagation.Therefore, the geometry of the vector field is equivalent to a domain with a hole as if a circular cylinder penetrated the initial and final region of excitation.Thus, the Betti number is actually one, and consequently, the harmonic field vector has at most a dimension of one.For the bunny in Fig. 11, we have the following error by the 2D continuous Helmholtz solver with moving frames: The magnitude of the vector Laplacian is not negligible in some regions of the domains, particularly near the boundaries of the atrium and non-smooth junctions of the bunny.However, the proposed scheme still yields the smooth harmonic potential U in this region.

Discussion
The proposed scheme provides an efficient high-order algorithm for computing covariant derivatives on any 2D curved surface with optimal accuracy.The scheme requires no metric tensor or Christoffel symbol from the analytically-known curved axis.This paper verifies that the covariant derivative using the tangent vectors is the same as that using orthonormal basis vectors on the sphere of the unit radius.Therefore, this scheme can be directly applied to various numerical schemes in more complex 2D realistic geometries.However, the drawbacks are as follows.The proposed scheme can have poor accuracy for a curved element with a highly-varying curvature and Jacobian.However, a modern adaptive mesh generator usually takes care of this issue by constructing a smaller highly curved element in a region of higher curvature.Moreover, the computational efficiency of the proposed scheme is almost the same as the covariant formulation with the Christoffel symbols.
Future studies should consider the difference between the Riemann curvature tensors using the tangent vector and orthonormal basis vectors.The Riemann curvature tensor means the Riemann curvature tensor of a wave propagation such that the first tangent vector of moving frames is aligned along with the propagation.The Riemann curvature tensor is closed related to the relative acceleration of propagation.Consider a biological signal propagation, such as cardiac action potential or neural spike propagation.The relative acceleration in a biological signal propagation is known to provide the stopping condition of the propagation as a representation of sink and source ratio [15].Consequently, we investigate how the relative acceleration of propagation can be computed by orthonormal moving frames, not by using the metric tensor and the Christoffel symbols.permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 1 Fig. 2
Fig. 1 Illustration to display the continuity of moving frames for spherical and LOC moving frames

Algorithm 1
Construct LOC moving frames for a triangular element while For each edge k do Compute the tangent vector of each edge by Eqs.(72), (73), and (74).Normalize the vectors and denote it by e Ek Compute the covariant divergence ∇ • e Ek end while Compare all e Ek and assign one with the smallest divergence to e 1 . .Compute e 2 by e 3 × e 1 with the unit surface normal vector e 3 return {e 1 , e 2 , e 3 } . .

Fig. 5
Fig.5 Distribution of the velocity vector for various .The contour represents the magnitude of the velocity vector.Some elements around the poles are removed to avoid the geometry singularities of the spherical coordinate system

26 ×Fig. 6
Fig. 6 Distribution of the velocity vector and the divergence ∇ ⋅ by LOC moving frames for p = 5 .The p-convergence of divergence on the sphere of the unit radius is shown for h = 0.344 5 .Covariant: covariant formulation, SPH: spherical moving frames, LOC: local moving frames On a 2D surface, Eqs.(94) and (95) contain covariant divergence, i.e., ∇ ⋅ and ∇ ⋅ J , which can be computed by Eq. (88).Inaccurate computation on the right-hand side of Eqs.

Fig. 8
Fig. 8 HHD of a curl-less vector, spherical moving frames multiplied by sin .a The irrotational component of the vector with the potential u, b the potential of the incompressible flow v, and c the potential of the harmonic flow U such as = ∇U , h=0.2, p = 10

Fig. 9
Fig. 9 HHD of a divergence-less vector, the Rossby-Haurwitz velocity vector.a The potential u of the irrotational component, b the potential v of the incompressible flow, and c the potential U of the harmonic flow such as = ∇U , h=0.2, p=10

Fig. 10
Fig. 10 HHD by the aligned moving frames along the propagational direction on an atrium.a Irrotational vector with the potential u, b incompressible vector with the potential v, and c harmonic vector with the potential U

Table 2
Mesh error for both meshes to compute the location errors of grid points in the L 2

Table 3
Covariant error of the vector shown in Eqs.(79) and (80) along the and axes.Error is computed in the L 2 norm on a projection mesh for = 0 and = π∕2 .SPH spherical moving frames, LOCSPH LOC-SPH moving frames

Table 4
Error between Eqs. (84) and (85) for the unit sphere.h = 0.4 and 498 elements.Some elements surrounding the poles are removed due to the singularities of the spherical coordinate axis