The Virtual Element Method in 50 lines of MATLAB

We present a 50-line MATLAB implementation of the lowest order virtual element method for the two-dimensional Poisson problem on general polygonal meshes. The matrix formulation of the method is discussed, along with the structure of the overall algorithm for computing with a virtual element method. The purpose of this software is primarily educational, to demonstrate how the key components of the method can be translated into code.


Introduction
The virtual element method, introduced in [6], is a generalisation of the standard conforming finite element method for the approximation of solutions to partial differential equations. The method is designed in such a way as to enable the construction of high order approximation spaces which may include an arbitrary degree of global regularity [9] on meshes consisting of very general polygonal (or polyhedral) elements. This cocktail of desirable features has attracted the method a lot of attention (see, for example, [1,10,4,13,11,3,18,26]), and is made possible through the virtual element space of trial and test functions, which is constructed on each mesh element from functions which are implicitly defined through local PDE problems. These local problems are designed in such a way that the virtual element space includes a subspace of polynomials of some prescribed degree (referred to as the degree of the method) alongside other, typically unknown, virtual functions. In this respect, and like many other approaches to polygonal meshes such as the Polygonal Finite Element Method [24,15] or BEM-based FEM [19], the virtual element method falls within the broad class of Generalised Finite Element Methods [22]. What sets the virtual element method apart from these other approaches, however, is that the extra nonpolynomial virtual functions never need to be determined or evaluated in practice. Instead, they are understood and used solely through certain defining properties of the virtual element space and through their degrees of freedom, which, along with the discrete bilinear form, are carefully selected to ensure that the final stiffness matrix can be directly and exactly computed.
There is a history of short, simple, codes being used to demonstrate the practical implementation details of various finite element methods. We refer, for instance, to the "Remarks around 50 lines of MATLAB" paper [2] which presented a simple and transparent MATLAB implementation of the conforming finite element method, the 99-line topology optimisation code presented in [20], or the mixed-FEM implementations presented by [5]. Unlike the large workhorse finite element libraries or commercially available 'black-box' finite element software packages, these codes are primarily designed to be educational, demonstrating how the theoretical concepts can be distilled into something practical. This is a tradition which we continue here, presenting a 50-line MAT-LAB implementation of the lowest order virtual element method for the Poisson problem in 2D on general polygonal meshes. To the best of our knowledge, this is the first publicly available implementation of the virtual element method, although there are various references containing details of the matrix formulation of the method, see for instance [7,13,14]. In particular, [7] contains detailed explanations of the formulation of the terms in the matrix equations for the high order virtual element method applied to a reaction-diffusion problem. In many ways, the work here can be seen as a spiritual successor to [7] in the sense that while we restrict ourselves to just the linear virtual element method for the Poisson problem, we take the process one step further and provide a clear, useable MATLAB implementation of the method.
The remainder of this work is structured as follows. In Section 2 we present the model problem of the Poisson problem. A very brief introduction to the virtual element framework is presented in Section 3, with a discussion of the discrete function spaces and bilinear forms. The details of the implementation of this method are given in Section 4, where we derive the matrix form of the discrete problem and show how this may be computed in practice. Section 5 contains a brief explanation of how to run the code using MATLAB. Finally, we offer some concluding remarks and ideas for possible extensions to the code in Section 6. The full code of the implementation is shown in Appendix A, and is available to download from: www2.le.ac.uk/departments/mathematics/research/virtual-element-methods-1/software alongside several examples of polygonal meshes.

Model problem
Let Ω ⊂ R 2 be a polygonal domain and consider the Poisson problem −∆u = f in Ω, with f ∈ L 2 (Ω) and g ∈ H 1/2 (∂Ω). This problem can be written in variational form as: find where (·, ·) denotes the standard L 2 (Ω) inner product. This variational problem possesses a unique solution by the Lax-Milgram lemma.

The Virtual Element Method
Let T h be a family of partitions of the domain Ω into non-overlapping polygonal elements with maximum diameter h whose boundaries are not self-intersecting. A theoretical analysis of the method requires certain additional assumptions on the regularity of the mesh, although since we do not present any analysis of the method we do not include these here. The resulting set of requirements is general enough, however, to include polygonal elements consisting of an arbitrary (but uniformly bounded) number of edges, which may also be non-convex. However, to simplify the implementation we restrict the mesh to include only elements which contain their own centroid, as defined in (4.4). We note that this class of elements includes those with co-planar edges, as commonly found in locally refined meshes with hanging nodes, and even non-convex elements.
For a polygonal element E with N E edges, we denote its vertices by ν i for i = 1, . . . , N E , and we adopt the convention that the edge e i connects ν i and ν i+1 , where the indices are understood to wrap within the range 1 to N E .
3.1. Virtual element function spaces. The discrete function space is defined to be h on the element E can be understood through the following three properties: • V E h includes the space P E of physical-frame polynomials on E with total degree ≤ 1 as a subspace.
• Any function in V E h can be uniquely identified by its values at the vertices of E, which are taken to be the degrees of freedom of the space. We note that this implies that the dimension of the space V E h is equal to N E . • Every function in V E h is a linear polynomial on each edge of E. The subspace of linear polynomials provides the approximation power of the virtual element space, and is responsible for the accuracy of the method. On triangular elements the space consists entirely of these linear polynomials, and thus the method reduces to the standard linear finite element method. However, on more general shaped polygonal elements the space will also include other, implicitly defined, 'virtual' functions, cf. (3.3). The method is designed in such a way that these will never need to be explicitly computed or evaluated, and are instead understood solely through their values at the vertices of E, which we take to be the degrees of freedom of the space In this respect, the virtual element space can be seen as a straightforward generalisation of the standard linear conforming finite element space on triangles to more general shaped elements.
The first observation we make about this space is that just the properties outlined above allow us to compute the Ritz projection Π E : V E h → P E of any function in the local virtual element space V E h onto the subspace of linear polynomials. This projection is defined for denotes the average value of w h at the vertices of E. This second condition is necessary to fix the constant part of Π E v h , and is clearly computable for any v h ∈ V E h from just its degrees of freedom.
From (3.1), the divergence theorem, and the fact that the Laplacian of a linear function is zero we have that, for any v h ∈ V E h and p ∈ P E , where n e denotes the unit normal vector to the edge e directed out of the element E. The final expression on the right hand side here can be exactly evaluated since v h is a linear polynomial on each edge of E, entirely determined by its values at the vertices, while the gradient of the linear polynomial p is a known constant. By picking a basis for the polynomial space P E , equation (3.2) can be written as a matrix problem which can be solved to find the coefficients of Π E v h with respect to this polynomial basis. We will come back to this in Section 4, although for now we just rely on the fact that this projection is computable.
The actual definition of the virtual element space which we use here is the lowest order space introduced in [6], given by , v| e ∈ P e for each e ∈ ∂E}, where P e denotes the space of linear polynomials on the edge e. The fact that the vertex values can be used as degrees of freedom to describe this space is proven in [6]. The global degrees of freedom for V h are simply taken to be the value of the function at each vertex in the mesh, thus imposing the continuity of the ambient space. The degrees of freedom at vertices on the domain boundary are fixed in accordance with the boundary condition. The dimension of the global virtual element space V h shall be denoted by N .
3.2. Discrete bilinear form. Define the bilinear form a E : H 1 (E) × H 1 (E) → R to be the restriction of a to the element E, i.e. a E (v, w) := (∇v, ∇w) 0,E for any v, w ∈ H 1 (E). Following [6], we introduce the discrete counterpart a E h : where dof r (v h ) denotes the value of the r th local degree of freedom defining v h in V E h with respect to some arbitrary (but fixed) ordering 1 . This means that S E is simply the Euclidean inner product between vectors of degrees of freedom. Finally, we define to be the discrete counterpart of a.
Crucial to the method is the observation that the local discrete bilinear forms satisfy the following two properties [6]: • Stability: there exists a constant C stab independent of h and E such that h . The requirement of polynomial consistency implies that the method satisfies the patch test commonly used in the engineering literature, expressing the fact that the method is exact when the solution is a piecewise linear polynomial with respect to the mesh T h , and provides the accuracy of the method. The stability property, on the other hand, ensures that the discrete bilinear form inherits the continuity and coercivity of the original variational form a, as proven in [6]. In the final matrix formulation of the problem, this property can be viewed as ensuring that the problem stiffness matrix is of the correct rank.
In light of these two properties, the two terms of a E h are referred to as the consistency and stabilising terms respectively since only the first term is non-zero when either v h or w h is a polynomial, and thus single-handedly ensures that the polynomial consistency property is satisfied, while the second term ensures that the stability property is satisfied even when v h or w h are in the kernel of Π E . For a proof that this choice of stabilising term S E does indeed satisfy the stability property, we refer to [13].
Moreover, both terms of a E h in (3.4) are computable using just the degrees of freedom of the virtual element space (to compute the projector Π E and to evaluate the stabilising term), and knowledge of the polynomial subspace P E (to evaluate the consistency term of a E h , which is made of integrals of polynomials, just like in a standard finite element method). This will be further demonstrated in Section 4, where it will also become apparent that this particular virtual element method can be implemented without requiring any quadrature to compute the stiffness matrix.
Still following [6], the linear form on the right-hand side of the variational problem The discrete problem which we solve can then be written as:

Implementation
As with a typical finite element method, we start by introducing the Lagrangian basis {ϕ i } N i=1 of V h with respect to the global set of degrees of freedom, defined by the property that ϕ i (ν j ) = δ ij , where δ ij is the Kronecker delta. We also introduce the Lagrangian basis of the local virtual , defined by the local equivalent of the same property.
We will also need a basis for the space P E of local physical frame linear polynomials on each element E. Many choices are possible here, although in keeping with the convention commonly adopted with virtual element methods, we choose the set of scaled monomials of degree 1. These are defined on the element E as where x E and y E respectively denote the x and y coordinates of the centroid of the element in the standard Cartesian coordinate system, and h E is the diameter of the element E. We denote by N P = 3 the cardinality of this basis and therefore the dimension of P E In the hope of avoiding confusion, we adopt the convention of indexing coefficients and basis functions in the basis of V E h using Latin letters, while those of P E will be indexed using Greek letters.
With these two bases at our disposal, we can now derive the matrix form of the discrete problem (3.6). Expanding the virtual element solution u h as for i, j = 1, . . . , N . Since both these terms are defined through sums over elements, the obvious way to compute the entries of K and F is by computing the non-zero local contributions from each element E in the form of the local stiffness matrix K E ∈ R N E ×N E and local forcing vector . . , N E , and adding them into the corresponding entries of K and F . This, of course, dictates that the overall structure of a virtual element method implementation will be much the same as for a standard finite element method, as outlined in Algorithm 1. The key point of departure from the standard finite element method is in how the element stiffness matrices should be calculated. Firstly, the computation of the local stiffness matrices relies on first computing the local Ritz projector Π E on each element. Secondly, where the implementation of a conventional finite element might rely on a mapping to a reference element, no such equivalent process is possible here because the mesh elements are allowed to be general (possibly non-convex) polygons. The full code of the implementation is given in Appendix A. In the remainder of this section we dissect the code to highlight how the various steps outlined in Algorithm 1 can be expressed in a form which can be implemented in code. Much of the matrix formulation presented in this section is similar to that in [7], although here we focus specifically on the case of the lowest order   Figure 1 are several examples of such meshes, distributed alongside the code. This information can also be generated in the same format using the Voronoi mesh generator PolyMesher [25], also written in MATLAB.

4.2.
Initialisation. The initialisation step of the code, shown in Listing 4.1, simply sets up various variables which will be useful to us later. In the interests of efficiency, we use a sparse matrix K to represent the stiffness matrix. In this step we also define the cell array linear polynomials containing three pairs of numbers indicating the degree of the associated polynomial in the x and y directions. Thus, the index of a specific polynomial in this array is taken to be the index of the polynomial in the basis M E . We note that the ordering imposed in the code coincides with the ordering in (4.1), although this choice is arbitrary. Some extra element-specific initialisation also takes place within the loop over all mesh element to compute various geometric properties of the element. The vector vert ids contains the global indices of the vertices forming the element E with an anti-clockwise ordering. As well as providing us with a means of looking up the coordinates of the vertices using the mesh structure, this also provides us with a very simple way to identify the global index of a particular local degree of freedom. This is possible because the indices of the degrees of freedom of the global virtual element space can be taken to be just the global indices of their associated vertices. Similarly, the local indices of the vertices of the element E dictate the local index of the associated degree of freedom. Therefore, the i th entry in the vector vert ids provides the global index of the i th local vertex and therefore also the global index of its associated local degree of freedom. Having access to this 'local to global' mapping is absolutely crucial when trying to assemble the local stiffness matrix and forcing vector into their global counterparts. When implementing more complex methods this sort of bookkeeping can quickly become very cumbersome, although here we are able to exploit the very simple arrangement of the degrees of freedom. Because of this property, we use the variable n sides to represent both the number of sides of E and equivalently the number of local degrees of freedom of V E h . The variable area denotes |E| and is computed using the formula where (x i , y i ) are the coordinates of the vertex ν i , and the indexing is understood to wrap within the range 1 to N E . The centroid (x E , y E ) of the element is stored in the vector centroid and calculated using the usual formula: where the indices are again to be understood to wrap within the range 1 to N E . In the code, we are able to combine some of the calculations which are necessary to find the area and centroid by storing the terms of the sum (4.3) in the vector area components.  An illustration of the labelling of the various geometric attributes on each element. Vertices are labelled as ν j , the edge connecting ν j to ν j+1 is denoted by e j , and e i denotes the lines segment connecting ν j−1 and ν j+1 . The outward unit normal in each case is denoted by n with the appropriate subscript.
for the projection Π E ϕ E i , either in the basis of V E h or in that of P E : Recalling (3.2) and using the fact that ∇m β is a constant vector for the linear polynomial m β , we find that for any m β ∈ M E , where in the last equality we have used the Lagrangian property of the basis functions ϕ E i at the vertices of E to determine that ϕ E i is only non-zero on the two edges e i and e i−1 which meet at the vertex ν i . As shown in [12], this can be further simplified to where we have denoted by e i the line segment connecting the vertices ν i−1 and ν i+1 , and n ei is the unit normal to e i such that n ei · n ej ≥ 0 for j = i, i − 1.
In view of this, we introduce the matrix G ∈ R N P ×N P and the vector B i ∈ R N P such that to encode the conditions above. The problem here is that the first row (and column) of G and B i are zero, since the gradient of a constant function is zero, and therefore G is rank deficient. This is overcome by using the second condition in the definition (3.1) of Π E , defining Therefore, the coefficients a i,α can be calculated by solving the (full rank) matrix equation Defining the matrix B ∈ R N P ×N E such that its i th column is the vector B i ∈ R N P , we obtain the matrix equation GΠ = B, where Π ∈ R N P ×N E is the matrix representation of the Ritz projector Π E , taking a vector of coefficients of a function expressed in terms of the basis of V E h to a vector of coefficients of the basis of P E , and has a i as its i th column.
We also introduce the matrix D ∈ R N E ×N P with D i,α := dof i (m α ) as a one-way 'change of basis' matrix for re-expressing polynomials in terms of the basis of V E h . Looking again at equation (3.2), it is apparent that we can use D as a shortcut to compute G and G, since (4.7) Finally, we can compute the local stiffness matrix as (4.8) where I ∈ R N E ×N E denotes the identity matrix. The first term of this sum corresponds to the consistency term of the discrete bilinear form, and the second term corresponds to the stabilising term. The code which computes the matrix form of the local Ritz projector and the local contribution of the stiffness matrix is given in Listing 4.2. The first task (lines 23-24 of Listing 4.2) is to initialise the two matrices D and B representing their namesakes D ∈ R N E ×N P and B ∈ R N P ×N E . For each of these matrices, we can immediately calculate the elements corresponding to the constant polynomial basis function m 1 . In the case of D, every element of the first column contains the value 1 since the constant function is 1 everywhere, while from (4.6) it may be observed that every element in the first row of B is equal to N −1 E . However, the remaining elements of D and B must be computed separately for each of the basis polynomials with total degree equal to 1, and for each local degree of freedom. Computing the entries of D is a straightforward task, since it just involves evaluating the basis monomials at the vertices of the polygon.
For the entries of B, however, we must evaluate the second expression in (4.6). The quantity | e i |n ei is simple to calculate, since it is just the vector | e i |n ei = (y i+1 − y i−1 , x i−1 − x i+1 ) due to the anti-clockwise orientation of the vertices around E. In the code, this result is stored in the variable vertex normal, so-called because it can be interpreted as a weighted normal vector at the vertex ν i . Again, the indices here are understood to wrap within the range 1 to N E , and in the code this is accomplished using the utility function mod wrap which modifies the standard mod function to produce output in the range 1 to N E rather than the range 0 to N E − 1. This modification to the mod function is necessary because arrays in MATLAB start at index 1, not 0.
To compute the entries of B we also need to evaluate ∇m β . Since m β is a linear polynomial, its gradient is simply a constant vector, and from the definition of M E it is clear that , and hence by representing the polynomial degree of m β in the x and y directions using a vector with one entry 1 and one entry 0, as in the cell array linear polynomials, the gradient can be very simply calculated by just dividing by h E .
With the matrices D and B computed, we are in a position to use (4.7) to calculate the matrix Π representing the projector Π E , as shown on line 37 of Listing 4.2. Consequently, with D, B and Π at our disposal, the local stiffness matrix can be computed as in (4.8).
The final step of this section is to add the elements of the local stiffness matrix to the positions in the global stiffness matrix for the corresponding global degrees of freedom. This is accomplished on line 41 through the local to global mapping discussed in Section 4.2.

4.4.
The local forcing vector. To calculate the local forcing vector given in (4.2), we must first compute the projection Π E 0 f . By definition, this satisfies which, because we are projecting to constants, can be simplified to where in the last relation we have used the barycentric quadrature rule on the polygon to approximate the integral. Since we are only considering the linear virtual element method, this is sufficiently accurate to ensure the optimal order of convergence in the H 1 (Ω) norm. It is the use of this quadrature which produces the requirement in Section 3 that the element must contain its own centroid. Clearly more general integration methods are possible (by triangulating the element, for example, or using a more advanced method such as [23,21,16,17], although for the sake of simplicity this is not something we pursue here. Since each basis function of V E h is defined to be 1 at a single vertex and 0 at the others, we can express The code to compute this and store the result in the appropriate positions in the global forcing vector is given in Listing 4.3. Applying the boundary conditions. The final step involves condensing the degrees of freedom associated with the boundary of the domain out of the linear system using the boundary condition, solving the resulting matrix equation, and re-applying the boundary data to the computed solution. This part of the procedure is exactly the same as for a standard finite element method, but for completeness we briefly review the process here.
Using the subscript B to denote the indices of degrees of freedom on ∂Ω, and I to denote those in the interior of Ω, the matrix problem KU = F can be expressed as where K IB = K BI by the symmetry of the bilinear form. Therefore, since U B is known, we must find U I by solving the problem  where r and θ are the standard polar coordinates centred at the origin, allows us to implement the standard example problem on an L-shaped domain. Then, running vem('meshes/L-domain.mat'); will produce the plot shown in Figure 3(b).

Conclusions and extensions
We have presented a 50-line MATLAB implementation of the linear virtual element method introduced in [6] for solving the Poisson problem on polygonal meshes in 2 spatial dimensions, available to download from: www2.le.ac.uk/departments/mathematics/research/virtual-element-methods-1/software alongside several example polygonal meshes. To the best of our knowledge, this is the first publicly available implementation of the virtual element method. It is clear from the literature surrounding the method that its capabilities extend far beyond what is presented here, although the intention behind this work is to exemplify how the method can be implemented in practice, in the simplest possible setting. The ideas we present here can, however, be generalised to much more  complicated situations by applying similar processes to compute the various required terms. The possible extensions of this code are endless: the implementation of higher order methods, more general elliptic operators including lower order terms and non-constant coefficients [13,8], basis functions with higher global regularity properties [9], mesh adaptation driven by a posteriori error indicators [11], or the consideration of time dependent problems [26] to name but a few.