Virtual Element Method for Quasilinear Elliptic Problems

We present a Virtual Element Method (VEM) for the solution of Dirichlet problems for the quasilinear equation $-\text{div} (k(u)\text{grad} u)=f$ with essential boundary conditions. Within the VEM the nonlinear coefficient is evaluated with the piecewise polynomial projection of the virtual element ansatz. Well posedness of the discrete problem and optimal order a priori error estimates in the $H^1$ and $L^2$ norms are proven. In addition, the convergence of fixed point iterations for the solution of the resulting nonlinear system is established. Numerical examples confirm the convergence analysis.


Introduction
In this work we present an arbitrary-order conforming Virtual Element Method (VEM) for the numerical treatment of quasilinear diffusion problems. Both two and three dimensional problems are considered and the method is analysed under the same mesh regularity assumption used in the linear setting [6,17], allowing for very general polygonal and polyhedral meshes.
Virtual element methods for linear elliptic problems are now well-established, see eg. [6,11,1,10,5,17,13] and [26] for a simple implementation. See also [7] for an extension to meshes with arbitrarily small edges and [16] where the mesh generality is exploited within an adaptive algorithm driven by rigorous a posteriori error estimates. While the VEM framework has been concurrently extended to a number of different problems and applications, the literature on VEM for nonlinear problems is scarce, the same being true for other approaches to polygonal and polyhedral meshes also. The Cahn-Hilliard problem is considered in [2], the stationary Navier-Stokes problem in [8], and inelastic problems in [4]. However, the first two problems are semilinear, while for the (quasilinear) latter problem no analysis is provided. The related nodal Mimetic Finite Difference method is analysed in [3] for elliptic quasilinear problems whereby the nonlinear coefficient depends on the gradient of the solution, however only low-order discretisations are considered. We also mention the arbitrary order Hybrid High-Order method on polygonal meshes for the general class of Leray-Lions elliptic equations [20], including the problems considered here. The HHO method belongs to the class of nonconforming/discontinuous discretisations and is, in fact, related to the Hybrid Mixed Mimetic approach and to the nonconforming VEM [23,19]. In [20], the convergence of HHO is proven under minimal regularity assumptions, but the rate of convergence of the method is not analysed.
The VEM presented here is based on the C 0 -conforming virtual element spaces of [1] whereby the local L 2 -projection of virtual element functions onto polynomials is available and the VEM proposed in [17] for the discretisation of linear elliptic problems with non-constant coefficients. In particular, to obtain a practical (computable) formulation, the nonlinear diffusion coefficient is evaluated with the element-wise polynomial projection of the virtual element ansatz. This results in nonlinear inconsistency errors which have to be additionally controlled.
We present an a priori analysis of the VEM which builds upon and extends the classical framework introduced by Douglas and Dupont [21] for standard conforming finite element methods. The analysis relies on the assumption that the nonlinear diffusion coefficient is bounded and Lipschitz continuous and is based on a bootstrapping argument: 1. existence of solutions for the numerical scheme is shown by a fixed point argument, 2. the H 1 -norm error is bounded by optimal order terms plus the L 2 -norm error, 3. using a standard duality argument and assuming that the discretisation parameter is small enough, the L 2 -norm error is bounded by optimal order terms plus potentially higher-order terms, 4. based on the existence result, L 2 -convergence is shown by a compactness argument, and now H 1 -convergence follows from step 2. Within this approach, we also obtain optimal order a priori error estimates in the H 1 -and L 2 -norms, albeit under the (higher) regularity assumptions needed by the duality argument. To the best of our knowledge, this work provides the first optimal order error estimate for a conforming discretisation of quasilinear problems on general polygonal and polyhedral meshes.
To simplify the presentation, we consider homogeneous Dirichlet boundary value problems only. To this end, we introduce the model quasilinear elliptic problem where Ω ⊂ R d is a convex polygonal or polyhedral domain for d = 2 or d = 3, respectively. The diffusion coefficient is a twice differentiable function κ : R → [κ * , κ * ] such that 0 < κ * ≤ κ * < +∞, and with bounded derivatives up to second order. Therefore κ is Lipschitz continuous, namely there exists a positive constant L such that for a.e t, s ∈ R.
The remainder of this work is structured as follows. We introduce the virtual element method in Section 2. The method is then analysed in Section 3, where the well-posedness and a priori analysis are presented. In Section 4 we establish the convergence of fixed point iterations for the solution of the nonlinear system resulting from the VEM discretisation. We present a numerical test in Section 5 and, finally, we provide some conclusions in Section 6.
We use standard notation for the relevant function spaces. For a Lipschitz domain ω ⊂ R d , d = 2, 3, we denote by |ω| its d-dimensional Hausdorff measure. Further, we denote by H s (ω) the Hilbert space of index s ≥ 0 of real-valued functions defined on ω, endowed with the seminorm | · | s,ω and norm · s,ω ; further (·, ·) ω stands for the standard L 2 -inner-product. The domain of definition will be omitted when this coincides with Ω, eg. | · | s := | · | s,Ω and so on. Finally, for ∈ N ∪ {0}, we denote by P (ω) the space of all polynomials of degree up to .

The Virtual Element Method
We introduce the virtual element method for the discretisation of problem (1.3), using general polygonal and polyhedral decompositions of Ω in two and three dimensions, respectively. We start by recalling the definition of the virtual element spaces from [1,17].  Let ω ⊂ R d , 1 ≤ d ≤ 3, be a d-dimensional polytope, that is, a line segment, polygon, or polyhedron, respectively. For any regular enough function v on ω, we define the following sets of degrees of freedom: where α is a multi-index with |α| := α 1 + · · · + α d and x α := x α1 1 . . . Let {T h } h be a sequence of decompositions of Ω into non-overlapping and not self-intersecting polygonal/polyhedral elements such that the diameter of any E ∈ T h is bounded by h.
On T h , we introduce element-wise projectors as follows. We denote by P h ≡ P ,E h : L 2 (E) → P (E), ∈ N, the standard L 2 (E)-orthogonal projection onto the polynomial space P (E). With slight abuse of notation, the symbol P h will also be used to denote the global operator obtained from the piecewise projections. Similarly, by P h ≡ P ,E h , ∈ N, we denote the orthogonal projection of (L 2 (E)) d onto the space P (E) = (P (E)) d , obtained by applying P ,E h componentwise. Further, we consider the projection R h ≡ R ,E h : H 1 (E) → P (E), for ∈ N, associating any v ∈ H 1 (E) with the element in P (E) such that with, in order to uniquely determine R h , the additional condition: Let k ≥ 1 be given, characterising the order of the method. We follow the construction of the corresponding C 0 -conforming VEM space presented in [1] to ensure that all of the above projectors, to be utilised in the definition of the method, are computable.
We first introduce the local spaces on each element E of T h , for d = 2. Let B 2 k (∂E) be the space defined on the boundary of E in the following way In [1] it is shown that the following degrees of freedom (DoF) uniquely determine the elements of V E h : with degrees of freedom given in agreement with the local degrees of freedom (2.3).
The construction of the space for d = 3 is similar, although now we define the boundary space to be where V f h is the two-dimensional conforming virtual element space of the same degree k on the face f . The local virtual element space is defined to be Finally, the global space and the set of global degrees of freedom for d = 3 are constructed from these in the obvious way, completely analogously to the case for d = 2.
The following are well established properties of the virtual element spaces introduced above [6, 1, 17]: • For each E ∈ T h , we have P k (E) ⊂ V E h as a subspace; ∇v are computable just by accessing the local DoFs of v given by (2.3) and (2.4) in the two and three dimensional case, respectively.
• The global virtual element space V h ⊂ H 1 0 (Ω) as a finite dimensional subspace.

2.2.
Virtual element method. The virtual element method of order k ≥ 1 for the discretisation of (1.1) reads: find u h ∈ V h such that where a h (·; ·, ·) is any bilinear form on V h defined as the sum of elementwise contributions a E h (·; ·, ·) satisfying the following assumption [6].
is bilinear and symmetric in its second and third arguments and satisfies the following properties: Remark 2.1. The above defining conditions are essentially those introduced in the linear setting [6,11,1,10,17] with, crucially, the nonlinear diffusion coefficient κ evaluated with the polynomial projection of the argument. We note also that the symmetry and stability assumptions imply the continuity in V h of the form a h (z; ·, ·), for z ∈ V h . Remark 2.2. The particular choice of local bilinear forms used in the numerical tests is given below in Section 5. We remark, however, that the following error analysis is valid whenever the assumption above is satisfied.

Error Analysis
We recall that k ≥ 1 is a fixed natural number representing the order of accuracy of the method (2.5).
The convergence and a priori error analysis of the VEM relies on the availability of the following best approximation results.
3.1. Approximation Properties. We recall the optimal approximation properties of the VEM space V h introduced above. These where established in a series of papers [6,1,16] under the following assumption on the regularity of the decomposition T h . Assumption 3.1. (Mesh Regularity). We assume the existence of a constant ρ > 0 such that • for every element E of T h and every edge/face e of E, h e ≥ ρh E • every element E of T h is star-shaped with respect to a ball of radius ρh E • for d = 3, every face e ∈ E h is star-shaped with respect to a ball of radius ρh e , where h e is the diameter of the edge/face e of E and h E is the diameter of E.
The above star-shapedness assumption can be relaxed by including elements which are union of star-shaped domains [6]. In particular, the following polynomial approximation result [14] is extended to more general shaped elements in [24] and the interpolation error bound below can be generalised by modifying the proof in [16], see also [25]. Theorem 3.2 (Approximation using polynomials). Suppose that Assumption 3.1 is satisfied and let s be a positive integer such that 1 ≤ s ≤ k + 1. Then, for any w ∈ H s (E) there exists a polynomial w π ∈ P k (E) such that Moreover, we have In the above bounds, C are positive constants depending only on k and on ρ.
The approximation properties of the virtual element space are characterised by the following interpolation error bound, whose proof can be found in [16]. Theorem 3.3 (Approximation using virtual element functions). Suppose that Assumption 3.1 is satisfied and let s be a positive integer such that 1 ≤ s ≤ k + 1. Then, for any w ∈ H s (Ω), there exists an element w I ∈ V h such that where C is a positive constant which depends only on k and ρ.
Then, using the fact that P k−1 h f is the L 2 projection on P k−1 (E), we can show the following lemma.

3.2.
Existence. We first show the existence of a solution u h of (2.5) using a fixed point argument.
To this end, for M > 0, we let Theorem 3.4. Let f ∈ L 2 (Ω) be given and assume that ( Proof. We devise a fixed point iteration for (2.5): for a fixed f ∈ L 2 (Ω), consider an iteration map It is easy to see that there exists h M > 0, such that for h < h M , T h v h is well defined, see for example [17]. For v h ∈ B M and w h = T h v h , in view of the stability assumption (2.7) and (3.3), we have Thus, choosing M sufficiently large, so that f ≤ M c , we get Therefore, the operator T h maps the ball v h ∈ B M into itself. By the Brouwer fixed point theorem, we know that T h has a fixed point, which implies that (2.5) has a solution u h ∈ B M .

Error bounds.
In our a priori error analysis, we follow a similar-in-spirit approach to the classical work of Douglas and Dupont [21] where standard conforming finite element methods were analysed in the same context.
We start with the following preliminary H 1 -norm error bound.
Proof. From Theorem 3.3, there exists a function u I ∈ V h , such that u − u I is bounded as desired. Thus, to show (3.6) it suffices to bound ∇(u h − u I ) . Let ψ = u h − u I , then using the stability Assumption 2.2 with c * = κ * α * , we have where u π is, on every element E ∈ T h , the polynomial approximation of u given by Theorem 3.2. Next, we will bound the various terms I i , i = 1, . . . , 5. We start with I 1 . Using Lemma 3.1, and the fact that r ≤ s, we have To bound I 2 , in view of (1.2), we get Also, using the fact that κ is bounded along with Theorem 3.2, we obtain (3.10) Using the fact that ∇u π ∈ P k−1 (E) and Assumption 2.2, we have thus, in view of the stability of P h , the fact that κ is Lipschitz continuous, u ∈ W 1 ∞ (Ω), Theorem 3.2 and the hypothesis κ(u) ∈ W r−1 ∞ (Ω), we deduce (3.11) Finally, we easily get (3.12) |I 5 | ≤ C( u − u π + u − u I ) ∇ψ ≤ Ch r u r ∇ψ .
Therefore, combining the above estimates (3.8)-(3.12) with (3.7) we obtain Then, in view of Theorem 3.2 and the stability of P h in L 2 -norm, we obtain the estimate Next, we shall demonstrate the following preliminary L 2 -norm, error bound.
Proof. We use a duality argument. Consider the (linear) auxiliary problem: find φ ∈ H 1 0 (Ω) such that −div(κ(u)∇φ) + κ u (u)∇u · ∇φ = u − u h . Noting that this equates to κ(u)∆φ = u − u h and given Ω is convex, we have φ ∈ H 2 (Ω) and (3.14) In variational form, the above problem reads In the sequel we will show Lemma 3.2, which in view of (3.14), gives For II in (3.16), using the Hölder inequality and the fact thatκ u ,κ uu are bounded uniformly on R, we get Next, in view of the Gagliardo-Nirenberg-Sobolev inequality, (3.21) v L3 ≤ C v 1/2 ∇v 1/2 , the Sobolev Imbedding Theorem and the elliptic regularity (3.14), we have Combining the previous estimates for terms I and II, we get the desired bound for h sufficiently small.
To complete the proof of Theorem 3.6, it remains to show that the consistency error bound (3.19) holds true. We do so through the following lemmas.
Proof. Let φ I ∈ V h be the approximation of φ given by Theorem 3.3 and using (1.3) and (2.5) we split the difference a(u; u, φ) − a(u h ; u h , φ) as Then, in view of (3.17), we rewrite term I as Employing Theorem 3.3 and (3.14), we obtain As for term II, using Lemma 3.1 we get In view of bounding term III, we write with φ 1 π | E ∈ P 1 (E) and u π | E ∈ P k (E), for any E ∈ T h given by Theorem 3.2. Using Theorems 3.2 and 3.3 we bound the term III 1 in (3.24) as Next, to estimate III 2 , we split this term as a summation over each E ∈ T h and use the polynomial consistency (2.6) and the definition ofκ u , given by (3.17), to get Then, following the steps for the estimation of I 4 in (3.11), using Theorems 3.2 and 3.3 along with (3.14), we can see that To bound III 2 2 , we first note, in view of (3.20), that (3.26) Further, using the stability property of P h , namely P h φ I L3(E) ≤C φ I L3(E) , and the Gagliardo-Nirenberg-Sobolev inequality (3.21), we obtain Using this in (3.26) and summing this new bound of (3.26) and (3.25) over all E ∈ T h and using Theorems 3.2 and 3.3, it follows that Finally, as a consequence of Lemma 3.3 below, we have Combining this with (3.23), the bounds for III 1 , and III 2 , the desired bound follows.
where φ 1 π ∈ P 1 (E) for all E ∈ T h , is given by Theorem 3.2, and r = min{s, k + 1}. Proof. Using polynomial consistency (2.6), the fact that P h ∇u π = ∇u π , with u π ∈ P k (E) given by Theorem 3.2 and the definition ofκ u given by (3.17) Let I = E I E , then we easily get . Using Theorem 3.2, we have ∇φ 1 π L6 ≤ C ∇φ W 1,6 and, hence, using a Sobolev imbedding, (3.28) ∇φ 1 π L6 ≤ C|φ| 2 . Now, using Theorem 3.2 once again, we get To bound II E , we rewrite this term as Then for II = E II E , using Theorem 3.2, it immediately follows that Next, we consider the term III E , which can be rewritten as Then using Hölder inequality (3.20), and we obtain for III 1 = E III E,1 . Hence, following the steps in the proof of Lemma 3.2, and using (3.28), we get Next, in view of the fact that ∇u π · ∇φ 1 π ∈ P k (E), we have Thus, for III 2 = E III E,2 , we get |III 2 | ≤ Ch u h − P h u h L3 ∇φ 1 π L6 ∇u π . Therefore, Theorem 3.2, and the Sobolev inequalities (3.21), (3.28), give Collecting the above bounds, yields for III = III 1 + III 2 Therefore from which the desired bound follows using once again Theorem 3.2.
Having concluded the proof of Theorem 3.6, in order to show optimal convergence rate of the error in H 1 and L 2 -norms, it remains to demonstrate that u h converge to u. Theorem 3.7. Under the same assumptions as in Theorems 3.5 and 3.6, the VEM solution u h converges to the exact solution u in H 1 0 (Ω). Proof. From Theorem 3.4 it follows that ∇u h is bounded from above. Therefore, we can choose a subsequence u h k such that for some z ∈ H 1 0 (Ω), u h k → z, weakly in H 1 0 (Ω), as h k → 0 and, thus, strongly in L 2 (Ω). Also, for arbitrary then z is the weak solution of (1.1). To show (3.31), we rewrite its left-hand side as Using the fact that u h k → z, and v h k → v, we see that (3.31) holds. Hence a(z; z, v) = (f, v), and thus u = z, since u is the unique solution of (1.1). Then, it follows that u h → u in L 2 (Ω). Hence, u − u h → 0 and the result follows from Theorems 3.6, and 3.5.
In view of Theorems 3.5, 3.6 and 3.7, the following a priori error estimates now readily follows.

Iteration method
In this section we show that, given a virtual element space V h , the sequence of solutions we obtain using fixed point iterations to solve the VEM problem (2.5) converges to the true solution u h ∈ V h of (2.5).
Starting with a given u 0 h ∈ V h we construct a sequence u n h , n ≥ 0, such that The convergence in H 1 of the sequence u n h as n → ∞ to a fixed point of (4.1), and hence a solution of (2.5), is an immediate consequence of the following result.
Proof. In view of Assumption 2.2 and the fact that a h (u n h ; ·, ·) is symmetric, we have with c = κ * α * . Then using (4.1), we obtain . Therefore, F(u n h ) is a decreasing sequence and, in view of the fact that ) → 0, as n → ∞, which completes the proof.

Numerical results
In order to test the VEM proposed in Section 2 we need to specify a bilinear form satisfying Assumption 2.2. We fix a E h as follows: with the VEM stabilising form S E given by here, I denotes the identity operator, − → v h is the vector with entries the degrees of freedom of v h ∈ V E h , and − → v h · − → w h is the euclidean scalar product of the degrees of freedom of v h , w h ∈ V E h . The above definition of the local bilinear form extends to the nonlinear setting the one considered in [17] and, similarly to the linear case, it is straightforward to show that it satisfies the stability condition (2.7). Following [6] instead, the projector R h can be used in place of P h in the stabilising term. The practical implementation of these projector operators and VEM assembly are discussed in [9,17].
In the examples below, approximation errors are measured by comparing the piecewise polynomial quantities P k h u h and P k−1 h ∇u h with the exact solution u and solution's gradient ∇u, respectively.
The tests are performed using the VEM implementation within the Distributed and Unified Numerics Environment (DUNE) library [12], presented in [15].
A representative example is shown in Figure 1. The polygonal mesh was generated using [27]. We use fixed point iterations analysed in Section 4 to solve the nonlinear system resulting from the VEM discretisation. This is compared below with Newton-Raphson iterations, defined as follows. Given an initial iterate u 0 h ∈ V h , we construct a sequence u n+1 h = u n h + δ n , n ≥ 0, by solving at each iteration the linearised problem: find δ n ∈ V h such that Here, the extra terms stemming from the linearisation of both the consistency and stability terms in a h are collected in the global form b h := E∈T h b E h , with the local form b E h , E ∈ T h , given by b E h (u n h ; δ n , v h ) = E κ u (P h u n h )P h δ k (P h ∇u n h ) · (P h ∇v h ) d x Numerical test. We consider the following test problem from [18]. We solve (1.1) on Ω = [0, 1] 2 with κ(u) = 1/(1 + u) 2 and the function f chosen such that the exact solution is u = (x − x 2 )(y − y 2 ). Note that, although the diffusion coefficient is not even bounded on the whole of R, it is smooth in a neighbourhood of the range of u. As initial guess for the nonlinear solve we use the constant zero function and the conjugate-gradient method is used to solve the linear system at each iteration. The relative errors for the approximation of u and its gradient as a function of the mesh size h are shown in Table 5 for k = 1 and a sequence of polygonal meshes, cf. the right-most plot in Figure 1. The numerical results confirm the theoretical rate of convergence. The table also displays the number of fixed point and Newton-Raphson iterations performed until the indicated stopping criteria is reached.  Table 1. Errors and empirical order of convergence (EOC) on a sequence of polygonal meshes. The Fixed Point and Newton-Raphson iterations needed to reach the tolerance 10 −10 are reported in the right-most columns.
The convergence history with respect to all meshes in Figure 1 are reported in the loglog plots of Figure 2 showing that the performance is similar in all cases. Note that, as k = 1, in the case of the sequence of triangular meshes, the VEM coincides with the standard linear finite element method.

Conlusions
With this paper, we show that the Virtual Element Method can be extended to nonlinear problems. In particular, we consider elliptic quasilinear problems with Lipschitz continuous diffusion in two and three dimensions and show that it suffices to evaluate the diffusion coefficient with the component of the VEM solution which is readily accessible. We prove optimal order a priori error estimates under the same mesh assumptions used in the linear setting.