In this paper, we study the discretization by an H 2 -conforming virtual element method (VEM) of the simplest elliptic equations in nondivergence form
where Ω ⊂ R d is a bounded convex polytopal domain, A ∈ [L ∞ (Ω)] d×d is a field of uniformly elliptic symmetric positive definite matrices, f ∈ L 2 (Ω), and g ∈ H 2 (Ω).
Here, ∇ 2 denotes the Hessian and the colon stands for the Frobenius inner product between matrices.
The value function of a stochastic differential equation satisfies a linear PDE such as (1.1), typically with low-order terms. 18,21 They also appear in the linearization of fully nonlinear equations, such as Monge-Ampère and Hamilton-Jacobi-Bellman, which are relevant to a number of applications, including differential geometry, optimal transport and stochastic control, fluid mechanics and meteorology, image processing. 24,28,29 The structure of (1.1) is deceivingly simple. For example, (1.1) with forcing f = 0 and discontinuous coefficient A given by
x |x| admits two solutions in the unit ball B 1 (0) centered at 0, namely u(x) = |x| α -1 and u(x) = 0, which happen to be of class H 2 (Ω) provided d > 2(2 -α) for any 0 < α < 1. For 0 ≤ µ < 1, the (pointwise) Cordes condition
gives sufficient conditions on A, possibly discontinuous, to guarantee a unique solution u ∈ H 2 (Ω) of (1.1) for Ω convex and f ∈ L 2 (Ω); this is a consequence of a perturbation theory for the Laplacian. 16,36 We refer to Ref. 39 for a brief review on the different notions of solutions and corresponding numerical methods discussed in the literature.
Finite difference methods have been proposed for (1.1) upon discretizing the Hessian ∇ 2 u directly and assuming continuity of A. These methods are known to require wide stencils to enforce monotonicity and consistency, the required size of those stencils depending on the ellipticity constant of A. 35,37 We refer to the semi-Lagrangian methods for linear and nonlinear elliptic problems by Debrabant and Jakobsen 25 and by Camilli and Falcone, 15 which are two-scale methods. We also refer to the methods by Bonnans et al. 10 for d = 2 and Bonnans et al. 9 for d > 2, which feature more compact stencils, to the extent permitted by the ellipticity of A, at the cost of requiring a Cartesian structure of the discretization grid.
Conforming VEM for Nondivergence Form Linear Elliptic Equations 77
On the other hand, a two-scale Galerkin method has been proposed by Nochetto and Zhang. 39 All these methods are monotone and are known to converge in the max-norm.
An essential difficulty for Galerkin discretizations of (1.1) is the lack of a natural variational formulation associated to this problem. If the matrix field A is differentiable, we clearly have that problem (1.1) is equivalent to the divergence form problem
which leads to the usual variational formulation. However, this approach is not applicable when either A is not differentiable or A is differentiable but rough, thereby giving rise to an advection-dominated diffusion problem. Moreover, in typical applications arising from fully nonlinear PDEs, the coefficients cannot be expected to be continuous everywhere and the location of discontinuities cannot be assumed to be known a priori. Thus, the discretization of (1.1) without any further assumption on the regularity of the coefficients is of great importance. An approach introduced in Ref. 41 for solving numerically problem (1.1), and that is applicable when the matrix field A is rough and satisfies the Cordes condition, consists in discretizing the following variational formulation: find u ∈ V g such that Ω (γA :
where
and where γ ∈ L ∞ (Ω; R * + ) is a suitable scaling function. Since the range of the Laplace operator ∆ : H 2 (Ω) → L 2 (Ω) is the whole of L 2 (Ω), because of the (crucial) convexity assumption on Ω, and the boundary condition is just of Dirichlet type for u, the variational formulation (1.2) is equivalent to the strong form second-order equation (1.1).
It can be shown that (1.2) is well-posed in V in the setting of the Cordes conditions. The proof hinges upon the Miranda-Talenti estimate; see Ref. 41 and Theorem 2.1. Hence, the variational formulation (1.2) allows for the direct finite element discretization of strong solutions via H 2 -conforming elements. 30,41 For instance, the C 1 -conforming Bogner-Fox-Schmit (BFS) is shown by Gallistl in Ref. 30 to produce second-order accurate solutions over rectangular partitions for d = 2. But, in general, the need for H 2 -conforming elements, which are notoriously cumbersome to construct and implement, makes the straightforward discretization of (1.2) challenging or even impractical when d ≥ 3.
Earlier, Smears and Süli introduced in Ref. 41 the alternative idea of relaxing the conformity requirements by considering nonconforming discrete variational formulations including stabilization terms. Their formulation is designed so that stability can be shown based on establishing a discrete analogue of the Miranda-Talenti estimate whose proof is highly technical. Later approaches essentially follow Ref. 41, proposing nonconforming stabilized discretizations relying on discrete Miranda-Talenti estimates; see e.g. Refs. 27, 30, 31, 33, 38, and 42. In contrast, in this paper, we propose and study the discretization of (1.2) using an H 2 -conforming VEM. Originally, the VEM was proposed in Ref. 6 as a generalization of the H 1 -conforming FEM to polygonal/polyhedral meshes. However, it was soon realized that the flexibility of its design could be exploited to construct entirely new discrete spaces satisfying specific properties, such as a given type of conformity. An H 2 -conforming VEM was already proposed by Brezzi and Marini in Ref. 14 for the solution of plate bending problems and has since been generalized by several authors including for the solution of polyharmonic problems. 2-5, 7, 8, 13, 19, 26 Most recently, based on a hierarchical construction in the space dimension, a complete family of virtual element spaces satisfying any degree of conformity in arbitrary dimensions has been presented by Chen et al. in Ref. 20.
In terms of both their design and implementation, VEMs with higher degree of conformity are conceptually similar to the basic H 1 -conforming VEM. This stands in stark contrast with their finite elements counterparts, which are characterized by ad hoc constructions leading to spaces of relatively large dimensions that are not affine equivalent and, hence, cumbersome to implement resorting on a reference element as standard. 22,34 For this reason, VEMs, in which elements are also constructed on the physical domain directly, without the use of reference-to-physical mappings, are competitive in terms of implementation difficulty and cost against finite elements with high degrees of conformity, while also providing some benefits such as a low number of degrees of freedom and a consistent construction of the relevant virtual element spaces regardless of the dimension and the required polynomial consistency order.
Virtual elements are of generalized finite element type in that they include both polynomial and nonpolynomial functions. For instance, the H 2 -conforming virtual element space of (optimal) order m -1 in the H 2 -seminorm considered herein contains the space P m of polynomials of degree up to m but no polynomial of higher degree. The space is completed by nonpolynomial functions in view of the conformity requirements. The dimension of the resulting element is inherently consistent with the required level of conformity and approximation power. 20 As far as we are aware, this is the first application of the VEM to elliptic problems in nondivergence form. Our approach yields a rather compact sparsity pattern, compared with two-scale methods, but at the expense of monotonicity. Therefore, our notions of stability and convergence are not expressed in L ∞ -norm but rather in the H 2 -norm.
We show that the H 2 -conforming VEM discretization of the prototype elliptic equation in nondivergence form (1.2) is well-posed by directly exploiting the continuous Miranda-Talenti estimate. Furthermore, we prove stability and optimal rates Conforming VEM for Nondivergence Form Linear Elliptic Equations 79 of convergence in the H 2 -norm, including the effect of quadrature. The method and analysis apply to both standard and polygonal/polyhedral meshes, under a common shape-regularity assumption, namely uniformly star-shapedness of both the elements and their interfaces. This assumption reduces to classical shape-regularity in the case of standard meshes made of simplices.
Crucially, the analysis holds under no additional regularity on the data (A, f, g), hence in exactly the same basic setting for existence of (1.1). This is in contrast to error estimates for VEMs applied to elliptic problems in divergence form, which typically involve regularity of the coefficients, see e.g. Ref. 17, as is also customary for FEMs. 22 The reason is that we are able to exploit the strong form (1.1) in the analysis; see Ref. 33 for a similar argument in the setting of discontinuous Galerkin methods. In view of showing optimal convergence under minimal regularity, we also construct and analyze a Scott-Zhang 40 type interpolation operator for the H 2conforming VEM spaces of Ref. 20. Furthermore, the analysis is extended to include the effect of quadrature under some minimal continuity assumptions required for pointwise evaluation of the data. These results are verified in practice with a set of numerical examples carefully chosen to examine the effect of the regularity of coefficients and solution on the convergence rate. For instance, they confirm that optimal convergence rates are preserved when the coefficients are rough but the solution is smooth.
H 2 -conforming VEMs are conceptually simple and greatly simplify the analysis upon circumventing the need of a discrete Miranda-Talenti estimate. In addition, they may have other potential advantages. For instance, the fact that (high-order) approximations of the gradient are readily available from the numerical solutions is of particular interest in optimal control applications. In contrast, state-of-the-art approaches based on finite difference methods typically provide low-order gradient approximations. 9,10,15,25 We further highlight that the degrees of freedom of the lowest order H 2 -conforming VEM, namely the value and gradient at the mesh vertices, are especially simple. This is true in all dimensions, making the extension of the method to moderately higher-dimensional problems plausible. Finally, the fact that we deal directly with the strong form (1.1) entails immediate access to residual error estimates, as in Ref. 30. Such estimates may be used to drive automatic mesh adaptivity, thereby exploiting the mesh flexibility associated with VEMs.
The rest of this paper is organized as follows. In Sec. 2, we discuss the Cordes condition and show unique solvability of (1.1). In Sec. 3, we introduce an H 2conforming VEM framework in arbitrary dimensions for which we show existence of a unique discrete solution, and prove a quasi-optimal error estimate in the H 2 (Ω)norm. In Sec. 4, we prove a variant of the error estimate taking into account the effect of numerical quadrature on the scheme. In Sec. 5 we recall the construction of H 2 -conforming VEMs satisfying the framework of our analysis. For simplicity, we limit ourselves to the case d = 2, 3 and analyze the corresponding quasiinterpolation operator. We use the latter to derive optimal convergence rates in Sec. 6. We finally conclude in Sec. 7 with a set of numerical experiments exploring the interplay of regularity of solution and data with the accuracy of the VEM.
In this section, we review some elements of the analysis of H 2 solutions to Eq. (1.1). We start by defining the Cordes condition in Sec. 2.1, and then in Sec. 2.2 we discuss the variational formulation (1.2) and its well-posedness. The proof of the well-posedness of (1.2) under the Cordes condition is standard; 36,41 we choose to reproduce it for completeness, and in order to make it easy to compare with the analysis of our discrete variational formulation in Sec. 3.2.
We equip R d×d with the Frobenius norm that we denote by | • |.
The matrix field A ∈ [L ∞ (Ω)] d×d is said to satisfy the µ-Cordes condition if A(x) satisfies the µ-Cordes condition for almost every x ∈ Ω.
Remark 2.1. The Cordes condition is often written as |A| 2 /(tr A) 2 ≤ 1/(d-1+ε) for some 0 < ε ≤ 1. This is equivalent to the definition above, with µ = √ 1 -ε. Proof. For any γ ∈ R, one has
In particular, if γ = tr A/|A| 2 , then |γA -
2 , from which the equivalence between (i) and (iii) follows. It is obvious that (iii) implies (ii), and to prove that (ii) implies (iii) it suffices to prove that tr A/|A| 2 = argmin γ>0 |γA -I d | 2 , which follows from using (2.1) and then writing the first-order optimality condition for the minimization problem. Conforming VEM for Nondivergence Form Linear Elliptic Equations 81
The above characterization of the Cordes condition motivates the following definition: Definition 2.2. (Admissible scaling) Assume that the matrix field A ∈ [L ∞ (Ω)] d×d satisfies the µ-Cordes condition for some 0 ≤ µ < 1. The function γ ∈ L ∞ (Ω; R * + ) is said to be a µ-admissible scaling of A if |γ(x)A(x) -I d | ≤ µ for almost every x ∈ Ω. By Proposition 2.1, if A satisfies the µ-Cordes condition for some 0 ≤ µ < 1, then there exists a µ-admissible scaling of A, and a possible choice is given by the formula γ(x) = tr A(x)/|A(x)| 2 (note that in this case the fact that γ ∈ L ∞ (Ω; R * + ) follows from the assumption that A is uniformly elliptic and bounded).
The following proposition will be useful in order to simplify our estimates in the following sections. Proposition 2.2. Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1 and that γ is a µ-admissible scaling of A. Then, for any measurable set
Proof. For almost every x ∈ K,
which concludes the proof.
Let us define the bilinear form a : H 2 (Ω) × H 2 (Ω) → R and the linear form
so that the variational problem (1.2) can be written as follows: find u ∈ V g such that
We emphasize that the variational formulation (2.3) is consistent with problem (1.1).
Proposition 2.3. (Consistency) For any u ∈ V g , the following are equivalent :
(i) A : ∇ 2 u = f almost everywhere in Ω.
(ii) γA : ∇ 2 u = γf almost everywhere in Ω.
(iii) u is solution to (2.3).
Proof. The equivalence between (i) and (ii) follows from the fact that γ > 0 almost everywhere in Ω. Since γA : ∇ 2 u ∈ L 2 (Ω) and γf ∈ L 2 (Ω), (ii) is satisfied if and only if
Using that ∆ is a bijection from V 0 to L 2 (Ω), we deduce that (ii) is equivalent to (iii).
Let us now explain how one can use the Lax-Milgram lemma in order to prove the well-posedness of (2.3), provided that the Cordes condition is satisfied and that the scaling function γ is chosen appropriately. We start by stating the continuity of a and L.
with hidden constants only depending on d.
Proof. This is immediately verified.
The proof of the coercivity of a relies on the Miranda-Talenti estimate, which we recall below.
Theorem 2.1. (Miranda-Talenti estimate) For any u ∈ V 0 , one has ∇ 2 u 0 ≤ ∆u 0 .
Proof. See Refs. 36 and 41. Note that the proof uses the facts that Ω is convex and that u| ∂Ω = 0. Proposition 2.5. (Coercivity) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1 and that γ is a µ-admissible scaling of A. Then, for any u ∈ V 0 , one has a(u, u)
Proof. One has
where we used that γ is a µ-admissible scaling of A (see Definition 2.2) for the second inequality and the Miranda-Talenti estimate (Theorem 2.1) for the last inequality. We conclude by invoking once again the Miranda-Talenti estimate.
Remark 2.2. (Poincaré-type inequality) The H 2 seminorm | • | 2 is a norm on V 0 , and for any u ∈ V 0 , one has u 2 |u| 2 , with a hidden constant depending only on d and Ω. One way to prove this inequality is using the a priori estimate for the Poisson equation, according to which u 2 ∆u 0 , and then the easily verified
Corollary 2.1. (Well-posedness) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1 and that γ is a µ-admissible scaling of A. Then there exists a unique solution u ∈ V g to (2.3).
Proof. Let u g ∈ V g . Then u is solution (2.3) if and only if u = u g + u 0 , where u 0 is solution to the following: find u 0 ∈ V 0 such that
The existence of a unique solution u 0 follows from Propositions 2.4 and 2.5 and the Lax-Milgram lemma.
In this section, we introduce and discuss the main properties of an H 2 -conforming VEM based on the variational formulation (2.3). After defining the scheme in Sec. 3.1, we prove its well-posedness in Sec. 3.2 and derive an error estimate in Sec. 3.3.
Let us first introduce the setting of the discretization. We assume that the following parameters are given as follows:
• A partition T h of Ω in finitely many nonoverlapping polytopes.
• A polynomial consistency order m ≥ 2 for the numerical scheme.
Moreover, we assume that T h satisfies the following standard shape-regularity properties:
• There exists ρ > 0 such that any K ∈ T h and any d -dimensional facet F of K, with 1 ≤ d ≤ d, are star-shaped with respect to respectively a d-dimensional ball of radius ρh K and a d -dimensional ball of radius ρh F . • There exists η > 0 such that any d -dimensional facet F of any
Here h ω denotes the diameter of ω ⊂ Ω. We further define h := max K∈T h h K . From now on, for any two quantities a, b ∈ R, we write a b if a ≤ Cb with some constant C > 0 depending only on d, Ω, ρ, η, and m.
For k ∈ N and K ∈ T h , we denote by P k (K) the space of polynomials of degree k on K, and by P k (T h ) the set of functions p : Ω → R whose restrictions p| K to each K ∈ T h belong to P k (K). We denote by Π 0 k : L 2 (Ω) → P k (T h ) (respectively L 2 (K) → P k (K)) the L 2 projection operator onto P k (T h ) (respectively onto P k (K), when applied to a function defined only on some K ∈ T h ). We naturally extend the definition of the projection operator Π 0 k to matrix-valued functions. It is also useful to define as follows, for any K ⊂ T h , the local bilinear and linear forms a K :
There are several variants of the construction of H 2 -conforming VEMs. For completeness, we describe a possible construction in Sec. 5. For now, we abstract over some aspects of the construction by assuming that, on any cell K ∈ T h , we are given the following (remember that m ≥ 2 denotes an arbitrary polynomial consistency order):
with a constant c π ≥ 1 not depending on u and K.
with constants 0 < c * ≤ 1 ≤ c * independent of u, v, and K.
(Note that inequalities c π ≥ 1, c * ≤ 1, and c * ≥ 1 are not restrictive; the reason we assume them is to make the expression of our error estimates simpler.) We shall show that the above virtual elements ansatz is sufficient to complete the analysis, postponing to Sec. 5 a complete description of a possible realization of the above framework.
Following (2.2) and the usual approach for constructing VEMs, we define, for any K ∈ T h , the bilinear and linear forms a h,K :
Conforming VEM for Nondivergence Form Linear Elliptic Equations 85
We define global virtual element spaces
and
where g I ∈ V h is a given function. In order for V g I h to be a good approximation of the space V g , the boundary values of g I should be chosen as some interpolation of those of g; an example of a suitable interpolation operator is described in Sec. 5.5.
We define the global counterparts a h :
where the projection operator Π * m :
The virtual element scheme that we study is the following: find
In order to show the well-posedness of scheme (3.5), we adapt to the discrete setting the arguments used in Sec. 2.2 in the case of the continuous problem (2.3): Corollary 3.1, Proposition 3.2, and Corollary 3.2 are counterparts in the discrete setting to, respectively, Proposition 2.4, Proposition 2.5, and Corollary 2.1.
Proof. This is immediately verified.
where c * is from (3.2).
Proof. This follows from Proposition 3.1 and the Cauchy-Schwarz inequality.
Proposition 3.2. (Coercivity) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1 and that γ is a µ-admissible scaling of A. Then for any u ∈
, where c * is from (3.2).
Proof. Recall that
One has
, where we used that γ is a µ-admissible scaling of A (see Definition 2.2) for the third inequality and the Miranda-Talenti estimate (Theorem 2.1) for the last inequality. On the other hand,
, where for the second inequality we used that ∆Π * m u ∈ P m-2 (K) and that Π 0 m-2 ∆u is the L 2 projection of ∆u onto P m-2 (K). Summing up the above, one has
We conclude using the Miranda-Talenti estimate.
Note that the above proves the coercivity of a h on V 0 h × V 0 h only when c * > µ. In that case, we can apply the Lax-Milgram lemma to establish the well-posedness of the VEM.
Conforming VEM for Nondivergence Form Linear Elliptic Equations 87 Corollary 3.2. (Well-posedness of the scheme) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1, that γ is a µ-admissible scaling of A, and that c * > µ in (3.2). Then there exists a unique solution u h ∈ V g I h to (3.5).
The existence of a unique solution u 0 h follows from Corollary 3.1, Proposition 3.2, and the Lax-Milgram lemma.
Remark 3.1. Assumption (3.2) on the scaling of s h,K is standard in the virtual element setting. However, while in the setting of divergence form elliptic equations it is usually sufficient to assume that c * > 0, here we require c * > µ. This requirement is not restrictive. Admissible definitions of the stabilization form s h,K can always include some scaling factor, which can always be chosen large enough so that c * > µ. We refer to Sec. 5.4 for an example with more discussion.
We now prove our main error estimate for scheme (3.5).
Theorem 3.1. (Error estimate) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1, that γ is a µ-admissible scaling of A, and that c * > µ in (3.2). Let u ∈ V g and u h ∈ V g I h be respectively the unique solutions to (2.3) and (3.5). Then, for any
Proof. Using the triangle inequality and Remark 2.2, one has u -
For convenience, we let δ h := u I -u h . Using Proposition 3.2 and the fact that u and u h are solutions to respectively (2.3) and (3.5),
For any K ∈ T h , by Proposition 2.2 one had γA 0,∞,K 1 ≤ c * , thus, using Proposition 3.1 and the triangle inequality,
To estimate the term
By Proposition 2.3, one has f = A : ∇ 2 u almost everywhere in Ω, and thus
We deduce that
, and then, using Proposition 2.2, that
Collecting the above bounds we easily conclude.
Remark 3.2. The term |u -Π 0 m u| 2,h in (3.6) can be estimated using the classical Scott-Dupont theory, 12 depending on the regularity of u. Estimation of the term u -u I 2 requires the construction of a suitable interpolant u I ∈ V g I h of u. This depends on the specific definition of the local virtual element spaces V h,K , K ∈ T h . In Sec. 5.5 we provide an instance of such spaces for which we are able to prove optimal interpolation error bounds, and thus deduce optimal rate of convergence for the resulting method, cf. Theorem 6.1.
Remark 3.3. In the context of divergence form problems, error estimates for VEMs usually involve regularity of the coefficients and data, see for instance Theorem 6.2 in Ref. 17. In contrast, the error estimate (3.6) does not involve the functions A, f , and γ. This is thanks to Proposition 2.2 and, importantly, to the fact that no integration by parts was performed in order to construct the variational formulation (2.3), which allowed us to directly use the almost-everywhere equality A : ∇ 2 u = f in the proof of Theorem 3.1, rather than relying on the variational formulation (2.3) of the continuous problem. Observe that a similar argument was used in Ref. 33 for the analysis of discontinuous Galerkin and C 0 -interior penalty finite element methods. Remark 3.4. Typically in VEMs, the numerical solution u h ∈ V g I h to (3.5) is not accessible in practice. Instead, the projection Π * m u h is accessible (see Sec. 5 for more details). For this reason, it makes sense to consider the error u -Π * m u h 2 . This error is readily controlled combining the above estimate with the following Conforming VEM for Nondivergence Form Linear Elliptic Equations 89 result, which we isolate in a separate proposition since it is useful in general and can also be combined with Theorem 4.1 or Theorem 4.2, rather than Theorem 3.1:
Proof. By the triangle inequality and then (3.1), one has
By standard polynomial projection estimates,
, which concludes the proof.
As any finite element method, the assembly of a VEM requires some form of quadrature; see e.g. Ref. 17, where this issue is discussed in the case of elliptic problems in divergence form. Given that the case of irregular data is of paramount importance, it is particularly relevant to prove that the analysis detailed in the previous section extends to a variant of the method with numerical quadrature. Here, we perform such analysis under only slightly more restrictive assumptions. Since quadrature requires pointwise evaluation, we assume that, for any K ∈ T h , we are given specific Sobolev representatives of A K , f K , and γ K of A| K , f | K , and γ| K , so that A K (x), f K (x), and γ K (x) are well-defined at every x ∈ K. Note that we allow A K1 (x) = A K2 (x) at x ∈ ∂K 1 ∩ ∂K 2 , K 1 = K 2 , and likewise for f and γ; this is a natural way to handle data with jumps at cell interfaces. Definition 4.1. (Cordes condition everywhere) Let 0 ≤ µ < 1. We say that A satisfies the µ-Cordes condition everywhere if A K (x) satisfies the µ-Cordes condition for every K ∈ T h and x ∈ K. We say that γ is an everywhere µ-admissible scaling of A if γ K (x) > 0 and |γ K (x)A K (x) -I d | ≤ µ for every K ∈ T h and x ∈ K.
Let us now describe the family of quadrature rules that we allow in our analysis. Assumption 4.1. (Quadrature) On any cell K ∈ T h , we are given a finite set of quadrature points X K ⊂ K and corresponding nonnegative quadrature weights
For any K ∈ T h and u, v ∈ V h,K , we define
We then define the global counterparts a Q h and L Q h of a Q h,K and L Q h,K by summing over all K ∈ T h as usual, and we consider the following scheme with quadrature:
The two following results are easily proved by using the Cauchy-Schwarz inequality
, the exactness of quadrature for polynomials of degree 2m -4, and arguing as in Propositions 3.1 and 3.2. Proposition 4.1. (Continuity, with quadrature) For any K ∈ T h and any u, v ∈ V h,K , one has
Proposition 4.2. (Coercivity, with quadrature) Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1 and that γ is an everywhere µ-admissible scaling of A. Then for any u
, where c * is from (3.2).
The well-posedness of scheme (4.1) can then be established using the Lax-Milgram lemma similar to Corollary 3.2. Next, we derive an error estimate. We start with the following lemma. Lemma 4.1. Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1, that γ is an everywhere µ-admissible scaling of A, and that c * > µ in (3.2). Let u h ∈ V g I h be the unique solution to (4.1). Then, for any u ∈ H 2 (Ω) and
Proof. Let δ h := u I -u h . Arguing as in Theorem 3.1,
Conforming VEM for Nondivergence Form Linear Elliptic Equations 91
For any K ∈ T h , using Proposition 4.1 and the triangle inequality, and arguing as in Proposition 2.2 to control γ K A K , one has
. Finally, we estimate the remaining term
, from which we conclude as in Theorem 3.1.
Note that the above lemma holds for any u ∈ H 2 (Ω). We now want to choose u as the solution to (2.3) and use the equality f K = A K : ∇ 2 u almost everywhere in order to estimate the rightmost term in (4.2). The difficulty is that this equality is not guaranteed to hold at quadrature points. A sufficient condition for it to hold is that A K , f K , and ∇ 2 u are continuous at quadrature points. This is the setting of our next result.
Theorem 4.1. (Error estimate with quadrature, continuous setting) Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1, that γ is an everywhere µ-admissible scaling of A, and that c * > µ in (3.2). Let u ∈ V g and u h ∈ V g I h be respectively the unique solutions to (2.3) and (4.1). Assume that, for any K ∈ T h and x ∈ X K , A K , f K , and ∇ 2 u are continuous at x. Then, for any
Proof. Lemma 4.1 applies, so we only need to estimate the rightmost term in (4.2). Since the equality f K = A K : ∇ 2 u holds at every quadrature point, one has, for any
, where we argued as in Proposition 2.2 to control γ K A K .
It may happen that A or f are discontinuous at some quadrature points (see for instance our numerical experiments on odd square and cubic meshes in Secs. 7.1 and 7.2). In this case, we can still obtain some error estimate provided that the data and the solution u to (2.3) satisfy the following property:
there exists an open set E ⊂ Ω such that u ∈ W 2,∞ (E) and, for any quadrature point x ∈ X K , K ∈ T h , one has x ∈ E ∩ K, and
The motivation for introducing the open set E is to only assume some partial continuity of A and f at quadrature points. The condition u ∈ W 2,∞ (E) is the weaker condition on u that we were able to exploit in our analysis; while it holds for instance in all the cases we considered in our numerical experiments, we are unfortunately not aware of a general sufficient condition on the data that would guarantee this regularity of the solution.
Theorem 4.2. (Error estimate with quadrature, discontinuous setting) Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1, that γ is an everywhere µ-admissible scaling of A, and that c * > µ in (3.2). Let u ∈ V g and u h ∈ V g I h be respectively the unique solutions to (2.3) and (4.1). If there exists an open set E ⊂ Ω satisfying (4.4), then, for any
Proof. Lemma 4.1 applies, so we only need to estimate the rightmost term in (4.2). By Proposition 2.3, A : ∇ 2 u = f almost everywhere in Ω, and thus also in E. Let K ∈ T h , x ∈ X K , and let (x n ) n≥0 ⊂ E ∩ K be a sequence of points converging to x and such that A K (x n ) :
, where we argued as in Proposition 2.2 to control |γ K (x)A K (x)|. We deduce that, on any K ∈ T h ,
∞,E∩K , which concludes the proof. Conforming VEM for Nondivergence Form Linear Elliptic Equations 93 We refer to Sec. 6 for convergence rate results based on Theorem 4.2, and Secs. 7.1 and 7.2 for examples of situations in which Theorem 4.2 applies but Theorem 4.1 does not.
Remark 4.1. It may seem surprising that the regularity of A, f , and γ is not involved in the error estimates (4.3) and (4.5). This is thanks to the fact that we were able to rely on the strong form of the equation A : ∇ 2 u = f in the analysis, as in Theorem 3.1 (see Remark 3.3). Observe also that the analysis relies on the important assumption that, on each cell K ∈ T h , the same quadrature rule Q K is used in the definitions of both the bilinear form a Q h,K and the linear form L Q h,K .
In this section, we detail a realization of the H 2 -conforming virtual element discretization framework introduced in Sec. 3, and then discuss the properties of associated quasi-interpolation operators.
The development of H 2 -conforming VEMs has been considered in a number of publications. Virtual elements for plate bending problems were proposed in Ref. 14. Subsequently, virtual elements for fourth-order problems in two dimensions have been analyzed in Refs. 2, 3, 4, and 8 and extended to three dimensions in Refs. 7 and 20; see also Refs. 13 and 19 for the special case of tetrahedral elements.
A complete family of virtual element spaces for all space dimensions was recently presented in Ref. 20. They use the so-called enhancement technique introduced in Ref. 1, and their construction is hierarchical in the space dimension. We review this construction in Secs. 5.1-5.4, and show, by referring to Ref. 20 as appropriate, that it fits in the abstract discretization framework of Sec. 3. For simplicity, we detail here only the spaces obtained for d ≤ 3 and refer the reader to Ref. 20 for the general case.
In Sec. 5.5, we present and analyze Scott-Zhang type and Lagrange type interpolation operators. A slightly different estimate can be found in Ref. 20. However, this cannot be used directly in our context as it was designed for polyharmonic problems with homogeneous boundary conditions. Note also that our Lagrange type interpolation error estimate has a slightly unusual form, due to the fact that we want to be able to apply it in low regularity settings; we refer to Sec. 5.5.2 for more details.
We require some extra notation. We shall use the symbols v and e to indicate a generic vertex and edge of the partition, respectively, and the symbol F for a face when d = 3. Further, we assume that a complete set of unit normal vectors is given on every edge e and face F , denoted by n e,i , i = 1, d -1, and n F , respectively. Similarly, to every vertex v we associate a basis of unit vectors n v,i , i = 1, d, of R d . We do not assume orthogonality of the families n v,i and n e,i ; rather, we assume i=1 n e,i , n 2 . These nondegeneracy assumptions will be more convenient than orthogonality in Sec. 5.5.
For any K ∈ T h , we denote by V K the set of its vertices and by E K the set of its edges; further, when d = 3, we denote by F K the set of faces of K. For a smooth enough function v defined on K, being it scalar or vector valued, we set v := 1
we let h v be an appropriate local mesh size parameter associated to v, such as the average diameter of the elements sharing v as a vertex.
Given ω ⊂ R d and k ∈ N, we denote by P k (ω) the space of polynomials of degree k on ω and let M k (ω) denote a basis for such space. This notation is extended to k < 0 by fixing P k (ω) = {0} and M k (ω) = {0} in this case. The L 2 -projection onto P k (ω) will be denoted by Π 0,ω k . The starting point for the sequential construction of virtual element spaces in two-and three-dimensions is one-dimensional spaces used on the edges composing the skeleton of the partition. These are fixed as standard polynomial spaces as follows: for a segment e we consider the spaces V 0 h,e := P m-1 (e) and V h,e := P 3∨m (e) typically used to construct the classical C j -conforming, j = 0, 1, finite elements in one dimension, respectively.
Let d = 2 and K ∈ T h . We first introduce the enlarged virtual element space
h,e , ∀ e ∈ E K }.
Remark 5.1. It may be surprising that we impose ∆ 2 v ∈ P m (K) rather than ∆ 2 v ∈ P m-4 (K), but this is standard when using the virtual element enhancement technique. The reason is that imposing ∆ 2 v ∈ P m-4 (K) would be too restrictive to be compatible with the additional constraints in Eq. (5.3). We refer to Ref. 1 for more details.
We clearly have P m (K) ⊂ V h,K . Moreover, for any v ∈ V h,K , we have that v and ∇v are continuous on ∂K in consequence of the compatibility conditions implied by v ∈ H 2 (K), cf. Refs. 4, 20, and 32.
We define the local Hessian projection operator Π 2 m : H 2 (K) → P m (K) by
, then the definition of V h,K permits the evaluation of Π 2 m v directly in function of the following (incomplete for V h,K ) set of degrees of freedom.
Conforming VEM for Nondivergence Form Linear Elliptic Equations 95 Definition 5.1. (Two-dimensional degrees of freedom) On K ∈ T h , we define the local degrees of freedom:
∂n j e,i q dℓ, for all q ∈ M m+j-4 (e), j = 0, 1, i = 1, d -1, and e ∈ E K ;
• |K| -1 K vq dx, for all q ∈ M m-4 (K).
A depiction of the degrees of freedom corresponding to m = 2, 3, 4 is shown in Fig. 1.
we have that the evaluation on the right-hand side of (5.1) only requires the knowledge of the following moments: the internal moments of v of order up to m -4, the edge moments of v of order up to m -3, and the edge moments of ∂v/∂n e of order up to m -2. On the other hand, the evaluation on the right-hand side of (5.2) only requires the knowledge of v and its gradient at the vertices. The proof that Π 0 m-2 ∇ 2 v is also computable from the given degrees of freedom of v is similar; see e.g. Ref. 17 for an analogue result.
The local H 2 -conforming virtual element space may now be defined as Clearly, we still have P m (K) ⊂ V h,K , and it is easy to see that the dimension of the space V h,K equals the number of degrees of freedom listed above and hence these may be used as degrees of freedom for V h,K ; see Ref. 20 for the unisolvency proof. Crucially, given that Π 2 m v only depends on such degrees of freedom, V h,K is a proper linear subspace of V h,K . Moreover, from the equality Π 0
m-4 v, we deduce the computability of Π 0 m on V h,K as well. Remark 5.2. The enhancement technique leading to the definition of V h,K is required to achieve the computability of Π 0 m v. As such, we note that the enhancement technique is not strictly necessary for the discretization of our model problem in two dimensions as the computability of Π 0 m-2 ∇ 2 v suffices in this case, and this is guaranteed also without enhancement. However, we do require the enhanced space in view of defining the face values of the three-dimensional discrete functions. Moreover, the enhancement is required if lower-order terms are present in the model, cf. Ref. 11.
Let d = 3 and K ∈ T h . We start by fixing the traces of virtual functions over the boundary of K using two-dimensional virtual element spaces. In particular, the space V h,F , F ∈ F K , defined according to (5.3) is used for the value trace and the standard H 1 -conforming two-dimensional virtual element space V 0 h,F is used for the normal derivative. We introduce the latter here briefly, referring to e.g. Refs. 1, 17, and 20 for the details.
Let F ∈ F K . We first define an enlarged virtual element space as V 0 h,F := {v ∈ H 1 (F ) : ∆v ∈ P m-1 (F ) and v| e ∈ P m-1 (e), ∀ e ∈ Γ F }. We consider the local gradient projection Π 1 m-1 : H 1 (F ) → P m-1 (F ) given by
h,F , the projection Π 1 m-1 v is proven to be computable through the local degrees of freedom: (i) the value v(v) for all v ∈ V F , (ii) |e| -1 e vq ds for all q ∈ M m-3 (e) and e ∈ Γ F , and (iii) 1 |F | F vq dx, for all q ∈ M m-3 (F ). Then, the local H 1 -conforming space is defined as
having the above as degrees of freedom. Such space is characterized by the fact that P m-1 (F ) ⊂ V 0 h,F and the property of computability, beside Π 1 m-1 , of the L 2projection Π 0 m-1 . Having the virtual trace spaces sorted, we can now implement once more the enhancement technique to define an H 2 -conforming virtual element over K ∈ T h . First, the enlarged virtual element space is defined by
where Γ K := e∈E K e; observe that similarly, ∂K = F ∈F K F . Once again we have P m ⊂ V h,K and continuity on ∂K up to the gradient of all functions in V h,K , cf. Ref. 20. It is easy to prove once again Proposition 5.1, namely the computability of the Hessian projector (5.1), (5.2), given the following (incomplete for V h,K ) degrees of freedom.
Definition 5.2. (Three-dimensional degrees of freedom) On K ∈ T h , we consider the local degrees of freedom as those listed in Definition 5.1 plus
q for all q ∈ M m+j-4 (F ), j = 0, 1 and F ∈ F K .
The degrees of freedom corresponding to m = 2, 3, 4 are shown in Fig. 2. The local virtual element space V h,K for K ∈ T h when d = 3 can now be defined as in (5.3). The proof that the set of degrees of freedom in Definition 5.2 is unisolvent in V h,K and that they allow the evaluation of Π 0 m alongside Π 2 m can be found in Ref. 20.
In view of completing the construction of a VEM according to Sec. 3, it remains to fix a local projection operator Π * m satisfying (3.1) and a local stabilization form s h,K satisfying (3.2).
The local Hessian projection is used for the projection operator Π * m of Sec. 3; hence, we fix
It is easy to deduce from (4.11) in Ref. 20, the trace inequality, and classical polynomial projection estimates, 12 that Π 2 m (u -Π 0 m u) 2,K |u| 2,K , from which we deduce that (3.1) holds. For any u, v ∈ V h,K , we let
Π 0,e m+j-4
and, when d = 3, we also introduce
Note in particular that only the form s 0 h,K is nonzero when = 2. We thus fix the local stabilization form as
where β > 0 is some arbitrary scaling factor. It is shown in Ref. 20 (see Eq. (5.2) therein) that, for all v ∈ V h,K ,
Both inequalities in (3.2) immediately follow, with constants c * and c * proportional to β. This implies that this choice of the stabilization form is admissible. The inequality c * > µ, required in our proof of well-posedness of the scheme, holds provided that β is large enough. Computing some value of β for which the inequality c * > µ is theoretically guaranteed would require a detailed study of the constants in the proof of (5.4). Obtaining a sharp estimate of the critical value of β is likely difficult. In practice, we simply fixed β = 1 in our numerical experiments, which led to satisfactory results (see Sec. 7).
It is clear that the local virtual element spaces defined above allow for the construction of a global space V h as a subspace of H 2 (Ω) according to (3.3). The elements of V h may be identified in standard fashion by the collection of the local degrees of freedom in Definitions 5.1 and 5.2. Then, the space V 0 h may be defined according to (3.4), which amounts to fixing the relevant boundary degrees of freedom. Likewise, the space V g I h may be defined according to (3.4) once the function g I ∈ V h is chosen.
Let us now describe interpolation operators suitable for use with the virtual element spaces V h , V 0 h , and V g I h .
Conforming VEM for Nondivergence Form Linear Elliptic Equations 99
We aim to design an interpolation operator I SZ h : H 2 (Ω) → V h satisfying the following property:
This operator can be used in order to choose g I as g I := I SZ h g, and u I as u I := I SZ h u in the error estimates (3.6), (4.3), and (4.5). Property (5.5) guarantees that u I coincides with g I on and thus belongs to V g I h . Following the Scott-Zhang construction, we consider the set of sides Σ K associated to any K ∈ T h , i.e. Σ K := E K if d = 2 and Σ K := F K if d = 3. We let Σ h := K∈T h Σ K and denote by n σ a unit normal vector to each σ ∈ Σ h .
To every vertex v ∈ V h , we associate sides σ v,i ∈ Σ h , i = 0, d, such that v ∈ ∂σ v,i (conceptually, σ v,0 will be used to fix the value degree of freedom v(v) and σ v,1 , . . ., σ v,d will be used to fix the gradient degree of freedom ∇v(v); why these sides may need to be chosen differently will become apparent below when we explain how to enforce property (5.5)). Similarly, when d = 3, to every edge e ∈ E h we associate sides σ e,i ∈ Σ h , i = 0, d -1, such that e ⊂ ∂σ e,i .
In dimension d = 2, for any v ∈ H 2 (Ω) and K ∈ T h , we define the local interpolant
In dimension d = 3, we keep the three above conditions and add the following:
• For any e ∈ E K : Π 0,e m-4
It is easily verified that this uniquely determines each one of the degrees of freedom characterizing I SZ K v. We then define
To establish (5.5), we need to make some additional assumptions. First we require that for any boundary vertex v ∈ V h ∩ ∂Ω, the side σ v,0 is included in ∂Ω, and, if d = 3 that for any boundary edge e ∈ E h , e ⊂ ∂Ω, the side σ e,0 is included in ∂Ω. This is not enough: for instance, if v ∈ V h ∩ ∂Ω belongs to the boundary of a side of the polytope Ω, then the gradient degree of freedom (∇Iv)(v) is prescribed by (Iv)| ∂Ω , and thus should depend only on v| ∂Ω . If v belongs to the interior of a side of Ω, then only the tangential part of (∇Iv)(v) is prescribed by (Iv)| ∂Ω . Accordingly, we define δ v := d in the first case and δ v := d -1 in the second case, and we assume that for i = 1, δ v , the side σ v,i is included in ∂Ω and the vector n v,i is tangential to σ v,i . This can always be achieved upon choosing the family of vectors n v,i appropriately. If d = 3, then similarly for any boundary edge e ∈ E h , Left: on vertices of Ω, the value degree of freedom and both components of the gradient degrees of freedom need to be imposed by the boundary data. Right: on vertices of T h that belong to the interior of an edge of Ω, only the value degree of freedom and the tangential part of the gradient degrees of freedom (in blue, as opposed to the normal part in green) need to be imposed by the boundary data.
e ⊂ ∂Ω, we define δ e := d -1 if e is included in the boundary of a side of Ω, and δ e := d -2 if e is included in the interior of a side of Ω, and we assume that for i = 1, δ e , the side σ e,i is included in ∂Ω and the vector n e,i is tangential to σ e,i .
We refer to Fig. 3 for an illustration of the above construction in dimension two.
Let us now study the stability and the accuracy of the local interpolation operators I SZ K so defined. Our analysis is based on the following direct specialization of one of the inequalities in the norm equivalence result in Lemma 4.7 in Ref. 20.
|∂ nv,i v(v)|.
For any
We first prove the following stability result (that we call intermediate since we will be able to improve it later, see Corollary 5.1). Lemma 5.1. (Intermediate Scott-Zhang stability) For any v ∈ H 2 (Ω) and K ∈ T h , one has
Conforming VEM for Nondivergence Form Linear Elliptic Equations 101 Proof. We apply the bounds in Proposition 5.2 to I SZ K v and estimate in turns all terms on the right-hand side of such bounds. One has
where the equality follows from the definition of I SZ K v. Similarly, for any σ ∈ Σ K , using the inverse trace inequality,
For any v ∈ V K ,
where we used successive trace and polynomial inverse inequalities for the first inequality. Similarly, for i = 1, d,
Finally, if d = 3, we prove similarly that for e ∈ E K and i = 1, 2, Π 0,e m-4
This concludes the proof.
By a Bramble-Hilbert argument, we deduce the following interpolation error estimate.
Proof. By construction, one has I SZ K p = p for any p ∈ P m (ω K ), and in particular for
where in the last inequality we used Lemma 5.1, as well as the fact that v -Π 0,ω K m v 2,K is controlled by the term v -Π 0,ω K m v 2,ω K . We conclude using classical Scott-Dupont theory. 12 Let us mention that Lemma 5.1 can be improved as follows once we have Theorem 5.1: Corollary 5.1. (Scott-Zhang local stability) For any v ∈ H 2 (Ω) and K ∈ T h , one has
Proof. We use the triangle inequality
, and conclude by applying Theorem 5.1 with s = 1.
Let us now provide some discussion of Lagrange type interpolation. In particular, in view of justifying our numerical results in Sec. 7, we need to consider Lagrange interpolation in the minimal regularity case. This will lead to an error estimate involving both an H 2 and a W 1,∞ term, see Theorem 5.2. This is unusual, but natural: H 2 is required because we are estimating an error in the H 2 norm, and W 1,∞ (and even C 1 ) is required because Lagrange interpolation in our setting relies on pointwise evaluations of the gradient. Let us stress out that these considerations about minimal regularity are not specific to the VEM: in the simple setting of Lagrange interpolation for C 0 conforming finite elements, it would similarly be natural to estimate the interpolation error in a mix of the H 1 and L ∞ norms.
For any v ∈ H 2 (Ω) ∩ C 1 (Ω) and K ∈ T h , we define I L K v ∈ V h,K as the unique function in V h,K sharing with v all the degrees of freedom defined in Definition 5.1 (or Definition 5.2 if d = 3). We then define
Mimicking the analysis of Scott-Zhang interpolation, we start by proving the following lemma, and then we state our main interpolation estimate. Lemma 5.2. (Intermediate Lagrange stability) For any v ∈ H 2 (Ω) ∩ C 1 (Ω) and K ∈ T h , one has
Proof. We do not give the full details, since the proof follows the sketch of the one of Lemma 5.1. The main difference is that one cannot use L 2 -based trace inequalities, which would require more differentiability of v than we want to assume, so we rely instead on their L ∞ based counterparts.
Proof. By construction, one has I L K p = p for any p ∈ P m (K), and in particular for
where we used Lemma 5.1 in the last inequality. We conclude using classical Scott-Dupont theory. 12 Similar to the Scott-Zhang case, Lemma 5.2 can now be improved as follows.
Corollary 5.2. (Lagrange local stability) For any v ∈ H 2 (Ω)∩C 1 (Ω) and K ∈ T h , one has
Proof. We use the triangle inequality
, and conclude by applying Theorem 5.2 with s = 1 and r = 0.
Let us also state the following simplification of Lemma 5.2 when v ∈ W 2,∞ (K).
|v| s+1,∞,K .
In this section, we show how the error estimates derived in Sec. 3.3 and Sec. 4 can be combined with the interpolation results in Sec. 5.5 in order to derive rates of convergence for schemes (3.5) and (4.1). Throughout the section, we assume that we are working in the concrete VEM framework defined in Sec. 5. As explained in Remark 3.4, two different notions of error that make sense in the VEM setting are u -u h 2 and u -Π * m u h 2,h , so we state our results for the all-encompassing error u -u h 2 + u -Π * m u h 2,h . Combining Proposition 3.3 and Theorems 3.1 and 5.1, we can prove the following optimal error bound for scheme (3.5). Theorem 6.1. (Rate of convergence) Assume that A satisfies the µ-Cordes condition for some 0 ≤ µ < 1, that γ is a µ-admissible scaling of A, that µ > c * in (3.2), and that g h = I SZ h g. Let u ∈ V and u h ∈ V h denote respectively the unique solutions to (2.3) and (3.5)
Proof. This follows from using Proposition 3.3, applying the error estimate (3.6) with u I = I SZ h u, and then using Theorem 5.1 and classical Scott-Dupont theory 12 in order to estimate respectively the terms u -I SZ h u 2 and |u -Π 0 m u| 2,h .
In the case of scheme (4.1) involving numerical quadrature, one may similarly deduce the following error bound from Proposition 3.3, Theorem 4.2, and Corollary 5.3. Theorem 6.2. (Rate of convergence, with quadrature) Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1, that γ is an everywhere µ-admissible scaling of A, that µ > c * in (3.2), that g ∈ C 1 (Ω), and that g h = I L h g. Let u ∈ V and u h ∈ V h denote respectively the unique solutions to (2.3) and (3.5).
Proof. This follows from using Proposition 3.3, applying the error estimate (4.5) with u I = I L h u, and then using Corollary 5.3 and classical Scott-Dupont theory in order to estimate respectively the terms |u -I L h u| 2 and |u -Π 0 m u| 2,h , |u -Π 0 m u| 2,∞,E∩K .
Note that in the above theorem we had to assume u ∈ W s+1,∞ (T h ) rather than u ∈ H s+1 (T h ), due to the rightmost term in estimate (4.5). One benefit of this assumption is that it implies that u ∈ C 1 (Ω) and thus allows to use Lagrange rather than Scott-Zhang interpolation, which is why estimate (6.1) only involves the regularity of u in a broken norm. However, the assumption that u ∈ W s+1,∞ (T h ) with s > 1 may be too strong for some applications. It is possible to prove rates of convergence under weaker assumptions, but one may need to use ad hoc arguments. We give below an estimate that applies to the example in Sec. 7.1; see also Sec. 7.3 for an example in which u ∈ W 2,∞ (Ω).
In the estimate below, u is assumed to have regularity W s+1,∞ on some submesh T * h ⊂ T h , while also having W 2,∞ regularity globally. This estimate is intended to be applied when the submesh T * h is known to cover a large proportion of Ω; for instance, as in Sec. 7.1, T * h could be the set of all cells of T h except for those that intersect with a finite union of one-dimensional curves. Theorem 6.3. (Rate of convergence, with quadrature and a rough solution) Assume that A satisfies the µ-Cordes condition everywhere for some 0 ≤ µ < 1, that Conforming VEM for Nondivergence Form Linear Elliptic Equations 105 γ is an everywhere µ-admissible scaling of A, that µ > c * in (3.2), that g ∈ C 1 (Ω), and that g h = I L h g. Let u ∈ V and u h ∈ V h denote respectively the unique solutions to (2.3) and (3.5). Assume that there exists an open set E ⊂ Ω satisfying (4.4), that u ∈ W 2,∞ (Ω), and that u ∈ W s+1,∞ (T * h ) for some 1 ≤ s ≤ m and T * h ⊂ T h . Then
Proof. Let E h := u -u h 2 + u -Π * m u h 2,h . We use Proposition 3.3, apply the error estimate (4.5) with u I = I L h u, and deduce that
In this section, we apply the numerical scheme (4.1) to problems featuring various regularity properties. These problems originate from or are inspired by the ones in Ref. 41. Our implementation is based on the code vem++. 23 The sparse linear system resulting from the discretization is solved using a direct solver. We use virtual element spaces featuring the polynomial consistency order m = 2. We use the barycenter quadrature rule: Q K [ϕ] := |K|ϕ(x K ), where x K denotes the barycenter of the cell K ∈ T h . This is consistent with Sec. 4, in which, for the case m = 2, the quadrature rule Q K is only assumed to be exact for constants. We use a Lagrange interpolation of the boundary data: g I := I L h g. For each considered problem, we compute the broken H 2 -seminorm, H 1seminorm, and L 2 -norm errors between u and Π * 2 u h , where u and u h are solution to respectively the continuous problem (2.3) and its discretization (4.1) -since the definitions of those errors involve integration of the exact solution u, in practice we approximate them using high-accuracy (fourth-order) quadrature rules.
We compare the H 2 -seminorm errors with the predictions of the theory. Observe that we did not prove higher-order error estimates for the error in H 1 and L 2 norms; this would be nontrivial, since this would require elliptic regularity estimates for the dual equation to (1.1). We still observe experimentally consistently higher orders of convergence in those norms. It is easily verified that γA ∈ [L ∞ (Ω)] d×d , γf ∈ L ∞ (Ω), and u ∈ W 2,∞ (Ω).
We solve scheme (4.1) on uniform square or cubic meshes obtained by dividing each edge of Ω into N = 4, 8, 16, . . . or N = 5, 9, 17, . . . intervals, and on polytopal meshes obtained by employing few iterations of the Lloyds algorithm from a random seed, see Fig. 4. For even square and cubic meshes, one has u ∈ C ∞ (T h ), thus Theorem 6.2 applies and predicts convergence of order 1 with respect to h in the H 2 norm. For the other meshes, Theorem 6.3 applies with T * h = {K ∈ T h | ∀ x ∈ K, ∀ 1 ≤ i ≤ d, x i = 0} and predicts convergence of order 1/2 with respect to h in the H 2 norm, since K∈T h \T * h |K| h. The numerical results are consistent with those predictions, as illustrated in Fig. 5.
As in the first problem, we consider Eq. (1.1) on the domain Ω = (-1, 1) d in dimensions d ∈ {2, 3}, with coefficients given by (7.1) and g(x) := 0, but now we choose f such that the exact solution u is
where we used Corollary 5.3 and classical Scott-Dupont theory for the second inequality. Now assume that |x K | ≤ h K . Using that
together with Theorem 5.2 and classical Scott-Dupont theory, one has
We note that
8 , and, by a direct computation, that |u| 2 1,∞,K + |u| 2 2,K h 1.2 K ≤ h 4 K |x K | -2.8 , which concludes the proof. The above result yields the error estimate u -Π * 2 u h 2,h c π c * c * -µ K∈T h h 4 K |x K | -2.8 1/2
, which we may further estimate depending on the mesh at hand. We consider either meshes which are uniform or graded toward the origin. On uniform meshes, using that h |x K | for all K ∈ T h , one has
, where the sum inside the parentheses is similar to integrating |x| -2 over Ω = (0, 1) 2 ; accordingly, we expect convergence of order 0.6 with respect to h, up to a logarithm. The numerical results we obtain are consistent with this expectation, see Fig. 8. Inspired by Lemma 7.1 and following a principle of error equidistribution, we design graded meshes by recursive quadrisection toward the origin so that h 4 K |x K | -2.8 lies below some given threshold for all K ∈ T h , see Fig. 7. According to the virtual element philosophy, we interpret hanging nodes as additional vertices of the polygonal cells. We observe from Fig. 8 that graded meshes allow recovering optimal convergence in the H 2 norm. 110 G. Bonnet, A. Cangiani & R. H. Nochetto In the H 2 norm and on uniform meshes, we observe convergence at the reduced rate 0.6/d with respect to the number of degrees of freedom, due to the singularity of the problem. The optimal rate of convergence in the H 2 norm is recovered with meshes suitably graded toward the origin.
G. Bonnet, A. Cangiani & R. H. Nochetto
G. Bonnet, A. Cangiani & R. H. Nochetto
G. Bonnet, A. Cangiani & R. H. Nochetto
G. Bonnet, A. Cangiani & R. H. Nochetto