Robust nonconforming virtual element methods for general fourth order problems with varying coeﬃcients

We present a class of nonconforming virtual element methods for general fourth order partial diﬀerential equations in two dimensions. We develop a generic approach for constructing the necessary projection operators and virtual element spaces. Optimal error estimates in the energy norm are provided for general linear fourth order problems with varying coeﬃcients. We also discuss fourth order perturbation problems and present a novel nonconforming scheme which is uniformly convergent with respect to the perturbation parameter without requiring an enlargement of the space. Numerical tests are carried out to verify the theoretical results. We conclude with a brief discussion on how our approach can easily be applied to nonlinear fourth order problems.


Introduction
In recent years the discretization of partial differential equations via the virtual element method (VEM) has seen a rapid increase. Introduced in [9] VEM began as an extension and generalization of both finite element and mimetic finite difference methods as discussed in [14]. In [9] the appropriate local and global VEM spaces are constructed and the approximation properties analysed for the Laplace equation. Another discretization of the Laplace problem was suggested in [1] while a nonconforming approach can be found in [7]. An extension to general, nonlinear second order elliptic PDEs for both conforming and nonconforming spaces is discussed in [21]. Similarly, another approach is taken in [11] for diffusion-convection-reaction problems.
The versatility of VEM has been showcased through the wide variety of problems it has been applied to over recent years. This has led to the construction of H(div) and H(curl)-conforming virtual element spaces in [12], conforming virtual elements for polyharmonic problems in [6], and the construction of methods for Stokes flow in [15,20,23], to name but a few. Especially the ease with which VEM spaces can be constructed to enforce desirable properties of the discrete functions even on general polygonal meshes makes the approach very interesting for a wide range of problems. An example of this is the construction of divergence free vector spaces in [15]. A further example is the construction of discrete spaces with higher order continuity conditions. The construction of even a lowest order C 1 conforming space is not straightforward within the standard finite element setting and higher order nonconforming spaces suitable for fourth order problems are also not readily available. Consequently, many software packages provide a large number of spaces for second order problems but often only provide the lowest order Morley element [32] for discretizing fourth order problems without requiring the use of splitting methods.
To construct conforming elements for fourth order problems, C 1 continuity is required which makes the methods highly complex. It is known that using traditional finite element methods, polynomials of at least degree five are needed to construct C 1 approximations which are piecewise polynomials. In contrast, it is shown in [19] that the virtual element construction of C 1 approximations to fourth order plate bending problems is much simpler and arguably more elegant. Additionally, the conforming virtual element method for polyharmonic problems, ∆ p u = f for p ≥ 1, has been addressed in [6] where the global VEM space consists of C p−1 functions. As well as this, the study of linear elliptic fourth order problems in three dimensions is considered in the conforming case in [13]. Another example of C 1 conforming elements can be seen in the application of the lowest order VEM space to the Cahn-Hilliard equation, investigated in [3]. Further studies of the application of virtual elements to the Cahn-Hilliard equation can be found in [30,31].
In this work we focus on studying nonconforming virtual element methods but for a wide range of problems. Although we focus on nonconforming VEM for fourth order problems, we highlight that due to the general framework we present, only minor modifications are needed to also include the study of C 1 conforming elements for these problems. Existing works which study nonconforming fourth order problems include the nonconforming approximation of the biharmonic plate bending problem, which is considered in [5,39,40]. A mixture of spaces have been suggested, some fully nonconforming [5,40] and others which include some level of continuity [39] though not the full C 1 continuity you would see in a fully conforming space. More recently, we see a C 0 conforming approach to fourth order perturbation problems being considered in [38]. There are other general approaches that can be employed such as hybrid high-order (HHO) methods [28] which, like VEM, can easily handle polyhedral meshes. The connection between VEM and HHO has been analysed and discussed for the Poisson problem in [29]. Hybrid high-order methods have also been extended to fourth order problems in [16] where a novel HHO method for the Kirchhoff-Love plate bending problem is presented. To our knowledge, the application of both higher order VEM and HHO methods to more general nonconstant coefficients and nonlinear fourth order problems is not available at the time of writing.
Arguably the most important ingredient of VEM is the construction of projection operators. In the available literature on fourth order problems, projection operators are constructed based on the underlying variational problem. The main idea of this approach is to construct only one projection which depends on the local contribution to the bilinear form. In [21,1] a different approach was taken for discretizing second order problems, which makes it straightforward to apply the method to nonlinear models. In this paper we generalize this approach and demonstrate how it can be applied to a wide range of fourth order problems. A major advantage of this approach is that it can be included more easily into existing software frameworks. A central building block for implementing Galerkin type schemes is the evaluation of nodal basis functions and their derivatives at given quadrature points.
To extend this to our VEM setting, these methods have to be replaced with the evaluation of projection operators defined on each element. We implemented this approach within the DUNE [8,26] software framework, requiring little change to the existing code base. From the user perspective switching between a finite element to a virtual element discretization is seamless especially within the available Python frontend [27,25].
We refer throughout this paper to two well known VEM spaces, the nonconforming space discussed in both [5,40] and the C 0 conforming space discussed in [39,38] by demonstrating how they fit into our generalized framework. The main contributions of this paper are detailed as follows.
(i) We present a general approach for constructing nonconforming VEMs for any order of accuracy, for solving a general fourth order PDE problem with nonconstant coefficients. We extend the work in [21] on second order elliptic problems and using the same techniques, we prove optimal error estimates in the energy norm (Section 5).
(ii) In Section 4, we introduce a general approach for constructing VEM projection operators based on constraint linear least square formulations that ensures they are fully computable from the available degrees of freedom.
(iii) The fourth order perturbation problem is studied in Section 6. We discuss how two standard spaces mentioned previously fit into our framework although note that our approach for constructing the projection operators differs from the approach taken in [5,40,38]. We furthermore present a novel fully nonconforming scheme which unlike modifications seen in [38], does not require an enlargement of the standard nonconforming VEM space.
(iv) Lastly, we carry out numerical experiments in Section 7 to confirm the a priori error analysis. We conclude by showing how our construction of projection operators allows us to straightforwardly solve complex nonlinear fourth order problems.

The continuous problem
Throughout this paper, we adopt the standard notation for Sobolev spaces H s (D) for nonnegative integers s, and for a domain D. We denote the norm and seminorm by · s,D and | · |s,D respectively. If D = Ω then the subscript shall be omitted. The notation (·, ·)D will be used to denote the L 2 (D) inner product. For a nonnegative integer l, let P l (D) denote the set of all polynomials up to degree l over D. We use the convention that P−1(D) = {0}. We denote the standard L 2 (D) orthogonal projection onto the polynomial space P l (D) by P l D . The tensor of all derivatives of a given order |µ| is denoted with D |µ| ϕ. Let ∂nϕ = ∇ϕ · n denote the normal derivative of a function ϕ over ∂D and let ∂sϕ = ∇ϕ · τ denote its tangential derivative where we use τ to denote a tangential vector.
Consider a general linear fourth order problem defined on a polygonal domain Ω ⊂ R 2 described by a bilinear form We make the minimal assumptions that κ, β, γ ∈ L ∞ (Ω). In later sections we impose further conditions on the coefficients. For now we assume that the coefficients satisfy κ ≥ κ0 > 0, for a constant κ0 and β, γ ≥ 0. We define β0 = min β(x), and γ0 = min γ(x). Note that we could also consider an even more general setting, e.g., take β ∈ L ∞ (Ω) 2×2 as in [21]. The results in this paper can be easily extended to cover this case but to keep the presentation simple we only consider scalar coefficients.
The variational problem for a given f ∈ L 2 (Ω) reads as follows: find u ∈ H 2 0 (Ω) such that Since the bilinear form is symmetric, we can define an energy norm |||·||| by |||v||| 2 = a(v, v). It follows easily that the bilinear form is coercive and continuous with respect to the energy norm, Hence it follows from the Lax-Milgram Lemma that (2.2) has a unique solution.
Note that if we were considering constant coefficients, taking κ, β, γ ∈ R, as in [13], then the strong form would reduce to the PDE studied there

The discrete problem and an abstract convergence result
In this section we provide some general ingredients needed for the discretization of our problem and present a Strang-type abstract error estimate. Let T h denote a tessellation of the computational domain Ω and denote the set of all edges in T h by E h . We split this set into boundary edges, E bdry , which again is made up of interior and boundary vertices.
For an integer s > 0, define the broken Sobolev space H s (T h ) by and on this space define the broken H s seminorm For a function v ∈ H 2 (T h ) we define the jump operator [·] across an edge e ∈ E h as follows. For an internal edge, For boundary edges, e ∈ E bdry h , let [v] := v|e. We denote with P k (K) the space of polynomials over a grid element K and define the piecewise polynomial space P k (T h ) for any k ∈ N with We now make the following basic assumptions. In particular, we stress that throughout the paper the polynomial order l is fixed.
Assumption 3.1. Assume the following holds for any fixed h > 0 and for a fixed l ≥ 2.
(A1) The mesh T h consists only of simple polygons. A simple polygon refers to the criteria that the boundary of each element must not intersect itself and is made up of a finite number of straight line segments.
(A3) There exists f h ∈ V h,l , which approximates the right hand side of our variational problem (2.2).
(A4) There exists a discrete bilinear form a h : V h,l × V h,l → R, such that for any u h , v h ∈ V h,l , The bilinear form a K h : V h,l |K × V h,l |K → R is the restriction of a h to an element K. We denote the restriction of the VEM space V h,l to an element K by V K h,l := V h,l |K .
(A5) Stability property: assume that there exists two constants α * , α * such that for all v h ∈ V K h,l . The criteria in the stability property (A5) is required to show that the discrete bilinear form is coercive and continuous.
Lemma 3.2. The broken Sobolev norm | · | 2,h is a norm on the spaces H 2 0 (Ω) and H 2,nc where the element wise contributions are given by Then, we have that |||·||| h is a norm on H 2,nc l (T h ). Therefore, under Assumption (A2), it follows that both | · | 2,h and |||·||| h are a norm on V h,l .
The following is now a direct consequence of the stability assumption (A5): Theorem 3.3 (Existence and uniqueness of solutions to the discrete problem). Under Assumption 3.1 the discrete problem: find u h ∈ V h,l such that admits a unique solution.
We now have the following Strang-type error bound, the proof of which is standard and identical to the method from e.g. [21].
where α * and α * are from the stability property (A5). The nonconformity error is given by We finish this section by collecting the remaining technicalities needed for the rest of the paper. In particular, we make the following regularity conditions on the mesh T h which are standard in the virtual element framework, see e.g., [9]. Assumption 3.5 (Mesh assumptions). Assume there exists some ρ > 0 such that the following hold.
(A6) For every element K ∈ T h and every edge e ⊂ ∂K, he ≥ ρhK where he = |e| and hK is the diameter of K.
(A7) Assume that each element is star shaped with respect to a ball of radius ρhK .
Finally, we recall some standard results for the L 2 projection operator. Definition 3.6. For any K ∈ T h define the L 2 (K) orthogonal projection onto the polynomial space P l (K), that is P l K : L 2 (K) → P l (K) by, and for any edge e ⊂ ∂K define the L 2 (e) orthogonal projection onto P l (e), P l e : L 2 (e) → P l (e) by, (P l e v, p)e = (v, p)e for all p ∈ P l (e). A proof of the following error estimates can be obtained using for example the theory in either [17,22].

The virtual element spaces
We dedicate this section to the virtual element discretization. We specify the chosen degrees of freedom by a dof tuple which allows us to easily encode a number of different local VEM spaces. A major part of the VEM method is the construction of projection operators. We detail a new construction of projection operators suitable for general VEM discretization of a wide range of fourth order problems with nonconstant coefficients. Throughout this section we provide examples of nonconforming VEM spaces for fourth order problems (C 1 nonconforming spaces). In particular, we use as examples the original nonconforming space detailed in [5,40] with V h,l ⊂ H 1 0 (Ω). As well as this space, we look at the C 0 conforming space detailed in [39,38], such that V h,l ⊂ H 1 0 (Ω). We conclude this section with defining the global VEM spaces and the bilinear forms.

Degrees of freedom tuple
We begin this section by introducing the concept of a degrees of freedom (dof) tuple, used to generically describe the degrees of freedom relating to a VEM space on each element of the grid. So let K ∈ T h be fixed in the following.
Definition 4.1. For the C 1 virtual element spaces, let the degrees of freedom tuple, M ∈ Z 5 , be defined as The entries correspond to the number of moments used in the definition of our degrees of freedom, with d v j , for j = 0, 1, encoding the information for the vertex moments, d e j for j = 0, 1, for the edge moments, and d i 0 for the inner moments. The subscript j = 0 corresponds to moments of the function values and j = 1 to derivative moments on the vertices and edges.
From the dof tuple, we are able to define the corresponding degrees of freedom. (D3) The moments of v h up to order d i Remark 4.3. We use the convention that D 0 v h = v h and ∂ 0 n v h = v h . If any of the entries in the dof tuple M are zero, this implies that we take constant moments, if an entry is less than or equal to −1, then this corresponds to no moments. So for example d v 0 = −1 implies that no vertex values are used. The case d v 0 = −1 is for example relevant for the C 0 nonconforming space presented in [21] but will not be considered in the discussion here.
Note that for the C 1 nonconforming spaces, we always have d v We only begin to prescribe the vertex derivative values when we consider C 1 conforming spaces or even higher order conforming spaces see for example [6,13].
We now give some examples of degrees of freedom tuples relating to some common VEM spaces to illustrate the idea.
. This dof tuple also corresponds to the space considered in [40]. A visualization of these dofs on triangles can be seen in Figure 1.
For l ≥ 2, consider the C 1 space introduced in [39,38] which is C 0 conforming. This space is described by the dof tuple . The dofs for this space are shown in Figure 2.
For l ≥ 3, consider the C 1 conforming space given in [6]. The dof tuple describing the local dofs for this space is given by : Degrees of freedom for polynomial orders l = 2, 3, 4, 5 on triangles for the C 1 -C 0 conforming space [38]. Circles at vertices represent vertex dofs, arrows represent edge normal dofs, circles on edges represent edge value moments and interior squares represent inner dofs.
Remark 4.5. As mentioned already, the concept of defining a dof tuple extends to describing the dofs relating to C 0 VEM spaces. Consider the simplest of the C 0 conforming serendipity spaces discussed in [10]. The dofs for this space can be described by the dof tuple

Local spaces and projection operators
We now focus on the crucial aspect of defining the VEM spaces and projection operators. Following the approach in [1,21] we introduce an enlarged local virtual element space to ensure that our projection operators satisfy certain L 2 projection properties which are stated at the end of this section.
Definition 4.6. The enlarged virtual element space V K h,l on K ∈ T h is defined as follows.
where n K e denotes the number of edges in the polygon K.
The enlarged space V K h,l is characterized by the following extended degrees of freedom tuple, denoted by Note that the number of extended dofs is equal to N dof s = n K e (2l − 1) + 1 2 (l + 1)(l + 2) so that we have N K V = N dof s . We denote with Λ K the set of extended degrees of freedom described by M K as given by Definition 4.2. We show next that this set of degrees of freedom is unisolvent in V K h,l .
Lemma 4.7. The set of extended degrees of freedom Λ K is unisolvent over the space V K h,l .
Proof. We show that if all the degrees of freedom vanish for v h ∈ V K h,l then v h ≡ 0. Using Green's formula [22], for a function v h ∈ V K h,l it follows that By using integration by parts on the last term and setting all degrees of freedom Λ K to zero, we see that Applying integration by parts on the second term we arrive at the following Since v h ∈ V K h,l , it follows from (4.2) that ∆ 2 v h ∈ P l (K) and ∆v h ∈ P l−2 (e). We therefore see that |v h | 2 2,K = 0. As in [5], this implies that v h = 0 due to the boundary conditions. Remark 4.8. In this paper we focus on the construction for C 1 nonconforming spaces. However, note that the following discussion on defining projection operators also covers some conforming spaces suggested in the literature for solving second order problems, e.g., the spaces from [10]. As well as this, as discussed above a conforming C 1 space can also be described in this framework (using the dof tuple (0, 0, l − 4, l − 3, l − 4)). The nonconforming C 0 space from [21] also fits the framework using the dof tuple (−1, −1, l − 1, −1, l − 2). Note that both of these spaces require a different enlarged space. Since we are not going to analyse these two spaces, we do not discuss these choices further but we would like to note that our definition of projection operators cover these cases as well with only minimum modifications.
We now turn our attention to the local VEM space V K h,l , which we define as a subspace of the enlarged virtual element space V K h,l . In order to define V K h,l , we first construct the following projections: an interior value projection, , and an edge normal projection Π e 1 : V K h,l → P l−1 (e). These projections have to be computable from the degrees of freedom Λ we can then define the VEM space V K h,l , the gradient projection, Π K 1 : V K h,l → (P l−1 (K)) 2 , and finally the hessian projection Π K 2 : V K h,l → (P l−2 (K)) 2×2 , which satisfy certain L 2 projection properties discussed in the following. These projections are then used to define the discrete bilinear form.
To begin with, we assume that we have a degrees of freedom tuple Definition 4.9. We say that a value projection Π K 0 , an edge projection Π e 0 , and an edge normal projection Π e 1 are dof compatible if for v h ∈ V K h,l they are a linear combination of the dofs Λ K M (v h ), and satisfy the following additional properties. and and Π e 0 q = q|e for all q ∈ P l (K).
and Π e 1 q = ∂nq|e for all q ∈ P l (K).
Note that there are multiple choices for defining the value, edge, and edge normal projections such that they are dof compatible. We provide an example for defining these projections based on constraint least squares problems after defining the gradient and hessian projections.
and the hessian projection Here n, τ denote the unit normal and tangent vectors of e, respectively.
One possible dof compatible choice for the value and two edge projections is shown in the following example.
Example 4.11. We consider projection operators obtained from constraint least squares problems. Consider the dof tuple (0, −1, d e 0 , d e 1 , d i 0 ). • We define the value projection Π K 0 v h ∈ P l (K) as the solution to the problem Minimize: subject to: Since Π K 0 is defined by a linear least squares problem with equality constraints, where the right hand side is given by dofs, it follows that Π K 0 is a linear combination of the dofs. From this definition it is clear that (4.3) holds.
• If we choose the edge projection to be the unique solution in P l (e) of then Π e 0 v h ∈ P l (e) and (4.4) is satisfied. To see this, we note that Π K 0 is linear and the other equations only involve linear systems of equations.
• Finally, if we define the normal edge projection Π e 1 v h ∈ P l−1 (e) to be the unique solution of then we satisfy (4.5). Note that we could replace the final constraint with . This is what we use in our implementation. Note that these definitions also cover the case of the C 0 nonconforming space with d e 1 = −1. The gradient projection in this case is identical to the one given in [21] but we get a projection for the hessian as well.
We can finally use given dof compatible projections Π K 0 , Π e 0 , and Π e 1 to define the local virtual element space on K.
Definition 4.12. The local virtual element space V K h,l is given as the following subset of the enlarged space.
We now show that the subset of dofs Λ K M are unisolvent for our local VEM space V K h,l .
Proof. Let v h ∈ V K h,l , and set the dofs Λ K M to zero. Then, noting that both the value and edge projection are computed using the dofs, it follows that and for j = 0, 1, i.e. the inner moments of order l, . . . , l − (d i 0 + 1) and the edge moments of order l − 2, . . . , l − (d e 0 + 1) of v h are also zero. Similarly, the higher order edge normal moments are also zero. So viewing v h as a function in the enlarged space, v h ∈ V K h,l ⊂ V K h,l it follows that v h ≡ 0 as a function in V K h,l with the extended degrees of freedom set to zero, using Lemma 4.7. Therefore the dofs are unisolvent for V K h,l .
This leads us to a final crucial result of this subsection. We see in the next Lemma that by construction, our value, gradient, and hessian projections possess important L 2 projection properties. In particular, due to construction of the VEM space, the value projection is identical to the L 2 projection into the space of polynomials of degree l for all VEM functions. Lemma 4.14. If the value, edge, and edge normal projections are dof compatible, then for any v h ∈ V K h,l it follows that, for s = 0, 1, 2, Proof. From Definition 4.9 it is clear that Combining the two we obtain the stated L 2 projection property of Π K 0 . Similar arguments are used for the other projections: due to the polynomial exactness of Π e 0 and the definition of the extended space V h,l we have Π e 0 v h = v h |e and therefore for a polynomial p ∈ P l−1 (K), where we used the L 2 projection property of Π K 0 . Applying integration by parts to the RHS completes the proof. Finally we have for the hessian projection using the results already proven and the dof compatibility of the edge normal projection, Π e 1 , so that the result follows using integration by parts.
Due to the dof compatibility, and the fact that P l (K) ⊂ V K h,l , it holds that P l (K) ⊂ V K h,l . This implies that Π K s p = D s p for p ∈ P l (K) for each s = 0, 1, 2.

Global spaces and the discrete bilinear form
We conclude with the definition of the global virtual element space and the discrete bilinear form. We keep the presentation brief, since it follows the general construction of VEM found in the literature. We define the global space V h,l by We set the local degrees of freedom which correspond to boundary vertices and boundary edges, e ∈ E bdry h , to zero. Note that the global degrees of freedom are unisolvent -this follows from the unisolvency of the local degrees of freedom and the definition of the local spaces.
Lastly, we define the discrete bilinear form.
Definition 4.17. For any u h , v h ∈ V K h,l , define the local discrete bilinear form a K h as We take the stabilization term S K (·, ·) to be a symmetric, positive definite bilinear form satisfying for constants c * , c * independent of h and K.
It is clear that the choice of a K h in Definition 4.17 with S K satisfying (4.10), will satisfy the stability property from (A5) see e.g., [9,21]. We take the standard choice for S K based on the scalar product of the local degrees of freedom scaled appropriately with some element wise constant approximations of the coefficientsκK ,βK , andγK . In our implementation we use the average of each coefficient κ, β, and γ evaluated at the vertices of an element K to calculateκK ,βK , andγK . Then S K is given by: Then, as in [21], it follows that this choice of S K satisfies (4.10) and therefore (A5) holds in Assumption 3.1.
Remark 4.18. Furthermore, notice that due to Lemma 4.14, it holds that Π K 0 p = p for all p ∈ P l (K) and therefore the stabilization term vanishes when one of u h or v h is a polynomial of degree l. The discrete bilinear form therefore possesses the standard polynomial consistency property.

Error analysis
In this section we collect all the building blocks needed to prove a general convergence result in the energy norm, which is presented in Theorem 5.7. Throughout this section we now assume that our dof tuple is H 2,nc l (T h ) conforming (recall Definition 4.15) and so (A2) is satisfied in Assumption 3.1, meaning that Theorem 3.4 holds. Proving convergence involves bounding all of the terms in the Strang-type estimate presented in Theorem 3.4. In the following we consider each term from (3.2) in the order in which they appear, starting with the approximation error.

Interpolation error estimate
The following estimate for the interpolation error, the first term in equation (3.2), is simply a consequence of standard scaling arguments.
for a constant C independent of h.

Load term
The next term to appear in the a priori bound (3.2) is an error estimate for the load term. To treat this term, following the methods in [5,21], we define f h to be the piecewise L 2 projection of f on T h , f h := P l K f . Using Lemma 4.14 we observe that for w h ∈ V K h,l , Consequently, the right hand side of the discrete variational problem (3.1) is computable. The following estimate now follows easily using for example, the method in [5].
Lemma 5.2. For l ≥ 2, let s be an integer with f ∈ H s (Ω) and define r := min(s − 1, l). Then, the following estimate holds for a constant C independent of h.
Proof. We can show the following using the bounds of Theorem 3.7

Nonconformity error
We now turn our attention to the next term in the a priori bound (3.2), the nonconformity error N (u, w h ). From now on we assume that l ≥ 2. Recall that our dof tuple is H 2,nc l (T h ) conforming and as per Definition 4.15, the values in our degrees of freedom tuple satisfy the following Recall that as a consequence of this for w h ∈ V h,l , and for any edge e in the grid This follows since the jump is zero at the vertices on the edges and ∂sq ∈ P l−2−1 (e) = P l−3 (e) which means that we can apply (5.2).
Assume that the solution u to (2.2) satisfies u ∈ H 4 (Ω). Then, for w h ∈ V h,l the nonconformity error N (u, w h ) satisfies

(5.5)
Proof. Using integration by parts, we can express the hessian and gradient terms as follows Since we assume u ∈ H 4 (Ω) we use (2.3) and an application of integration by parts to see that Therefore the nonconformity error is equal to We also use the following identities for rewriting the boundary terms in a way that is useful later on in Theorem 5.5. The next corollary looks at how (5.5) simplifies when our VEM space is C 0 conforming. We now bound each term in (5.5) to achieve an error estimate for the nonconformity error. This essentially involves the jump properties of the VEM space (5.2) -(5.4) as well as standard interpolation estimates detailed in Theorem 3.7.
Theorem 5.5 (Nonconformity error bound). Let (A1) -( A7) hold. Assume that the solution u to (2.2) satisfies u ∈ H l+1 (Ω) and assume that the coefficients satisfy κ ∈ W l−1,∞ (Ω) and β ∈ W l−2,∞ (Ω). Then, for l ≥ 3 the nonconformity error satisfies the following estimate for any w h ∈ V h,l , for a constant C independent of h. For l = 2, the nonconformity error satisfies for a constant C independent of h.
Using arguments as in the proof of Lemma 5.3 it follows that E2 = I1 + I2. These terms can be bounded as in (5.14) and (5.15) respectively. Therefore, it follows that Note that when the VEM space is conforming, V h,l ⊂ H 1 0 (Ω), it follows that w h is a viable test function in H 1 0 (Ω) and therefore we can take v = w h in (5.17). Hence, the nonconformity error reduces to the following To treat E1 in the case that V h,l ⊂ H 1 0 (Ω) we now introduce the interpolation of w h into the lowest order conforming VEM space, Π 1 w h ∈ H 1 0 (Ω). Then, it holds that Using standard estimates, we see that and finally, This concludes the proof.

Energy norm estimate
Now that we have successfully bounded the nonconformity error, we look at the final term in Theorem 3.4 and prove convergence in the energy norm.
To ease the presentation of the next Theorems and Corollary, we observe that the following inequality holds . Notice that the inequality in (5.18) holds due to standard discrete Poincaré inequalities, and for example can be shown similarly as in [5].
Theorem 5.6. Assume that (A1) -( A7) hold, defined in Assumptions 3.1 and 3.5. Let u ∈ H l+1 (Ω) be the solution to continuous problem (2.2). Assume that the coefficients satisfy κ ∈ W l−1,∞ (Ω), β ∈ W l,∞ (Ω), and γ ∈ W l+1,∞ (Ω). Then it holds that for a constant C independent of h. We define the remaining constants c0, c1, and c2 as follows; let c2 = κ Proof. To show that (5.19) holds, we note that Note that the first term is bounded easily by standard interpolation estimates. For the other term, we recall Remark 4.18, which implies that the stabilization part of the discrete bilinear form vanishes. Therefore, Then, denoting the coefficients as α2 := κ, α1 := β and α0 := γ, we can show that To see that (5.20) holds, we use the results of Lemma 4.14 to express our projections as L 2 projections. Hence, it follows that The result now follows.
We now have the following convergence theorem which is a result of Theorems 5.1, 5.5, 5.6 and Lemma 5.2.

Perturbation problem
We turn our attention to the following fourth order perturbation problem. For a polygonal domain Ω ⊂ R 2 the perturbation problem reads as follows We make the minimal assumptions that f ∈ L 2 (Ω) and ∈ R such that 0 < ≤ 1. Taking κ(x) = 2 , β(x) = 1, and γ(x) = 0 we can examine the error analysis from the previous section, with the energy norm now becoming |||v||| 2 h = 2 |v| 2 2,h + |v| 2 1,h . It is well known that for example the lowest order C 1 nonconforming space on triangles (the Morley element, [32]) does not lead to a scheme that is robust with respect to → 0 (see for example [35,33]). There have been a range of modifications suggested to the original Morley element, for example in [34,37,38]. In [33] a modification is suggested which on triangles corresponds to our C 1 -C 0 conforming space in the lowest order setting and convergence of the method is proven in this case. We give error estimates for the higher order version of those two spaces in the following. In addition we study a new modified C 1 nonconforming discretization, C 1 -mod, which has the same degrees of freedom as the original C 1 nonconforming space but is stable with respect to the perturbation paramter . This is achieved by a modification to the gradient projection, Π K 1 , given next.
Definition 6.1. We define the modified gradient projection to be the following for any v h ∈ V K h,l . We denote with Π l−1 the interpolation into a H 1 conforming VEM space of order l − 1. We use the "lazy" version of the serendipity spaces discussed in [10]. The dofs for this space were mentioned in Remark 4.5 and for order l − 1 are described by the dof tuple (0, −1, l − 3, −1, l − 4). Therefore the dofs are a subset of the dofs defining the nonconforming C 1 space. The vertex values and the l − 3 moments on the edges uniquely define Π l−1 v|e ∈ P l−1 (e) so that the gradient projection given above is computable using the dofs for the C 1 nonconforming space.
The following property is now obtained for the modified gradient projection. Lemma 6.2. For the modified gradient projection detailed in Definition 6.1 it holds that Our discrete bilinear form a h (·, ·) is now chosen as in Definition 4.17 with coefficients κ(x) = 2 , β(x) = 1, and γ(x) = 0 but using the modified gradient projection operator as detailed in Definition 6.1. So for any u h , v h ∈ V K h,l , define the local discrete bilinear form a K h as where we choose the same stabilization term S K presented previously.
To prove convergence in the energy norm for this modified scheme we use the ideas seen in [37,36]. We consider a modified bilinear form b(·, ·) by changing the lower order contribution.
for all w, v ∈ H 2 0 (Ω). Using b(·, ·) we can prove a a Strang-type lemma similar to Theorem 3.4. Theorem 6.3 (Abstract a priori error bound). Let l ≥ 1 be an integer. Let u be the solution to problem (2.2) with coefficients κ(x) = 2 , β(x) = 1, and γ(x) = 0. Suppose that u h ∈ V h,l is the solution to (3.1) using the modified local bilinear form (6.3). Under Assumptions 3.1, it holds that As in the previous section we can bound each term in the above abstract estimate to obtain an error estimate in powers of h and ε. Theorem 6.4. Assume that (A1) -( A7) hold, defined in Assumptions 3.1 and 3.5. Let l ≥ 2 and suppose that u ∈ H l+1 (Ω) is the solution to the continuous problem (2.2) with coefficients κ(x) = 2 , β(x) = 1, and γ(x) = 0. Suppose that u h ∈ V h,l is the solution to the discrete problem (3.1) using the modified gradient projection (Definition 6.1) in the discrete bilinear form. Assume that f ∈ H l−1 (Ω). Then, under these assumptions it follows that for a constant C independent of h.
Proof. The proof follows similar arguments used in Section 5 and therefore we keep the presentation here brief.
It is straightforward to show the following two bounds We focus on the remaining two terms involving the modified bilinear form b(·, ·).
For the nonconformity type error term, it holds that We show this by multiplying the strong equation (6.1) by Π l−1 w h ∈ H 1 0 (Ω) and integrating over Ω Notice that J2 can be bounded as before in the nonconformity error proof (Theorem 5.5). Then for J1 using Cauchy-Schwarz and standard estimates, Now, for J3 we bound this term using the optimal interpolation properties of the lower order VEM space as well as stability of the interpolation operator.
Finally, for the last term in (6.4) we can show that This follows from Due to the L 2 properties in Lemma 4.14 the hessian terms cancel. Secondly, using the property in Lemma 6.2 it holds that Π K Therefore, the result follows from combining all of the intermediary results.
The next corollary details all error estimates for the following spaces applied to the perturbation problem (6.1). We analyse the fully nonconforming space, C 1 -nc, with dof tuple (0, −1, l − 3, l − 2, l − 4), the C 1 -C 0 conforming space with dof tuple (0, −1, l − 2, l − 2, l − 4) (see Example 4.4), and finally the new modified scheme C 1 -mod with dof tuple (0, −1, l − 3, l − 2, l − 4). For completeness, we repeat the result obtained in Theorem 6.4 for the modified scheme C 1 -mod. Recall that this space has the same dof tuple as the C 1 -nc space however we define the discrete bilinear form for this problem using the modified gradient projection (Definition 6.1). Due to keeping track of the coefficients in the previous error analysis section, the following corollary is simply a consequence of Theorem 5.7 and Theorem 6.4. Corollary 6.5. Under the same assumptions as Theorem 5.7 with f ∈ H l−2 (Ω) the approximate solution in the C 1 nonconforming space C 1 -nc of order l satisfies The solution in the C 1 -C 0 conforming space, assuming that f ∈ H l−1 (Ω), satisfies the following error estimate Finally, for the modified C 1 nonconforming scheme defined by the modified gradient projection in Definition 6.1 and assuming that f ∈ H l−2 (Ω), it holds that In all three cases, the constant C remains independent of h and .
Remark 6.6. Note that for = O(1) all methods converge in the energy norm with order l − 1 while for = O(h) the convergence rate of the original C 1 -nc method reduces to O(h l−2 ) while the modified scheme retains the order h l−1 . The H 1 conformity of the C 1 -C 0 space leads to an improvement of the order to O(h l ) but requires more degrees of freedom.

Numerical results
In this section we show some numerical results to verify the a priori bounds from the previous sections for three virtual element spaces.
• C 1 -nc: the nonconforming space from [5,40]  using the modified gradient projection (see Section 6). • C 1 -C 0 : the continuous space from [39,38] defined by the dof tuple (0, −1, l − 2, l − 2, l − 4) (see Example 4.4) and using the default projection operators (see Example 4.11). We focus on results for l = 2, 3, 4 on both structured triangular grids and to demonstrate the flexibility of the virtual element method, on grids consisting of Voronoi cells. Figure 3 contains the images of the first few mesh refinements for this choice of grid. Grid data for all considered spaces and mesh refinements is shown in Table 1.
We first study convergence of the methods for a linear problem of the general form (2.2) with varying coefficients. As a second example we study the simple perturbation problem (6.1) with constant coefficients, and we show results with varying h, . For both examples, we solve the problem on the domain Ω = (0, 1) × (0, 1). The approximation errors are measured in the relative energy norm.
Finally we highlight how our approach based on general projection operators allows us to handle nonlinear problems by showing results for the Cahn-Hilliard equation and the Willmore flow of graphs.
The code used to perform the simulations is based on the DUNE software framework [8]. We implemented our VEM approach within the module DUNE-FEM [26]. This is an extension module for DUNE that provides interfaces for the implementation of general grid based numerical schemes on general unstructured grids. It is open source software implemented in C++. In our implementation we construct quadrature rules by subdividing polygons into triangles and applying a quadrature of sufficient order on each triangle. The required order of the quadrature is the same as required by standard finite elements. The consistency error due to quadrature was studied in [21] for the second order method and similar arguments can be applied here. Other possibly more efficient methods for constructing quadrature rules are available for general polygons, for example, [4]. Recently a Python based frontend was added to DUNE [27,25]. The domain specific language UFL [2] can be used to describe mathematical models. A detailed tutorial including some VEM examples, e.g., for linear elasticity and also the Willmore flow example described here, showcases the flexibility of the approach [25].

Linear varying coefficient problem
We start by studying the more general linear fourth order problem with varying coefficients.
We choose the forcing so that u(x, y) = (sin(2πx) sin(2πy)) 2 is the exact solution on the domain Ω =(0, 1) 2 . We show results on both the structured simplex grid and the Voronoi cells for polynomial orders l = 2, 3, 4 in Table 2.
All methods converge with the expected order of l − 1 and produce very similar errors on a given grid with the C 1 -C 0 space requiring more degrees of freedom.  Table 1: Grid data and total number of degrees of freedom of the sequence of simplex meshes (top) and Voronoi meshes (below) for all three spaces C 1 -nc, C 1 -mod, and C 1 -C 0 . The second and third columns correspond to the number of grid elements and the mesh size, h, respectively. The final columns show the total number of degrees of freedom for the three VEM spaces, for each order l = 2, 3, and 4. Notice that the C 1 -nc and C 1 -mod spaces have the same dof tuple and hence share the same number of dofs for all polynomial orders and grids.
Simplex meshes grid data grid size mesh size l = 2 l = 3 l = 4

Perturbation problem
We now study the perturbation problem considered in Section 6 for various values of ∈ (0, 1] and mesh size h. The aim is to verify the convergence orders discussed in Section 6 especially the improved order of the new scheme C 1 -mod compared to the original C 1 -nc space. We are only considering problems without boundary layers, a discussion of problem (6.1) with boundary layers can be found in [38].
Example 7.2. Consider problem (6.1) with right hand side given by f = 2 ∆ 2 u − ∆u. As exact solution we use u(x, y) = (sin(2πx) sin(2πy)) 2 which satisfies the Dirichlet boundary conditions on the domain Ω = (0, 1) 2 . The results on both the structured triangular grid and the Voronoi grid for the values = 10 −2 , = 10 −8 and for polynomial orders l = 2, 3 are summarized in Table 3, the results for l = 4 are in line with expectations.

Nonlinear problems
We conclude this section with some preliminary results, which demonstrate that the method discussed in this paper is well suited to solve complex nonlinear fourth order problems. We choose two problems which have been studied in the literature, both energy minimization problems solved by a gradient descent algorithm. In both cases the mathematical models are time dependent fourth order problems. We use a Rothe approach in which first the problem is discretized in time. The resulting spatial problems are stationary fourth order problems with linearizations of the form studied here.
Note that no special linearization is required, we use a standard Newton solver to handle the fourth order nonlinear problems arising from the implicit time discretization. Please note that this final subsection is only a first investigation into applying our approach to more complicated settings and it is beyond the scope of this paper to do a detailed analysis. We therefore restrict the presentation to the C 1 -nc space with l = 3. We show results on triangular grids on Ω = (0, 1) 2 . Example 7.3. For our first problem we solve the Cahn-Hilliard equation using a fully implicit backward Euler method to discretize the problem in time. A virtual element method for this problem was studied in [3] where a different approach for defining the projection operators, requiring a special linearization, restricted the method to l = 2.
Let ψ : R → R be defined as ψ(x) = (1−x 2 ) 2 Table 3: Example 7.2 with = 10 −2 (l = 2, 3 correspond to first two rows of tables) and = 10 −8 (l = 2, 3 correspond to last two rows of tables). Relative energy errors and eocs on simplex grids (left column) and Voronoi grids (right column) are shown. The observed convergence rates of l − 1 for h large are in accordance with our convergence results summarized in Corollary 6.5 and confirm results from Example 7.1. Note that the errors on a given grid are larger for the standard C 1 -nc method compared to the other two due to lower convergence on the coarser grids (see two top tables). In the case of small shown in the lower two rows of tables, the convergence of the three methods differs considerably as expected from Corollary 6.5. The C 1 -nc method converges with order l − 2, the C 1 -mod methods converges with order l − 1 while the C 1 -C 0 method converges with order l but requires more degrees of freedom (see Table 1).  Note that this problem requires slightly different boundary conditions compared to the problems studied so far. Some snapshots from a simulation with γ = 0.02 on a 60 × 60 grid and time step τ = 10 −3 are displayed in Figure  4. The initial conditions were a small perturbation in the centre of the domain. The first snapshot is taken at a point in time where the phase separation is already well developed. The following figures then show the usual coarsening ending with the red phase concentrated in approximately a circle in the centre of the domain. Example 7.4. In our final example we study the minimization of the Willmore energy of a surface, in the case where the surface is given by a graph over a flat domain Ω. The corresponding Euler-Lagrange equations can be rewritten as a fourth order problem for a function u defined over Ω. We use a second order two stage implicit Runge-Kutta method as suggested in [24]. The resulting problem is a system of two nonlinear fourth order partial differential equations. The corresponding linearized problem is a system of equations for the two Runge-Kutta stages which is of the form (2.1).
The time stepping scheme [24] produces an approximation u n h which can be constructed in the following way: given u n h and time step τn, solve the following coupled system for the two functions u n,1 h , u n,2 h ,

Conclusion
We have presented a general approach for constructing nonconforming VEM projection operators suitable for solving a wide range of fourth order problems. We have analysed the resulting methods and shown that the virtual element approximation in each of the considered spaces converges to the true solution in the energy norm with optimal convergence rates. We also introduced a novel modified nonconforming scheme for solving the fourth order perturbation problem which remained convergent as → 0. Unlike modifications seen in the literature, our change did not require an enlargement of the space and was obtained via an adjustment to the gradient projection. These results were verified with numerical experiments on a variety of polygonal meshes. We introduced a new concept of describing the degrees of freedom via a dof tuple which allowed us to easily encode a variety of different VEM spaces. We followed the approach taken in [1,21] and introduced dof compatible projections which were constructed in a way that made them computable entirely from the available degrees of freedom. In particular, we gave examples of how dof compatible projections could be constructed based on constraint least squares problems. We were then able to construct the VEM spaces in a way that ensured that the value, gradient, and hessian projections were all identical to L 2 projections.
The ease with which our approach can extend to the application of nonlinear fourth order problems was also demonstrated with some additional numerical experiments. However, note that the theory behind the numerical results displayed is out of the scope of this paper and is future work.