hp–Version discontinuous Galerkin methods on polygonal and polyhedral meshes

An hp–version interior penalty discontinuous Galerkin method (DGFEM) for the numerical solution of second–order elliptic partial differential equations on general computational meshes consisting of polygonal/polyhedral elements is presented and analysed. Utilizing a bounding box concept, the method employs elemental polynomial bases of total degree p (Pp-basis) defined on the physical space, without the need to map from a given reference or canonical frame. This, together with a new specific choice of the interior penalty parameter which allows for face-degeneration, ensures that optimal a priori bounds may be established, for general meshes including polygonal elements with degenerating edges in two dimensions and polyhedral elements with degenerating faces and/or edges in three dimensions. Numerical experiments highlighting the performance of the proposed method are presented. Moreover, the competitiveness of the p-version DGFEM employing a Pp-basis in comparison to the conforming p-version finite element method on tensor-product elements is studied numerically for a simple test problem.

1. Introduction. The original articles concerned with the construction and mathematical analysis of Discontinuous Galerkin (DG) methods date back over 50 years ago. For hyperbolic partial differential equations, in 1973 Reed & Hill, cf. [29], developed the first DG discretization of the neutron transport equation. Independently, DG methods were constructed for elliptic problems based on weakly enforcing Dirichlet boundary conditions; see, for example, [7,8,25,28]. In particular, we highlight the works of Nitsche [28] and Baker [9], which form the basis of the class of interior penalty DG methods, cf. also [5,36]. Since the very early work, DG methods were partially abandoned, in part due to the increase in the number of degrees of freedom compared, for instance, with their conforming counterparts. However, in the last two decades there has been a renewed interest in the field of discontinuous discretizations both from a theoretical and computational viewpoint, cf. [17,24,30,18], for example. This resurgence is due to the inherent advantages offered by DG schemes, such as, for example, the limited interelement communication, which is restricted only to neighbouring elements, the local conservativity property, the simplicity in treating meshes with hanging nodes, and the development of efficient hp-adaptivity refinement strategies. Moreover, recently in [10,11,12,15] it has been shown that the underlying DG polynomial bases may be efficiently constructed in the physical frame, without needing to map local polynomial spaces defined in a given reference/canonical frame. In this way, DG methods can easily deal with general-shaped elements, including polygonal/polyhedral elements, cf. [1,10,15]. The flexibility of DG methods in handling general meshes has no immediate counterpart in the conforming framework, where the design of suitable finite element spaces for meshes of polygons/polyhedra is far from being a trivial task. Several examples include the Composite Finite Element Method [23,22], the Polygonal Finite Element Method [33,34], the Extended Finite Element Method [19] and the most recent Virtual Element Method [13].
At present, the design of solvers and preconditioners for DG discretizations on nonstandard grids lends itself to huge developments in the field of numerical analysis. Indeed, to the best of our knowledge, the only study regarding solution techniques for this class of problems is reported in [2], where a non-overlapping Schwarz preconditioner for composite DG finite element methods on complicated domains is analyzed. In the current article we exploit the theoretical framework developed in [15] to study the performance of a two-level and W-cycle multigrid solver. The possibility to employ general-shaped elements in the physical framework makes the choice of multilevel schemes natural. The flexibility afforded by this approach allows us to define the set of grids needed in the multigrid algorithm by agglomeration; thereby, the definition of the associated subspaces is straightforward, since inter-element continuity is not required. This property overcomes the usual difficulties encountered in the construction of agglomeration multigrid schemes in the conforming framework, where the agglomeration strategy must be followed by a proper definition of the conforming subspaces. In [16], for example, the sublevels are obtained by combining a graph based agglomeration algorithm and re-triangulations, thus resulting in a set of non-nested grids, while the associated nested subspaces are defined by introducing suitable interpolation operators. The resulting V-cycle multigrid algorithm converges uniformly with respect to the meshsize h provided that the number of levels is kept fixed.
In this paper we analyze the convergence of a two-level scheme and W-cycle multigrid method for the solution of the linear system of equations arising from the hp-version of the interior penalty DG scheme on polygonal/polyhedral meshes [15], thereby, extending the theoretical framework developed in [4] for standard quasi-uniform meshes. Our analysis is based on the smoothing and approximation properties associated with the proposed method: the former corresponds to a Richardson iteration, whose study requires a result concerning the spectral properties of the stiffness matrix, while the latter is inherent to the interior penalty DG scheme itself and exploits the error estimates derived in [15]. We show that, under suitable assumptions on the agglomerated coarse grid, the two-level method converges uniformly with respect to the granularity of the underlying partition and the polynomial approximation degree p, provided that the number of smoothing steps is chosen of order C(p)p 2 , where C(p) is a constant, which in general depends on p and the geometric properties of the grids. This implies that the generation of good quality agglomerated meshes is fundamental for the performance of the solver. We prove an analogous result for the case of W-cycle multigrid scheme. However, we remark that, in addition to the geometric assumptions on the agglomerated grids, we also require that the maximum number of element faces remains bounded. Due to the agglomeration process, the latter assumption is critical in the case of multilevel methods, which implies that our analysis is only valid if the number of levels is reasonably small. Throughout the analysis, we also track the dependence of the error reduction factor of the two solvers on the polynomial approximation order p, thereby recovering a similar result to the case when standard quasi-uniform meshes are employed, as well as the geometrical properties of the agglomerated grids.
The rest of this paper is organized as follows. In Section 2 we introduce the interior penalty DG scheme for the discretization of second-order elliptic problems on general meshes consisting of polygonal/polyhedral elements. Then in Section 3, we recall some preliminary analytical results concerning this class of schemes. In Section 4 we define the multilevel framework and introduce several technical results. We then focus first on the analysis of two-level method, followed by the extension to the W-cycle multigrid solver, where we assume that the number of levels obtained by agglomeration is kept limited. The main theoretical results are investigated through a series of numerical experiments presented in Section 5. In particular, we show that, in general, the limitation on the number of levels employed in the W-cycle multigrid solver does not seem to be restrictive in practice.
2. Model problem and discretization. We consider the weak formulation of the Poisson problem, subject to a homogeneous Dirichlet boundary condition: with Ω ∈ R d , d = 2, 3, a convex polygonal/polyhedral domain with Lipschitz boundary and f a given function in L 2 (Ω).
For the sake of brevity, throughout this article, we write x y and x y in lieu of x ≤ Cy and x ≥ Cy, respectively, for a positive constant C independent of the discretization parameters. Moreover, x ≈ y means that there exist constants C 1 , C 2 > 0 such that C 1 y ≤ x ≤ C 2 y. When required, the constants will be written explicitly.
In view of the forthcoming multigrid analysis, we denote by {T j } J j=1 a sequence of partitions of the domain Ω, each of which consists of disjoint open polygonal/polyhedral elements κ of diameter h κ , such that Ω = κ∈Tjκ , j = 1, . . . , J. We denote the mesh size of T j , j = 1, . . . , J, by h j = max κ∈Tj h κ . To each T j , j = 1, . . . , J, we associate the corresponding discontinuous finite element space V j , j = 1, . . . , J, defined as where P pj (κ) denotes the space of polynomials of total degree at most p j ≥ 1 on κ ∈ T j . A suitable choice of the sequences {T j } J j=1 and {V j } J j=1 leads to the so-called h-, p-, and hp-multigrid schemes. In particular, the h−multigrid method is based on employing a constant polynomial approximation degree for each j, j = 1, . . . , J, (i.e., p j = p), on a set of nested partitions {T j } J j=1 , such that the coarse level T j−1 , j = 2, . . . , J, is obtained by agglomeration from T j in such a way that

2)
i.e., we consider a bounded variation hypothesis between subsequent levels. In the p-multigrid method, the partition is kept fixed for any j, j = 1, . . . , J, while we assume that the polynomial degrees vary moderately from one level to another, i.e., 3) The hp-multigrid method is obtained by combining these two strategies. Note that with the above choices we obtain nested finite element spaces V j , j = 1, . . . , J, i.e., V 1 ⊆ V 2 ⊆ · · · ⊆ V J .

Grid assumptions.
In this section, we introduce some additional notation from [15], and outline some key definitions and assumptions. For any T j , j = 1, . . . , J, when no hanging nodes/edges are included in the partition, we define the interfaces of the mesh T j as the set of (d − 1)-dimensional facets of the elements κ ∈ T j . The presence of hanging nodes/edges, on the other hand, can be handled by defining the interfaces of T j as the intersection of the (d − 1)-dimensional facets of neighboring elements. This implies that, for d = 2, an interface will always consist of a line segment. For d = 3, we assume that for each interface of an element κ ∈ T j , a sub-triangulation into co-planar triangles is given. We then denote by "face" a (d − 1)-dimensional simplex (i.e., a line segment for d = 2), which is part of the boundary of κ ∈ T j . As a consequence, in the two dimensional case, the face and interface of an element κ ∈ T j coincide. With this notation in mind, we denote by F j the set of all mesh interfaces if d = 2 and the set of all open triangles belonging to the sub-partition of all mesh interfaces if d = 3. Moreover, we have that F j = F I j ∪ F B j , where F I j is the set of interior element faces of T j , such that F ⊆ ∂κ + ∩ ∂κ − for any F ∈ F I j , where κ ± are two adjacent elements in T j . The set F B j contains the boundary element faces, i.e., F ⊂ ∂Ω for F ∈ F B j . With this notation, we introduce the following assumptions on the partitions T j , j = 1, . . . , J; in the case of the h-and hp-multigrid schemes, these assumptions must be satisfied for the meshes generated by the underlying agglomeration process.
Assumption 2.1. The number of faces N κ of any κ ∈ T j , j = 1, . . . , J, is uniformly bounded, i.e., there exists a constant C F such that Assumption 2.1 is critical in our multilevel framework, because the number of faces N κ grows with the number of levels due to the agglomeration process. As a consequence, this assumption is only satisfied if the number of levels is kept limited. However, it will be demonstrated in Section 5, that this assumption does not seem to be a limitation in practice. Assumption 2.2. For any κ ∈ T j , j = 1, . . . , J, we assume that . . , J, denote a covering of Ω consisting of shape-regular d-dimensional simplices K. We assume that, for any κ ∈ T j , j = 1, . . . , J, there exists K ∈ T j such that κ ⊂ K and card κ ∈ T j : κ ∩ K = ∅, K ∈ T j such that κ ⊂ K 1.
Consequently, for each pair κ, K ∈ T j , with κ ⊂ K, To keep the notation as simple as possible we will assume that our decompositions are quasiuniform.
Assumption 2.4. We assume that the mesh size h j , j = 1, . . . , J, is such that We remark that the above assumption can be weakened and, according to [15], only a local bounded variation property is needed for our theoretical analysis; see Remark 4.9 below for details. We will also need the following definitions.
Definition 2.5. For each κ ∈ T j , j = 1, . . . , J, we denote by F κ the set of all possible dsimplices contained in κ and having at least one face in common with κ. Moreover, we denote by κ F , an element in F κ sharing a face F with κ ∈ T j . Definition 2.6. For any κ ∈ T j , j = 1, . . . , J, we denote by T κ the family of all possible subtessellations T κ of κ consisting of d-simplices τ , such thatκ = τ ∈Tκτ , and write h τ to denote the diameter of τ ∈ T κ . 4 2.2. DG formulation. The definition of the proceeding DG method is based on employing suitable jump and average operators. To this end, for (sufficiently smooth) vector-valued and scalar functions τ and v, respectively, we define jumps and averages across F ∈ F j , j = 1, . . . , J, as follows: where v ± and τ ± denote the traces of v and τ on F taken from the interior of κ ± , respectively, and n ± the outward unit normal vectors to ∂κ ± , respectively, cf. [6].
On any level j, j = 1, . . . , J, we consider the bilinear form A j (·, ·) : V j × V j → R, corresponding to the symmetric interior penalty DG method, defined by where σ j ∈ L ∞ (F j ) denotes the interior penalty stabilization function. To define the stabilization function σ j , j = 1, . . . , J, we first introduce an inverse inequality for polygonal/polyhedral elements; to this end, we recall the following definition. Definition 2.7. Let T j , j = 1, . . . , J, denote the subset of elements κ ∈ T j , such that each polygonal/polyhedral element κ ∈ T j can be covered by at most m Tj shape-regular d-simplices K i , i = 1, . . . , m Tj , such that and |K i | c as |κ|, The following inverse inequality for general-shaped elements is derived in [15,Lemma 4.4]. Lemma 2.8. Let κ ∈ T j , j = 1, . . . , J, be a polygonal/polyhedral element, and let F ⊂ ∂κ be one of its faces, and T j be defined according to Definition 2.7. Then, for each v ∈ P pj (κ), we have and κ F ∈ F κ as in Definition 2.5. The positive constant C inv is independent of |κ|/ sup κ F ∈κ |κ F |, p j and v.
The interior penalty stabilization function σ j : F j → R + is then given by with C j σ > 0 independent of p j , |F | and |κ|, cf. [15]. In this article, we develop two-level and W-cycle multigrid schemes to compute the solution of the following problem on the finest level J: find u J ∈ V J such that (2.6) 3. Preliminary results. We first endow the finite element spaces V j , j = 1, . . . , J, with the following DG norm: The well-posed of the DG formulation is established in the following lemma. Lemma 3.1. The following continuity and coercivity bounds, respectively, hold where C cont and C coerc are positive constants, independent of the discretization parameters, provided that C j σ > C F , j = 1, . . . , J. Proof. See [15]. The proceeding error estimates are based on the following approximation result, which is a simplified version of the analogous bound presented in [ with s = min{p j + 1, k}, and Analogously, We point that, as for Lemma 3.2 stated above, any bound derived under the validity of Assumption 2.1 will necessarily lead to a dependence on C F in the resulting constant. Next, we state error bounds for the underlying interior penalty DG scheme in terms of both the DG and L 2 (Ω)-norms, cf. [15]. Theorem 3.3. Assume Assumptions 2.1, 2.2, 2.3, and 2.4 hold. We denote by u j ∈ V j , j = 1, . . . , J, the DG solution of problem (2.6) posed on level j, i.e., If the solution u of (2.1) satisfies u| κ ∈ H k (κ), k > 1 + d/2, such that Eu| K ∈ H k (K), for each κ ∈ T j , j = 1, . . . , J, where κ ⊂ K, K ∈ T j , then the following results hold Proof. The error bound (3.3) stems from the general result derived in [15,Theorem 5.2]. We now proceed with the proof of the bound on the L 2 (Ω)-norm of the error, cf. (3.4). To this end, we employ a standard duality argument: let w ∈ V , be the solution of the problem Exploiting a standard elliptic regularity assumption, we note that According to Galerkin orthogonality, we immediately obtain which together with (3.3) gives the desired result. We also need to introduce an appropriate inverse inequality; to this end, we first recall the following result, cf. [20,Lemma 3.7].
Lemma 3.4. Let K be a shape-regular simplex. Then for any v ∈ P p (K) there exists a simplex κ ⊂ K, with the same shape as K and faces parallel to the faces of K, with dist(∂κ, ∂K) < C as diam(K)/p 2 , for some constant C as > 0, independent of v, K and p, such that With the above result we now prove the following lemma. Lemma 3.5. For any v ∈ V j , j = 1, . . . , J, the following inverse estimate holds Proof. For κ ∈ T j \T j , we recall the family T κ of sub-tessellations T κ of κ defined in Definition 2.6. We then have, by standard inverse estimates on simplicial elements, that In order to obtain a sharp bound, we take the supremum over all T κ ∈ T κ , namely, For κ ∈T j , we consider the covering of κ by shape-regular simplices K i , i = 1, . . . , m Tj , such that see Definition 2.7. We recall the following inverse estimates on the simplex K i , cf. [31,20], for any u ∈ P pj (K i ). By exploiting the covering of κ ⊂ m T j i=1 K i and the bounds (3.7) and (3.8), we obtain We now defineκ i ⊂ K i to denote the simplex relative to K i as outlined in Lemma 3.4. Hence, utilizing Lemma 3.4 and Definition 2.7, gives sinceκ i ⊂ κ, and henceκ i ⊂ K i ∩ κ ⊂ K i , cf. [15]. Substituting (3.10) into (3.9), and employing inequality (3.6) gives where, given (3.6), we assume that h Ki h κ . We then take the minimum between (3.5) and (3.11) to deduce the desired result. The inverse estimate presented in Lemma 3.5 is fundamental to the proof of the following upper bound on the maximum eigenvalue of A j (·, ·). We remind that the analogous result on standard grids can be found in [3].
Theorem 3.6. Given that Assumptions 2.1, 2.2 and 2.4 hold, then for any u ∈ V j , j = 1, . . . , J, we have that Proof. Given the continuity of the bilinear forms A j (·, ·) stated in Lemma 3.1, we restrict ourselves to estimate the two terms involved in the DG norm. The local contributions of the H 1 seminorm can be simply bounded by applying Lemma 3.5 and the quasi-uniformity of the partition, i.e., κ∈Tj |u| 2 For the norm of the jump across F ∈ F I j , with F ⊂ ∂κ + ∩ ∂κ − , we employ the inverse inequality (2.5); thereby, Summing over the internal faces and employing Assumptions 2.1, 2.2, and 2.4 gives An analogous result also holds for boundary faces; the statement of the theorem now follows immediately.
The theoretical results derived in this section form the basis of the analysis of the proposed multigrid algorithms presented in the following section.
4. Two-level and W-cycle multigrid algorithms. The forthcoming analysis is based on the classical multigrid theoretical framework already employed in [4] for high-order DG schemes on standard quasi-uniform meshes. The two key ingredients in the construction of our proposed multigrid schemes are the inter-grid transfer operators and the smoothing scheme. The prolongation operator connecting the space V j−1 to V j , j = 2, . . . , J, is denoted by I j j−1 : V j−1 → V j , while its adjoint with respect to the L 2 (Ω)-inner product (·, ·) is the restriction operator I j−1 j : V j → V j−1 : As a smoothing scheme, we choose a Richardson iteration, whose operator is defined as: For the definition of the solvers, we first address the two-level method. Given the following problem in Algorithm 1 we outline the two-level cycle, where MG 2lvl (z 0 , m 1 , m 2 ) denotes the approximate solution obtained after one iteration, with initial guess z 0 and m 1 , m 2 pre-and post-smoothing steps, respectively. As a multilevel extension of Algorithm 1, we consider a standard W-cycle scheme. On level j, we consider for a given g ∈ V j . The approximate solution obtained by applying the j-th level iteration to the above linear system, with initial guess z 0 and m 1 , m 2 number of pre-and post-smoothing steps, respectively, is denoted by MG W (j, g, z 0 , m 1 , m 2 ). On the coarsest level j = 1, the corresponding subproblem is solved based on employing a direct method, i.e., MG W (1, g, z 0 , m 1 , m 2 ) = A −1 1 g, while for j > 1 we apply the recursive procedure outlined in Algorithm 2. We observe that Algorithm 1 can be considered as a special case of Algorithm 2, corresponding to J = 2.

Convergence analysis of the two-level method.
We first define the following norms based on the operator A j , j = 1, . . . , J, Hence, For the proceeding analysis, we need to introduce some additional hypotheses on the agglomerated meshes employed both within the two-level method studied in this section, as well as the W-cycle multigrid algorithm analyzed in Section 4.2. To this end, for any F ∈ F j ∩ F j−1 , j = 2, . . . , J, we denote by κ ± j and κ ± j−1 the neighboring elements sharing the face F in T j and T j−1 , respectively. It is trivial to see that κ ± j ⊂ κ ± j−1 , since the grids are nested. We then assume that, there exists Θ > 0 such that and which implies, together with (2.3), that for any F ∈ F j ∩ F j−1 , j = 2, . . . , J. Remark 4.1. The above assumption is satisfied if the agglomeration algorithm preserves the shape-regularity of the elements. In Figure 4.1, we show two examples of possible macroelements: the agglomerate on the left is not suitable to guarantee assumption (4.4) due to the presence of a dominant dimension, while the element on the right can be considered appropriate. Moreover, we note that the fulfilment of the above geometric assumptions (4.3) and (4.4) can be considered a good criterion in evaluating the quality of the agglomerated grids employed in the multigrid algorithm.
In order to undertake the convergence analysis of the two-level solver outlined in Algorithm 1, we follow the approach developed in [4]. We then provide an estimate based on the error propagation operator, which is defined as We now study separately the smoothing property and the approximation property. Before proceeding with the analysis, we first observe that by Theorem 3.6, we can bound Λ j , j = 1, . . . , J, in (4.1) as follows The last result is employed to prove the smoothing property in the next lemma; see [4,Lemma 4.3] for the proof.
for 0 ≤ t < s ≤ 2 and m ∈ N \ {0}.   Proof. For any v ∈ V J , we consider the following equality . (4.9) Next, we consider the solution η of the following problem for φ ∈ L 2 (Ω), and let η J ∈ V J and η J−1 ∈ V J−1 be the corresponding DG approximations in V J and V J−1 , respectively, given by (4.10)

By Theorem 3.3 and the hypotheses (2.2) and (2.3), we deduce that
and from a standard elliptic regularity assumption, it follows that (4.14) Substituting (4.14) into (4.9) gives the desired result.
The convergence result for the two-level method, involving the error propagation operator E 2lvl for any v ∈ V J , with . Therefore, the two-level method converges uniformly provided the number of pre-and post-smoothing steps satisfy for a positive constant χ > C 2lvl . Proof. The statement of the theorem follows in a straightforward manner by applying the smoothing property (4.7) twice, the approximation property (4.8) and exploiting the bounded variation assumptions (2.2) and (2.3).
We observe that the rate of convergence is independent of the mesh size, but depends on p J , as in the case of standard quasi-uniform meshes. Moreover, a dependence on p J is also hidden inC(p J ), which also involves the geometric properties of the partitions. As a consequence, a good quality agglomerated coarse grid is fundamental to guarantee the uniformity of the solver.

4.2.
Convergence of the W-cycle multigrid algorithm. As mentioned at the beginning of Section 2, we recall that in the multilevel case, Assumption 2.1 represents a critical issue. In the following analysis, we then assume that the number of levels is limited, in such a way that the number of interfaces on each level can be bounded by a constant C F that does not lead to an excessive over-penalization due to the penalization parameter C j σ , j = 1, . . . , J.
To proceed, we first need to establish the equivalence between DG norms on subsequent grid levels. We point out that, in contrast to the case of standard quasi-uniform grids presented in [4], such an equivalence result does not follow in a straightforward manner; indeed, here we need to exploit the hypotheses (4.3) and (4.4) introduced in the previous section. Under these assumptions, the proof of the following result follows immediately. where C equiv = C equiv (Θ), in general, depends on the quality of the agglomerated grids. Lemma 4.5 is essential to deduce the stability of the operators I j j−1 and P j−1 j , j = 2, . . . , J. In particular, we state the following bounds.
Lemma 4.6. There exists a positive constant C stab , independent of the mesh size, the polynomial approximation degree and the level j, j = 2, . . . , J, such that The proof of Lemma 4.6 is based on employing inequality (4.16); for details, see [4,Lemma 4.6].
Remark 4.7. We stress that the constant C stab depends on C equiv in (4.16), which means that the quality of the agglomerated meshes plays a crucial role in keeping this constant bounded, thus resulting in the uniformity with respect to the mesh size and the number of levels as shown in Theorem 4.8 below.
The error propagation operator associated to Algorithm 2 is defined as where G j = Id j − B −1 j A j and P j−1 j is defined analogously to (4.6), cf. [21,14]. Then the convergence estimate for the W-cycle multigrid scheme follows from Theorem 4.4 and the stability estimates (4.17) and (4.18).
Theorem 4.8. Let C 2lvl and C stab be defined as in Theorem 4.4 and Lemma 4.6, respectively. Then, there exists a constant C > C 2lvl , independent of the mesh size, the polynomial approximation degree and the level j, j = 1, . . . , J, such that, if the number of pre-and post-smoothing steps satisfy otherwise, (4.20) then |||E j,m1,m2 v||| 1,j ≤ CΣ j |||v||| 1,j ∀v ∈ V j , (4.21) with Proof. The proof follows the derivation of the analogous result presented in [4,Theorem 4.7]. For j = 1, the statement of the theorem trivially holds. For j > 1, by an induction hypothesis, we assume that (4.21) holds for j − 1. By the definition of the error propagation operator E j,m1,m2 v in (4.19), it follows that The first term corresponds to a two-level method between level j and j − 1. We now observe that the smoothing property of Lemma 4.2 and the approximation property of Lemma 4.3 can be extended to any level V j , j = 2, . . . , J, and we therefore have, by Theorem 4.4, that The bound on the second term is obtained by applying the smoothing property (4.7) for j = 2, . . . , J, the stability estimates (4.17) and (4.18) and the induction hypothesis; thereby, we get We then obtain We now bound Σ j−1 with Σ j as follows: ifC(p j−1 ) ≤C(p j ), then Otherwise, ForC(p j−1 ) ≤C(p j ) by (4.23), we have that We then observe that if m 1 and m 2 are such that it follows that ForC(p j−1 ) >C(p j ), starting from (4.24) and following the same steps, we deduce the statement of the theorem, provided m 1 and m 2 are such that As in the two-level case, inequality (4.21) implies that the convergence of the method is guaranteed if the number of smoothing steps is chosen sufficiently large, cf. (4.20). Moreover, compared to the case of standard quasi-uniform grids, cf. [4], the bound (4.20) on the number of smoothing steps involves a strong dependence on the geometrical properties of the underlying agglomerated meshes, which in principle, could lead to restrictive conditions on the hierarchy of grids employed. However, we remark that, in practice, the numerical simulations indicate that the proposed multigrid algorithms converge uniformly, even when low quality agglomerated grids are employed; moreover, an increase in the polynomial order does not seem to require a higher number of smoothing steps to obtain a convergent iteration, cf. Section 5 for details. Remark 4.9. Whenever the agglomerated grids are not quasi-uniform, i.e., Assumption 2.4 is not satisfied, Theorem 4.4 and Theorem 4.8 still hold. More precisely, we need to introduce the ratio θ j between the maximum and minimum element size on level j Moreover, we assume that there exists a constant C mesh , independent of the granularity of the mesh, such that θ j ≤ C mesh , j = 1, . . . , J.
Then the results in Theorem 4.4 and Theorem 4.8 hold with cf. (4.22), and the bound (4.20) is modified as follows (4.25) Remark 4.10. We recall that in Theorem 4.8 and Remark 4.9, in order to guarantee the convergence of the method, we require a lower bound on the number of smoothing steps, cf. (4.20) and (4.25). In fact, forC(p j−1 ) ≤C(p j ), we obtain An analogous result can be obtained forC(p j−1 ) >C(p j ). Moreover, we note that we have considered the general case of (4.25), since (4.20) can be regarded as a particular case. Each initial grid is then subsequently coarsened in order to obtain a sequence of nested partitions by employing the software package MGridGen [26,27]. Before testing the performance of the twolevel and W-cycle multigrid solvers presented in Algorithm 1 and Algorithm 2, respectively, we first address the issue of the choice of the penalization coefficient C j σ in (2.4). According to Lemma 3.1, the bilinear form A j (·, ·) is coercive provided that C j σ > C F , with C F an upper bound for the maximum number of element interfaces in the partition T j ; see Assumption 2.1. In Table 5.1, we report the coercivity constant C coerc of (3.1) for a fixed value of C j σ ≡ C σ = 10 for j = 1, . . . , 4. We observe that, despite the increase in the value of C F from grid G1 to grid G4, the bilinear form is uniformly coercive for a constant value of the penalization coefficient, which in general does not satisfy the theoretical assumption. As a consequence, in the following, we set C j σ ≡ C σ = 10 for j = 1, . . . , 4. We now consider the grids in Set 1, and numerically evaluate the constant C 2lvl Σ J , J = 2, in Theorem 4.4 and the constant CΣ 3 in Theorem 4.8, for the h-version of the two solvers, based on selecting m 1 = m 2 = m = 2p 2 , cf. Figure 5.2. Here, we observe that C 2lvl Σ 2 and CΣ 3 are roughly (asymptotically) constant, as the polynomial degree p increases; thereby, this implies thatC(p J ), J = 2, 3, respectively, is approximately O(1), as p increases. Next, we investigate the performance of the two-level and W-cycle multigrid schemes in terms of the convergence factor where N denotes the number of iterations required to attain convergence up to a (relative) tolerance of 10 −8 and r N and r 0 are the final and initial residual vectors, respectively. In Table 5.2, we report the iteration counts and the convergence factor (in parenthesis), needed to attain convergence of the h-version of the two-level (TL) method and W-cycle multigrid scheme (with 3 and 4 levels), as a function of the number of elements (given by the choice of different grid sets), and the number of smoothing steps (m 1 = m 2 = m). Here, we have fixed the polynomial approximation order on each level p j ≡ p = 1. We first observe that, although the agglomerated grids, in general, do not necessarily strictly satisfy  Table 5.1: Value of the coercivity constant C coerc for the sets of grids considered in Figure 5.1 with C j σ = C σ = 10, j = 1, . . .  significantly increase with the number of elements in the underlying mesh; moreover, for the Wcycle solver, the number of iterations remains bounded with the number of levels. As expected, the convergence is faster for larger values of m and the solvers are convergent provided the number of smoothing steps is sufficiently large. For each grid, we have also reported the iteration counts N CG iter for the Conjugate Gradient (CG) method, which shows that the two proposed solvers outperform the CG scheme in terms of the number of iterations required to attain convergence, even when a small number of smoothing steps are employed. Table 5.3 presents analogous results for the first three sets of meshes, in the case when p = 3. Here, we observe that, as expected, the convergence factor increases, but the increase in p does not require an increase in the minimal number of smoothing steps needed to ensure that the underlying multilevel solvers are convergent.
A more exhaustive investigation of the effect of increasing p is reported in Table 5.4, where we consider a fine grid of 1024 elements and the corresponding agglomerated meshes (Set 2 in Figure 5.1). We observe that, even though both multilevel solvers converge for a fixed value of m, with increasing p, the number of iterations required to attain convergence increases as p grows. However, the twolevel and W-cycle multigrid solvers still employ less iterations, than the number required by the CG method. In Table 5.5, we report the number of iterations of the hp-version of the two solvers as a function of the number of smoothing steps and the number of levels for varying p. Here, p = p J denotes the polynomial approximation degree on the finest level V J , while, because of the   Table 5.3: Iteration counts of the h-version of the two-level and W-cycle solvers and the CG method as a function of m and the number of levels (C j σ ≡ C σ = 10, p = 3).
hp-approach, the polynomial order is decreased from the finest level to the coarser ones in such a way that p j−1 = p j − 1. We observe that the introduction of the hp-multigrid is detrimental for the convergence of the method, since the minimum number of smoothing steps needed to obtain a convergent method increases with respect to the h-version. Then, we can conclude that, as the pversion of the method does not exhibit uniform convergence with respect to the polynomial order, the hp-approach turns out to be not very effective for problems resulting from high-order discretizations.
As a second numerical test, we consider the same problem discretized by the interior penalty DG method of Section 2 on an initial mesh of triangles, thus reproducing a more common scenario in the  Table 5.4: Iteration counts of the h-version of the two-level and W-cycle solvers and the CG method as a function of p and the number of levels (C j σ ≡ C σ = 10, m = 5).  Table 5.5: Iteration counts of the hp-version of the two-level and W-cycle solvers and the CG method as a function of m and the number of levels (C j σ ≡ C σ = 10).
framework of finite element discretizations. In analogy to the case of polygonal meshes, we consider four sets of nested grids obtained by agglomeration, cf. Figure 5.3. The sets considered derive from initial meshes of 528 (Triangle Set 1), 1086 (Triangle Set 2), 2198 (Triangle Set 3) and 4318 (Triangle Set 4) elements. In Table 5.6, we show the iteration counts needed to attain convergence with respect to a fixed tolerance of 10 −8 as a function of the set (i.e., the number of elements) and the number of smoothing steps of the h-version of the two-level and W-cycle multigrid solvers, with p j = p = 1. We recall that, as in the previous numerical test, we have considered C j σ ≡ C σ = 10, for each j. The results are similar to the case of initial polygonal meshes, with uniform convergence with respect to the granularity of the mesh, and in the case of the W-cycle solver, also with respect to the number of levels. We again attain improved performance, compared to the standard CG method, in terms of the number of iterations required to attain convergence.  Table 5.6: Iteration counts and converge factor (in parenthesis) of the h-version of the two-level and W-cycle solvers and iteration counts of the CG method as a function of m (C j σ ≡ C σ = 10, p = 1). Starting mesh of triangles.