Andrea Cangiani

Emmanuil Zhaonan Dong
email: zhaonan.dong@inria.fr

Emmanuil H Georgoulis
email: emmanuil.georgoulis@le.ac.uk

hp-VERSION DISCONTINUOUS GALERKIN METHODS ON ESSENTIALLY ARBITRARILY-SHAPED ELEMENTS

Introduction

Recent years have witnessed a coordinated effort to generalize mesh concepts in the context of Galerkin/finite element methods. A key argument has been that more general-shaped elements/cells can potentially lead to computational complexity reduction. This effort has given rise to a number of recent approaches: mimetic finite difference methods [START_REF] Beirão Da Veiga | The mimetic finite difference method for elliptic problems[END_REF], virtual element methods [START_REF] Beirão Da Veiga | Basic principles of virtual element methods[END_REF][START_REF] Brenner | Virtual element methods on meshes with small edges or faces[END_REF], various discontinuous Galerkin approaches such as interior penalty [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], hybridized DG [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF] and the related hybrid high-order methods [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF]. Earlier approaches involving non-polynomial approximation spaces, such as polygonal and other generalized finite element methods [START_REF] Sukumar | Conforming polygonal finite elements[END_REF][START_REF] Fries | The extended/generalized finite element method: an overview of the method and its applications[END_REF], have also been developed and used by the engineering community. All the above numerical frameworks allow for polygonal/polyhedral element shapes (henceforth, collectively termed as polytopic) of varying levels of generality.

Simultaneously, various classes of fitted and unfitted grid methods for interface or transmission problems exploit generalized concepts of mesh elements in an effort to provide accurate representations of internal interfaces. Several unfitted finite element methods have been proposed in recent years: unfitted finite element methods [START_REF] Barrett | Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces[END_REF], immersed finite element methods [START_REF] Guo | A group of immersed finite-element spaces for elliptic interface problems[END_REF][START_REF] Gong | Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions[END_REF], virtual element methods [START_REF] Chen | An interface-fitted mesh generator and virtual element methods for elliptic interface problems[END_REF], unfitted penalty methods [START_REF] Bastian | An unfitted finite element method using discontinuous Galerkin[END_REF][START_REF] Massjung | An unfitted discontinuous Galerkin method applied to elliptic interface problems[END_REF][START_REF] Huang | An unfitted interface penalty finite element method for elliptic interface problems[END_REF][START_REF] Wu | An unfitted hp-interface penalty finite element method for elliptic interface problems[END_REF][START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF], see also [START_REF] Johansson | A high order discontinuous Galerkin Nitsche method for elliptic problems with fictitious boundary[END_REF] for unfitted discretization of the boundary, cutCell/cutFEM [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method[END_REF][START_REF] Burman | CutFEM: discretizing geometry and partial differential equations[END_REF], and unfitted hybrid high-order methods [START_REF] Burman | An unfitted hybrid high-order method for elliptic interface problems[END_REF], to name just the few closer to the developments we shall be concerned with below. A central idea in the majority of these methods is the weak imposition of interface conditions in conjunction with some form of penalization, see, e.g., [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF], an idea going back to [START_REF] Babuška | The finite element method for elliptic equations with discontinuous coefficients[END_REF]. These approaches are often combined with level set concepts [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] to describe the interfaces accurately. Nonetheless, in their practical implementation, the interface is typically represented via piecewise smooth polynomial approximations to the level sets.

The interior penalty discontinuous Galerkin (IP-dG) approach appears to allow for extreme generality with regard to element shapes/geometries. Indeed, in contrast to aforementioned families of general mesh methods, IP-dG can handle arbitrary number of faces per element with solid theoretical backing involving provable stability and convergence results; see [START_REF] Cangiani | hp-version space-time discontinuous Galerkin methods for parabolic problems on prismatic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] for details. This property becomes extremely relevant upon realising that IP-dG (as well as other classical dG methods, such as LDG) associate local numerical degrees of freedom to the elements only, and not to other geometrical entities such as faces or vertices. As such, the nature and dimension of the local discretization space is independent of the number of vertices/faces. The latter observation implies also naturally a form of complexity reduction: classical total degree ('P -type') local polynomial spaces in physical coordinates are admissible on box-type or highly complex element shapes [START_REF] Cangiani | hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF][START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. We refer to our monograph [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] for details on the admissible polygonal/polyhedral element shapes for which the IP-dG method is, provably, both stable and convergent. The mild element shape assumptions in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] are such to ensure the validity of crucial generalizations of standard approximation results, such as inverse estimates, best approximation estimates, and extension theorems. Thus, the developments presented for IP-dG in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] can be potentially ported also to other classical dG approaches within the unified framework of [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]; we also note the recent static condensation approach presented in [START_REF] Lozinski | A primal discontinuous Galerkin method with static condensation on very general meshes[END_REF] in this context.

The question, therefore, of further extending rigorously the applicability of hpversion IP-dG methods to meshes consisting of curved polygonal/polyhedral elements arises naturally. Indeed, such a development is expected to provide multifaceted advantages compared to current approaches, including, but not restricted to, the treatment of curved interfaces as done, e.g., in [START_REF] Massjung | An unfitted discontinuous Galerkin method applied to elliptic interface problems[END_REF][START_REF] Huang | An unfitted interface penalty finite element method for elliptic interface problems[END_REF][START_REF] Wu | An unfitted hp-interface penalty finite element method for elliptic interface problems[END_REF][START_REF] Burman | An unfitted hybrid high-order method for elliptic interface problems[END_REF][START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF]. For instance, allowing for extremely general curved elements enables the exact representation of curved computational domains, e.g., arising directly from Computer Aided Design programs. Allowing also for arbitrary local polynomial degrees, provides the possibility of achieving required accuracy via local (polynomial) basis enrichment (p-version Galerkin approaches) without increased mesh-granularity. If, nonetheless, local mesh refinement is also required/desired, IP-dG methods can be immediately applied on refined curved elements without local re-parametrizations of the local Galerkin spaces. This is in contrast to the need to perform costly re-parametrisations upon mesh refinement in other approaches, e.g., Isogeometric Analysis [START_REF] Hennig | Adaptive mesh refinement strategies in isogeometric analysis-a computational comparison[END_REF] or, indeed, even to keep track of the domain-approximation variational crimes of standard finite element discretizations. Exact geometry representation can also be highly relevant in representing locally discontinuous/sharply changing PDE coefficients, e.g., in permeability pressure computations in porous media, coefficients defined via level-sets of smooth functions, or shape/topology optimization applications.

Furthermore, exact geometry representation is relevant in the p-version Galerkin context: to achieve spectral/exponential convergence for smooth PDE problems posed on general curved domains, we are required to use isoparametrically mapped elements. This is both cumbersome to implement and costly as the polynomial degree increases [START_REF] Murti | Numerical inverse isoparametric mapping in remeshing and nodal quantity contouring[END_REF][START_REF] Murti | Numerical inverse isoparametric mapping in 3d fem[END_REF]. A successful alternative to isoparametric maps is the use of non-linear maps on element patches [START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF][START_REF] Melenk | hp-finite element methods for singular perturbations[END_REF] to represent domain geometry. Nevertheless, if the elemental maps are not a priori provided, it is difficult to construct them in practice, especially in three dimensions.

Finally, curved element capabilities should ideally be developed in conjunction with the already developed highly general polytopic mesh IP-dG methods, allowing for instance elements with arbitrary number of faces. This is particularly pertinent in the contexts of adaptivity and multilevel solvers, which benefit from element agglomeration [START_REF] Antonietti | Multigrid algorithms for hp-version interior penalty discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF][START_REF] Antonietti | An agglomeration-based massively parallel non-overlapping additive schwarz preconditioner for high-order discontinuous galerkin methods on polytopic grids[END_REF] to achieve coarser representations. With regard to adaptivity, mesh coarsening is essential in keeping the computation sizes at bay, at least in the case of evolution problems. The extreme coarsening capabilities via element agglomeration, therefore, have the potential in retaining structure, e.g., possible coefficient heterogeneities at the discrete level for instance.

It is, therefore, desirable to design and analyze IP-dG and related methods posed on meshes comprising of elements with arbitrary number of curved faces, under as mild geometric assumptions as possible. To address this central, in our view, question, this work aims at rigorously extending the applicability of IP-dG methods on meshes comprising of essentially arbitrarily curve-shaped polytopic elements with arbitrary number of faces per element; this includes, in particular, curved elements not exactly representable by (iso-)parametric polynomial element mappings.

The theoretical developments presented below regarding stability and a-priori error analysis of IP-dG methods hinge on new, to the best of our knowledge, extensions of known inverse and trace inequalities. More specifically, we extend the hp-version trace inverse estimate presented in [START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF], allowing for more general curved element shapes; see also [START_REF] Huang | An unfitted interface penalty finite element method for elliptic interface problems[END_REF] for an earlier, related result. Trace inverse estimates are crucial in the proof of stability of IP-dG methods and, simultaneously, determine the so-called discontinuity-penalization parameter for a given mesh. This is crucial on meshes of such generality: insufficient penalization results in loss of stability, while excessive penalization typically results in accuracy loss. Also, we prove new hp-version L ∞ -L 2 and H 1 -L 2 inverse inequalities on extremely general curved domains. Particular care has been given so that these new inverse estimates are 'shape-robust', in the sense that there is no hidden dependence of the element shape in the constants. We believe that these extensions of known inverse estimates to be of independent interest, due to their frequent use in the analysis of finite element methods.

The new inverse estimates are combined with ideas from the analysis of polytopic dG methods [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], resulting in significant generalization of the results presented therein. More specifically, by relaxing certain earlier coverability assumptions, (postulating the ability to cover tightly general-shaped elements by unions of simplices of similar size, cf. [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Definition 10]) as well as by proving a new stability result for norms of polynomials under domain perturbations (Lemma 4.14), we prove stability and a new hp-version a priori error analysis for the IP-dG method on essentially arbitrary element shapes. The a priori error analysis follows closely the proof from [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF]: upon establishing an inf-sup stability result of the method in a streamlinediffusion-like norm, standard Strang-type arguments with hp-best approximation results lead to an error bound. The inf-sup result justifies also the good stability properties of the method in convection-dominated problems. The theoretical tools presented may also be of interest in Nitsche-type formulations of unfitted grid interface methods. To emphasize the mesh-generality of the proposed approach, we shall refer to the framework presented below as discontinuous Galerkin method on essentially arbitrarily-shaped elements (dG-EASE).

The remainder of this work is organised as follows. Upon describing the advectiondiffusion-reaction model problem in Section 2, we introduce the hp-version interior penalty discontinuous Galerkin method in Section 3. We prove new inverse estimates in Section 4, along with the necessary hp-approximation results. In Section 5, we present the stability and a-priori error analysis. Finally, the performance of the dG methods is assessed in practice through a series of numerical experiments presented in Section 6.

Model problem

To highlight the versatility of dG-EASE, we consider the class of second-order partial differential equations with nonnegative characteristic form over an open bounded Lipschitz domain Ω in R d , d ≥ 1, with boundary ∂Ω. This class includes general advection-diffusion-reaction problems possibly of changing type, see, e.g., [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. The model problem reads: find u ∈ V such that

(2.1) -∇ • (a∇u) + ∇ • (bu) + cu =f in Ω,
for some suitable solution space V, and a = {a ij } d i,j=1 , symmetric with a ij ∈ L ∞ (Ω), so that at each x in Ω, we have

(2.2) d i,j=1 a ij (x)ξ i ξ j ≥ 0, for any ξ = (ξ 1 , . . . , ξ d ) T ∈ R d ; also b = (b 1 , . . . , b d ) T ∈ W 1 ∞ (Ω) d , c ∈ L ∞ (Ω) and f ∈ L 2 (Ω).
To supplement (2.1) with suitable boundary conditions, following [START_REF] Oleinik | Second Order Equations with Nonnegative Characteristic Form[END_REF], we first subdivide the boundary ∂Ω into ∂ 0 Ω = x ∈ ∂Ω : d i,j=1 a ij (x)n i n j > 0 , and ∂Ω\∂ 0 Ω with n = (n 1 , . . . , n d ) T denoting the unit outward normal vector to ∂Ω. Loosely speaking, we may think of ∂ 0 Ω as being the 'elliptic' portion of the boundary ∂Ω. We further split the 'hyperbolic' portion of the boundary ∂Ω\∂ 0 Ω, into inflow and outflow boundaries ∂ -Ω and ∂ + Ω, respectively, by

∂ -Ω = {x ∈ ∂Ω\∂ 0 Ω : b(x) • n(x) < 0} , ∂ + Ω = {x ∈ ∂Ω\∂ 0 Ω : b(x) • n(x) ≥ 0} .
If ∂ 0 Ω is nonempty, we shall further divide it into two disjoint subsets ∂Ω D and ∂Ω N , with ∂Ω D nonempty and relatively open in ∂Ω. It is evident from these definitions that

∂Ω = ∂Ω D ∪ ∂Ω N ∪ ∂ -Ω ∪ ∂ + Ω.
It is physically reasonable to assume that b • n ≥ 0 on ∂Ω N , whenever ∂Ω N is nonempty; then, we impose the boundary conditions:

(2.3) u = g D on ∂Ω D ∪ ∂ -Ω, n • (a∇u) = g N on ∂Ω N ;
For an extension, allowing also for b • n < 0 on ∂Ω N , we refer to [START_REF] Cangiani | Discontinuous Galerkin methods for mass transfer through semipermeable membranes[END_REF]. Additionally, we assume that there exists a positive constant γ 0 such that

(2.4) c 0 (x) := c(x) + 1 2 ∇ • b(x) 1/2 ≥ γ 0 a.e. x ∈ Ω.
For a proof of the well-posedness of (2.1), (2.3), subject to (2.4), we refer to [START_REF] Oleinik | Second Order Equations with Nonnegative Characteristic Form[END_REF][START_REF] Houston | Stabilised hp-finite element approximation of partial differential equations with nonnegative characteristic form[END_REF].

Discontinuous Galerkin method

We shall now define the interior penalty discontinuous Galerkin (dG) method posed on essentially arbitrarily-shaped elements. A key attribute of the method is the use of physical frame basis functions, i.e., the elemental bases consist of polynomials on the elements themselves, rather than mapped from a reference element. The implementation challenges arising from this non-standard choice with regard to construction of the resulting linear system will be discussed below.

3.1. The mesh. Let T = {K} be a subdivision of Ω into non-overlapping subsets (elements) K ∈ T with, possibly curved, Lipschitz boundaries and let h K := diam(K). The mesh skeleton Γ := ∪ K∈T ∂K is subdivided into the internal part Γ int := Γ\∂Ω and boundary part ∂Ω. We further explicitly assume that the (d -1)dimensional Hausdorff measure of Γ is globally finite, thereby, not allowing for fractal-shaped elements.

We note immediately that we allow mesh elements K ∈ T which are essentially arbitrarily-shaped and with very general interfaces with neighbouring elements.

For instance, two elements may interface at a collection of (d -1)-dimensional (possibly curved) faces, as those shown in Figure 3.1. The precise assumptions on the admissible element shapes are given in Section 4 below.

K K K • • • • • • • • • • • • • • • • • • Figure 1. Curved elements K, K , K for d = 2
with possibly many curved faces; • denotes a vertex.

3.2. Discontinuous Galerkin method. We define the hp-version discontinuous finite element space S p T , subordinate to the mesh T = {K} and a polynomial degree vector p := {p K }, possibly different for each element K, by

(3.1) S p T = {v ∈ L 2 (Ω) : v| K ∈ P p K (K), K ∈ T }.
For any elemental face F ⊂ Γ int , let K and K be the two elements such that F ⊂ ∂K ∩ ∂K . The outward unit normal vectors on F of ∂K and ∂K are denoted by n K and n K , respectively. For a function v : Ω → R that may be discontinuous across Γ, we define the jump v and the average {v} of v across F by

(3.2) v = v| K n K + v| K n K , {v} = 1 2 (v| K + v| K ) .
Similarly, for a vector valued function w, piecewise smooth on T , we define

w = w| K • n K + w| K • n K , {w} = 1 2 (w| K + w| K ) .
When F ⊂ ∂Ω, we set {v} = v, v = vn and w = w • n with n denoting the outward unit normal to the boundary ∂Ω.

For any element K ∈ T , we define inflow and outflow parts of ∂K by

∂ -K = {x ∈ ∂K : b(x) • n K (x) < 0}, ∂ + K = {x ∈ ∂K : b(x) • n K (x) ≥ 0},
respectively, with n K (x) denoting the unit outward normal vector to ∂K at x ∈ ∂K. Further, we define the upwind jump of the (scalar-valued) function v across the inflow boundary

∂ -K of K ∈ T by v (x) := lim →0 + v(x + b(x) ) -v(x -b(x) ) , when x ∈ ∂ -K\∂Ω. Finally, we define the broken gradient ∇ T v of a function v ∈ L 2 (Ω) with v| K ∈ H 1 (K), for all K ∈ T , element-wise by ∇ T v | K := ∇(v| K ).
The discontinuous Galerkin method on essentially arbitrarily-shaped elements (dG-EASE for short) reads: find

u h ∈ S p T such that (3.3) B(u h , v h ) = (v h ) for all v h ∈ S p T , with the bilinear form B(•, •) : S p T × S p T → R defined as B(u, v) := B ar (u, v) + B d (u, v),
where B ar (•, •) accounts for the advection and reaction terms:

(3.4) B ar (u, v) := Ω ∇ T (b • u) + cu v dx - K∈T ∂-K\∂Ω (b • n) u v ds - K∈T ∂-K∩(∂ΩD∪∂-Ω) (b • n)uv ds,
B d (u, v) := Ω a∇ T u • ∇ T v dx + Γint∪∂ΩD σ[[u]] • [[v]] ds - Γint∪∂ΩD { {a∇u} } • [[v]] + { {a∇v} } • [[u]] ds,
while the linear functional : S p T → R is defined by

(3.6) (v) := Ω f v dx - K∈T ∂-K∩(∂ΩD∪∂-Ω) (b • n)g D v ds - ∂ΩD g D (a∇v) • n -σv ds + ∂ΩN g N v ds.
The nonnegative function σ ∈ L ∞ (Γ int ∪ ∂Ω D ) appearing in (3.5) and (3.6) is the discontinuity-penalization function, whose precise definition, which depends on the diffusion tensor a and the discretization parameters, will be given below. We note that a 'good' choice of discontinuity penalization is instrumental for the stability of the method, while simultaneously not affecting the approximation properties in the general mesh setting considered herein.

For simplicity of presentation, we shall assume that the entries of the diffusion tensor a are element-wise constants on each element K ∈ T , i.e., (3.7) a ∈ [S 0 T ] d×d sym , Our results can be applied to the case of general a ∈ [H 1/2 (Ω)] d×d sym by slightly modifying the bilinear form above as proposed originally in [START_REF] Georgoulis | A note on the design of hp-version interior penalty discontinuous Galerkin finite element methods for degenerate problems[END_REF] and extended to polytopic meshes in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. In the following, √ a denotes the (positive semidefinite) square-root of the symmetric matrix a; further, āK := | √ a| 2 2 | K , where | • | 2 denotes the matrix-2-norm. Also, in the interest of accessibility, we shall not consider problems with high contrast diffusion tensors, with the usual weighted averaging modification of the method [START_REF] Burman | A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems[END_REF][START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Di Pietro | Discontinuous Galerkin methods for anisotropic semidefinite diffusion with advection[END_REF]; the extension to that setting is completely analogous to the analysis presented below.

Remark 3.1. The parameter σ is typically selected to be face-wise constant in the definition and implementation of IP-dG methods. To ensure that only physically correct penalization takes place, σ is chosen below to be proportional to the quantity n T an; see [START_REF] Georgoulis | A note on the design of hp-version interior penalty discontinuous Galerkin finite element methods for degenerate problems[END_REF] for details. As such, σ will vary along a curved element face even for element-wise constant diffusion a, thereby justifying the terminology "penalization function" as opposed to the standard terminology "penalization function" from the literature. Further, the theory presented below can also be extended with minor modifications to curved faces F , such that n T an > 0 only on a strict subset of that face and n T an = 0 on the remaining part. That way one can reduce or even remove unphysical penalization on the hypersurfaces where the PDE may change type.

Inverse and approximation estimates

A key challenge in the error analysis presented below is the availability of inverse estimation and approximation results with uniform /explicit constants with respect to the shape of the elements in a given mesh.

A trace type inverse estimate for elements with one curved face has been recently proven in [START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF] under a shape-regularity assumption; see (4.1) below. Results in this direction have also appeared under various geometric assumptions in [START_REF] Wu | An unfitted hp-interface penalty finite element method for elliptic interface problems[END_REF][START_REF] Massjung | An unfitted discontinuous Galerkin method applied to elliptic interface problems[END_REF][START_REF] Burman | An unfitted hybrid high-order method for elliptic interface problems[END_REF], among others. Here, we extend these results by proving trace-inverse estimates for elements that are locally star-shaped, Lipschitz domains (see Assumption 4.1 below). Moreover, given the importance of trace-inverse estimates for the stability of interior penalty dG methods, the new estimate constant is expressed via explicit and practically verifiable, geometric quantities (Lemma 4.4 below).

In the same vein, we also extend the classical (Markov-type) H 1 -L 2 inverse estimate to elements with piecewise C 1 , locally star-shaped boundaries (see Assumptions 4.1 and 4.3 below). The proof builds upon and extends on earlier ideas from [START_REF] Kroó | On Bernstein-Markov-type inequalities for multivariate polynomials in Lq-norm[END_REF] and [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. Here, we are particularly concerned with explicit quantification of the respective constant for a given element geometry. We note that H 1 -L 2 inverse estimates are also relevant in the determination of penalty parameters in IP-dG methods for biharmonic operators [START_REF] Dong | Discontinuous galerkin methods for the biharmonic problem on polygonal and polyhedral meshes[END_REF].

Also, we revisit a key stability argument that enabled the use of 'degenerate' polytopic element shapes, i.e., ones containing very small/degenerating faces/edges compared to the element diameter, first proposed in [START_REF] Cangiani | hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]; see also [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] for improvements. This result is crucial in offering a practical choice of the discontinuitypenalization parameter for general polytopic meshes. The stability argument was based on two ingredients: 1) control of integral norms of polynomials with respect to domain perturbations using [START_REF]Inverse-type estimates on hp-finite element spaces and applications[END_REF]Lemma 3.7], and 2) an L ∞ -L 2 inverse estimate. To retain this capability in the current setting, we prove an extension of [START_REF]Inverse-type estimates on hp-finite element spaces and applications[END_REF]Lemma 3.7] (see also [START_REF] Kroó | On Bernstein-Markov-type inequalities for multivariate polynomials in Lq-norm[END_REF]Lemma 6] for a related result) for generalized/curved prismatic elements; see Lemma 4.14 below. Moreover, we also prove an extension of the classical L ∞ -L 2 inverse estimate for generalized/curved prismatic elements. The latter two new estimates, in conjunction with a revised concept of coverability (compared to [START_REF] Cangiani | hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]) are enough to provide extensions to previously known stability results for IP-dG within the present level of mesh generality. Assumption 4.1. For each element K ∈ T , we assume that K is a Lipschitz domain, and that we can subdivide ∂K into mutually exclusive subsets {F i } n K i=1 satisfying the following property: there exist respective sub-elements

• • • K Fi F i x 0 i • • • n K Fi F i m i x 0 i
K Fi ≡ K Fi (x 0 i ) ⊂ K with d planar faces meeting at one vertex x 0 i ∈ K, with F i ⊂ ∂K Fi , such that, for i = 1, . . . , n K , (a) K Fi is star-shaped with respect to x 0 i . We refer to Figure 2(left) for an illustration for d = 2; (b) m i (x) • n(x) > 0 for m i (x) := x -x 0 i , x ∈ K Fi ,
and n(x) the respective unit outward normal vector to F i at x ∈ F i . (We refer to Figure 2(right) for an illustration for d = 2.) Remark 4.2. Some remarks on the above (very mild) mesh assumption are in order:

(i) The sub-domains {F i } n K i=1 are not required to coincide with the faces of the element K: each F i may be part of a face or may include one or more faces of K. Also, there is no requirement for {n K } K∈T to remain uniformly bounded across the mesh. (ii) We can make Assumption 4.1(b) stronger by further postulating that: it is possible to fix the point x 0 i such that there exists a global constant c sh > 0,

K K K • • • • • • • • • • • • • • • • • • Figure 3. Curved elements K, K
m i (x) • n(x) ≥ c sh h K F i ;
this is the case, of course, for straight-faced polytopic elements, cf., [START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF][START_REF] Wu | An unfitted hp-interface penalty finite element method for elliptic interface problems[END_REF]. Note that (4.1) does not imply shape-regularity of the K Fi 's; in particular K Fi 's with 'small' F i compared to the remaining (straight) faces of K Fi are acceptable. Such anisotropic sub-elements K Fi 's may be necessary to ensure that each K Fi remains star-shaped when an element boundary's curvature is locally large; see, e.g., K Fi in Figure 2 and a collection of both 'shape-regular' and 'anisotropic' K Fi 's in Figure 3. (iii) On certain geometrically extreme cases, satisfying Assumption 4.1 may require a small number of refinements of the elements K ∈ T of a given initial mesh. (iv) F i is not required to be connected. However, by splitting F i to its connected subsets, re-indexing the F i 's to correspond to unique K Fi , we can allocate one K Fi to each F i ; we shall take the latter point of view in what follows to avoid further notational complexity. Lemma 4.4. Let element K ∈ T be a Lipschitz domain satisfying Assumption 4.1. Then, for each F i ⊂ ∂K, i = 1, . . . , n K , and for each v ∈ P p (K), we have the inverse estimate:

(4.2) v 2 Fi ≤ (p + 1)(p + d) min x∈Fi (m i • n) v 2 K F i .
Proof. We partition F i into r (d -1)-dimensional curved simplices denoted by F j i , j = 1, . . . , r, which are subordinate to the vertices possibly contained in F i ; r is large enough to accommodate this requirement. Further, we construct a partition of K Fi into (curved) sub-elements K j i , by considering the simplices with one (curved) face F j i and the remaining vertex being x 0 i ; this is possible due to the star-shapedness of K Fi with respect to x 0 i as per Assumption 4.1(a). We refer to Figure 4 for an illustration when d = 2. Notice that each F j i may include at most one constituent curved face of F i , or part thereof.

• • • K Fi K j i F j i x 0 i K j i F j i F j i F j i x 0 i Figure 4. Partitioned curved sub-element K Fi ⊂ K ∈ T ; • de- notes a vertex of K (left)
; detail with K j i and related faces F j i and F j i (right).

Let now F j i denote the straight/planar related face defined by the d -1 vertices of F j i . Let also K j i be the largest straight-faced simplex contained in K j i with face F j i parallel to F j i and the remaining faces being subsets of the straight faces of K j i . The Divergence Theorem implies

K j i \K j i ∇ • (v 2 m i ) dx = ∂(K j i \K j i ) v 2 m i • n ∂(K j i \K j i ) ds = F j i v 2 m i • n ds + F j i v 2 m i • n F j i ds,
with n ω denoting the outward normal vector of a domain ω ⊂ R d and m i as in Assumption 4.1(b), upon observing that

m i •n ∂(K j i \K j i ) = 0 on ∂(K j i \K j i )\(F j i ∪F j i ). Now, denoting by |•| 2 the Euclidean distance in R d ,
the product rule and elementary estimates imply

K j i \K j i ∇ • (v 2 m i ) dx ≤ 2 max K j i |m i | 2 v∇v L∞(K j i \K j i ) + d v 2 L∞(K j i \K j i ) |K j i \K j i |, noting that ∇ • m i = d.
The right-hand side of the above inequality converges to zero as |K j i \K j i | → 0, which, in turn, is achieved as r → ∞. Thus, Assumption 4.1(b) gives min

x∈F j i (m i • n) v 2 F j i ≤ F j i v 2 m i • n ds ≤ F j i v 2 m i • n F j i ds + , for some = O(|K j i \K j i |) as r → ∞.
Each of the finite F i 's is, in turn, image of a finite number of Lipschitz functions locally. Let L be the Lipschitz constant of a parametrisation of F j i with respect to F j i , giving

|F j i | ≤ L|F j i |. At the same time, we have |K j i \ K j i | ≤ Lh F j i |F j i |, as the maximum Euclidean distance between F j i and F j i is bounded from above by Lh F j i . Hence the area |K j i \ K j i | converges to zero faster than |F j i | by an order of h F j i .
At the same time, a standard trace-inverse estimate on simplices, [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF], yields

F j i v 2 m i • n F j i ds ≤ max x∈F j i (m i • n F j i ) (p + 1)(p + d)|F j i | d|K j i | v 2 K j i .
Combining the above, we have that, for any δ > 0, there exists an r large enough such that

v 2 F j i ≤ max x∈F j i (m i • n F j i ) min x∈F j i (m i • n) (p + 1)(p + d)|F j i | d|K j i | v 2 K j i + min x∈F j i (m i • n) ≤ (1 + δ) (p + 1)(p + d)|F j i | d|K j i | v 2 K j i ≤ (1 + δ) (p + 1)(p + d)|F j i | d|K j i | v 2 K j i ,
as the first ratio on the first estimate tends to 1 as r → ∞ . In the last inequality we used the bound

|F j i | ≤ |F j i | and that K j i ⊂ K j i .
Another application of the Divergence Theorem and elementary calculations give

d|K j i | = K j i ∇ • m i dx = F j i m i • n ds ≥ min x∈F j i (m i • n)|F j i |,
or

d|K j i | + d|K j i \K j i | ≥ min x∈F j i (m i • n)|F j i |,
or

|F j i | |K j i | ≤ d min x∈F j i (m i • n) 1 + |K j i \K j i | |K j i | ≤ (1 + δ)d min x∈F j i (m i • n) ,
for any δ > 0 when r is sufficiently large. Combining the above, we deduce

v 2 Fi = r j=1 v 2 F j i ≤ r j=1 (1 + δ) 2 (p + 1)(p + d) min x∈F j i (m i • n) v 2 K j i ≤ (1 + δ) 2 (p + 1)(p + d) min x∈F K i (m i • n) v 2 Ki .
Taking, finally, r → ∞, allows for δ → 0 and the result (4.2) follows.

Remark 4.5. It is important to stress that the right-hand side of (4.2) is a function of x 0 i defining K Fi . Since the closure of the original (curved) element K is compact in R d , it is possible to minimise the right-hand side of (4.2) by selecting an 'optimal' x 0 i . Moreover, upon making the stronger assumption (4.1), we arrive at the familiar trace inverse estimate for star-shaped, shape-regular elements with piecewise smooth boundaries:

v 2 ∂K ≤ Cp 2 h -1 K v 2 K .
Example 4.6. Let K = B(0, R) ⊂ R d be the ball of radius R centred at the origin. Then, selecting F 1 = ∂K =: S(0, r), we have

v 2 S(0,R) ≤ (p+1)(p+d)R -1 v 2 B(0,R) .
Within this geometric setting, we can specify the constants of the classical trace inequality for H 1 -functions. The result below is a mild extension of [START_REF] Cangiani | Adaptive discontinuous Galerkin methods for elliptic interface problems[END_REF]Lemma 4.1], (cf. also [START_REF] Wu | An unfitted hp-interface penalty finite element method for elliptic interface problems[END_REF]) following closely the classical proof from [START_REF] Agmon | Lectures on elliptic boundary value problems[END_REF]. Lemma 4.7. Let K ∈ T be a Lipschitz domain satisfying Assumption 4.1. Then, for all ζ > 0, we have the estimate

(4.3) v 2 Fi ≤ d + ζ min x∈Fi (m i • n) v 2 K F i + max x∈Fi |m i | 2 2 ζ min x∈Fi (m i • n) ∇v 2 K F i , for all v ∈ H 1 (Ω) and i = 1, . . . , n K .
Proof. The Divergence Theorem and the fact that

m i • n = 0 on ∂K Fi \ F i imply Fi v 2 m i • n ds = K F i ∇ • (v 2 m i ) dx ≤ d v 2 K F i + 2 max x∈Fi |m i | 2 v K F i ∇v K F i ,
from which the result already follows.

Remark 4.8. Summing over i = 1, . . . , n K , assuming (4.1) and that h

K F i ∼ h K , (4.3) gives the classical trace estimate v 2 ∂K ≤ C h -1 K v 2 K + h K ∇v 2 K . F F 0 ρ Figure 5. A reference generalized prism K for d = 2. 4.2. Basic H 1 -L 2 inverse estimate. H 1 -L 2 inverse
estimates for polynomials on d-dimensional simplicial or box-like domains are proven via directional arguments, if explicit dependence on the polynomial degree is desired, see, e.g., [START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF].

Generalizations of these estimates on convex domains use an analogous method of proof [START_REF] Kroó | On Bernstein-Markov-type inequalities for multivariate polynomials in Lq-norm[END_REF]. Here, in the same spirit, we extend further the domain generality in H 1 -L 2 inverse estimates, by also employing directional arguments on curved prismatic subdomains; the general case then follows by covering general Lipschitz domains by these curved prisms.

Definition 4.9.

Let F 0 := [0, 1] d-1 ⊂ R d and φ : F 0 → R a Lipschitz continuous scalar function. A reference generalized prism is a domain K ≡ Kφ ⊂ R d given by K ≡ Kφ := {x ∈ R d : 0 ≤ x i ≤ 1, i = 1, . . . , d -1, 0 ≤ x d ≤ φ(x 1 , . . . , x d-1 )},
with the properties: 1) [0, 1] d ⊂ K, and 2) the straight line connecting any pair (x, y) ∈ F 0 × F lies fully in K.

F := {x ∈ R d : 0 ≤ x i ≤ 1, i = 1, . . . , d -1, x d = φ(x 1 , . . . , x d-1 )}.
Also, we set ρ := sup{ρ ≥ 1 : F 0 × [0, ρ] ⊂ K} and r := max x∈ F 0 φ(x) + 1.

We refer to Figure 5 for an illustration.

Remark 4.10. A sufficient but, crucially, not necessary condition for K ≡ Kφ to be a reference generalized prism is that φ is a contraction. Since, however, Kφ will be used in conjunction with affine maps below, it will become possible to consider φ with Lipschitz constants greater than one.

Remark 4.11. The 'height' r is a measure of anisotropy of the reference generalized prism. Note that we can take ρ = 1 without essential loss of generality. Indeed, if ρ > 1, the change of variables x d → x d /ρ implies a modification of the Lipschitz function φ, reducing its Lipschitz constant. Star-shapedness with respect to F 0 is also ensured (cf., Remark 4.10).

In light of the above remark, we consider the case ρ = 1 only, in what follows.

Lemma 4.12. Let v ∈ P p ( K), p ∈ N, with K ⊂ R d a reference generalized prism. Then, we have the inverse estimate

(4.4) ∇v 2 K ≤ C B inv p 4 v 2 K , with C B inv ≡ C B inv (d, r) := 288(d -1)r 2 + 12d.
Proof. We begin by introducing some notation. Let F be a hyperplanar region in R d and let vector v ∈ R d . We define a zone Z(F, v) ⊂ R d , to be the geometric locus given by Z(F, v) := {z + βv : z ∈ F, β ∈ R}. (Thus, for instance, the domain [0, 1] × R = Z([0, 1], (0, α)) for any α ∈ R.) Using this notation, we now construct 2(d -1) suitable subsets { Kxj ± } d-1 j=1 , so that the union of K0 := [0, 1] d together with

{ Kxj ± } d-1 j=1 cover K. R x1 - R x1 + Kx1 - Kx1 + Figure 6. R = [0, 1] × [0,
r] and K, with r = 2, and their respective truncated prisms R x1 ± and Kx1 ± .

We first present the construction for d = 2 for accessibility. Set F x1 0,-:= [(2r) -1 , 1] and F x1 0,+ := [0, 1 -(2r) -1 ]. Then, elementary geometric arguments reveal that the rectangle R := [0, 1]×[0, r] can be covered by the union of R 0 := [0, 1] 2 and the 2 truncated prisms defined as:

R x1 ± := R ∩ Z F x1 0,± , v x1 ± , with v x1 ± := (±1/2, r) T |(1/2, r)| ,
with | • | denoting the standard Euclidean distance; we refer to Figure 6 for an illustration with r = 2. Correspondingly, for d = 3, R := [0, 1] 2 × [0, r] can be covered by R 0 := [0, 1] 3 together with 2 'x 1 -direction tilted', truncated prisms:

R x1 ± := R ∩ Z F x1 0,± , v x1 ± , with v x1 ± := (±1/2, 0, r) T |(1/2, 0, r)| , with F x1 0,-:= [(2r) -1 , 1] × [0, 1] and F x1 0,+ := [0, 1 -(2r) -1 ] × [0, 1]
the respective prism bases. At the same time, R can be also covered by R 0 together with the 2 'x 2 -direction tilted', truncated prisms:

R x2 ± := R ∩ Z F x2 0,± , v x2 ± , with v x2 ± := (0, ±1/2, r) |(0, 1/2, r) T | , with F x2 0,-:= [0, 1] × [(2r) -1 , 1] and F x2 0,+ := [0, 1] × [0, 1 -(2r) -1
] the respective prism bases. (We note the 'overloading' of notation with respect to dimension.) The construction for d ≥ 4 follows in a completely analogous fashion by considering 2 'x j -direction tilted', truncated prisms for each j = 1, . . . , d -1.

Since K ⊂ [0, 1] d-1 × [0, r], we consider the sets K0 = R 0 together with

Kxj ± := K ∩ R xj ± ,
for each fixed j = 1, . . . , d-1; see Figure 6 for an illustration for d = 2 and r = 2.

First, we observe the estimates

(4.5) |v xj | 2 ≤ 8r 2 |v x d | 2 + 2(4r 2 + 1)|v xj ± • ∇v| 2 , for j = 1, . . . , d -1.
Using the latter, we have, respectively, (4.6)

∇v 2 K -v x d 2 K - d-1 j=1 v xj 2 K0 ≤ d-1 j=1 v xj 2 Kx j + + v xj 2 Kx j - ≤ d-1 j=1 8r 2 v x d 2 K + 2(4r 2 + 1) v xj + • ∇v 2 Kx j + + v xj -• ∇v 2 Kx j - .
We now estimate each term on the right-hand side of (4.6). For e d := (0, . . . , 0, 1) T , let x := K ∩ {x + αe d : α ∈ R}, i.e., the vertical line contained in K and passing through a point x. Then, we have

(4.7) v x d 2 K = F 0 x v 2 x d dx d dx ≤ 12p 4 v 2 K ,
from Fubini's Theorem and an one-dimensional inverse estimate, see, e.g., [START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF]Theorem 3.91]. We set

xj x,± := K ∩ {x + αv xj ± : α ∈ R}. Then, v xj ± • ∇v 2 Kx j ± = F x j 0,± x j x,± (v xj ± • ∇v) 2 dζ dx ≤ 48r 2 p 4 4r 2 + 1 v 2 K ,
for j = 1, . . . , d, upon noticing that the length of xj x,± is bounded from below by r2 + 1/4/r (the length of the portion of xj x,± contained in R 0 ). Thus, (4.6) implies

∇v 2 K -v x d 2 K - d-1 j=1 v xj 2 K0 ≤ 288p 4 (d -1)r 2 v 2 K .
The result already follows by combining the last estimate with (4.7) and the corresponding inverse estimates for v xj 2 [0,1] d . Notice that (4.4) retrieves the known constant for d = 1. If K is cuspoidal, (i.e., not Lipschitz) (4.4) does not hold in general; we refer to [START_REF]Sharp Lp Markov type inequality for cuspidal domains in R d[END_REF] for a counterexample. Remark 4.13. Careful inspection of the above proof shows that, in fact, we have proven the sharper inverse estimate

∇v 2 K ≤ C B inv p 4 ρ-2 v 2 K .
In view of Remark 4.11, however, a linear scaling results into a modified K (and, possibly modified r) for which (4.4) is sharp.

4.3.

Stability with respect to domain perturbation. We now prove a stability result with respect to domain perturbation in the spirit of [START_REF]Inverse-type estimates on hp-finite element spaces and applications[END_REF]Lemma 3.7] (see also [START_REF] Kroó | On Bernstein-Markov-type inequalities for multivariate polynomials in Lq-norm[END_REF]Lemma 6]).

Lemma 4.14. Let K a reference generalized prism and consider its subset K := K ∩ ( K -e d ); here A + z := {x + z, x ∈ A}, for A ⊂ R d and z ∈ R d . Then, for all v ∈ P p ( K), p ∈ N, and for any 0 < ≤ (8p) -2 , we have the estimate

(4.8) 1 2 v 2 K ≤ v 2 K .
Proof. Set x, := K ∩ ( x -e d ). Then, we have, respectively, (4.9)

v 2 K\ K = F 0 x\ x, v 2 dx d dx ≤ F 0 v 2 L∞( x\ x, ) dx ≤ F 0 v 2 L∞( x) dx ≤ 32 p 2 v 2 K ,
by Markov's inequality (see, e.g., [START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF]Theorem 3.92]) since the length of x is bounded from below by one. Selecting now 0 < ≤ (8p) -2 , the result follows, by simply observing that

v 2 K -v 2 K = v 2 K\ K . Example 4.15. Consider K with φ(x) = 1 + x • (1 -x) + (16p) -2 sin(απ1 • x), for some α ≥ 1, with 1 := (1, 1, . . . , 1) T ∈ R d-1 . This is chosen so that φ(x) ∈ [1 + x • (1 -x) -, 1 + x • (1 -x) + ],
for within the range required for the statement of Lemma 4.14 to hold. For sufficiently large α, K is not star-shaped with respect to F 0 = [0, 1] d-1 . Nevertheless, φ is sufficiently approximated by φ 2 (x) := 1 + x • (1 -x), which is star-shaped with respect to F 0 = [0, 1] d-1 . Thus, for v ∈ P p ( K), from Lemma 4.12, we have

∇v 2 K ≤ 2 ∇v 2 Kφ 2 ≤ 4C B inv (d, 2)p 4 v 2 K .
4.4. L ∞ -L 2 inverse estimate. We continue by proving an L ∞ -L 2 inverse estimate for reference generalized prisms.

Lemma 4.16. Let K ≡ Kφ ⊂ R d a reference generalized prism. Then, the inverse estimate

(4.10) v 2 L∞( K) ≤ (32p 2 ) d v 2 K ,
holds for all v ∈ P p ( K), p ∈ N.

Proof. Let x max ∈ K such that v L∞( K) = |v(x max )|. Then, either x max ∈ K0 = [0, 1] d ⊂ K
)| 2 = v 2 L∞( K) = v 2 L∞( K(xmax))
. Since x max is a vertex, we employ an one-dimensional trace inverse estimate [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF] iteratively with respect to dimension, to deduce |v(x max )| ≤ d j=1 (p + j) v K(xmax) . Here we have used the fact that the dimensions of K(x max ) are grater than 1. Combining the two cases and taking the maximum constant, the result already follows.

Remark 4.17. An inspection of the proof of Lemma 4.16 show that the only geometric assumption needed is that the curved face F is given as the graph of a Lipschitz function φ. So Lemma 4.16 holds without assuming property 2) in Definition 4.9. 4.5. Inverse estimates on general domains. We now extend the above inverse estimates to general curved polytopic elements K ∈ T . To that end, we shall relax the concept of p-coverability of polytopic elements introduced in [START_REF] Cangiani | hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], (see also [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]) from simplicial coverings of general-shaped elements K ∈ T , to coverings involving affinely mapped generalized prisms. Definition 4.18. An element K ∈ T is said to be p-coverable with respect to p ∈ N, if there exists a set of m K ∈ N generalized prisms Kj and corresponding affine maps Φ j , such that the mapped generalized prisms K j := Φ j ( Kj ), j = 1, . . . , m K , form a, possibly overlapping, covering of K with the additional properties The motivation for the above definition is the stability result for polynomials with respect to domain perturbation given in Lemma 4.14 above. If K is p-coverable, (4.11) implies that there exists a covering of affinely mapped generalized prisms K j and respective sub-prisms K j := K j, , 0 < ≤ h Kj (8p) -2 , j = 1, . . . , m K , such that K j ⊂ K. Then, we have, for any v ∈ P p (K),

(4.13) 1 2 v 2 Kj ≤ v 2 K j ≤ v 2 K .
We now show that (4.11) is implied by Assumptions 4.1 and 4.3. Therefore, p-coverability for an element satisfying Assumptions 4.1 and 4.3 is ensured under the validity of (4.12) only.

Lemma 4.19. Let K ∈ T satisfying Assumptions 4.1 and 4.3. Then, there exists a set of m K ∈ N generalized prisms Kj and corresponding affine maps Φ j , such that the mapped domains K j := Φ j ( Kj ), j = 1, . . . , m K , form a cover of K with the property

(4.14) dist(∂K j , K) ≤ h Kj (8p) -2
for all j = 1, . . . , m K , with the notation of Definition 4.18.

Proof. From Assumption 4.3, ∂K is comprised of a finite number of closed (codimension one) C 1 surfaces (∂K) j , j = 1, . . . , zK , for some zK ∈ N. By possibly further subdividing the (∂K) j 's into subsets, say, (∂K) j , j = 1, . . . , z K , Assumption 4.1, ensures that for each of (∂K) j there exist a point x 0 j ∈ K such that the curved simplex K (∂K)j ≡ K (∂K)j (x 0 j ) is star-shaped with respect to x 0 j and that g j (x) := (x -x 0 j ) • n(x) > 0 for any x ∈ (∂K) j . (More than one (∂K) j are allowed to share the same x 0 j .) Since (∂K) j is C 1 , g j is continuous on (∂K) j and, thus, there exists a positive number δ j , such that g j (x) ≥ δ j . Now, for any x0 j ∈ R d , with |x 0 j -x0 j | < δ j , we have

|g j (x) -(x -x0 j ) • n(x)| ≤ |x 0 j -x0 j | < δ j
, and, therefore, gj (x) := (x -x0 j ) • n(x) > 0 for any x ∈ (∂K) j . Hence, (∂K) j is star-shaped with respect to B(x 0 j , δ j ) in K; that is any line connecting any point of (∂K) j with a point of B(x 0 j , δ j ) lies wholly in K. This implies that K (∂K)j is starshaped in K with respect to any subset of B(x 0 j , δ j ) and, in particular, with respect to any (d-1)-hypercube passing through x 0 j and contained in B(x 0 j , δ j ). In general, however, B(x 0 j , δ j ) ⊂ K, but we have dist(∂B(x 0 j , δ j ), K) ≤ |x 0 j -x0 j | < δ j . For the boundary pieces (∂K) j with B(x 0 j , δ j ) ⊂ K, we fix δ j to its largest possible value ensuring B(x 0 j , δ j ) ⊂ K. If, however, B(x 0 j , δ j ) ⊂ K, we select δ j small enough, so that dist(∂B(x 0 j , δ j ), K) ≤ h Kj (8p) -2 . On the other hand, the C 1 smoothness of (∂K) j ensures that there exists a finite tessellation comprising of diagonally scaled and rotated (d -1)-hypercubes approximating (∂K) j to a desired accuracy, say δ j /2. Consider now the truncated prisms intersecting (∂K) j defined uniquely by the 2 d-1 vertices of each element of the tessellation and the 2 d-1 vertices of a second (d-1)-hypercubical base contained in B(x 0 j , δ j ) and passing through x 0 j . The union of the latter generalized prisms covers K (∂K)j within a distance δ j . Considering the corresponding construction for all j, we conclude the construction of a finite cover of K by affinely mapped generalized prisms such that (4.14) holds.

Remark 4.20. The purpose of the construction in the proof of Lemma 4.19 is to assert the existence of at least one covering with the required properties, and not to construct the 'optimal' one.

4.15) v 2 Fi ≤ C INV (p, K, F i ) (p + 1)(p + d)|F i | |K| v 2 K , with (4.16) C INV (p, K, F i ) := min C reg (K, F i ), 2c -1 as 32 d p 2(d-1) , K p-coverable, C reg (K, F i ), otherwise, with C reg (K, F i ) := |K|/ |F i | sup x 0 i ∈K min x∈Fi (m i •n)
, and c as > 0 as in Def. 4.18.

Proof. If K is not p-coverable, using (4.2), we simply have

(4.17) v 2 Fi ≤ (p + 1)(p + d) min x∈Fi (m i • n) v 2 K F i ≤ (p + 1)(p + d)|F i | |K| C reg (K, F i ) v 2 K .
If, on the other hand, K is p-coverable, then K Fi ⊂ K ⊂ ∪ j=1,...,m K K j and, thus,

v 2 Fi ≤ |F i | v 2 L∞(Fi) ≤ |F i | v 2 L∞(K F i ) ≤ |F i | max j=1,...,m K v 2 L∞(Kj ) .
Now, Lemma 4.16 (together with an elementary scaling argument), along with (4.11) and (4.13), imply

v 2 L∞(Kj ) ≤ 32 d p 2d |K j | -1 v 2 Kj ≤ 2(32p d ) 2 c -1 as |K| -1 v 2 K .
Combining the last two estimates, taking the supremum over x 0 i ∈ K, the inverse estimate constant is then given by the minimum of the two estimates.

The above result generalizes both [20, Lemma 11] and [23, Lemma 4.9] in a number of ways. The coverings are now allowed to consist of curved domains. Also, elements with arbitrary number of (curved) faces are now admissible and an earlier hypothesis on uniform boundedness of m K across the mesh has now been removed by a more careful analysis. Note that, when K ∈ T is a polytopic element with straight faces, Lemma 4.21 collapses to [20, Lemma 11] with improved constants.

Remark 4.22. The sub-division {F i } n K i=1 of the (curved) element boundary ∂K is typically not unique. We can seek to minimize the coefficient (4.16) by considering different candidates for {F i } n K i=1 . However, such optimization would be practically beneficial only for rather "exotic" element shapes as, in most cases, we can simply resort to (4.1). Of course, extremely general curved "exotic" element shapes must be used only when deemed beneficial for the particular problem at hand. In such cases, a basic geometric study for improving the constant (4.16) (and, therefore, as we shall see below, the dG discontinuity-penalization function, cf. Remark 5.3 below) may be in order. In any case, Lemma 4.21 is sharp for each given subdivision {F i } n K i=1 and directly generalizes the inverse estimates in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. Next, we present an H 1 -L 2 -inverse inequality for polynomials on a general curved element K ∈ T . Lemma 4.23. Let K ∈ T satisfy Assumptions 4.1 and 4.3. Then, for each v ∈ P p (K), the inverse estimate

(4.18) ∇v 2 K ≤ C B INV (p, K) p 4 ρ 2 K v 2 K ,
holds, with ρ ω denoting the radius of the largest inscribed circle of a domain ω ⊂ R d , and

(4.19) C B INV (p, K) := 4 min ρ cov (K), ρ p-cov (p, K) , if K p-coverable 4ρ cov (K), otherwise, with (4.20) ρ cov (K) := m K j=1 C B inv (d, rj ) rj ρ K ρ Kj 2 , and
(4.21) ρ p-cov (p, K) := c -1 as (32) d (p -1) 2d max 1≤ ≤m K C B inv (d, r ) r ρ K ρ K 2 ,
for K , = 1, . . . , m K , cover of K consisting of affinely mapped generalised prisms.

Proof. From Lemma 4.19, there exists a cover of K, consisting of affinely mapped generalized prisms K j , j = 1, . . . , m K . Thus, for v ∈ P p (K), Lemma 4.12, (with a standard affine scaling) and (4.13) imply:

∇v 2 K ≤ m K j=1 ∇v 2 Kj ≤ m K j=1 C B inv (d, rj ) 2r 2 j p 4 ρ 2 Kj v 2 Kj ≤ 4 m K j=1 C B inv (d, rj ) r2 j p 4 ρ 2 Kj v 2 K ,
with rj denoting the r of K j as per Definition 4.9. Thus, we have

(4.22) ∇v 2 K ≤ ρ cov (K) p 4 ρ 2 K v 2 K .
Note that ρ cov (K) grows with ρ K / min j=1,...,m K ρ Kj growing.

On the other hand, if K is p-coverable, there exists ∈ {1, . . . , m K } such that

∇v 2 K ≤ |K| ∇v 2 L∞(K) ≤ |K| ∇v 2 L∞(K )
. Using now Lemmata 4.16 and 4.12 (with scaling), as well as (4.12), we deduce, respectively, (4.23)

∇v 2 L∞(K ) ≤ 32 d (p -1) 2d |K | -1 ∇v 2 K ≤ 32 d (p -1) 2d p 4 C B inv (d, r )2r 2 ρ -2 K |K | -1 v 2 K .
The result already follows by combining (4.22) and (4.23).

In the special case of an element K being star-shaped with respect to a contained ball, we can have a more precise statement in terms of the constants involved. Corollary 4.24. Let K ⊂ R d domain which is star-shaped with respect to a ball B(x, ρ K ) ⊂ K, x ∈ K. Then, for any v ∈ P p (K), we have the inverse estimate

∇v 2 K ≤ C(d) h K ρ K d+5 p 4 ρ 2 K v 2 K ,
for some universal constant C(d) > 0 that can be estimated explicitly. Thus, if additionally, K is shape-regular, i.e., h K ∼ ρ K , we retrieve the classical inverse estimate ∇v 2 K ≤ Cp 4 /h 2 K v 2 K , with C now also dependent on the shape-regularity constant.

Proof. We have B(x, ρ K ) ⊂ K ⊂ B(x, h K ). A comparison of the area of the largest (d-1)-hypercube contained in B(x, ρ K ), given by ρ d-1 K /2 (d-1)/2 , with the surface of B(x, h K ), shows that we can cover K using z := 2 (d+1)/2 π d/2 h d-1 K /(Γ(d/2)ρ d-1 K ) + 1 mapped right generalized prisms K j , j = 1, . . . , z, whose bases are given by rotations of the largest (d -1)-hypercube contained in B(x, ρ K ). So, we have

∇v 2 K ≤ z j=1 ∇v 2 Kj ≤ z j=1 C B inv (d, rj ) 2r 2 j p 4 ρ 2 Kj v 2
Kj ;

here we have used scaling via Φ j : Kj → K j affine mapping [0, 1] d to a rotation of the largest d-hypercube contained in B(x, ρ K ). Since each K j is right, we have rj ≤ h K /ρ Kj + 1, for all j = 1, . . . , z. Also, from the star-shapedness with respect to B(x, ρ K ), we have ρ Kj ≥ ρ K /2. Combining the above, we deduce

∇v 2 K ≤ 8C B inv ( h K /ρ Kj + 1, d)( h K /ρ Kj + 1) 2 p 4 ρ 2 K z j=1 v 2 Kj .
Using the (pessimistic) bound z j=1 v 2 Kj ≤ z v 2 K , and combining the numerical constants, the result follows.

The last result holds under weaker domain assumptions compared to [45, Theorem 1] and, in contrast to the main result in [START_REF]Sharp Lp Markov type inequality for cuspidal domains in R d[END_REF], it offers explicit dependence on the domain size in the case of piecewise C 1 domains. We also note [START_REF] Kroó | On Bernstein-Markov-type inequalities for multivariate polynomials in Lq-norm[END_REF]Theorem 3], which provides a similar bound for the special case of K being a d-ellipsoid, in conjunction with John's Ellipsoid Theorem. It is interesting to investigate the extension of the above inverse estimates with explicit constants to cuspoidal domains in the spirit of [START_REF]Sharp Lp Markov type inequality for cuspidal domains in R d[END_REF]; this will be considered elsewhere.

Example 4.25. We revisit Example 4.6 for d = 2, with K = B(0, R) a circular element with radius R.

Let K 1 = B(0, R) ∩ {(x 1 , x 2 ) : x 1 ∈ R, - √ 2R/2 ≤ x 2 ≤ √ 2R/2} and K 2 = B(0, R) ∩ {(x 1 , x 2 ) : - √ 2R/2 ≤ x 1 ≤ √ 2R/2, x 2 ∈ R}, so that K = K 1 ∪ K 2 .
We further subdivide each K i in half to form prisms with flat base; for each of these, we can select r = 2. Thus, Lemma 4.23 implies

∇v 2 K ≤ C B INV p 4 R 2 v 2 K , with C B INV ≤ 2304.
The constant in this special case can be improved considerably upon taking advantage of the circle's symmetries. 7. Denoting by r the length of each of the (equal length) n small faces and with h K its diameter, we consider the case r h K . If r < h K /p 2 , we can cover K by one triangle , namely, the smallest simplex containing K. Then K is p-coverable and C B INV (p, K) remains bounded, independently of n. Hence, when the two geometric scales h K and r are significantly different, K is essentially a simplex in this context.

On the other hand, for p large enough and fixed r and n, we have r > h K /p 2 and, hence, we cannot cover K as before. Instead, we consider a family of n nonoverlapping simplices K j ⊂ K, each defined by one 'small' face of length r and the vertex x 0 . Then, we have c as = n -1 in Definition 4.18 and ρ Kj ∼ h K /n. Since also ρ K ∼ h K and q = 1, we compute

C B INV (p, K) ∼ n -1
. This is reasonable, as sufficiently high polynomial degree p basis functions can resolve the scale of the 'sawtooth' face ensemble. 4.6. Best approximation estimates. We now turn to hp-version polynomial approximation bounds over general domains. The setting here remains essentially unchanged compared to the case of just polytopic elements presented in [START_REF] Cangiani | hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. More specifically, under a mild set of covering assumptions and upon postulating the existence of so-called function space domain extension operators, we are able to apply hp-version best approximation results in various seminorms. Definition 4.27. Given a mesh T , we define a covering T = {K} of T to be a set of open shape-regular d-simplices K, such that for each K ∈ T , there exists a K ∈ T with K ⊂ K. For a given T , we define the covering domain Ω := ∪ K∈T K.

For illustration, in Figure 8 we show a single two-dimensional curved element K ∈ T , along with a covering simplex K ∈ T with K ⊂ K.

K K • • • • • Figure 8. A simplex K ∈ T covering an element K ∈ T .
Assumption 4.28. For a given mesh T , we postulate the existence of a covering T , and of a (global) constant O Ω ∈ N, independent of the mesh parameters, such that max

K∈T card K ∈ T : K ∩ K = ∅, K ∈ T such that K ⊂ K ≤ O Ω .
For such T , we further assume that h K := diam(K) ≤ C diam h K , for all pairs K ∈ T , K ∈ T , with K ⊂ K, for a (global) constant C diam > 0, uniformly with respect to the mesh size h K .

Remark 4.29. Assumption 4.28 ensures the shape-regularity of the mesh covering T only. Shape-regularity of the mesh T is not assumed. We refer to Figure 3.6 in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] for an example on how these two concepts may differ considerably.

The validity of Assumption 4.28 allows for the use known hp-version approximation estimates on simplicial elements [START_REF] Babuška | The h-p version of the finite element method with quasi-uniform meshes[END_REF][START_REF]The optimal convergence rate of the p-version of the finite element method[END_REF][START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF], on each K and, subsequently restrict the error over K ⊂ K. However, it requires to extend the exact solution u into Ω in a stable fashion. To that end, we shall use the following classical result. Theorem 4.30 ([57]). Let Ω be a domain with a Lipschitz boundary. Then there exists a linear extension operator E :

H s (Ω) → H s (R d ), s ∈ N 0 , such that Ev| Ω = v and Ev H s (R d ) ≤ C E v H s (Ω)
, where the constant C E > 0 depends only on s, Ω.

Subsequent refinements of the dependence of the constant C E on the domain shape in Theorem 4.30, have been presented for instance in [START_REF] Sauter | Extension operators and approximation on domains containing small geometric details[END_REF][START_REF] Carstensen | A posteriori error analysis for elliptic PDEs on domains with complicated structures[END_REF].

For the estimation of the best approximation error on the mesh skeleton Γ int ∪∂Ω, we require the trace estimate on curved domains from Lemma 4.7.

We now have all the ingredients to assert the validity of the following hpapproximation error bounds. Lemma 4.31. Let K ∈ T satisfy Assumptions 4.1 and 4.28, and let K ∈ T be the corresponding simplex with K ⊂ K as per Definition 4.27. Suppose that v ∈ L 2 (Ω) is such that Ev| K ∈ H l K (K), for some l K ≥ 0, and that Assumption 4.28 is satisfied. Then, there exists an operator π p : H l K (K) → P p (K), such that

(4.24) v -π p v H q (K) ≤ C 1 h s K -q K p l K -q Ev H l K (K) , l K ≥ 0,
for 0 ≤ q ≤ l K , and

(4.25) v -π p v Fi ≤ C 1/2 ap (p, K, F i )|F i | 1/2 h s K -d/2 K p l K -1/2 Ev H l K (K) , l K > d/2, with C ap (p, K, F i ) := C 2 min h d K |F i | sup x 0 i ∈K min x∈Fi (m i • n) -1 , p d-1 ,
s K = min{p + 1, l K }, and C 1 , C 2 > 0 constants depending only on the shaperegularity of K, q, l K , on C diam (from Assumption 4.28) and on the domain Ω.

Proof. Let Π p : H l (K) → P p (K) be a known optimal hp-version approximation operator on simplices, see, e.g., [START_REF] Babuška | The h-p version of the finite element method with quasi-uniform meshes[END_REF][START_REF]The optimal convergence rate of the p-version of the finite element method[END_REF][START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF]. We define π p : H l (K) → P p (K) by π p v := Π p (Ev). To prove (4.24), we begin by observing that

v -π p v H q (K) = Ev -Π p (Ev) H q (K) ≤ Ev -Π p (Ev) H q (K) .
Thus, Assumption 4.28 and standard hp-approximation estimates on simplices (e.g. [START_REF] Babuška | The h-p version of the finite element method with quasi-uniform meshes[END_REF][START_REF]The optimal convergence rate of the p-version of the finite element method[END_REF][START_REF] Schwab | p-and hp-Finite element methods: Theory and applications in solid and fluid mechanics[END_REF] yield the desired bound; we refer to the proof of [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Lemma 3.7] for a similar argument for polytopic elements.

To prove (4.25), we use the trace inequality (4.3) with ζ = p to deduce

(4.26) v -π p v 2 Fi ≤ C h K sup x 0 i ∈K min x∈Fi (m i • n) h 2s K -1 K p 2l K -1 Ev 2 H l K (K) , noting that max x∈Fi |m i | 2 2 ≤ h 2 K . On the other hand, we observe that v -π p v 2 Fi ≤ |F i | v -π p v 2 L∞(K F i ) ≤ Ev -Π p (Ev) 2 L∞(K)
. Hence, employing a classical hp-approximation estimate for the maximum norm error from [START_REF] Babuška | The h-p version of the finite element method with quasi-uniform meshes[END_REF][START_REF]The optimal convergence rate of the p-version of the finite element method[END_REF], (cf. also [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Lemma 20] we arrive at

(4.27) v -π p v 2 Fi ≤ C|F i | h 2s K -d K p 2l K -d Ev 2 H l K (K) , for l K > d/2.
The result follows by taking the minimum between the bound in (4.26) and the bound in (4.27).

A priori error analysis

We are now ready to briefly discuss a priori error bounds for sufficiently smooth exact solutions, thereby generalizing the respective results presented in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] to the case of curved polytopic meshes. The line of argument is similar to the case of straight polytopic meshes presented in detail in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF].

A crucial ingredient of the analysis for the proof of stability of the dG-EASE method is the precise definition of the discontinuity-penalization function σ appearing in the method (3.3). It is important to define σ sufficiently large for stability, while at the same time not substantially larger than what is required, as it could potentially cause loss of accuracy and/or conditioning issues. Additionally, following [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], we provide a stronger inf-sup stability result with respect to a 'steamlinediffusion'-type augmented norm, when the wind b is non-zero. The size of the 'steamline-diffusion' coefficient depends crucially on Lemma 4.23, whose constant provides information on the stabilization capabilities of the method.

The dG norm for which we seek to prove a priori error bounds is given by

| v| := | v| 2 ar + | v| 2 d 1/2
, where

| v| 2 ar = c 0 v 2 + 1 2 K∈T |b • n| v 2 ∂-K + |b • n|v 2 ∂+K∩∂Ω ,
with c 0 given in (2.4), and

| v| 2 d = √ a∇ T v 2 + √ σ[[v]] 2 Γint∪∂ΩD .
Definition 5.1. For a mesh T , we define the set F int of interfaces F ⊂ Γ int by

F int := {F ⊂ Γ int : there exist K, K ∈ T with F = ∂K ∩ ∂K };
correspondingly, we set F D := {F ⊂ ∂Ω D : there exists K ∈ T with F = ∂K ∩ ∂Ω D }. For notational compactness, we also define F int,D := F int ∪ F D . (Note that F may comprise of one or more faces of K, K .) Moreover, each interface F may be contained in one or more F i 's of the elements K, K as per Assumption 4.1. Thus, there exists a subset

{F K i } i∈I K F with index set I K F ⊂ {1, . . . , n K }, such that F ⊂ ∪ i∈I F K F K i and, correspondingly, a set I K F ⊂ {1, . . . , n K } such that F ⊂ ∪ i∈I F K F K i .
For technical reasons (cf. [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] and the references therein), we shall make use of the following extensions Bd : (H

1 (Ω) + S p T ) × (H 1 (Ω) + S p T ) → R and ˜ : (H 1 (Ω) + S p T ) → R
σ : Γ int ∪ ∂Ω D → R for every interface F ∈ F int,D , F = ∂K ∩ ∂K , by (5.1) σ| F := 2 max K∈{K,K } |I K F | max i∈I K F C INV (p K , K, F K i )|F K i | āK (p K + 1)(p K + d) |K| ;
when F ∈ F D we set K = K . Then, for all w, v ∈ H 1 (Ω) + S p T , we have

(5.2) Bd (w, w) ≥ 1 2 | w| 2 d and Bd (w, v) ≤ 2| w| d | v| d .
Proof. The idea of proof is standard and makes use of the trace inverse estimate developed above. The novel attribute here is the choice of σ which requires some care since the star-shapedness of each interface F may correspond to different boundary segments F i in either side of the interface. To that end, for w ∈ H 1 (Ω) + S p T , we have Bd (w, w)

≥ | w| 2 d -2 Γint∪∂ΩD { {aΠ∇w} } • [[w]
] ds. Therefore, Lemma 4.21 and the stability of the orthogonal L 2 -projection give

Π √ a∇w 2 F ∩∂K ≤ |I K F | max i∈I K F C INV (p K , K, F K i )|F K i | (p K + 1)(p K + d) |K| √ a∇w 2 K ,
Coercivity already follows by a Young's inequality. The proof of continuity is standard and, therefore, omitted for brevity.

Remark 5.3. The stability of the dG-EASE method is guaranteed under extremely general mesh assumptions thanks to the judicious choice of the penalization function (5.1). As discussed also in Remark 4.22, the latter ultimately depends on the choice of subdivisions {F i } n K i=1 of ∂K appearing in Assumption 4.1. Of course, whenever possible, by simply following the recipe in Remark 4.2(ii), we can easily arrive at a practical value of the penalization function for general curved elements.

We shall additionally assume for simplicity of the presentation that

(5.3) b • ∇ξ ∈ S p T ,
for all ξ ∈ S p T , as is a standard in this context, cf. [START_REF]Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF] and also [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Chapter 5]. Assumption (5.3) can be further relaxed at the expense of an additional mild suboptimality with respect to the polynomial degree p; see [START_REF]Discontinuous hp-finite element methods for advection-diffusion-reaction problems[END_REF]Remark 3.13] and Remark 5.8 below.

Theorem 5.4. Let T = {K} a subdivision of Ω ⊂ R d , consisting of, possibly curved, elements satisfying Assumptions 4.1, 4.3 and 4.28. Then, assuming that (5.3) holds and that the discontinuity-penalization function σ is given by (5.1), there exists a constant Λ s > 0, independent of h and of p, such that:

(5.4) inf

w∈S p T \{0} sup v∈S p T \{0} B(w, v) | w| s | v| s ≥ Λ s , with | v| s := | v| 2 + K∈T λ K b • ∇v 2 K 1/2
, whereby

λ K := min ρ K C B INV (p K , K) , F ⊂∂K i∈I K F C INV (p K , K, F K i ) |F K i | |K| -1 max{ b L∞(K) , σ K }(p K + 1)(p K + d) , for K ∈ T , p K ≥ 1, σ K := max{σ a K , σ b K }, with σ a K := max F ⊂∂K σ| F , and
σ b K := 2 max F ⊂∂K max K∈{K,K } C B INV (p K , K) āK (p K + 1)(p K + d) ρ K . (5.5)
Proof. The proof follows in a completely analogous fashion to the proof of [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Theorem 5.2] and is, therefore, largely omitted for brevity: the key idea is to set

v ≡ v(w) := w + αw s , for w ∈ S p T with w s | K := λ K b • ∇w, K ∈ T , with R α > 0. Then, it is sufficient to prove that | v| s ≤ C * | w| s , and B(w, v) ≥ C * | w| 2
s , and then to set Λ s = C * /C * , for some C * , C * > 0 constants independent of the discretization parameters. The last two conditions are proven by using the inverse estimates above along with a judicious use of α.

Remark 5.5. Theorem 5.4 extends respective results on polytopic meshes from [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], to general meshes consisting of polytopic and/or curved elements with arbitrary number of faces. Moreover, this choice removes a dependence of the inf-sup constant Λ s on the inverse inequality constants C INV and C B INV ; cf. [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF][START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]. Theorem 5.6. Let T = {K} be a subdivision of Ω ⊂ R d , consisting of general curved elements satisfying Assumptions 4.1, 4.3 and 4.28. Let also T = {K} an associated covering of T consisting of shape-regular simplices as per Definition 4.27. Assume that (5.3) holds. Assume that u ∈ H 1 (Ω) the exact solution to (2.1), (2.3), is such that u| K ∈ H l K (K), l K > 1 + d/2, for each K ∈ T . Let u h ∈ S p T , with p K ≥ 1, K ∈ T , be the solution of (3.3), with σ as in (5.1). Then, we have

| u -u h | 2 s ≤ C K∈T h 2s K K p 2l K K (D K (F, C ap , λ K , p K ) + G K (F, C INV , C ap , p K )) Eu 2 H l K (K) , with s K = min{p K + 1, l K }, D K (F, C ap , λ K , p K ) = c 0 2 L∞(K) + ζ 2 K + λ -1 K + λ K β 2 K p 2 K h -2 K + āK p 2 K h -2 K + β K h -d K p K F ⊂∂K |I K F | max i∈I K F C ap (p K , K, F K i )|F K i | , and
G K (F, C INV , C ap , p K ) = ā2 K p 3 K h -d-2 K F ⊂∂K∩(Γint∪∂ΩD) σ -1 |I K F | max i∈I K F C ap (p K , K, F K i )|F K i | + ā2 K p 4 K h -2 K |K| -1 F ⊂∂K∩(Γint∪∂ΩD) σ -1 |I K F | max i∈I K F C INV (p K , K, F K i )|F K i | + h -d K p K F ⊂∂K∩(Γint∪∂ΩD) σ|I K F | max i∈I K F C ap (p K , K, F K i )|F K i | , s K = min{p K + 1, l K }, ζ K := c/c 0 L∞(K) , c 0 is in (2.4), β K := b L∞(K) ,
and C is a positive constant, which depends on the shape-regularity of T , but is independent of the discretization parameters. In the special case in which the coefficient a is strictly positive definite a.e. in Ω while b = 0, Assumption 4.3 can be removed from the hypotheses.

Proof. The proof follows on very similar lines to the respective one for polytopic meshes and can be found in [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF]Section 5.2].

The above hp-version a priori error bounds hold without any assumptions on the relative size of the faces F , F ⊂ ∂K, of a given curved element K ∈ T . To aid the understanding of the rates of convergence resulting from the above results, we set p K = p ≥ 1, h = max K∈T h K , s K = s, s = min{p + 1, l}, and l > 1 + d/2, and assume that diam(F ) ∼ h K , for all faces F ⊂ ∂K, K ∈ T , so that |F | ∼ h (d-1) K . Then, Theorem 5.6 reduces to

| u -u h | d ≤ C h s-1 p l-3 2 u H l (Ω) ,
i.e., it proves optimal convergence in h and suboptimal in p by p 1/2 . At the other end of the spectrum, consider the case of transport equation, i.e., when a ≡ 0. In this case, the dG norm | •| degenerates to | •| ar ; note that, then we have λ K = O(h K /p 2 K ), and the a priori error bound in Theorem 5.6 reduces to

| u -u h | ar ≤ h s-1 2 p l-1 u H l (Ω) .
This bound is, again, optimal in h and suboptimal in p by p 1/2 and completely generalizes the error estimate derived in our previous work [START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF] to essentially arbitrarilyshaped meshes under the same assumption (5.3).

Remark 5.7. We remark on typical cases which result to simplified formulas for

λ K . Assuming |K| ∼ h d K , ρ K ∼ h K , and |F | ∼ h (d-1) K
for an element K ∈ T and for its immediate neighbours, both constants C B INV and C INV will be defined by the first term in the maxima in (4.16) and (4.19), respectively. Then, we deduce λ K ∼ h K /p 2 K for the important case of advection-dominated problems. Remark 5.8. For general advection fields b, the proof of the inf-sup condition needs to be modified by using a slightly different norm involving Π(b • ∇ T ) instead of (b•∇ T ) in the s-norm, yielding an error bound which is optimal in h but suboptimal in p by p 3/2 for the purely hyperbolic problem. Of course, if we modify the method by including a streamline-diffusion stabilization term as done in [START_REF] Houston | Stabilized hp-finite element methods for first-order hyperbolic problems[END_REF], then an hpoptimal bound can be derived without enforcing (5.3).

Numerical examples

We test the dG-EASE approach through a series of numerical experiments using curved elements, ranging from basic domain approximation to highly complex element shapes arising from random element agglomeration of a fine background triangulation.

In the case of the agglomeration-constructed elements, the background (curved) triangulation is also used for the assembly step. In particular, the discontinuitypenalisation function σ| F is fixed following the recipe in (5.1) with the subdivisions {F i } n K i=1 of ∂K, K ∈ T , appearing in Assumption 4.1, given by unions of faces of the background triangulation. Moreover, for simplicity the background triangulation is also used for integration, exploiting parallellization of the quadrature process [START_REF] Dong | GPU-accelerated discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF], see also [START_REF] Cangiani | hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] for a more detailed discussion of implementation of such methods. Nonetheless, very often it is possible to use substantially coarser subdivisions than the background triangulation the elements have been constructed from, e.g., a subdivision with one simplex per straight face.

For curved elements, the current implementation performs quadrature by a sufficiently fine sub-triangulation approximating the curved element, exactly as in the agglomerated-element case. We stress, however, that in this case the subtriangulation is only used to generate the quadrature rules. These calculations are fully parallelizable : in [START_REF] Dong | GPU-accelerated discontinuous Galerkin methods on polygonal and polyhedral meshes[END_REF] it is shown that quadrature cost becomes irrelevant if modern GPU architectures are used in the implementation of assembly. Of course, this is not the only possibility. For instance, domain-exact quadrature algorithms for many curved domains exist, see, e.g., [START_REF] Artioli | Algebraic cubature on polygonal elements with a circular edge[END_REF] and the references therein for such algorithms.

6.1. Example 1: curved elements. We begin by testing the method on triangular elements with (non-parametric) curved faces. Specifically, we consider a two-dimensional diffusion problem with a = I 2×2 , I 2×2 denoting the 2 × 2-identity matrix, b = (0, 0) , c = 0 and f so that u(x 1 , x 2 ) = sin(πx 1 ) sin(πx 2 ) in (2.1). We solve this problem on an irregular annular domain constructed as the unit disc centred at origin, with a circular hole centred at (0.25, 0.25) and radius 0.4 removed; cf. Figure 9 for an illustration. We construct a sequence of domain-fitted curvilinear meshes as follows. First, using the mesh generator from [START_REF] Persson | A simple mesh generator in MATLAB[END_REF], we construct a sequence of meshes approximating the domain Ω comprising of 65, 109, 527, 2266, and 9411 quasi-uniform triangular elements, respectively. The 65-element mesh is shown in Figure 9. Then, exploiting the knowledge of the level-set function of ∂Ω, elements containing straight faces approximating the curved boundary are marked. Finally, all marked elements are treated as curved triangular elements with two straight faces and one curved face described by the domain level-set function, thus capturing the domain exactly.

In Figure 9 (top row), we present the convergence history of | u -u h | and u -u h Ω against √ Dof, with Dof the number of degrees of freedom on the aforementioned curvilinear meshes with 65, 109, 527, 2266, and 9411 elements, for p = 1, 2, 3, 4, respectively. We clearly observe that, for each fixed p, all errors converge to zero at the optimal rates O(h p ) and O(h p+1 ), respectively, as the mesh size h tends to zero. In the two bottom plots in Figure 9, we investigate the convergence history of the dG-EASE solution under p-refinement, using the two meshes with 65 and 109 curved elements, respectively, in linear-log scale. Here, we observe exponential convergence of all errors against √ Dof.

6.2. Example 2: convergence study. We now investigate the convergence of dG-EASE on a highly complex mesh comprising of elements arising from agglomeration of a very fine background mesh, which also contains curved boundary-touching elements. Set a = I 2×2 and = 0.01, b = (1 -x 2 , 1 -x 1 ) , c = 2 and f so that u(x 1 , x 2 ) = sin(πx 1 ) sin(πx 2 ) in (2.1) for d = 2, on a domain Ω ≈ (0, 1) 2 enclosed by a piecewise curved sinusoidal boundary; we refer to Figure 10 for an illustration. We impose non-homogeneous Dirichlet boundary conditions on ∂Ω. The mesh is constructed as follows. An initial curved mesh, fitted to the sinusoidal boundary via the level set approach described above, is subdivided into a very fine background subdivision consisting of approximately 500K sub-elements. The latter is, in turn, agglomerated into 30, 132, 555, 2151, 8337 curved/polygonal elements using a standard mesh partitioning software. The parameters chosen in the partitioning software have been selected to yield a high-frequency 'sawtooth' vertical boundary for many of the agglomerated elements. We refer to Figure 10 for an illustration of the resulting meshes with 30 and 132 agglomerated elements.

In Figure 11, the convergence history for p = 1, 2, 3, 4 for the errors | u -u h | and u -u h Ω against √ Dof is presented for the aforementioned agglomerated meshes with 30, 132, 555, 2151, 8337 elements. Here, we clearly observe that, for each fixed p, the dG-and L 2 (Ω)-norm errors converge to zero at the optimal rates O(h p ) and O(h p+1 ), respectively, as the mesh size h tends to zero. Further, we report also the error in the stronger 'streamline-diffusion' norm | u -u h | s in Figure 11; here we have chosen λ K = O(ρ K /p 2 K ). For each fixed p, the errors converge to zero at the optimal rates O(h p ), as the mesh size h tends to zero.

Finally, in Figure 11 (bottom-right), we also investigate the convergence history of the dG-EASE solution under p-refinement, using the mesh with 132 elements shown in Figure 10(right). Here, we observe exponential convergence of the three norm errors against √ Dof. Interestingly, we observe that the difference between the errors | u -u h | and | u -u h | s is insignificant. 6.3. Example 3: stability study. We continue by assessing the stability of the dG-EASE method for convection-diffusion problems in the presence of unresolved lower-dimensional sharp solution layers. To this end, for d = 2, we set a = I 2×2 and = 10 -4 , b = (1, 1) , c = 0 and f = 1 in (2.1). We solve this problem on a variant of the domain Ω from Example 2 above, in which circular internal pieces of the domain of various radii have been removed; we refer to Figure 12(left) for an illustration of the domain and sample mesh of essentially arbitrarily-shaped elements obtained using a completely analogous construction to that used in Example 2. We close the problem by prescribing homogeneous Dirichlet boundary conditions on ∂Ω (i.e., including the internal boundaries at the holes). We expect strong exponential boundary layers on the top and right portions of the curved boundary, as well as variable intensity layers at the outflow portions of the internal hole boundaries.

In Figure 12 (right), we provide the dG-EASE solution using p = 1 and the mesh of 531 elements shown on the left plot. This mesh is not fine enough to resolve the singularly perturbed behaviour in the vicinity of the outflow portions of the boundary. Nevertheless, the dG-EASE method provides a stable discretization with very localized, expected, oscillatory behaviour at the vicinity of the outflow boundary. The stable behaviour of dG-EASE with respect to the size of the Péclet number P e := b / is expected due to the upwind flux used in B ar (•, •); nonetheless, to the best of our knowledge, its performance in the context of elements with such geometrical shape generality has not been tested before in the literature. To highlight the behaviour of the method on different meshes, we report the dG-EASE solution, obtained with meshes composed of 129 and 2048 linear elements in Figure 13 (top). In both cases the mesh not sufficiently fine to resolve the exponential boundary layer behaviour, while the finer mesh with 2048 linear elements sufficiently resolves the parabolic layers initiated at the holes.

Finally, we test the hyperbolic limit case by setting = 0. The DG-EASE solution, shown in Figure 13 (bottom), remains stable and there is no oscillation around the outflow boundaries, as expected by a stabilised method. 6.4. Example 4: changing type PDE across a curved interface. To highlight a number of attractive features of the dG-EASE approach, we consider a coupled parabolic-hyperbolic partial differential equation, whose type changes across a sinusoidal interface Γ. Let Ω = Ω 1 ∪ Ω 2 ∪ Γ with

Ω 1 = {(x 1 , x 2 ) : -1 ≤ x 1 ≤ 1, A sin(ωπx 1 ) ≤ x 2 ≤ 1}, Ω 2 = {(x 1 , x 2 ) : -1 ≤ x 1 ≤ 1, -1 ≤ x 2 ≤ A sin(ωπx 1 )},
for A, ω > 0 whose precise values will be given below; we refer to Figure 14 for an illustration. On this geometrical setting, we consider the problem: This problem is hyperbolic when x 2 ≤ A sin(ωπx 1 ), x 1 ∈ (-1, 1), and parabolic otherwise. The normal flux of the exact solution is continuous across the interface Γ with equation x 2 = A sin(ωπx 1 ), while the solution itself has a discontinuity across the interface. This problem is a variant of an example from [START_REF] Georgoulis | Discontinuous Galerkin methods on shape-regular and anisotropic meshes[END_REF][START_REF] Cangiani | hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes[END_REF]. As such, there is no discontinuity penalisation imposed at the interface Γ. Moreover, we point out that b • n = at the interface in this example.

Our aim is to highlight the performance of dG-EASE of arbitrary order, when the mesh is fitted with respect to a discontinuity of the exact solution. To that end, we focus on p-version convergence, using 64 rectangular elements with curved faces exactly fitting the interface; we refer to Figure 14 for an illustration with A = 0.025, and ω = 8 and ω = 16, respectively.

Interestingly, the mesh is not aligned with the inflow and outflow parts of the boundary ∂Ω. This is due to the oscillating coefficient of the first order term. In Figure 14, the inflow parts of the boundary are marked in red; on these parts, inflow boundary conditions are imposed. Correspondingly, this pattern continues in the internal element faces in which the inflow parts of ∂ -K are also not aligned with the faces. As such, face integral terms in the dG method may be computed only on parts of a face of a rectangular element. Nonetheless, the method is able to cope unaltered with this complication. The quadrature is implemented in the composite fashion described in Example 1 above.

We begin by setting A = 0.025 and ω = 8. In Figure 15, we record the convergence history of the dG-EASE solution under p-refinement, using the mesh shown in Figure 14 (left) and p = 1, . . . , 13. Although the elements are perfectly aligned with the interface Γ, the mesh is still coarse: each element includes roughly one full oscillation of the wind b. Still we observe exponential convergence of | u -u h | and u -u h Ω errors against √ Dof under p-refinement. This result reinforces the claim that dG-EASE on perfectly aligned meshes with appropriate quadrature rules can lead to spectral accuracy. In contrast, if the mesh is not aligned exactly with the solution's discontinuity, the error is only expected to decay at an algebraic rate, according to standard best approximation results.

Next, we set ω = 16 and we record the convergence history under p-refinement, for p = 1, . . . , 17, for the fixed mesh from Figure 14(right). Here 64 elements constitute a very coarse mesh as, at the interface, there are now two full oscillations of the wind b per element. Again, we observe exponential convergence of | u -u h | and u -u h Ω against √ Dof under p-refinement, after an initial plateau for p ≤ 5. This is expected as the dG-EASE approach is not designed as a multiscale framework. Nevertheless, for p ≥ 6, exponential convergence is observed.

Acknowledgements

We are grateful to the anonymous referees and to the editor for their constructive comments which helped to improve this work substantially. AC gratefully acknowledges support from the MRC (MR/T017988/1), ZD from IACM-FORTH, Greece, and EHG from The Leverhulme Trust (RPG-2015-306). This work was supported by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the "First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant" (Proj. no. 3270).

  and B d (•, •) corresponds to the diffusion part:(3.5)

Figure 2 .

 2 Figure 2. Elements K ∈ T tr are assumed to satisfy Assumption 4.1 (a) (left) and (b) (right); • denotes a vertex.

  Figure 3. Curved elements K, K with, respectively, 8 and 4 sub-elements satisfying Assumption 4.1.

Assumption 4 . 3 .

 43 We assume that the boundary ∂K of each element K ∈ T is the union of a finite (yet, arbitrarily large!) number of closed C 1 surfaces. Assumption 4.1 is sufficient for the proof of the trace estimates presented below. Requiring both Assumptions 4.1 and 4.3 is sufficient for the validity of the H 1 -L 2 inverse estimate presented below. 4.1. Basic trace estimates. We now discuss the new trace-inverse estimate and a version of the standard Sobolev trace estimate for Lipschitz domains satisfying Assumption 4.1.

( 4

 4 .11) dist(∂K j , K) ≤ h Kj (8p) -2 and (4.12) |K j | ≥ c as |K|, for all j = 1, . . . , m K , where h Kj := sup x∈ F 0 |Φ j ( x,j )| and c as is a positive constant, independent of K and of T , with dist(∂K j , K) := sup x∈∂K inf y∈K |x -y| the onesided Hausdorff distance of ∂K j from K, and x,j := Kj ∩ {x + αe d : α ∈ R}.

Lemma 4 . 21 .

 421 Let K ∈ T Lipschitz satisfying Assumption 4.1. Then, for each v ∈ P p (K), we have the inverse inequality

(

K x 0 Figure 7 .

 07 Figure 7. Example 4.26. K ∈ T with 'multiscale' boundary behaviour.

Example 4 . 26 .

 426 Let d = 2, and consider the polygonal element K ∈ T with 'multiscale' boundary behaviour depicted in Figure

Remark 4 . 32 .

 432 We note the correspondence between C INV (p, K, F i ) from Lemma 4.21 and C ap (p, K, F i ), in the typical case h d K ∼ |K|. The key attribute of both expressions is that they remain bounded for degenerating |F i |, allowing for the estimates (4.15) and (4.25) to remain finite as |F i | → 0. Remark 4.33. If the constant C ap (p, K, F i ) in (4.25) is taken with the first term, then the approximation result (4.25) holds with l K > 1/2.

  of the bilinear and linear forms (3.5) and (3.6), which are given replacing { {a∇w} } and { {a∇v} } with { {aΠ∇w} } and { {aΠ∇v} } in B d and , respectively, where Π : [L 2 (Ω)] d → [S p T ] d denotes the orthogonal L 2 -projection operator onto the (vectorial) finite element space. Observe that Bd (w, v) = B d (w, v) and ˜ (v) = (v) when w, v ∈ S p T . Similarly, we define B(w, v) := Bd (w, v) + B ar (w, v). Next, we discuss the coercivity and continuity of Bd . Lemma 5.2. Let (3.7) hold and consider a mesh T satisfying Assumption 4.1. With the notation introduced in Definition 5.1, define the discontinuity-penalization function

Figure 9 .

 9 Figure 9. Example 1. Top: convergence history for p = 1, 2, 3, 4 in the | •| and • Ω norms against √ Dof for the curved triangular mesh exemplified in the bottom (right) plot with 65, 109, 527, 2266, and 9411 elements, respectively. Bottom (left): Convergence history for p = 1, 2, . . . , 7 in the | •| and • Ω against √ Dof for the meshes with 65 and 109 elements.

Figure 10 .

 10 Figure 10. Example 2. The computational domain Ω and two meshes with 30 and 132 elements, respectively.

Figure 11 .

 11 Figure 11. Example 2. Convergence history for p = 1, 2, 3, 4, in | •| (top-left), • Ω (top-right) and | •| s (bottom-left) against √ Dof for the meshes exemplified in Figure 10 with 30, 132, 555, 2151, 8337 elements, respectively. Bottom right: convergence history for p = 1, 2, . . . , 7 in the three norms against √ Dof for the mesh with 132 elements shown in Figure 10 (right).

Figure 12 .

 12 Figure 12. Example 3. Domain with holes and 531-element mesh (left). Corresponding solution obtained with p = 1 (right).

-x 2 1 u x 2 x 2 + u x 1 + 1 , u x 1 +π 2 x 3 1 12 )

 2111312 Aωπ cos(ωπx 1 )u x 2 + u = 0, in Ω Aωπ cos(ωπx 1 )u x 2 + u = 0, in Ω 2 ,coupled with inflow and Dirichlet boundary conditions, so that the analytical solution is given byu(x 1 , x 2 ) =    sin( π 2 (1 + x 2 -A sin(ωπx 1 ))) exp(-(x 1 + ), in Ω 1 , sin( π 2 (1 + x 2 -A sin(ωπx 1 ))) exp(-x 1 ). in Ω 2 .

Figure 13 .

 13 Figure 13. Example 3. Solutions computed for = 10 -4 (top) and = 0 (bottom) using linear elements. The meshes are composed of 129 (left) and 2048 (right) elements with compex shapes.

Figure 14 .

 14 Figure 14. Example 4. The domain Ω with the 64 element mesh fitted to the interface (blue) for A = 0.025, and ω = 8 (left) and ω = 16 (right). The inflow parts of the boundary are drawn in red.

Figure 15 .

 15 Figure 15. Example 4. Convergence history for p = 1, 2, . . . in the | •| and • Ω against √ Dof for the curved triangular mesh with 64 elements shown in Figure 14, for A = 0.025, ω = 8, 16.

  or x max ∈ K\ K0 . Let now K(x max ) be the pyramid with vertex x max and baseF 0 = [0, 1] d-1 . If x max ∈ K0 , then, we have v 2 L∞( K) ≤ 32 d p 2d v 2 K (see, e.g., [56, eq. (3.6.4)],) whereas if x max ∈ K\ K0 , we have, respectively, |v(x max