Virtual Elements on polyhedra with a curved face

We revisit classical Virtual Element approximations on polygonal and polyhedral decompositions. We also recall the treatment proposed for dealing with decompositions into polygons with curved edges. In the second part of the paper we introduce a couple of new ideas for the construction of VEM-approximations on domains with curved boundary, both in two and three dimensions. The new approach looks promising, although sound numerical tests should be made to validate the efficiency of the method.


Introduction
The aim of this paper is twofold. On the one hand, it is an attempt to give an idea of the nature of the Virtual Element Methods (VEM) that were introduced around ten years ago for the numerical solution of PDEs. This part would be essentially aimed at researchers that are not specialists in Numerical Analysis and Scientific Computing. There, on a polygonal domain Ω in 2d (polyhedral in 3d) we will recall the classical VEM approach. In the second part we will first recall the approximation presented in [5] for dealing with 2d problems on a domain with a curved boundary. Then we will introduce new ideas which can simplify the previous approach, and, at the same time, allow the extension to 3-dimensional problems where the computational domain has a boundary (or part of it) described by a curved surface.
Note that, in many cases, an approximation of the boundary by a polyhedral surface would be too rough for the applicative point of view, unless the mesh is very fine (which, in three dimensions, could be outrageously expensive). Moreover, the approximation of the curved boundary with straight lines in 2d (flat polygons in 3d) might produce a loss of convergence. This occurs with Finite Element approximations, and is known as Babǔska paradox (see [4]). Since on triangular decompositions Virtual Elements of degree k = 1 coincide with piecewise linear Finite Elements, it is highly possible that a similar situation might happen with VEM, even in more general cases. Thus the interest of devising VEM-approximations able to deal with curved boundaries (see, e.g., [13], [7], [2], [5]).
The possibility of dealing with curved boundary is one of the main reasons that ensured the enormous success of the so-called IsoGeometric Analysis (IGA; see e.g. [17]), where the solution is approximated by splines of the same type used (in almost 100% of industrial design) to describe the boundary of the domain.
Here, however, contrary to IGA, we do not need to correlate the two ways of describing the boundary and discretizing the problem. In some sense, the type of description of the boundary and the type of discretization of the problem are, in essence, totally independent one of the other and the two choices do not need to be correlated. This, obviously, does not mean that you can choose the two discretizations independently, but that you can choose the type of discretization (the decomposition and the degree to be used in each polyhedron), essentially, without considering the type of description of the boundary: by splines, trigonometric functions, or just (brutally) by a given analytic expression.
We will use the Poisson problem as super-simple example in order to explain the construction of the numerical method (passing from the PDE to a system of linear equations that can then be solved on the computer).
The problem of constructing three-dimensional Virtual Element spaces on general polyhedra having one or more curved faces is still open (at the fully general level). Here we propose several possible ways of tackling the problem.
We point out that the "discretizations of order k" described here will satisfy the so-called patchtest of order k (very popular among Engineers) that, very roughly, says that: whenever the exact solution is a polynomial of degree k, then the solution of the discretized problem will coincide with the exact solution. Although the patch-test is only a necessary, but not a sufficient condition for convergence, it makes you feel more confortable to have it satisfied.
From a certain point of view, insisting on the correct treatment of curved boundaries, the approach presented here could constitute a viable alternative to the classical Iso-Geometric Analysis (already mentioned), allowing a wider choice of the domain decomposition, payed with a lower regularity of the subspace (basically H 1 instead of the H p regularity of IGA). Note that, compared with Finite Elements, the methods proposed here would indeed allow a much easier use of more regular subspaces (say, like H 2 or H 3 ) although using more degrees of freedom than IGA or, more generally, than splines.
An outline of the paper is as follows.
After a short reminder of classical notation, we will start recalling classical VEM spaces on polygons/polyhedra, and their use for the approximate solution of PDE problems.
This will include both the definition of the VEM spaces, and a hint on the classical ways of using them to solve PDE problems on the computer. We anticipate that, within each element of the decomposition, the VEM spaces, in general, contain (together with all polynomials of a certain degree) also functions that are not polynomials, and whose point values are not known. Typically, in each element they will be solutions of suitable boundary value problems whose boundary data and right-hand side are identified by a certain number of parameters (the degrees of freedom). The solution of such element-wise boundary value problems is clearly out of the question, and each element of the Virtual Element space will be known only through its degrees of freedom. This implies that its manipulation will not be as easy as it happens when using just polynomials (or splines), and we will briefly indicate how all this can be dealt with on the computer.
Then we will discuss the use of polygons with a curved edge, recalling what was done in [5] and proposing some new interesting alternatives. At this point we will be ready to discuss the treatment of polyhedra with a curved face. For simplicity we will deal only with the case in which each element has just one curved face, that is part of the boundary of the computational domain, but the extension to more general cases should not be too difficult.
We believe that the approaches that we suggest might open new perspectives on the treatment of curved edges/faces, and lead to new more suitable methods. However, sound numerical tests will have to be performed to validate the efficiency and, most important, the accuracy of the methods, but this goes beyond the scopes of the present paper.

Notation
Throughout the paper we will use the common notation for functional spaces L 2 (D) and H s (D) on a domain D, with scalar product and norm (·, ·) 0 , · 0 ( and (·, ·) s , · s , respectively); H 1/2 (Γ) will denote the space of traces of functions of H 1 (D) on the boundary Γ (see [21]).
When convenient, we will use the notation ∆ 2 and ∆ 3 to indicate the two-dimensional or the three-dimensional Laplacian, respectively.
For k ≥ 0, P k (D) will denote the space of polynomials of degree ≤ k on a domain D. As common, P −1 = {0}.
Throughout the paper, C will represent, as usual, a fixed positive constant, not necessarily the same from one occurrence to another.
We also point out that, given a domain Ω and a subset P of it, by saying that "P is internal to Ω" we mean that the closure of P is a subset of • Ω . Given a finite dimensional linear space V and an integer number S, a linear mapping is said to be a set of generators for V if it is surjective. If moreover G is also injective, then it is one-to-one, and its inverse D : V → R S could be represented as a set of S linear operators (δ 1 , ..., δ S ), each in L(V, R), that are then called degrees of freedom: Note that, even when G is not injective (but only surjective) we could always construct a right inverse D : that however, this time, will not be unique.
In what follows, a set of generators for a space V will be typically denoted by G V .

Projectors in VEM spaces
As we shall see, most of the Virtual Element spaces that we are going to construct will be made of functions that are solutions of local (i.e., element by element) boundary value problems for PDEs. Typically, in two dimensions, we will know explicitly only their values at the interelement boundaries, and possibly some of their moments inside each element. The computation of other quantities, as for instance their pointwise value inside the elements, will be totally out of reach. Consequently, we will often (actually: almost always) use in their place some suitable projection on polynomial spaces. There too, from the information that we have, some projections will be computable, and others will not. For the moment we just briefly anticipate the classical structure of VEM spaces of order k on a polygonal element P , in order to see some projectors that are actually computable out of the degrees of freedom.
The general structure of a VEM space with order of accuracy k ≥ 1 on a polygon P is: Typically the degrees of freedom for the spaces in (3) where the degrees of freedom for the spaces in (4) should provide directly • the values on each edge of P , and • the moments up to the order k L in P .
Clearly the spaces (3) would be obtained for k L = k − 2.
Here below we will see some computable projections that we are going to use in the sequel.
* The H 1 0 (P )-projection. Given a polygonal element P , for each v ∈ H 1 (P ) we define its projection Π ∇ k v onto the space P k as the solution, in P k (P ), of Actually, (5) identifies Π ∇ k v only up to a constant, that can easily be fixed, for instance, with the additional requirement that The left-hand side of (5) is a product of polynomials, and is obviously computable. Integrating the right-hand side by part we have and we will be able to perform our computation as far as we know explicitly the moments of v up to the order k − 2 on P , the moments of v up to the order k − 1 on each edge of P .
As we anticipated, normally on each edge we actually know directly the whole (polynomial) values of the functions of our discretized spaces, and the integrals on the edges can be computed exactly. Moreover, the first term on the right-hand side of (7) can be easily computed out of the degrees of freedom. Now let us see some other important projectors.
* The L 2 projection of the gradient. With an almost identical argument we can compute the Indeed, integrating by parts the right-hand side we have and both terms are immediately computed as in (7). Note that, referring to (4), we can actually compute the L 2 projection of ∇v on (P s ) 2 whenever v ∈ V k,k L with k L ≥ s − 1.
* The dofi-dofi projector. It is also important to note that there are other projectors from V k (P ) to P k that are computable out of the degrees of freedom Dv. The simplest one, that we call Π D , would be defined for each v ∈ V k (P ) as the (unique) solution in P k of where obviously (· , ·) R S is the usual Euclidean scalar product in R S , and S is the number of degrees of freedom of V k (P ). Note that this can be done even if in (9) the mapping D is not representing the degrees of freedom, but just any other identifier of the type (2), provided that D is injective from P k to R S . Indeed, as far as S is equal to the number of degrees of freedom, (9) is just a pompous way to say that D(Π D v) = Dv.
* The L 2 projection. In general, the L 2 -projection of an element v on P s (P ) can be computed only when we know the moments of v of order up to s. If we know the integrals P vm j for j = 1, ...S, where the m j are a basis for P s , and S is the dimension of P s , then we can orthonormalize the m j . Denoting by m j the orthonormal basis, i.e., such that * The Serendipity-like projectors. A particularly relevant class of possible alternative projectors is given by the Serendipity-like projectors. The basic idea (or, so to speak, the seed of it) is the following: For an element P we denote by η(P ) the minimum number of straight lines necessary to cover the whole boundary ∂P . Note that in our setting two edges (consecutive or not) may lie on the same straight line (see an example in Fig. 1). Then we observe that for every k < η(P ) a polynomial of degree k will be completely identified by its value on the boundary ∂P . Then, always for k < η(P ) we can define a projector Π η k from V k (P ) to P k (P ) defined by which clearly has a unique solution. For k ≥ η(P ) we will have to add some internal information. Typically, we may choose an integer r with k − η(P ) ≤ r ≤ k − 2, add to (10) the condition and consider the pair of equations (10)- (11), possibly to be solved in the least squares sense, mimicking what is done in the construction of Serendipity VEM spaces (see [11]). An important remark is in order, concerning the L 2 − projection of the gradient, in the context of the stabilization procedures that will be discussed in the next pages.

Remark 2.
In several cases, it would be convenient (and we daresay very convenient) to construct polynomial projections Π(∇v) of the gradient that satisfy Π∇v = 0 => ∇v = 0 (12) (or, in oher words, that are injective). Clearly, if the polygonal element that we are considering has many edges (and consequently, even for a low k, many degrees of freedom), condition (12) would require that we project onto a big polynomial space, and will make the whole idea too expensive. However, as we shall see in a while, the advantages of (12) are relevant in many cases, and the decision on whether to enforce it or not could be delicate in several circumstances.

The simplest model problem
Let Ω ⊂ R d , d = 2, 3, be a domain with boundary Γ. We assume that Γ is the union of two (for simplicity, connected) parts, denoted by Γ D and Γ N . As usual we will assume thatΓ We consider the simple model problem Setting, for ϕ ∈ H 1/2 (Γ D ), the variational formulation of (13) can be written as where a(u, v) is the bilinear form defined by and and the boundary integral could be replaced by a suitable duality for a less regular g N . For the sake of simplicity, here we will discuss, separately, only the two cases: Full Dirichlet (when Γ D ≡ Γ; actually, very simple) and Full Neumann (when Γ N ≡ Γ).
In both cases, we will assume that we are given a sequence of decompositions {T h } h of Ω in polytopes P with diameter h P , and we indicate with |h| the maximum of h P for P ∈ T h . We will also make the usual assumptions that each polytope (of each T h ) is star-shaped with respect to a ball of radius ρ P ≥ Ch P , and each edge/face of a polytope P has length ≥ Ch P .

Virtual Elements for 2D polygons 2.1 The local VEM spaces
For the sake of simplicity, we begin by recalling the original plain vanilla VEM spaces on twodimensional polygons with straight edges, as presented in [8] and already anticipated in (3). Let P be a polygon, and let k be an integer ≥ 1.
We define first and then It is clear that all polynomials of degree ≤ k belong to V k (P ), and that (as anticipated in the previous section) an element v of V k (P ) is uniquely determined by: the values of v at the vertices of P, For reasons that will be clear soon, the degrees of freedom (19)- (21) should be suitably scaled so that they all scale in the same way. We do not enter the details here, for which we refer for instance to [10].

The two-dimensional "full Dirichlet" case
Let Ω ⊂ R 2 be a polygonal domain with boundary Γ. Let moreover g and f be smooth-enough functions defined on Γ and in Ω, respectively. Our model problem (14) becomes Since now Γ D ≡ Γ, we can simply write H 1 g (Ω) instead of H 1 g,Γ D (Ω), and (14) becomes For a generic ϕ defined on Γ (here either ϕ = g or ϕ = 0), and for an element P ∈ T h , we define the local spaces as with their global counterpart The degrees of freedom in V ϕ h (Ω), on top of the value of ϕ on Γ, will obviously be: • The values at the internal vertices of the decomposition, • (for k ≥ 2) The moments on each internal edge up to the degree k − 2, • (for k ≥ 2) The moments inside each element up to the degree k − 2.
It is immediate to check that (as seen in the subsection 1.2) we can follow (5) and (6) and define, in each element P , the projector Π ∇ k , that will be computable using the above degrees of freedom. Next, always on each element P , we consider the restriction a P (u, v) of the bilinear form a (u, v) to P , and then on each V ϕ k (P ) (as defined in (24)) we set The first term in the right-hand side of (26) is the consistency term which is in general singular, and the second term is needed for stability. We can take as S P (u, v) any bilinear form that scales as a P (u, v) and is positive on the kernel of Π ∇ k . In particular we would require that there exist two positive constants α * and α * such that The simplest example would be (setting S = number of dofs in P ): where dof i is the i th degree of freedom properly scaled in such a way that S P (u, v) scales as a P (u, v). This is the reason why the degrees of freedom need to be properly scaled, as briefly anticipated at the end of Sect. 2.1. But other choices can be convenient, as where u t and v t are the tangential derivatives of u and v, respectively. Needless to say, the above expressions for S P (u, v) can be multiplied by a fixed constant (say, 5, or 1/5, or whatever fixed constant you might like, as far as it is independent of h).

Remark 3.
Instead of using the Π ∇ k operator in the definition of the consistency part of the discrete bilinear form (26), we might use the L 2 -projection of the gradient, and define Actually, this choice is preferable for more general problems, in particular in presence of variable coefficients, since there the Π ∇ k operator might produce a loss of order of convergence for high values of k (see [12]).

Remark 4.
For the case that we are considering here, most choices of the stabilizing form S P will give quite satisfactory results. However, for more complex problems the choice of the stabilizing term (or, for a given term, the choice of the coefficient to be used in front of it) might become quite delicate: a poor choice would not jeopardize the convergence of the method (in the limit for h going to 0), but would give poor results for the reasonably affordable decompositions.

Remark 5.
In some sense, we might say that, in several applications (just to take an example: for the extension of the present approach to fourth order problems), the choice of the right stabilizing term (with the right coefficient) can become the most delicate point (see [14], [15] For a detailed treatment of the right-hand side we refer to [8].
Here we just recall the definition of f h on each polygon. With Π 0 k := L 2 −projection operator onto P k , we set The discretised problem reads now: Remark 6. We point out that, if the boundary value g is not the trace of a polynomial, we will loose the property that every local space V g k (P ) contains all the polynomials of degree ≤ k on P . If, for some reason, we are interested in such property, then we should take first a piecewise polynomial approximation g h of g, and then use V g h k (P ) and V g h h (Ω) in place of V g k (P ) and V g h (Ω), respectively.

A tiny historical note
As explicitly claimed from the very beginning (see e.g. [8]) the Virtual Element Methods are a direct follow-up (and initially almost a simple re-formulation) of Mimetic Finite Differences (MFD, see e.g. [6] and the references therein). In MFD the unknowns are just the degrees of freedom (as, indeed, it always happens in the final computer code) and the formulation is done directly in terms of them (typically, nodal values or averages); in other words, with MFD we do not have a discrete functional space playing the role that here (and generally in VEMs) is played by V h . There too, as in VEM, the contribution of each element to the final stiffness matrix, for instance for a problem as (13), is the sum of a consistency part and of a stabilizing part. The first part is obtained considering first the polynomial subspace of dofs, meaning the subspace whose elements are obtained, each, as the values at nodes and the averages of a single polynomial of the prescribed degree k. So, for k = 1 we will have the subset (of dimension 3) generated by the (say) nodal values of 1, x, y. The consistency part of the MFD stiffness matrix will then act only on this subspace, and there it will reproduce the corresponding expected behaviour of differential operators and/or integrals that appear in the continuous bilinear form (as (15)). The stabilization part, instead, will vanish on the above subspace, and will be symmetric-positive-definite on its complement (more precisely, in the complement of its kernel). In other words, we have the perfect analogue of (26)(that does actually coincide with the one used in MFD for a judicious choice of the respective stabilizing parts).
The two methods, MFD and VEM, are just part of a pack of different approaches to the treatment of PDEs using decompositions of the computational domain into elements whose shape is more general than just simplexes and boxes. Among the most ancient ones (together with MFD) another building block is surely made by Discontinuous Galerkin Methods (see, e.g. [3], or the most recent [19], [18]). There the discrete spaces are just made of polynomials, and suitable tricks are needed in order to enforce some kind of weak continuity from one element to another. Many variants are related with them, from Hybridizable Discontinuous Galerkin methods (where, roughly speaking, a weak continuity is enforced via Lagrange Multipliers), or Hybrid-High-Order and Weak Galerkin methods, where one works with two piecewise polynomial spaces: one at the boundary and one inside. See e.g. [16], [18], [22] and the references therein. See also [20] for similar ideas. Needless to say, if necessary one could always use one type of discretisation in some elements and another one in other elements, as, in the end, the computer will always solve "just a problem in R n " having as unknowns the degrees of freedom.

The two-dimensional "full Neumann" case
We consider now the full Neumann case, i.e., Γ N ≡ Γ. We recall that, in this case, g N and f must satisfy the compatibility condition Γ g N ds = Ω f dx due to the Gauss divergence theorem. We also recall that the solution u will be determined only up to an additive constant, that in the computer code can be fixed, as common, just by fixing its value to be equal to zero at some point (usually, a vertex of Ω). From now on, when discussing a full Neumann problem, we will implicitly assume that we have chosen, once and for all, a vertex of the discretization, and prescribed that the solution of the continuous problem u, and all the elements of the discrete subspaces, vanish there.
The discrete problem differs from the full Dirichlet case only in the right-hand side, upon performing a slight change in the definition of the local and global discrete spaces (24) and (25). By defining V k (P ) = {v ∈ H 1 (P ) s.t. : v |e ∈ P k (e) ∀ edge e, and ∆v ∈ P k−2 (P )}, the discrete problem reads: In (30) a h (u h , v h ) and (f h , v h ) are the same as in the previous subsection, and < g N , v h > Γ is a boundary integral that can be computed since the functions v h are polynomials on each edge, completely known.

Recalling 3D VEM for classical polyhedra
The extension of the above construction to the three dimensional case is, in some sense, both immediate and tricky. Indeed, the first attempt that comes to mind, given a polyhedron P and an integer k, is to extend slavishly what we did in the two-dimensional case. Let us first examine the full Dirichlet case.
To begin with, we define V h (Ω) := {v ∈ H 1 (Ω) ∩ C 0 (Ω) such that: v |e ∈ P k (e) ∀ internal edge e, We remark that V h (Ω), as defined in (31), is infinite dimensional (we are not making requirements on the values of v on ∂Ω). Then, given a (smooth enough) function ϕ defined on ∂Ω we can mimic (24)-(25) and restrict our space to that, now, will be finite-dimensional. In particular, for a fixed given ϕ, the degrees of freedom in (32) will be: The values at the internal vertices of the decomposition, Then, apparently, one can follow, for a three-dimensional problem like (13), the same path that we used for the two-dimensional case, both for the full Dirichlet and the full Neumann case.
Given a decomposition T h of Ω into polyhedra P , an integer k ≥ 1, and a smooth-enough function g defined on Γ, we can define the finite dimensional subspaces V g h and V 0 h of V h (Ω). Still following slavishly the 2-dimensional path, for each P ∈ T h and for each virtual element function v we can define its H 1 0 (P )-projection Π ∇ k v as the unique solution (up to a constant that can be easily fixed) in P k of Ooops! When we attempt to compute the right-hand side of (34) using the degrees of freedom (33), we have Now, the second term in the right-hand side of (35) does not cause any trouble: for p k in P k we have that ∆p k ∈ P k−2 and the term can be computed using the degrees of freedom (33). But for the first term in the right-hand side of (35) we would need to know the moments of v on each face up to the order k − 1 (the degree of ∂p k ∂n ), while in (33) we have the moments only up to k − 2. The way-out, as presented first in [1], is: where Π ∇,F k v is the two-dimensional projection of v onto P k (F), as defined in (5)-(6), whose computation, in turn, on each face F requires the moments of v on F only up to the order k − 2. One might consider this as a typical use of the approach described in subsection 1.2.
For the Neumann case, also the treatment of the right-hand side needs a more careful approach than in the 2D case. The integrals on Γ N might be computed as in (36): Finally, the term (f, v h ) can be treated as in (28). For more details we refer for instance to [9].

Remark 7.
As pointed out in subsection 1.2, we might prefer other projectors (instead of Π ∇ k ) and use the "projected v" in place of v for other quantities that cannot be computed directly from the degrees of freedom of v.
Once Π ∇ k has been defined, we can follow step-by-step (with obvious minor modifications) the path of the two-dimensional case (as suggested in [1]) both for the full Dirichlet and the full Neumann case.

Remark 8. Similarly to what was discussed and suggested in Subsection 2.3, here too we could
consider the possibility of using other degrees of freedom, in particular on the "Neumann boundary faces". This could be seen as: As usual the unknowns in the computer are the degrees of freedom. To these degrees of freedom, in most elements (with few exceptions, as we shall see) we could associate a space of local functions, some explicitly computable (e.g. polynomials), others not explicitly computable but well defined, as solutions of local PDE problems (that we will not solve, and use instead their projections on polynomial spaces). Other degrees of freedom might be used directly to compute projected polynomials, more in the spirit of Mimetic Finite Differences. What we will carefully preserve will always be the Patch test, ensuring that the method will be exact whenever the exact solution is a polynomial.

2D VEM for polygons with a curved edge
Assume now that we have a "polygon" P with one curved edge, that we call η, belonging to ∂Ω. The case of polygons having two or more curved edges (always on ∂Ω) could be treated in a very similar manner.

The full Dirichlet case for 2D problems with curved boundary
For the full Dirichlet problem, we can follow slavishly what has been done in (24)-(25), as well as in (26). Essentially, since the values of the VEM spaces on the boundary are assigned (both for test and trial functions) the fact that the boundary is curved is just a minor nuisance, requiring the computation of integrals on the curve. As we shall see, this will not be the case for the Neumann boundaries.

The full Neumann case for 2D problems with curved boundary
We recall in this subsection the procedure introduced in [5], which we refer to for more details. In a natural way, on the curved edge η we introduce Remark 9. Definition (37), simple as it may seem, opens the door towards a more delicate discussion. Indeed, for a curved edge η the dimension of the space P k (η) might be, in certain cases, far from obvious. In facts, the dimension could change between a minimum of k + 1 (when the edge η is straight), to a maximum of k(k + 1)/2 when η is curved enough (here meaning that the only element in P k (R 2 ) that vanishes identically on η is the polynomial ≡ 0). If we want to be allowed to treat systematically η as curved we must then consider k(k + 1)/2 generators (typically, the values at k(k + 1)/2 points as in Figure 3), but the actual dimension might be as few as k + 1. We shall come back to this problem (and related troubles) pretty soon.
For the moment we can keep going on the same track as before, setting and v |η ∈ P k (η)}, and then V η k (P ) = {v ∈ H 1 (P ) such that v |∂P ∈ B η k (∂P ) and ∆v ∈ P k−2 (P )}.
So far so good. Now we observe that an element v ∈ V η k (P ) will be uniquely determined by: • the values of v at the vertices of P, (39) It is crucial to point out that, in (41), we did not forget to say "for k ≥ 2". Indeed, as we already discussed in Remark 9, for k = 1, on a straight edge e we do have clearly that all polynomials of degree ≤ 1 vanishing at the two endpoints of e would be identcally zero on the whole edge. But this is not anymore true for a curved edge, where a polynomial of degree 1 in R 2 , vanishing at two distinct points would (indeed!) vanish identically on the straight line connecting the two points, but might very well be different from zero, for instance, on a circular arch passing through the same two points. And, in fact, hic sunt leones (that is: here come the difficulties). Actually, when analysing and coding a numerical method, the basic step is often to define the finite dimensional space where we are looking for the solution, and the basis to be used for it (the degrees of freedom). The difficulty here is in finding the dimension of such a space, or, in the terminology of Numerical Analysis, the correct number of degrees of freedom. In particular, to forbid curved edges that are nearly flat would be quite cumbersome in applications: the crucial point would be the definition of nearly flat: resonably clear for the human mind, but a nasty source of troubles in a computer code, taking into account that there the definition must be quantitative, and a minor change in it can produce a significant change in the numerical solution. And, sticking for simplicity to the case k = 1 (but the difficulty pops up for every k), on a general curved edge η the space P k (η) has dimension 3, but when the edge η is almost a straight line the number of parameters necessary to identify an element of P k (η) tends to 2 (in a sense difficult to be made precise, but surely prospecting troubles). Referring to Figure 2, for k = 1 the dimension of V k (P ) depends on the shape of the curved edge. In particular, in exact arithmetic the dimension of V 1 (R) is clearly 4. But in the computer we should treat it carefully, since a simple minded treatment might lead to singular or nearly singular matrices depending on the number of digits that we are using.
The typical choice for VEMs (as done in [5]) is to treat every edge η that could be curved as curved, and to consider as unknowns associated with the edge the values (on that edge) of all polynomials of P k (R 2 ) that vanish at the two endpoints of the edge, that therefore would sum up to (k + 1)(k + 2)/2 − 2 parameters. Typically, for a "curved" edge η this is done considering first the straight segment Q connecting the two endpoints of η, then considering an equilateral triangle T η having Q as one of its edges, and finally considering on T η the traditional degrees of freedom that one would have for P k on that triangle: • the value at the vertex not belonging to η, • for k ≥ 2 the values at k − 1 equally spaced points on each edge, and • for k ≥ 3 the values at (k − 1)(k − 2)/2 internal nodes.
(These will be called generating points). Hence, on an edge η that has been declared as curved it will be simpler (in the computer code) to replace (41) with the value at the generating points.

Normal dof
Trace generator k=1 k=2 k=3 Normal dof Therefore, even when the so called "curved edge" is in fact straight, we will associate with it (k + 1)(k + 2)/2 − 2 parameters: the values at the generating points (or, alternatively, the moments on η against the traces (on η) of all P k polynomials vanishing at the two endpoints of η). Clearly this will be a set of generators (in the sense of (1)), but not a set of degrees of freedom, in the traditional sense.
All this, unfortunately, will require some additional paraphernalia. Indeed, the obvious path, once we defined the local spaces (38), is to collect them in the global space where T h is our decomposition of Ω into "polygonal" elements, possibly with a curved edge. Then we would like to construct an approximate problem of the type (29). But in the computer code, the unknowns of the discretized problem cannot be the elements of V h (that are functions), and must be, instead, their generators (that are elements of some suitable R N ). Remember that: given the generators, the function is uniquely determined, but not the other way around, as (for instance when the edge is straight) the same function on η could be generated by several different generators. Hence, once we have defined the global spaces in (43), we must consider the space M V h (M is for Midwife!) given by that, in general, might have a dimension bigger than that of V h . Hence we will have a linear mapping (C is for Children) that is, in other words: For every given element v ∈ M V h , C selects a unique element v = C(v) ∈ V h , generated by v, that will be the function that we are interested in.
We observe that, given a generator v ∈ M V h , the polynomial Π ∇ k (C(v)) can be computed as in (5)-(6) since the function v = C(v) is defined all over the boundary of the element. Needless to say, other projectors from V k (P ) to P k could be computed, if needed, for instance using the degrees of freedom as in (9). Then, with obvious notation, for every element P we can define the bilinear form where again S P is a bilinear form on R N that scales like a P (Π ∇ k (C(u)), Π ∇ k (C(v))) and is positive on the kernel of D(Π ∇ k C). Note that, even though D (a right inverse of C as in (2)) is not uniquely defined (unless C is injective), D(Π ∇ k (C(v))) can be uniquely defined since Π ∇ k (C(v)) is a polynomial. Then, as we did for polygons with straight edges, we can collect the definitions (44) setting The discretised problem reads now: The treatment of the term (f h , C(v)) is discussed in [5], which we refer to. Neumann boundary conditions were not dealt with in [5], but the computation of the second term in the right-hand side does not pose additional difficulties.

Remark 10.
For problems where curved edges occur inside the domain Ω, one has to be careful with the choice of generators. Indeed, an edge that is internal to Ω will naturally belong to the boundary of two different elements. But the functions v ∈ V h (as defined in (43)) must be single valued on the common edge, so that the set of generating points used must be the same for the two elements. Moreover, as pointed out in [5], it is better to stabilize the generators for the curved common edge only once: choosing, once and for all, one of the two elements having the edge in common, and then using them in the stabilizing term only when dealing with the chosen element.
The approach presented in the previous section does not extend easily to three-dimensional problems. In this subsection we introduce a couple of new ideas that could be used instead.

Using a subset of degrees of freedom.
Another way to set the values of VEM functions on curved edges of ∂Ω (in particular on the part of ∂Ω where natural boundary conditions are prescribed) would be to use the degrees of freedom on the straight edges plus the internal ones. To explain the procedure, let P be a polygon, with just one curved edge η belonging to ∂Ω. We would like to define a space of functions like (17)- (18), that is, a space of continuous functions, polynomials of degree k on each edge, and with Laplacian polynomial of degree k − 2 in P . The problem is that we do not know how to define these functions on the curved edge. The information that we have (candidates to become degrees of freedom) are: the values of v at the vertices of P, These conditions are not enough to individuate a VEM-function, but they are enough to individuate a polynomial of degree k. Indeed, to make the simplest possible example, let us consider a trianglelike element, i.e., with two straight edges and one, η, curved. For every integer k ≥ 1 we see that using the values of v on the two straight edges (amounting to 2k +1 dofs), and the internal moments up to the degree k − 2 (amounting to k(k − 1)/2 dofs) we can identify uniquely a polynomial p k in P k . Indeed, 2k + 1 + k(k − 1)/2 = (k + 1)(k + 2)/2 ≡ dim(P k ).
With a given VEM-function v we can associate a polynomial p * k ∈ P k , computed through the conditions It is trivial to check that (always in the case of a triangle-like element) conditions (46) are a set of unisolvent degrees of freedom for P k , so that p * k is unique. Then we can take Now conditions (45) plus (47) determine v uniquely, and we point out that, by construction, global continuity is guaranteed. We observe that, in this particular case of triangle-like elements, the function v is actually a polynomial, coinciding with p * k on the whole element. This result could suggest to have a decomposition made of triangle-like around the boundary. The approach would be simple and clean, but not really in the spirit of VEM, whose strong point is to allow polygons of arbitrary shape.
Let us then consider a general polygon, with a number of straight edges higher than 2. Then we would have more information than necessary, i.e., the number of conditions (45) would be bigger than the dimension of P k . In this case we should use a least-square solution, keeping of course fixed the values at the two endpoints of η to guarantee the global continuity. More precisely, if N is the number of conditions (45), ordered in such a way that the values at the two endpoints of η are the last two, we solve the following problem: find p * k ∈ P k (P ) such that p * k = v at the two endpoints of η, and Then, as before, we define v |η = p * k|η .
Once the function v is individuated by (45) and (49), we can define Π ∇ k v as in (5) (or Π 0 k−1 ∇v as in (8)) and write the discrete bilinear form as in (26) (or as in (27), respectively).

Remark 11.
This strategy cannot be used, as such, if the curved edge is inside Ω and common to two or more elements. In this case one might think of a master and slave approach (only one element sets the degrees of freedom to be used on the curved edge) or take a suitable combination of the effects of the elements sharing the curved edge. Once the functions are defined on the curved edges, everything goes along the same lines used for polygons with straight edges, and the discrete problem can be written exactly as in (30).

Remark 12.
A similar approach, although conceptually very different, would consist in using, associated with the curved edge, suitable degrees of freedom to be used in the construction of the various projectors, without defining a functional space within the element. In other words, in the element with a curved edge we would consider only degrees of freedom without introducing a functional space. Hence, associated with the element with the curved edge we will still have a local stiffness matrix (as we would have when using Mimetic Finite Differences). This would correspond, somehow, to use Finite Elements (or Virtual Elements) in the other elements, and Mimetic Finite Differences in the elements having a curved edge with natural boundary conditions (in the spirit, somehow, of Subsection 2.3). Although this might look as a monster from the theoretical point of view, the computer code will not suffer (in particular if you already have a VEM code and an MFD one) as, in any case, it will deal only with degrees of freedom and never with functions. The only drawbacks would appear only in the proof of error estimates, although, as a combination of two reliable methods, it is reasonable to expect the usual level of accuracy. The above procedure would be in this case the following: we use again a least-square approach to compute a polynomial p * k ∈ P k , this time without any need for fixing its values at the two endpoints. Then, to individuate a VEM-function v we add to conditions (45) the moments of p * k of order k − 1 as degrees of freedom on η: Conditions (45) and (50) allow us to compute the projection Π ∇ k v as in (5) (or Π 0 k−1 ∇v as in (8)), and we can write the discrete bilinear form as in (26) (or as in (27), respectively).

Coupling FEM and VEM: the superimposed polygon
Another possible approach that might be considered (always on the part of the boundary where Neumann boundary conditions are imposed) is the use of a superimposed polygon P including Ω, such that the two boundaries (∂ P and ∂Ω) are never closer to each other more than Ch for some fixed C > 0. Let then T ( P ) be a decomposition of P in polygons (for simplicity, in convex polygons), that naturally produces a decomposition of Ω consisting of normal polygons (corresponding to the elements of T ( P ) that are all contained in Ω) and polygons with a curved edge (corresponding to the restrictions to Ω of the polygons in T ( P ) that contain parts of Ω). For simplicity we assume that there are no polygons in T ( P ) that do not contain at least a part of Ω.
A simple way of realizing this is to construct a ribbon of quadrilaterals, as in Fig.4, around ∂Ω, that naturally defines two polygons: one containing Ω (that will be our P ), and one (say, Q), contained in Ω. The polygon Q can then be decomposed as we need, while the quads in the ribbon will be decomposed in triangles. Clearly, the union of the decomposition of Q and the strip of triangles forms a decomposition of P .
Once P and its decomposition have been constructed, we can define the discrete space V h ( P ) as follows. On polygons belonging to Q we will use virtual elements of degree k, as we did in (18), while on the triangles in the ribbon we simply take usual finite elements of degree k. Notice that global continuity is ensured, since the degrees of freedom at the interelements are the same for Virtual and Finite elements. The unknowns of our problem will then be the elements of V h ( P ). The next step is to construct a suitable projection operator Π ∇ k in each element of the decomposition. For elements in Q which, according to the construction, will be just polygons as those that we considered before, we proceed as in the previous Section 2. The discrete bilinear form will be defined as in (26) (or in (27)) in the elements contained in Q, after computing the Π ∇ k -operator as in (5) (or Π 0 k−1 ∇v as in (8)). For the triangles in the ribbon we do not need any projection, since we already have polynomials. Hence, on the true element Ω T := Ω ∩ T Figure 4: Ω and the ribbon we can simply take Π ∇,Ω T k v = p k|Ω T . Then, the discrete bilinear form will be defined simply as (15) in each Ω T .

3d VEM for polyhedra with a curved face
The extension of the approach described in Subsection 4.2 to three-dimensional problems looks particularly hard. On the one hand, the number of generating points grows in a significant way with the degree k. On the other hand, most important, ensuring global continuity looks particularly complicated to realize. Instead, the new approaches indicated in Subsections 4.3 and 4.4 look more promising. Clearly, extensive numerical tests will need to be performed to validate them in terms of feasibility, efficiency, and accuracy. We briefly sketch here the extension of the two approaches. Here too, like we said in subsection 4.1, the full Dirichlet case presents no difficulties. We then concentrate on the treatment of the full Neumann case, or of the part of the boundary where Neumann conditions are imposed.

Using a subset of degrees of freedom.
The extension of the procedure highlighted in subsection 4.3 requires some care. Actually, on a polyhedron the number of degrees of freedom not positioned on ∂Ω is in general higher than the dimension of P k . Thus, a least-square solution would be necessary, but this is in contrast with the continuity requirement at the interfaces. In 2d the problem was circumvented easily just by fixing the values at the two endpoints of the curved edge. In 3d we should fix the value on the whole boundary of the curved face, which is what we are trying to define through the construction of the polynomial. Hence, in 3d we have only two possibilities, and not three like in 2d.
One possibility is to have a decomposition of tetrahedra-like around the boundary, that is, polyhedra with 3 flat faces and one, say σ, curved face on ∂Ω. We see that the values of a VEMfunction v on σ could be obtained through its dofs on the three flat faces plus the internal moments up to the order k − 3 only (and not k − 2). More precisely, for every integer k ≥ 1, using the values of v on the three straight edges (amounting to 3k + 1 dofs), the moments up to the degree k − 2 on the three flat faces (amounting to 3k(k − 1)/2 dofs), and the internal moments up to the degree k − 3 (amounting to k(k − 1)(k − 2)/6 dofs), we can identify uniquely a polynomial p * k ∈ P k . Indeed we have: (3k + 1) + 3k(k − 1) 2 + k(k − 1)(k − 2) 6 = k 3 + 6k 2 + 11k + 6 6 ≡ dim(P k ).
We can then compute a polynomial p * k ∈ P k (P ) through the conditions (It is not difficult to check that the dofs (51) are unisolvent for P k ). Once the polynomial p * k has been computed, we take v |σ = p * k|σ , and proceed as we did for polyhedra with flat faces. We point out that, by construction, the global continuity is guaranteed. We also observe that, contrary to what happens in 2D, the function v does not coincide with the polynomial, since for computing p * k we did not use all the internal moments of v, but just the moments up to order k − 3. To summarize, a function v is completely determined by the following conditions: The values at the vertices of P, for k ≥ 2 the moments e v p k−2 de ∀p k−2 ∈ P k−2 (e) ∀ straight edge e, for k ≥ 2 the moments F v p k−2 dF ∀p k−2 ∈ P k−2 (F) ∀ flat face F, for k ≥ 2 the moments P v p k−2 dP ∀p k−2 ∈ P k−2 (P ), Conditions (52)-(53) determine v uniquely, allowing to compute Π ∇,F k v on each flat face F, and then Π ∇,P k v as in Section 3 (see (34)-(36)), so that the discrete bilinear form (26) is computable. On a generic polyhedron, with more than three flat faces, the number of conditions (51) is bigger than the dimension of P k . In such a case we might use a least-square procedure to compute a polynomial p * k , and then use its moments on σ as degrees of freedom on the curved face. However, in order to guarantee the global continuity, before computing p * k we have to take care of the flat faces.
•1 st step: on each flat face F with one curved edge η we use least-squares to compute a polynomial p F k ∈ P k (F) as in (48), and set v |η = p F k|η , thus ensuring continuity of v on ∂σ. •2 nd step: we compute a polynomial p * k ∈ P k (P ) as the least-square solution of (51) •3 rd step: we set σ v ∇p k · n dσ = σ p * k ∇p k · n dσ, ∀p k ∈ P k (σ).
This procedure, together with conditions (52), determines uniquely v, and allows to compute Π ∇,F k v on each flat face F, and then Π ∇,P k v as in Section 3 (see (34)-(36)), so that the discrete bilinear form (26) is computable. Contrary to the case of tetrahedra-like elements, here we do not assign v on the curved face. Instead we impose degrees of freedom, more in the spirit of Mimetic Finite Differences (see always Remark 12).

Coupling FEM and VEM: the superimposed polyhedron
Let us sketch how to extend to the 3-D case what we did in Subsection 4.4, extending to the skin what we did in the two-dimensional case for the ribbon. For this, we assume that we are given a polyhedron P that contains strictly our three-dimensional domain Ω, and that P is (suitably) decomposed in polyhedra. More precisely, we will have a skin (that we manage to have made of tetrahedra), while the interior of P will be decomposed into polyhedra as needed. In particular here too we will assume that all the elements of the decomposition of P that are internal to P are also internal to Ω, as we had in the 2-dimensional case. Then we can proceed, mutatis mutandis, as we did in the two-dimensional case. Once P and its decomposition have been constructed, we can define the discrete space V h ( P ): we will use virtual elements of degree k, as we did in Section 3, in the internal polyhedra, while on the tetrahedra in the skin we simply take finite elements of degree k. The elements of V h ( P ) will be our unknowns. The discrete bilinear form will be defined as in (26) (or in (27)) in the elements contained in Ω, after computing the Π ∇ k -operator as in Section 3. In the tetrahedra, having already polynomials, we do not need any projection. Hence, like in the 2d case, in the elements Ω T := Ω ∩ T with one curved face we take Π ∇,Ω T k v h = p k|Ω T . Then, in each Ω T the discrete bilinear form will simply be (15) as for Finite Elements.