Arbitrary-order intrinsic virtual element method for elliptic equations on surfaces

We develop a geometrically intrinsic formulation of the arbitrary-order Virtual Element Method (VEM) on polygonal cells for the numerical solution of elliptic surface partial differential equations (PDEs). The PDE is first written in covariant form using an appropriate local reference system. The knowledge of the local parametrization allows us to consider the two-dimensional VEM scheme, without any explicit approximation of the surface geometry. The theoretical properties of the classical VEM are extended to our framework by taking into consideration the highly anisotropic character of the final discretization. These properties are extensively tested on triangular and polygonal meshes using a manufactured solution. The limitations of the scheme are verified as functions of the regularity of the surface and its approximation.


Introduction
Surface partial differential equations of elliptic and parabolic types are often used for the simulation of diverse phenomena in many fields of applications, for example in biology, atmospheric dynamics, and image processing [53,55]. One of the main motivation that prompted this work is related to the modeling of gravity-driven flows in earth-sciences, such as, e.g., flood forecasting, landslide and debris flow dynamics, avalanche simulations [6,20,39,41]. The numerical solution of surface PDEs has seen wide-spread interest in the last few years with different approaches being proposed, including continuous Finite Element Methods (FEM), discontinuous Galerkin (DG), finite volumes, trace-FEMs, etc. [4,6,7,33,36,41,54,55]. A recent survey of surface FEM was published in [37], where both steady and moving surfaces are considered, the latter finding a unifying theory in [38].
Although FE-based approaches are very successful in the numerical treatment of surface PDEs, they share the limitation that an explicit form of the basis functions is required in the formulation of the method, and thus are restricted mostly to triangular/quadrilateral elements. This restriction is overcome by the Virtual Element Method (VEM) that was designed from the very beginning to work on generally shaped elements with high order of accuracy. In fact, in the VEM approach (i) it is possible to decompose the computational domain into very general polygonal elements; (ii) an explicit form of the basis functions is not required; (iii) approximation of arbitrary order and arbitrary regularity are straightforward in two and three dimensions. VEM was originally developed as a variational reformulation of the nodal mimetic finite difference (MFD) method [10,13,24,47] for solving diffusion problems on unstructured polygonal meshes. A survey on the MFD method can be found in the review paper [45] and the research monograph [12]. The scheme inherits the flexibility of the MFD method with respect to the admissible meshes and this feature is well reflected in the many significant applications that have been developed so far, see, for example, [3, 5, 8, 9, 14-19, 27, 29, 30, 32, 51, 56, 57, 64].Because of its origins, VEM is intimately connected with other FE-based approaches. The connection between the VEM and finite elements on polygonal/polyhedral meshes is thoroughly investigated in [26,35,46], between VEM and discontinuous skeletal gradient discretizations in [35], and between the VEM and the BEM-based FEM method in [28]. VEM was originally formulated in [11] as a conforming FEM for the Poisson problem, and was later extended to convection-reaction-diffusion problems with variable coefficients in [15]. However, the VEM technology has seen so far very few applications to surface PDEs, and only with first-order polynomial accuracy [42].
One of the major difficulties in the high-order numerical solution of surface PDEs is the achievement of a consistent approximation of both the geometry and the PDE. The work of [34] develops a general technique for high-order polygonal approximation of a smooth manifold, but this method requires the explicit knowledge of the distance function within the tubular neighborhood of the surface. The task of approximating this distance function to high order is still an open problem [50]. A recent approach based on this idea was presented in [4], where the authors study a high-order (up to four) DG scheme based on a piecewise polynomial approximation of the surface triangulation. However, extensions to polygonal grids with high polynomial orders have not yet been addressed. While VEM provides an ideal framework to work at high-order on generally shaped cells, according to [42] the main difficulty is the high-order approximation of the surface, limiting their current developments to polynomials of order one. The same authors suggest the use of the approach in [34] to extend their VEM scheme to higher order polynomials, without however eliminating the difficulty of the approximation to a consistent order of the distance function.
In this paper, we develop a novel VEM-based approach for the solution of elliptic surface PDEs that works at all polynomial orders. We avoid the difficulties related to high-order surface approximation by employing intrinsic geometry and following the approach described in [7] to adapt the virtual element technology to the surface PDE. Using this approach, we first rewrite the partial differential equation in covariant form in such a way that the geometric information, essentially the metric tensor, is completely encoded in the equation itself. As a consequence, the numerical scheme can be constructed directly on the two-dimensional local chart where the surface parametrization is defined, thus enabling the full exploitation of the VEM machinery. Here, we restrict our attention to the case when the surface is defined by a single chart, a case of great interest for example in gravity-driven flows on terrain surfaces, such as water flow and sediment transport in mountain areas [20,21,39,40,52]. In principle, our proposed approach can be applied to the more general situation of a surface defined by an atlas if the transition between charts is done with care by enforcing proper smoothness as described in [44]. This is shown by testing our approach on the sphere by using the well-known charts arising from the stereographical projection.
Our method starts from a partition of the surface into polygons with curvilinear edges, assuming that the parametrization of the surface is known at relevant quadrature points. Proceeding from the covariant PDE, we construct a high-order scheme exploiting the ability of the VEM approach to discretize problems that are anisotropic and with spatially variable coefficients [15]. In practice, we redefine the PDE on a local coordinate system using intrinsic geometric quantities and operators, which contain explicitly the metric information deriving from the surface. Then, all the VEM projection operators are calculated using this local coordinate system and the knowledge of the parametrization is used to define the needed quantities, thus incurring in no explicit geometric error. Hence the final scheme is defined on a planar two-dimensional domain (the surface chart) and all the available machinery to achieve high-order on polygonal cells can be exploited. The price we pay is that now the PDE contains the anisotropic metric tensor and all the coefficients vary in space as a function of the regularity of the surface. The virtual element method has proved its efficiency in handling these situations [49] and can be implemented directly in this two-dimensional setting. In addition, with our approach the convergence theory extends straight-forwardly to surface problems without additional efforts.
Our development of the intrinsic VEM proceeds as follows. In Sect. 2, we describe the local reference system of choice and the geometric setting. Next, we define the differential operators and the corresponding PDE in covariant form. In Sect. 3, we summarize the VEM adopted in this work and discuss the necessary adaptations to the problem at hand. The final Sect. 4 reports the results of extensive numerical experiments that assess the effectiveness, accuracy, and robustness of the proposed approach.

The surface partial differential equation and its Galerkin discretization
Notation Throughout the paper, we use the standard definition and notation of Sobolev spaces, norms and seminorms (see [1]), which can be directly extended to a compact manifold (see [63]). Given an open and bounded subset of ℝ d , d = 2, 3 , we denote with L p ( ) and W k,p ( ) the Lebesgue and Sobolev spaces, with W k,2 ( ) = H k ( ) the classical Hilbert space. Norms and seminorms in H k ( ) are denoted by || ⋅ || H k ( ) and | ⋅ | H k ( ) , respectively, and (⋅, ⋅) denotes the inner product in L 2 ( ) . We omit the subscript in the inner product notation when is the whole computational domain. In a few situations, for the sake of clarity, we may prefer to use the integral notation of the inner product. Consider a compact surface ⊂ ℝ 3 with boundary over which the following elliptic partial differential equation is defined: where the solution u ∶ → ℝ is a scalar function defined on the surface, ∶ → ℝ 2 is a given divergence-free velocity field tangent to the surface, the function ∶ → ℝ is a non-negative reaction coefficient. We denote by Δ G and ∇ G the Laplace-Beltrami and the tangential gradient operators, respectively, and by ⟨⋅, ⋅⟩ G the intrinsic scalar product. These operators will be given precise definitions depending on the chosen coordinate system. Classically, we assume f ∈ H −1 ( ) , ∈ [W 1,∞ ( )] 2 , and − 1 2 ∇ G ⋅ > 0 . Here we consider homogeneous Dirichlet problems with the more general boundary conditions described in [25].

Remark 1
The well-posedness of Problem 1 follows from the application of the Lax-Milgram theorem since the classical theory of elliptic equations can be extended to surface PDEs in a straight-forward manner. In particular, due to the coercivity of the bilinear form a(⋅, ⋅) , the continuity of the bilinear forms a(⋅, ⋅) , b(⋅, ⋅) , and c(⋅, ⋅) and the linear functional F(⋅) , and under the assumption that − 1 2 ∇ G ⋅ > 0 , the solution u exists and is unique and belongs to H 1 0 ( ) if f ∈ H −1 ( ) . For a detailed description of the properties, well-posedness, and regularity of the variational problem on manifolds we refer to [43] and to [62,63].
The discrete approximation of this problem reads as follows: These mathematical objects are defined and discussed in the next sections.

Geometrical setting
We assume that the surface is C m regular, i.e.: is a homeomorphism with its image (i.e., there exists an open neighborhood of , V ⊆ ℝ 3 such that (U) = V ∩ ); iii The differential d ∶ ℝ 2 → ℝ 3 is injective in U (i.e., it has maximum rank, in our case 2).
The map is the local parametrization of centered in and we denote with −1 ∶ V ∩ → U its inverse map, called the local chart, in . The set (U) ⊂ is called a coordinate neighborhood, while (s 1 , s 2 ) are the local coordinates of any point ∈ (U).

Remark 2
Throughout the paper we assume that is contained in only one chart. This is not a limitation under the assumption of C ∞ regularity of (or C m with m sufficiently large) since we can always find compatible local parametrizations covering . Indeed, given two points and ∈ with local parametrizations and such that U ∩ U ≠ � , the transition map • −1 is a C ∞ (or C m ) diffeomorphism. Thus, it is always possible to find an atlas for formed by appropriate charts that maintains all the required continuity properties. A proper selection of these charts is fundamental to obtain a numerically well-conditioned reference system in our approach. For an example of a constructive methodology for the definition of smooth multi-charts see [44].
For simplicity, from now on we will drop the subscripts and in both the local coordinates = (s 1 , s 2 ) and the global Cartesian coordinates = (x 1 , x 2 , x 3 ) . In summary, we have the following explicit definitions of these transformations: We want to choose a coordinate system to give a workable meaning to the partial differential equation and related differential operators. We define the local reference system following the approach in [6,7]. To this aim, we compute the pair of tangent vectors {̂ 1 ( ),̂ 2 ( )} on the tangent plane T : This pair is orthogonalized via Gram-Schmidt, yielding the orthogonal frame { 1 , 2 } . The ensuing metric tensor is given by: The associated scalar product between two vectors and is given by where " ⋅ " is the canonical ℝ 2 scalar product. Tensor G represents the realization of the first fundamental form with respect to the chosen reference system (chart). For a C m -regular surface (see definition 1), the determinant detG of the metric tensor is a well-defined and bounded function, and the metric tensor itself is coercive, i.e., it is symmetric and positive-definite and has a symmetric and positive-definite inverse. In other words, we can find constants g * and g * such that [31]: We can now write the intrinsic differential operators with respect to the local coordinate system, and we collect the appropriate definitions in the following proposition, which we state without proof: Proposition 1 (Intrinsic Differential Operators) Given f ∶ → ℝ a scalar differentiable function on and denoting with ∇ and ∇⋅ the gradient and divergence operators in ℝ 2 , the intrinsic differential operators expressed in the local coordinate system are given by the following expressions: We would like to recall that our reference frame is covariant and thus scalar products must act on vectors written in contravariant components. This applies both to velocity vector and to the divergence operator as well.

Remark 3
We can use the orthonormal reference frame 1 , 2 as a base for ℝ 2 with which we can express differential operators and vector quantities. In this case we need to take into consideration the Jacobian matrix = [ 1 , 2 ] of the parametrization, and recall that G = T . This applies in particular to the velocity vector field , which can be written as: where ̂ = [w (1) , w (2) ] T is the velocity vector written with respect to 1 , 2 .
In this setting, we can give the definition of the integral of a function over a surface as follows: Definition 2 Let f ∶ → ℝ be a continuous function defined on a regular surface , which we assume contained in the image of a local parametrization ∶ U → . The f on is given by We can relate any function f ∶ → ℝ to a specific coordinate system using the above coordinate transformations, i.e.: In the following we will make use only of the local coordinate system and will write f ( ) omitting the hat symbol. The classical tools deriving from Stokes theorems hold with the intrinsic operators without any modification. In particular, the intrinsic Green formula can be stated as in the following lemma.
Lemma 1 (Intrinsic Green formula) Let ⊂ ℝ 3 be a surface with smooth boundary and given two functions u ∈ C 2 ( ) and v ∈ H 1 ( ) , then: where ∶ → ℝ 2 denotes the vector tangent to and normal to with components written with respect to the local reference frame (i.e. = 1 1 + 2 2 ).
In view of remark 3, we reformulate the bilinear forms a(⋅, ⋅) , b(⋅, ⋅) , and c(⋅, ⋅) the linear functional F(⋅) of the intrinsic variational formulation (2) on the chart −1 ( ) through: and Therefore, the intrinsic variational formulation (2) is equivalent to solving the advection-diffusion-reaction equation in variational form: where the equation coefficients are defined by The problem as above formulated is still well-posed and maintains all the properties listed in remark 1. Indeed, since our surface is assumed to be regular, there exist two positive constants c where * = c G g * and * = C G g * , and g * and g * are the constants introduced in (5). The coercivity of a(⋅, ⋅) with respect to the H 1 -norm follows immediately by noting that �u� H 1 ( ) ≤ ‖u‖ H 1 ( ) and applying the Poincaré inequality in H 1 0 ( ) . Moreover, there exist two positive constants w max and max such that

The virtual element method
In this section, we discuss the virtual element approximation of problem 2. The numerical method that we use in this work is based on refs. [2,11,15], which define optimal approximations of the finite dimensional spaces on polygonal meshes when the equation coefficients are variable in space. As already observed in remark 2, we work on a single coordinate neighborhood Ω = −1 ( ) and consider the (global) parametrization ∶ Ω → . We start from a partition of formed by surface polygonal elements E with edges denoted by e . Through the parametrization , we can associate the partition of with a partition of Ω formed by elements E and possibly curvilinear edges e . Because of the regularity assumption on the surface, every element E is in a one-to-one relation with one and only one polygonal element E in Ω . To avoid curvilinear edges, we use the surface vertices of E to define the vertices of the polygons E in Ω through the inverse parametrization and connect them with straight segments to define the partition of Ω . This procedure maintains the above one-to-one relationship between elements E in and E in Ω (see Fig. 1).
In addition, any function in E can be expressed in E by composition with the inverse parametrization. Thus, all the local functional spaces of interest can be defined indifferently on E or E . The definitions of the building blocks of the virtual element method is done in Ω using standard two-dimensional Cartesian coordinates. These constructions are needed to evaluate the surface bilinear forms and the righthand side linear functional of the weak formulation (2) by a careful use of the metric tensor.
The conforming virtual element space Let T = {Ω h } h be a set of decompositions Ω h of the computational domain Ω into a finite set of nonoverlapping polygonal elements E . The subindex label h is the maximum of the diameters of the mesh elements, i.e., h E = sup � , �� ∈E | � − �� | . Each element E has a nonintersecting boundary denoted by E formed by straight edges e , center of gravity E and area |E| . A few regularity assumptions are needed on the mesh family {Ω h } to prove the convergence of the VEM and derive the error estimates in the L 2 and H 1 norms. We present these Let k ≥ 1 be an integer number and E ∈ Ω h a generic mesh element. The conforming virtual element space V h k of order k ≥ 1 built on mesh Ω h is obtained by gluing together the local approximation spaces denoted by V h k (E): where ℙ k (E) and ℙ k (e) denote the polynomial spaces of degree at most k defined over an element E or an edge e , respectively. By definition, each space V h k (E) contains ℙ k (E) and the global space V h k is a conforming subspace of H 1 0 (Ω) . The definition of the virtual element bilinear forms a h (⋅, ⋅) , b h (⋅, ⋅) , and c h (⋅, ⋅) , and the forcing term F h (⋅) requires the definition of the elliptic and orthogonal projections operators.

Elliptic projection The elliptic projection operator
Equation (14) allows the removal of the kernel of the gradient operator. The elliptic projection operator Π ∇,E k is a polynomial-preserving operator, i.e., Π ∇,E k q = q for every q ∈ ℙ k (E) . One of its major property is that the projection [11], which are defined as follows.
The degrees of freedom of the virtual element function v h ∈ V h k (E) are given by the set of values: where M k−2 (E) is the set of scaled monomials that span the linear space of polynomials of degree up to k − 2 . These set of values are unisolvent in V h k (E) , cf. [11], and thus, every virtual element function is uniquely identified by them. The degrees of freedom of a virtual element function in the global space V h k are given by collecting the elemental degrees of freedom (D1)-(D3). Their unisolvence in V h k is an immediate consequence of their unisolvence in every elemental space V h k (E) .
Orthogonal projections From the degrees of freedom of a virtual element function v h ∈ V h k (E) we can also compute the orthogonal projections Π 0,E k v h and Π 0,E k−1 ∇v h , cf. [2]. In fact, the definition of the orthogonal projection Π 0,E k v h reads as The right-hand side is the integral of v h against the polynomial q , and is computable from the degrees of freedom (D3) of v h when q is a polynomial of degree up to k − 2 , and from the moments of Π ∇,E k v h when q is a polynomial of degree k − 1 and k, cf. (12). Clearly, the orthogonal projection Π 0,E k−1 v h is also computable. In turn, using the definition of the orthogonal projection Π 0,E k−1 ∇v h and integrating by parts, we find that The bilinear form S E h (⋅, ⋅) in the definition of a E h (⋅, ⋅) provides the stability term and can be any symmetric positive definite bilinear form defined on E for which there exist two positive constants c * and c * such that Note that S E h (⋅, ⋅) must scale like the restriction of a(⋅, ⋅) on the mesh element E . Also, the stabilization term in the definition of a E h (⋅, ⋅) gives a zero contribution if one of its two entries is a polynomial of degree (at most) k since Π ∇,E k is a projection on the polynomial space. In this work, we consider two possible implementations of the stability term: -The choice originally provided in [11], which is sometimes called the "dofi-dofi stabilization" in the virtual element literature, and reads as where DOF i (⋅) is the map between a virtual function and its degrees of freedom; -The formula proposed in [48], which is sometimes called the "D-recipe stabilization" in the virtual element literature, and reads as where A is the matrix resulting from the implementation of the first term in the bilinear form a E h (⋅, ⋅) : where i (and j ) are the "canonical" basis functions generating V h k (E) , i.e., the functions whose i − th (or j − th ) degree of freedom is equal to 1 and all other degrees of freedom are 0. We note that these basis function are unknown in the virtual element framework, but their projections Π 0,E k−1 ∇ i (and Π 0,E k−1 ∇ j ) are computable. The stabilization term, and, in particular, condition (22), is designed in order that a E h (⋅, ⋅) satisfies the two fundamental properties: k-consistency: for all v h ∈ V h k and for all q ∈ ℙ k (E) it holds -Stability: there exist two positive constants * , * , independent of h and E , such that The virtual element forcing term To approximate the right-hand side of (3), we split the term into the sum of elemental contributions and approximate every local linear functional by means of the orthogonal projection Π 0,E k v h : With these definitions the VEM scheme in problem 2 is completely determined.

Convergence properties
The numerical analysis of the scheme requires the following hypotheses on the mesh, typical of VEM methods.
Assumption 1 (Mesh regularity assumptions) There exists a positive constant independent of h (and, hence, of Ω h ) such that (i) Every element E of every mesh Ω h is star-shaped with respect to a disk with radius ≥ h E ; (ii) Every edge e ∈ E has length h e ≥ h E .
The star-shapedness property (i) implies that the polygonal elements are simply connected subsets of ℝ 2 . In turn, the scaling assumption (ii) implies that the number of edges in each elemental boundary is uniformly bounded over the whole mesh family {Ω h }.
The following theorem summarizes the results for the virtual element approximation in problem 2. The proof of these results found in [15] is easily extended to our setting. Indeed, we choose to write the theorem in terms of the chart Ω and its discretization Ω h , but it can be written equivalently in terms of the surface and its discretization h , since the norms of the parametrization and its inverse are uniformly bounded by hypothesis.
-The H 1 -error estimate holds: -The L 2 -error estimate holds: The constant C may depend on the coefficient bounds * , * , w max and max , the stability constants * and * , the mesh regularity constant , the size of the computational domain |Ω| , and the approximation degree k.
The approximate solution u h is not explicitly known inside the elements. Consequently, in the numerical experiments of Sect. 4, we approximate the error norms as follows: Here, Π 0 k u h is the global projector on the space of discontinuous polynomials of degree at most k built on mesh Ω h , and ||u − Π 0 k u h || H 1 (Ω h ) is the norm in the broken Sobolev space H 1 (Ω h ) that is defined by summing the H 1 (E)-norms of each element E . Operator Π 0 k u h is obtained by taking the elemental L 2 -orthogonal projections Π 0,E k u h in every mesh element E , which are computable from the degrees of freedom of u h , so that Π 0 k u h |E = Π 0,E k u h|E .

Numerical results
In this section we present numerical results on synthetic test cases to support the statements of the previous sections by means of experimental evidence. Our test cases are grouped into four main categories. The first two sets of experiments, Test Cases 1 and 2, are aimed at showing the correctness of our implementation and the order of convergence of the proposed VEM scheme up to fourth order of accuracy. In the third set of experiments, Test Case 3, we explore the limits of the VEM approach as the metric tensor G becomes more and more anisotropic (large condition numbers) as a function of the regularity of the surface. Finally, Test Case 4 considers the use of stereographic projection to build a two chart atlas for the sphere to show the applicability of our approach in a multi-chart case.
In the first three experiments we consider the surface provided by the graph of the following height function, a simple trigonometric perturbation of a portion of a sphere embedded in ℝ 3 : where r is the radius of the sphere, and a and k are the amplitude and the frequency of the cosine trigonometric perturbation. We use the Monge parametrization given x 2 )} and work on the single chart represented by the domain Ω = {(s 1 , s 2 ) ∶ s 1 , s 2 ≥ 0 and √ (s 1 ) 2 + (s 2 ) 2 ≤ 1} . For r → 1 the metric tensor G tends to become singular as one of the two tangent vectors increases indefinitely at the boundary of the surface, leading to large spectral condition numbers (G) . Analogously, the condition number of G increases when the frequency k and the amplitude a are increased. Fig. 2 shows the three-dimensional plot of the surface in the left panel, and the spatial distributions of g 11 , g 22 , and √ detG , respectively, in the next three columns. The rows are relative to the case r = 2 , with the sphere ( a = 0 ) shown in the top, while the trigonometric deformation of the sphere is shown in the middle row for the case a = 0.5 and k = 5 , and in the bottom row for the case a = 2 and k = 5 . In this latter cases we note that G has a sinusoidal behavior with g 11 ≈ 1 in the regions where g 22 ≈ 20(150) , leading to (G) ≈ 20(150) . Note that G does not enter the reaction term and has a small effect in the advection term, as it amounts to a rotation and a stretching of the advective field. On the other hand, it has a large effect on the coercivity of the diffusion bilinear form, and thus we concentrate on the latter. The condition number of the VEM stiffness matrix A can be bounded by [22,58]: where (A LAP ) is the condition number of the VEM stiffness matrix of Laplace equation. In our case, we have that: Note that, √ det(G) (G −1 ) ≥ 1 is smooth although possibly unbounded when r → 1 (or k and a are large), as mentioned above. We would like to remark that this is not a contraddiction of Eq. (10), but rather a consequence of the fact that we use a single Monge parametrization. The presence of the metric tensor in the equation always deteriorates, possibly drastically, the condition number of the system matrix.
In Test Case 4, we consider Laplace equation (i.e., = 0 and = 0 ) on = S 2 and use two parametrizations, one for the northern and one for the southern hemispheres, given by: and: We proceed by discretizing the unit disk as reference domain once and for all for both hemispheres, with a polygonal approximation of the boundary. We then use the appropriate charts to express the VEM linear and bilinear forms in the northern and southern cells. We connect the two domains together by means of a simple Jacobi domain decomposition approach, and to avoid iterations we use the manufactured solution as boundary condition at the domain interface. Note that the use of curved edges would allow to solve the problem without the need to decompose the computational domain, exploiting the fact that the transition map for the two charts is readily available and sufficiently regular.
In all the experiments, numerical errors are evaluated by defining a manufactured solution u(s 1 , s 2 ) ∶ Ω → ℝ and calculating the resulting forcing function f (s 1 , s 2 ) by substitution into the original equation. Using u(s 1 , s 2 ) = sin(2 s 1 ) sin(2 s 2 ) and taking into account the contributions of the metric G , the general form of f (s 1 , s 2 ) can be evaluated as: where ̂ = [w (1) , w (2) ] T .

Test case 1
We use two families of polygonal meshes discretizing the domain (see Fig. 3). To avoid geometric error in the refinement process we uniformly distribute 8 nodes on the curvilinear boundary of Ω and approximate it with linear interpolation. All the refined meshes are built on this geometry. The meshes in the first family are constrained Delaunay triangulations obtained using Triangle [59,60], dividing by a factor 4 the area target of the elements at each level. The second family of meshes is obtained by means of PolyMesher [61] by imposing approximately the same number of elements of the triangulation at each corresponding level. Note that the sides of the boundary elements may contain as extra node one of the fixed vertices used to define the curvilinear boundary, and are thus formed by more than one edge. Convergence is tested on four different grid refinement levels, in the cases of a = 0 and r = 1.1, 1.01, 1.001 . Correspondingly, the results are reported in Figs. 4 and 5 for the triangulations and the polygonal meshes, respectively. The experimental convergence rates are optimal for all the tested polynomial orders, as can be seen from the figures and from the slopes of the lines, which are obtained by

Test case 2
This test case is designed to verify the robustness of the scheme for increasingly accurate approximations of the curvilinear boundary. To this aim we look at the errors for a fixed mesh size and vary the number of vertices used to discretize the curvilinear boundary, using only the polygonal mesh. Two sets of meshes are defined: one with approximately 25 elements and the other with 100 elements.
In each set we consider 5 mesh levels characterized by different approximations of the curvilinear boundary. In the first level the boundary is discretized with 8 vertices, as in the previous test case. The subsequent levels are obtained by doubling each time the number of nodes located on the curvilinear boundary to arrive at the final level with 128 vertices (see, e.g., Fig. 6). Since the size of the cells remains approximately the same, the boundary sides of the boundary elements are formed by an increasing number of straight edges. While the mesh levels approximate the curvilinear boundary with increasing accuracy, the fact that the length of the edges of the boundary elements becomes unbalanced may lead to increased errors in the VEM solution. However, we know from the literature [23] that such unbalance may only affect the constant that appears in the error estimates. Such constant is increased by a factor proportional to the square root of log(1 + ) , being the ratio between the maximum and the minimum edge length. Hence, in our experiments we expect the errors to remain approximately constant as we refine the boundary. The results of the simulations are shown in Fig. 7, where we report the L 2 and H 1 errors (top and bottom, respectively) with respect to the manufactured solutions as a function of the number of points discretizing the curved boundary for the two set of meshes (left and right). It is evident that the proposed scheme is robust with respect the increasing unbalance of the edge lengths as the errors for each polynomial order remain approximately constant.

Test case 3
This test case is aimed at studying the robustness of the scheme to spatial variability and strong anisotropy of the diffusion bilinear form with spatially variable anisotropy ratios. We verify convergence of the proposed VEM scheme by solving our equations on the same families of meshes of Test Case 1 with r = 2 , k = 5 and two different values for a, a = 0.5 and a = 2 . Figure 2, second and third rows, shows the surface and the spatial distribution of g ii for k = 5 and a = 0.5 and 2, respectively. The spatial variability of the anisotropy ratio (ratio between g 11 and g 22 or its inverse) varies between 1 and 150 for a = 2 corresponding to large basis vectors for the tangent plane. This test case challenges the ability of the discretization scheme to handle large and spatially varying anisotropy ratios. Figures 8 and 9 show the numerical convergence of the L 2 and H 1 norms of the error as a function of h . We note that asymptotic behavior of the error is reached as soon as the mesh size is able to resolve the spatial scales of variation of the metric tensor. For this reason the convergence lines in the figures are obtained by interpolation of the last two point values for each polynomial order. Pre-asymptotic convergence is more evident for the higher order polynomials. We attribute this behavior to the smoothing effects of a lower order interpolation. Indeed, the first two or three point values for all polynomial orders display at least the same convergence order of linear polynomials.
The high degree of anisotropy of this test case causes a loss of convergence in the higher polynomial orders. This behavior is more evident for the polygonal meshes. To better quantify this convergence loss, Tables 1 and 2 report the convergence order sequence for the triangulations and polygonal meshes, respectively. We note that, for first and second order polynomials, almost optimal order of convergence is reached for both mesh sets and both a = 0.5 and a = 2 . In the higher orders, a clear loss of at least half an order is evident. This can be attributed to difficulties in resolving the large anisotropy ratios that are typical of this test case. The fact that this occurs only for the higher orders is due to the ill-conditioning of the resulting linear system that controls the ratio between the norms of the residuals and the errors in the linear system solver.

Test case 4
This test case shows the convergence properties of the proposed framework in a multiple charts setting. We discretize the full unit disk with 5 polygonal mesh levels following the strategy reported in Test Case 1, so that the level = 0 is characterized by 100 cells and the last one ( = 4 ) by 25600 cells. The VEM solution, reconstructed on = S 2 using nodal only degrees of freedom, is shown in Fig. 10 for the coarsest mesh.
The experimental convergence on these mesh levels is reported for the L 2 and H 1 norms of the error as a function of h in Fig. 11 and in Table 3. The numerical results show that the proposed approach is functioning as expected and that the use of two different charts doe not influence the optimal convergence of the scheme. Obviously, this is a very favorable case as the stereographical projection produces charts and, if necessary, transition maps that are sufficiently smooth to allow high order. In the future it will be important to study how to derive charts and transition maps with specified regularity for different surfaces, possibly starting from the work of [44].

Conclusions
We have developed an arbitrary-order virtual element method for the discretization of elliptic surface PDEs. The approach employs a local parametrization Table 1 TC3: experimental errors and convergence rates for the triangular mesh set  of the surface to properly re-define the PDE on the local chart. This allows the straight-forward definition of a two-dimensional VEM discretization at all polynomial orders, overcoming the difficult task of the consistent approximation of the surface and of the distance function of its tubular neighborhood. The drawback of the approach is that the geometrically intrinsic form of the PDE contains the metric information, which may be strongly non-isotropic and highly variable in space, depending the regularity of the surface. The choice of the VEM scheme  is motivated by the need to ensure robustness and high order of convergence for these anisotropic and spatially variable coefficients.
The developed scheme has been tested on several numerical examples showing varying degrees of regularity. In fact, optimal orders of convergence up to 5 has been reached for surfaces with relatively small curvatures. Only when curvatures and metric information become extremely large loss of convergence is noticed. This loss of convergence is related to the presence of strongly anisotropic diffusion tensors and strongly aligned advective fields due to the behavior of the metric tensor. Handling strong anisotropy is still a major challenge in the numerical solution of PDEs by the virtual element method, and is left for future research. This difficulty can be relaxed also by employing multiple charts that decrease the anisotropic characteristic of the metric tensor. However, proper regularity of the transition maps between the different charts must be ensured to achieve full order convergence. To verify the ability of our formulation to work with multiple charts we tested the proposed scheme on the full sphere by employing two charts arising from the stereographical projection. Future work will be addressed to define for general surfaces appropriate multiple charts with regular transition maps starting from the work of [44].
One of the major advantages of the developed VEM formulation is that can be used efficiently to minimize geometric errors of curvilinear boundaries. We have tested our approach on an hemispherical surface discretized by a fixed number of polygonal cells, where the boundary edges formed by an increasing number of nodes. The resulting errors were independent of the edge discretization, showing the robustness of the VEM scheme in this situation.
Funding Open access funding provided by Università degli Studi di Padova within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not 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/.