Finite element theory on curved domains with applications to discontinuous Galerkin finite element methods

In this paper we provide key estimates used in the stability and error analysis of discontinuous Galerkin finite element methods (DGFEMs) on domains with curved boundaries. In particular, we review trace estimates, inverse estimates, discrete Poincaré–Friedrichs' inequalities, and optimal interpolation estimates in noninteger Hilbert–Sobolev norms, that are well known in the case of polytopal domains. We also prove curvature bounds for curved simplices, which does not seem to be present in the existing literature, even in the polytopal setting, since polytopal domains have piecewise zero curvature. We demonstrate the value of these estimates, by analyzing the IPDG method for the Poisson problem, introduced by Douglas and Dupont, and by analyzing a variant of the hp‐DGFEM for the biharmonic problem introduced by Mozolevski and Süli. In both cases we prove stability estimates and optimal a priori error estimates. Numerical results are provided, validating the proven error estimates.


Abstract
In this paper we provide key estimates used in the stability and error analysis of discontinuous Galerkin finite element methods (DGFEMs) on domains with curved boundaries. In particular, we review trace estimates, inverse estimates, discrete Poincaré-Friedrichs' inequalities, and optimal interpolation estimates in noninteger Hilbert-Sobolev norms, that are well known in the case of polytopal domains. We also prove curvature bounds for curved simplices, which does not seem to be present in the existing literature, even in the polytopal setting, since polytopal domains have piecewise zero curvature. We demonstrate the value of these estimates, by analyzing the IPDG method for the Poisson problem, introduced by Douglas and Dupont, and by analyzing a variant of the hp-DGFEM for the biharmonic problem introduced by Mozolevski and Süli. In both cases we prove stability estimates and optimal a priori error estimates. Numerical results are provided, validating the proven error estimates.

INTRODUCTION
When modeling second-and fourth-order (as well as higher order) elliptic partial differential equations (PDEs), one may be required to consider a domain that cannot be expressed as a finite union of polytopes, for example, the unit ball, B 1 (0) ≔ {x ∈ R d :| x| <1} ⊂ R d . This necessity could be driven by the domain considered in the underlying application, where the domain is for example Lipschitz continuous, and piecewise C 1, , ∈ (0, 1), but not piecewise smooth, or for example the domain is C 1 , and thus not polytopal. Such domains arise naturally in the theory of PDEs, for example, a natural assumption for the Monge-Ampère equation [1][2][3][4][5][6] is that the domain is uniformly convex [2,3,5], and oblique boundary-value problems [7][8][9][10] in nondivergence form, with bounded and measurable coefficients, require a C 2 boundary assumption [9], both of which rule out the possibility of a polyhedral domain.
When it comes to finite element methods (FEMs), it is useful if the domain is polytopal, then since one can discretize the domain, Ω, exactly by polytopes, that is, there exists a family of shape-regular meshes ( h ) h>0 on Ω for which Ω = ∪ K∈ h K (the sets K are often d-simplices or parallelepipeds). If the boundary of Ω is curved, an exact mesh consisting of a finite set of polyhedrons cannot be obtained; one must instead use curved elements. In [11], the author introduces the concept of exact curved domain approximation by curved d-simplices, following [12,13], providing an optimal (with respect to the parameter h) finite element interpolant (interpolating with and without boundary conditions), with estimates in W m,p -norms, m ∈ N 0 , p ∈ [1, ∞].
We will see, however, that in order to design and analyze discontinuous Galerkin finite element methods (DGFEMs) for second-and fourth-order elliptic PDEs on domains with curved boundaries, one requires further estimates, in particular: inverse estimates; discrete Poincaré-Friedrichs' inequalities; simplicial curvature bounds; and optimal interpolation estimates in noninteger Sobolev norms.
One is often motivated to use DGFEMs, other nonconforming FEMs, and mixed FEMs over conforming FEMs, due to the structural and computational challenges that conforming FEMs impose. For conforming FEMs, it is required that the approximation space is a subset of the space of weak solutions to the PDE, examples of this being the spaces H 1 0 (Ω) and H 2 0 (Ω) for second-and fourth-order elliptic problems, such as the Poisson problem and biharmonic clamped plate problem, which we shall consider as our model second-and fourth-order problems. In the H 1 0 (Ω) case, this can be achieved by considering piecewise polynomials that are globally continuous, however, for H 2 0 (Ω), one must also enforce continuity of the gradient across neighboring elements. An example of this being the Argyris finite element [14], which can be rather expensive to implement, requiring polynomials of degree five on two dimensional simplicial polynomials. In contrast, nonconforming methods weakly enforce this regularity by penalizing jumps of the discrete functions, and their derivatives across the edges of neighboring elements, and as a result, the methods that we consider only require a polynomial degree greater than or equal to the number of derivatives in the weak formulation of the PDE. Furthermore, since curved domain approximations require the composition of piecewise polynomials with functions that are not piecewise polynomials (the details of this will be made clearer in Section 3), applications of the chain rule show that in general, the piecewise derivative no longer maps from the finite element space into itself (as is often seen in discontinuous Galerkin [DG] finite element spaces), complicating the derivation and structure of inverse estimates. For penalty FEMs for fourth-order problems, we will see that this leads to the necessity of discrete Poincaré-Friedrichs' inequalities.
When defining a finite element space, it is not always necessary to utilize the composition of polynomials with nonaffine maps. For example, in [15] an hp-DGFEM is proposed that allows the computational domain to be curved, instead defining the finite element functions to be piecewise polynomials on the (potentially) curved elements, instead of on the reference element. The authors prove that the scheme is well-posed and approximates the true solution optimally, provided the physical elements themselves satisfy mild curvature assumptions. It is mentioned in [15] that a given mesh may have to be further refined, so that such a mesh assumption can hold, however, this is akin to the validity of one of the necessary assumptions of this paper (see (3.24) in Assumption 3.21), which could require further refinement of the initial mesh. In contrast, in our case, such a refinement would only be necessary on the elements corresponding to a nonaffine mapping from the reference element. Furthermore, it is discussed in [15] the computational challenges that may arise from choosing functions that are polynomials on the physical curved element.
One can also consider inexact domain approximation. For example, the authors of [16] propose a virtual element method (VEM) for the Poisson problem on curved domains, where the element is curved in a polynomial manner, and thus the mesh does not always approximate the domain exactly. The authors prove that the error arising in approximating the boundary does not dominate the overall approximation of the true solution, and in fact the method approximates the true solution optimally (in terms of the mesh size). The VEM generally requires that any bilinear forms can be calculated exactly, using the degrees of freedom of the approximation space. It is not clear in the literature that such formulations are extendable to nondivergence form elliptic equations with L ∞ (Ω) coefficients (whereas the results of the current paper are shown to be of benefit to [17,18]); in particular, the discrete spaces are typically augmented by (theoretically) solving the PDE locally on the physical elements. This consideration is likely to be nontrivial in the setting of nondivergence form equations. Furthermore, exact domain approximation is useful in the design of numerical methods. For example, in the works [17][18][19] a discrete analogue of the Miranda-Talenti estimate [9] motivates the design of the methods, and due to this, the exact approximation of the domain is key to the stability of the methods.
The remainder of this paper is organized as follows: in Section 2 we shall discuss the existence and uniqueness of weak solutions to the Poisson and biharmonic equations, and discuss conforming FEMs, and DGFEMs (the latter of which falls into the category of nonconforming FEMs) on polytopal domains, with the goal of highlighting important features, such as the stability and consistency of such schemes. In Section 3 we review the key tools from finite element analysis that are well known in the polytopal case, in the context of curved simplicial finite elements. In Section 4 we will provide the numerical methods for the Poisson and biharmonic problems, and prove that they are stable, and in Section 5 we prove that the numerical solutions satisfy optimal a priori error estimates in H k -type norms. Finally, in Section 6 we provide numerical experiments that validate the error estimates of Section 5.
Let Ω ⊂ R d be Lipschitz continuous, and consider the following second-and fourth-order elliptic boundary-value problems, for k = 1, 2, find u k : Ω → R such that: where f ∈ L 2 (Ω). When k = 1, (2.1) is the well-known Poisson problem, and for k = 2, (2.1) is the biharmonic clamped plate problem. In particular, one can show that in each case, there exists a unique weak solution u k ∈ H k 0 (Ω). That is, u k satisfies where the bilinear forms a k ∶ H k 0 (Ω) × H k 0 (Ω) → R are given by Note that the existence of such functions follows from applying the Lax-Milgram theorem [20]; in particular one must show that the bilinear forms are coercive in the H k -norm. In the case that k = 1, this follows from the Poincaré inequality [20], and for k = 2, the following identity (see (4.1)) , which, coupled with the Poincaré inequality, also proves the coercivity of a 2 .
The derivation of the weak formulations (2.2) follows from the following integration by parts identities, valid for functions u, v ∈ C ∞ (K), where K ⊂ R d has a Lipschitz boundary, and extendable to u, v in suitable Sobolev spaces by density: and where n K is the unit outward normal to K. Taking K = Ω, the choice of u, v ∈ H k 0 (Ω) justifies the lack of the appearance of boundary integrals in (2.3) and (2.4) (however, for this we utilize the density of C ∞ c (Ω) in H k 0 (Ω)). For a conforming FEM, one assumes that the finite dimensional space V k,h ⊂ H k 0 (Ω), and so one may obtain a conforming FEM by directly substituting the finite element functions into the bilinear forms. That is, one seeks u k,h ∈ V k,h such that Indeed, since V k,h ⊂ H k 0 (Ω), the properties of the bilinear forms are still valid on V k,h × V k,h , and so the existence and uniqueness of a numerical solution follow in a similar manner to the existence and uniqueness of a weak solution. In particular, the bilinear form a k is coercive on V k,h × V k,h in the H k (Ω) norm, and so we obtain the stability estimate where C k is a positive constant independent of the approximation parameter h. Since the problem (2.7) is equivalent to solving a linear system of equations, the stability estimate implies uniqueness, which in turn implies invertibility of the matrix describing the corresponding linear system, which also yields existence. Furthermore, we see that the true solutions, u k ∈ H k 0 (Ω) satisfy and so A FEM that satisfies (2.9) is called consistent, and (2.10) is referred to as Galerkin orthogonality, which, when combined with the stability estimate (2.8), yields Cea's lemma: One obtains optimal error estimates, by noting that the infimum over V k,h is bounded above by any choice of z h ∈ V k, h . In particular, assuming that u k ∈ H s (Ω) ∩ H k 0 (Ω), s ≥ k, one may choose z h to coincide with a suitable interpolant, yielding where the constant C * is independent of h. Unlike conforming FEMs, where the approximating space V k,h is a subset of H k 0 (Ω), nonconforming FEMs involve approximating spaces for which this is not true; in the case of DGFEMs one only has V k, h ⊂ L 2 (Ω), and for the C 0 -interior penalty method proposed in [21], one has V 2,h ⊂ H 1 0 (Ω), which is nonconforming in the sense that V 2,h is not contained in H 2 0 (Ω). For DGFEMs, one also has analogues of stability, consistency, and optimal error estimates. However, since the finite element functions do not have sufficient global regularity, one cannot directly substitute u h , v h ∈ V k,h ≔ {v ∈ L 2 (Ω) ∶ v| K ∈ P p (K) ∀K ∈  h } (implicitly, we assume p ≥ k, and that ( h ) h>0 is a family of regular simplicial meshes on Ω) into the bilinear forms a k , k = 1, 2.
Such functions do, however, satisfy a property of piecewise regularity; since, } (in particular, piecewise polynomials are piecewise smooth) and so, assuming u k ∈ H k 0 (Ω) ∩ H 2k (Ω) are the weak solutions to the PDE, we can sum the integration by parts identities (2.5) and (2.6) over all K ∈  h , obtaining (see Definition 3.1, as well as (3.2) and (3.1) for the relevant notational conventions in present in the identities that follow): where n F denotes a fixed choice of unit normal to F. Let us define 18) we arrive at the following DGFEMs for the approximation of the solutions u k , k = 1, 2, to (2.1): find Furthermore, we see that but, the remaining terms present inÃ 1 andÃ 2 are not bounded quite as simply. If F is a face of K ∈  h , trace estimates yield for w k ∈ H 2k (K) , where C depends upon the shape-regularity constant of  h . Then, applying inverse estimates [22] of the form for w ∈ P p (K), gives us so long as w k ∈ P p (K) ⊂ H 2k (K), k = 1, 2. Then, utilizing (2.24)-(2.26), and the Cauchy-Schwarz inequality with a parameter, yields the following for any 1 > 0, and where the final inequality holds due to the fact that the number of elements that share a given face is bounded in terms of the dimension, d. Similarly, for any 2 > 0, and any v h ∈ V 2, h , we see that . (2.28) The above estimates lead one to supplement the bilinear formsÃ k , k = 1, 2, with additional bilinear forms S k , J k : V h × V h → R, where the bilinear forms J k penalize interface jumps of the inputs and their piecewise weak derivatives up to order 2k − 1 across interior faces, and up to order k − 1 on boundary faces, and the bilinear forms S k preserve the symmetry of the scheme. Clearly, the choice of J k and S k lead to different FEMs; in [23] the authors present and analyze nine DG methods from [24][25][26][27][28][29][30][31][32] for the Poisson problem (k = 1), and in [33] the first interior penalty discontinuous Galerkin method for the Biharmonic problem (k = 2) was proposed and analyzed. Furthermore, in [34] an hp-FEM is introduced for the Biharmonic problem (k = 2) with symmetric and nonsymmetric penalties. For other examples of nonconforming methods for second-and fourth-order elliptic problems see [17-19, 21, 35-38]. Thus, we may take where j F , j = 1, 2, 3, are positive parameter choices independent ofh F , that are chosen sufficiently large, in order to compensate for the jumps across F ∈  i,b h present in estimates (2.27) and (2.28), as well as the jump estimates resulting from the terms included for symmetry that are present in S 1 and S 2 (these terms are bounded in exactly the same manner as in the derivation of estimates (2.27) and (2.28)).
By (2.14) and (2.15), we see that J k (u k , v h ) = 0 for all v h ∈ V h , and so the bilinear forms are also consistent, that is, they satisfy (2.20); furthermore, they are symmetric. These particular choices of J k (and thus A k ) coincide with the interior penalty discontinuous Galerkin (IPDG) method of [26] (k = 1) and the h-version of the symmetric hp-DG method of [34], with the parameters Analogously to deriving (2.27) and (2.28), one can show the following stability estimates [23,34].
where the norms, ‖⋅‖ h, 1 , and ‖⋅‖ h, Δ are defined by 31) and the constants C *, 1 and C *, Δ depend only on the dimension, the domain Ω, the polynomial degree, and the shape-regularity constants. These estimates of course yield existence and uniqueness of u k, h satisfying for some polynomial , unless F K (and thus F −1 K ) is affine, that is, the mesh is polytopal. This leads one to obtain estimates of the form which would not directly lead to the derivation of the stability estimate (2.30) of A 2 (since we are no longer able to estimate in the ‖⋅‖ h, Δ -norm, as the Laplacian structure is no longer preserved). This leads us to define a new variant of A 2 with the goal of replacing the inner product leading to coercivity in the norm . In order to achieve such a stability estimate, one is required to prove a discrete Poincaré-Friedrichs' inequality, in order to bound the H 1 -terms of the right-hand side of (2.32) by H 2 terms, and factors that are present in J 2 (⋅, ⋅).
Finally, we discuss error estimates. Since the methods are consistent, one has and thus, for any z h ∈ V k, h , the triangle inequality yields where t k K ≔ max{p + 1, s k K }; in the case of quasiuniform meshes, the above becomes , and so the estimate is optimal with respect to the mesh size. For k = 1, the estimate is provided in [23] for the case that s 1 K = 2 for all K ∈  h , that is, the integer case, and for k = 2, the estimate (2.37) is provided in [34]. In the case of curved finite elements, the method for proving optimal error estimates is the same (except there are a few more terms that we must estimate), however, one still requires a suitable interpolate. In the context of (2.37), this means that there is an element z h ∈ V k, h , uniquely determined by a function w k ∈ H s k (Ω;  h ), such that for each K ∈  h , each integer 0 ≤ q ≤ min{p, 2k − 1}, and each multi-index , with 0 ≤ | | ≤ q, , (2.38) where C may depend upon the polynomial degree, Ω, and the shape regularity constant, but is independent of h K . A goal of the proceeding section will be to prove (2.38) in the curved case, which will yield optimal error estimates for both the schemes we propose, and, since the polytopal case can be viewed as a special case of the curved case, we will provide optimal estimates for the IPDG method of [26] for the Poisson problem in noninteger Sobolev norms. The first estimate of (2.38) is proven in [11] for the case that s k K is integer valued, we aim to provide such estimates in H s -norms, for noninteger s.

CURVED DOMAIN APPROXIMATION AND FINITE ELEMENT ESTIMATES
We will begin this section by providing the details of [11], which provides us with a notion of exact domain approximation, along with essential scaling arguments that allow us to prove the desired trace and inverse estimates. Such estimates will allow us to prove that our proposed FEMs are stable, yielding existence and uniqueness of numerical solutions. This requires the following notation.

Notation
Definition 3.1 (Face and vertex sets). Given a mesh  h , we denote by  i,b h , the set of faces of  h , by  i h the set of interior faces of  h , and by  b h , the set of boundary faces.

Definition 3.2 (Jump and average operators). For each face
, with corresponding unit normal vector n F which, for convention, is chosen so that it is the outward normal to K, we define the jump operator, ⟦⋅⟧, over F by and the average operator, ⟪⋅⟫, by For an element K, we define the inner product ⟨⋅, ⋅⟩ K by Any ambiguity in this notation will be resolved by the arguments of the bilinear form.

Definition 3.4 (≲ and ≈ symbols). Herein we write
and u, but otherwise possibly dependent on the polynomial degree, p, the shape-regularity constants of  h , C  , and d.

Curved simplices
The ability to define a nonaffine approximation of a domain, Ω ⊂ R d relies upon the Ω satisfying a notion of piecewise regularity, which motivates the following definition.
we may express the boundary of Ω, Ω, as a finite union where each Γ n ⊂ R d is of zero d-dimensional Lebesgue measure, and admits a local representation as the graph of a uniformly C k function. That is, for each n, and at each x ∈ Γ n there exists an open neighborhood V n of x in R d and an orthogonal coordinate system (y n 1 , … , y n ), such that V n = {(y n 1 , … , y n ) ∶ −a n j < y n j < a n j , 1 ≤ j ≤ }; as well as a uniformly C k function n defined on V ′ n = {(y n 1 , … , y n −1 ) ∶ −a n j < y n j < a n j , 1 ≤ j ≤ − 1} and such that | n (y n ′ )| ≤ a n ∕2 for every y n ′ = (y n 1 , … , y n −1 ) ∈ V n′ ,

Definition 3.6 (Curved d-simplex). An open set K ⊂ R d is called a curved d-simplex
if there exists a C 1 mapping F K that maps a straight reference d-simplexK onto K, and that is of the form is an invertible map and Φ K ∈ C 1 (K; R ) satisfies where ‖ ⋅ ‖ denotes the induced Euclidean norm on R d × d .
. Given a curved d-simplex K, with the associated straight reference d-simplexK, and map F K ∶K → K, with F K =F K + Φ K , we define the associated straight d-simplex: Remark 3.8 The associated d-simplex,K, is a straight d-simplex that "approximates" K.
Proof. See Remark 2.3 of [11]. ▪ Remark 3.10 (Affine mesh). In the case that the domain has a flat boundary, one employs an affine approximation of the domain, in which case, the corresponding functions Φ K in (3.5) are all zero.
Lemma 3.14 The mapping F K is a C 1 -diffeomorphism fromK onto K and satisfies

Lemma 3.15
Let us denote by c , 2 ≤ ≤ m, m ∈ N, the constants There exist constants c − , 2 ≤ ≤ m, depending continuously on c K , c 2 (K), …, c m (K), such that

Scaling arguments
Lemma 3. 16 Assume that K is a curved d-simplex of class C m . Let l be an integer, where the constants C depend continuously on c K , c 2 (K), …, c m (K).

Lemma 3.17
Assume that K is a curved d-simplex of class C m , and that F is a face of K; we denote byB F the restriction ofB K toF ≔ F −1 K (F). Let l be an integer, 1 ≤ l ≤ m, s ∈ [0, l − 1/2). Then, for any v ∈ H l (K), the function F (v) belongs to H s (F), and we have where the constant C depends continuously on c K , c 2 (K), …c m (K).
A key tool in the derivation of optimal interpolation estimates on affine meshes is the following scaling argument (see Theorem 3.1.2 of [14]): Here, we are considering the affine equivalent straight d-simplicesK andK, and an invertible affine One can see that (3.15) and (3.17) are similar. The main difference is the presence of the lower order seminorms on the right-hand side of (3.15).
Let us look at the particular example of the H 2 -seminorm when F K is not affine. The chain rule, and the multivariable change of variables formula yields A sufficient assumption that yields an estimate of the same order as (3.17) with l = p = 2 (in terms of ‖B K ‖), is to assume that c , given by (3.12), is uniformly bounded for = 2. This coupled with the fact that C K < 1 gives us sup Applying the above to (3.18) yields In order to appropriately bound the determinant term, one must note that Ultimately, this gives us This motivates the two following definitions, generalizing the prerequisite assumptions, allowing one to obtain analogous estimates in higher order seminorms.

Definition 3.18
The family ( h ) h of meshes is said to be regular if there exist two constants, and c, independent of h, such that, for each h, any K ∈  h satisfies where K is the diameter of the sphere inscribed inK. Furthermore, we have Remark 3.19 Condition (3.20) is referred to as nondegeneracy (e.g., in [22]).

Definition 3.20
The family ( h ) h of meshes is said to be regular of order m if it is regular and if, for each h, Assumption 3. 21 We assume that any two elements sharing a face have commensurate diameters, that is, there is a C  ≥ 1, independent of h, such that for any K and K ′ in  h that share a face. Finally, we assume that each F ∈  b h satisfies for some n ∈ {1, …, N}, with Γ n given as in (3.4). This implies that each boundary face is completely contained in a boundary portion Γ n , as well as ensuring that our approximation of the domain Ω is exact.
Remark 3.22 The assumptions on the mesh given by Assumption 3.21, in particular (3.23), show that if F is a face of K, then A final, necessary step, before providing optimal interpolation estimates and inverse estimates for (continuous and discontinuous) curved Lagrange finite element spaces, is to relate the estimates of this section to the local mesh size, h K . The general rule of thumb in this context is that ‖B K ‖ is of order This notion is made more concrete by the following theorem from [14].

Corollary 3.24
Assume that the family ( h ) h of meshes satisfies (3.20). Then, there exists a positive constant C depending only on , such that for any K ∈  h with an associated straight elementK, that Proof. See [39], Corollary 3.24. ▪

Lagrange finite element spaces
The finite element spaces we consider in this paper consist of discontinuous piecewise polynomial functions, which fall into the class of discontinuous (curved) Lagrange finite element spaces. In general, a finite element is a triple (K, P K , Σ K ) where K is a subset of R d , P K is a finite dimensional space on K, and Σ K is a set of continuous linear forms on P K , which we will call the degrees of freedom.
In the context of Lagrange finite element spaces, the continuous linear forms are given by (local) point evaluations. In the simplicial case, the placement of these points is naturally described using the barycentric coordinates of the simplex.
Definition 3.26 (Barycentric coordinates). Given a straight d-simplexK, with verticeŝ a 1 , …,â +1 ∈ R , we define the barycentric coordinates ofK, 1 , … , , +1 via the following (invertible) system and for any ∈ J(p), we associate the pointâ ∈K with barycentric coordinates i = i /p, i = 1, …, d + 1. Then, we call (K,P K ,Σ K ) a straight Lagrange finite element of type p, whereP , for f ∈P K , and we recall that P p (K) is the space of all polynomials with total degree less than or equal to p.
where (K,P K ,Σ K ) is a straight Lagrange finite element of type p.
Definition 3.29 (Discontinuous Galerkin finite element space). The discontinuous Galerkin finite element space V h, p is defined by Piecewise polynomial functions naturally satisfy a property of piecewise regularity. This is accurately captured by considering the notion of broken Sobolev spaces.
Higher order discrete derivatives are defined in a similar way. We define a norm on W s,r (Ω;  h ) by with the usual modification when r = ∞.
Definition 3. 32 We define the following for with the usual modification when r = ∞. Note that |⋅| W s,r * (Ω; h ) is a norm when s = 0, and a semi-norm when s ∈ N. We also define |⋅| H s * (K) and |⋅| H s * (Ω; h ) in the usual way.
Remark 3. 33 We can use these semi-norms to equivalently phrase estimates such as (3.19), which can now be written as We now provide trace and inverse estimates that we will be utilized frequently. In particular, the noninteger order trace estimate will be utilized in proving the second estimate of (2.38).

3.5
Trace and inverse estimates Lemma 3.34 Assume that  h is a regular mesh on Ω. Then, for any K ∈  h , we have that Proof. From the multivariable change of variables formula, we obtain Under a second change of variables, we obtain where F * is a (d − 1)-face of a fixed reference d-simplex, K * , and G K ∶ K * → K, G K (x * ) ≔ÃKx * +ãK, withÃK ∈ GL(R ),ãK ∈ R , andÃF is the restriction ofÃK to F * = G −1 K (F). Since the trace operator is continuous from H r (K * ) → L 2 (K * ) for r > 1/2 [40], we see that where 1 and 2 are positive, continuous functions that we will soon provide.
Recall the definition of the H r -semi norm: We note that sincex 1 ,x 2 ∈K, which, when applied to (3.44), gives us We apply the multivariable change of variables formula once more, obtaining Of course, we also have We obtain the functions 1 and 2 in a similar manner, except since G K is affine, the scaling argument is simpler, and we have that From the nondegeneracy condition (3.20), it follows (from the proof of Theorem 4.4.20 in [22]) that the collection of the invertible matrices given by the affine maps from Overall, we have obtained where the final inequality follows from (3.28 where the positive constant C I is independent of K and h K . Proof. We first note that (3.47) is trivial when m = 0, since then s = 0, and |⋅| W m,q * = |⋅| W s,q * = ‖⋅‖ L q , so we will assume that m ≥ 1. We will first prove (3.47) when s = 0. By (3.14), for j ∈ N, 1 ≤ j ≤ m, q ∈ {2, ∞}, and any K ∈  h , we have ) .
( 3.48) Now, for 0 ≤ r ≤ j, where the inequality is due to (3.28). Now, let K * be a fixed reference element, and takẽ GK ∶ K * →K, withGK(x * ) =ÃKx * +ãK, withÃK ∈ GL(R ) andãK ∈ R . As in the proof of Lemma 3.35, it follows thatÃK belongs to a compact subset BL of GL(R d ). Now, defining v * (x * ) =v(GK(x * )), it follows that v * ∈ P p (K * ), where P p (K * ) is of finite dimension, depending only on K * , d and p, thus by the equivalence of norms on finite dimensional spaces, we see that (3.49) Thus, applying the above inequality, (3.15) with l = 0, and (3.28), to (3.48), we obtain Since our choice of 1 ≤ j ≤ m was arbitrary, we may take 1 ≤ k ≤ m, and sum the above over 1 ≤ j ≤ k, obtaining We obtain (3.47) with s = 0, by setting k = m above. We will now prove (3.47) for 1 ≤ s ≤ m.
In this case we will argue by induction, and as our base case, we shall prove the result for s = 1. Take 1 ≤ j ≤ m, and let | | = j. Then we may write D v = D (D v) for some | | = j − 1, | | = 1. One must note that by the chain rule, where the components of (Dv•F −1 K )DF −1 K do not necessarily belong to P p (K). It is the case, however, that Dv ∈ P p (K) for any | | = 1. One can see that where we are denoting We also have that Applying (3.52), (3.53), and (3.54) to (3.51), and summing over all | | = j, we obtain Lastly, applying (3.9) and (3.28) to the above estimate, we obtain (noting that  h is regular of order m) Again, our choice of 1 ≤ j ≤ m was arbitrary, and so we can sum (3.55) over 1 ≤ j ≤ k for any 1 ≤ k ≤ m, obtaining To proceed to argue by induction, we will assume that for 1 ≤ s ≤ k ≤ m − 1 we have and we will use this to show that To this end, let us take s + 1 ≤ j ≤ m and let | | = j. Again we write D v = D (D v) for some | | = j − 1, and | | = 1, and so, analogous to our previous argument, we obtain Applying our inductive hypothesis (3.56) with k = j − 1 ≥ s, we obtain Applying the above to (3.57), we obtain Let us momentarily assume that, for any 1 ≤ s ≤ m − 1, , where | | = j, and s + 1 ≤ j ≤ m was arbitrary. Summing over all | | = j, and then all s It is also clear that ∑ min{1,s}≤| |≤s and so we obtain s + 1 ≤ k ≤ m, which concludes our inductive argument, and yields (3.47) for 1 ≤ s ≤ m, by taking k = m. It remains to show that (3.58) is in fact true. Let us recall the formula [11].
, (3.59) where E(r, i) is the set given by: and the c 's, ∈ E(m, r) are some given constants, bounded independently of h K . From this, we obtain where the final inequality follows from (3.13), and the fact that the mesh is regular of order m ≥ s + 1. Applying (3.28), and noting that by definition, if ∈ E(r, i), then | | = i and ∑ r l=1 l l = r, we obtain as desired. Note that the estimates we have derived are independent of the choice of K ∈  h . ▪

Interpolation estimates
The proofs of the following lemmas can both be found in [11], that is,

61)
where the constant C is independent of h K , u, and K.
Remark 3. 38 We note that the classical Lagrange interpolation operator can only be applied functions that have well defined point values. Even in two dimensions, it is not in general the case that functions in H 1 have well defined point values. This leads one to define other interpolation operators that require less regularity, in particular, we define a local interpolation operator that is well defined on L 2 functions (one of the first examples is due to Clemént [41], using local averaging; however the one we will define is provided in [11] and is slightly different).

Definition 3.39 (Local L 2 projection).
For v ∈ L 2 (Ω), and K ∈  h , we definêv to be the unique element of P p (K) that satisfies

Lemma 3.41 (H r -multipliers).
Assume that u ∈ H r (Ω), 0 < r < 1 and ∈ C 0,1 (Ω). Then, there exists a constant C depending only on d and r, such that Then, there exists a z k, h ∈ V h, p , p ≥ 2k − 1, and a constant C, independent of u k , and h K , but dependent on max K s k K , such that for each K ∈  h , each nonnegative integer q ≤ 2k − 1, and each multi-index with | | = q, we have Proof. We will first discuss how we will obtain the second bound of (3.65). Let k ∈ {1, 2}. We either have that q ≤ 2k − 1 is a nonnegative integer, and t k K −q > 2k−1∕2−q ≥ 1∕2, or we have that q ≤ 2 ≤ p, and t k K − q = min{p + 1, m + 1, s 1 K } − q > 5∕2 − 2 = 1∕2. Thus, under the hypotheses of the lemma, for k = 1, 2 we have that t k K − q > 1∕2. Since the family of triangulations is regular of order m, it follows that for any such that | | = q, and any v ∈ V h, p , that D (u − v) ∈ H t k K −q (K). In particular, t k K − q > 1∕2. Thus, we may apply the trace estimate (3.42) with Let us assume that there exists a z h ∈ V h, p satisfying the first estimate of (3.65). Then, setting v = z h above we obtain (3.66) Thus, to obtain both estimates of (3.65), it suffices to prove that the there exists a z h ∈ V h, p such that the first estimate of (3.65) holds, as well as the following: Since, applying the above estimate to (3.66), and noting the factor h r K K in the second inequality of (3.66), we obtain the second estimate of (3.65). Note that we already have such bounds in the case that s K is an integer, and as such, we shall assume from this point on that s K ∉ N.
We will now prove the first estimate of (3.65). Let satisfy | | = q, and let K * be a fixed reference simplex. Then, from (3.14) we obtain We take the function z h ∈ V h, p , defined as follows: z h | K = Π h u| K where Π h is the local interpolation operator, given by (3.63). Due to (3.62), this operator reproduces polynomials in P p (K), and so we may apply Theorem 5 of [42] in conjunction with Theorem 1.8 of [43] (applying Theorem 1.8 of [43] allows us to consider noninteger Sobolev spaces when applying the Bramble-Hilbert lemma), obtaining for min{1, q} ≤ j ≤ q where by assumption t k K > q (note that the final inequality follows from a scaling argument similar to the one used in estimate (3.43), noting that K * andK are affine equivalent, and the mesh is shape regular). We now decompose t K = K + r K , where K ≥ k is an integer, and r K ∈ (0, 1). We see that Recalling formula (3.59), and applying the triangle inequality, we obtain .
We now apply (3.64) to the above estimate, obtaining (noting thatK is contained in the unit ball, and thus diam(K) ≤ 2) By (3.9), and the fact that the triangulation is regular of order m ≥ k + 1 (and that N ∋ K < t K ≤ m + 1, so K ≤ m), we estimate T 1 as follows For the second term, we see that Since the triangulation is regular of order m ≥ k + 1 (and K + 1 ≤ m + 1), applying (3.9) above yields the following for ∈ E( K , i) Applying (3.73) to (3.72) gives us We now apply (3.71) and (3.74) to (3.70), obtaining Applying the change of variables formula in the L 2 -norms in (3.75), and the scaling argument (3.45) and (3.46) to the |D i u•F K | H r K (K) term for i = K (noting that this argument is valid for any r K ∈ (0, 1), as long as the function has H r K -regularity) in (3.75), in conjunction with (3.28), we obtain where the constant C is independent of h K and the choice of K ∈  h (note that we have utilized the continuous embedding H s K (K) ⊆ H t K (K), where the constant in the embedding only depends upon d and r K , due to Proposition 2.1 of [44]). We note, however, that the terms of the sum on the right-hand side of the final inequality of (3.76) are not present in the H t K -norm. Furthermore, for 1 ≤ i < K , we note the following: for any M ∈ [P 0 (K * )] dim(D i u) , where the first inequality follows from a scaling argument, and the fact that the mesh is regular, and the final equality holds due to the fact that constant functions are in the kernel of |⋅| H r . We now use the fact that the embedding H 1 (K * ) ⊆ H r K (K * ) is continuous, obtaining where the penultimate inequality follows from an application of Theorem 1.8 of [43], and the final inequality follows from the fact that the mesh is regular. Thus, we obtain Applying the above to (3.76) gives us Finally, applying (3.79), (3.69), and (3.28) to (3.68), we obtain which is the first estimate of (3.66). Estimate (3.67) is obtained in a similar manner, utilizing (3.13). ▪

Discrete Poincaré-Friedrichs' inequalities Lemma 3.43 (Discrete Poincaré-Friedrichs' inequality). Assume that { h } h is regular of order 2 family of triangulations, and let v ∈ V h, p . Then, the following inequality holds
where the positive constant, C, depends only on the shape-regularity constants of the mesh, d, and Ω.
Proof. Let K ∈  h , and take v ∈ V h, p . We see that , subtracting (1/2)∫ K v 2 from each side and multiplying by 2 yields Summing the above over all K ∈  h , and denoting n F to be a fixed choice of unit normal to F ∈  i,b h , we obtain for any > 0. We then apply the trace inequality (3.41), obtaining ∑ Applying the above estimate to (3.81), we obtain, for any > 0, .

Lemma 3.44 (Gradient Poincaré-Friedrichs' inequality). Assume that { h } h is regular
of order 2 family of triangulations, and let v ∈ V h, p . Then, the following inequality holds where the positive constant, C P , depends only on the shape regularity constants of the mesh, d, and Ω.

Tangential operators and curved simplex curvature bounds
In order to appropriately define and bound the bilinear forms that define our method, we need to be able to define tangential differential operators (i.e., operators that involve derivatives that are tangential to the faces of the curved simplices K of the mesh), and bound the curvature terms arising in the bilinear form (these curvature terms appear both on boundary faces, and on interior faces if the dimension d ≥ 3). Tangential differential operators. For F ∈  i, b , denote for s > 1/2 the space of H s -regular tangential vector fields on F by H s Below we define the tangential gradient ∇ T ∶ H s (F) → H s−1 T (F) and the tangential divergence div T ∶ H s T (F) → H s−1 (F), where 1 ≤ s ≤ 2 (note that in the case that Ω is piecewise C m , with m ≥ 2, we are able to consider 1 ≤ s ≤ m).
We see that F ⊂ K, for some K ∈  h . Since K is piecewise C 2 (see the proof of Lemma 3.46), for a.e. x ∈ K, there exists a neighborhood W x of x in K, sufficiently small to allow the existence of a family of C 2 curves that satisfy the following: a curve of each family passes through every point of W x , and the unit tangent vectors to these curves form an orthonormal system (assumed to be oriented with respect to n, where n is the unit outward normal to K) at every point of W x . We take the lengths s 1 , …, s d − 1 along each of these curves, respectively, to be the local coordinate system, and denote t 1 , …, t d − 1 to be the unit tangent vectors along each curve, respectively. In this notation, we have the following for v : For ∈ C 1 (K), and ∈ C 1 (K) , with | K = ∑ −1 j=1 j t j , we obtain which extend to ∈ H s (K), s > 3/2, by density and the construction of the trace operator. Furthermore, one can see that by rearranging the first identity of (3.84), that ∇ T = ∇ − n n (and thus div T ) is well defined a. e. on K, and is independent of the choice of normal n.
We approach (3.84) and (3.85) in the context of traces and Sobolev spaces, in the following lemma. In particular we are able to decompose the Laplacian, Δ, in terms of the tangential Laplacian Δ T ≔ div T ∇ T , the mean curvature of the face, and first-and second-order normal derivatives.

Lemma 3.45
Let Ω be a piecewise C 2 domain, and let { h } h>0 be a family of meshes on Ω that is regular of order 1 and satisfies Assumption 3.21. Then, for any h > 0, for each K ∈  h and each face F ⊂ K, the following identities hold: for all v ∈ H s (K), s > 5/2, where n F is a fixed choice of unit normal to F,  F ≔ ∇ T ⋅ n F is the mean curvature of the face F, and F is the trace operator from K to F.
Proof. Let us take U ∈ C 3 (K), and for F ∈  i,b h , let u = U| F . Then, as the family of meshes { h } h>0 is regular of order 1, it follows that F ⊂ K for some K ∈  h , where K is piecewise C 2 (see the proof of Theorem 3.46). Thus, we may extend (without relabeling) the unit normal to F, n F (note that this choice of unit normal is fixed, and that (3.86) and (3.87) are independent of this choice), by n F ∈ C 1 (K) (note that the extension may not be normal to the other faces of K, when restricted there), and so also define an extension of the tangential gradient, ∇ T ∶ C 3 (K) → C 1 (K) , by This can be rearranged to yield Upon restricting to F, we obtain Thus, by density and the construction of the trace operator, this extends to u ∈ H s (K), s > 3/2, giving us For the identity (3.87), we follow a similar approach to [45], in which the statement is essentially proven for d = 2, 3. Now, for x ∈ F let us take a local coordinate system s 1 , …, s d − 1 , on a neighborhood W x of x in F. Expressing F locally as the graph of a C 2 function , we see that Furthermore, let us assume that the coordinates have been chosen so that ∇ s ′ (0) = 0 (denoting s ′ = (s 1 , …, s d − 1 )), so that the local coordinates {s ′ , s d } = {s ′ , (s ′ )} are tangent to the hyperplane {s d = 0} at x = (0, (0)). Then, in W x , we have that where U j , U jk denote the first-and second-order partial derivatives in the j and j, k components of U, respectively. Thus, at x, that is, at s ′ = 0, we have This decomposition is valid at any x ∈ F, and so we obtain Thus, by density, applying (3.88), for u ∈ H s (K), s > 5/2, we obtain Proof. First, let us assume that F ∈  b h . Then, since Ω is piecewise C 2 , F ⊂ Γ n ⊂ Ω, where Γ n is a C 2 portion of Ω. It then follows that for a given F ∈  b h , n F is of class C 1 (F). For any two vector-valued functions 1 , 2 : F → R d tangent to F, it then follows that Thus, for an arbitrary where the constant above depends on Ω, as the portions Γ n are determined by Ω. If F ∈  i h , then we may express F locally as the graph of a function determined by one of the maps F K that make up the mesh  h ; we also have that F K ∈ C 2 , as the family of meshes is regular of order 1. That is, since F ⊂ K for some K ∈  h , there exists a (straight) reference faceF, such that F = F K (F). Furthermore, there exists a straight approximating faceF =F K (F) (F K is the affine part of F K ), which provides us with a local coordinate system. AsF is flat, after a suitable change of coordinates, one has thatF ⊂ {(x ′ , 0) ∶ x ′ ∈ R −1 }. Furthermore, without loss of generality, we may assume that F does not intersectF. Let us denoteF ′ = {x ′ ∈ R −1 ∶ (x ′ , 0) ∈F}. Now, defining F K ∶F ′ → R by are of class C 2 , and furthermore, for any K ∈  h , K may be expressed as the finite union of the closures of F ∈  i,b h , and thus for all K ∈  h , K is piecewise C 2 . Furthermore, expressing F as the zero level set of the function h F K (x ′ , x ) = x − F K (x ′ ), we see that , for any two vectors 1 , 2 : F → R d tangent to F (with components k 1 , … , k , k = 0, 1), and hence orthogonal to n F , it follows that ( 1 ) T ∇ T n T F 2 = ( 1 ) T ∇n T F 2 . Furthermore, denoting ij ≔ 1 − ij , where ij is the Kronecker-delta symbol, we see that and so One also has that where the final inequality follows from (3.22), and C int is independent of both h, and the choice of F, since the family of meshes is regular of order 1. Thus, defining

93)
for all v ∈ H s (K), s > 5/2, where F ⊂ K, and F is the trace operator from K to F.
Then, by definition, we see that Furthermore, let us take 1 , 2 ∈ R d , and decompose them in terms of their tangential and normal components, that is, k = ( k ) T + ( k n F )n F , k = 1, 2. Then, we see that on F where the penultimate inequality is due to (3.90), as ( k ) T are tangent vectors. Since this holds for all 1 , 2 ∈ R d , we deduce that sup x∈F |∇ T n T F (x)| ≤ C, which, combined with (3.94) yields which is (3.93). Finally, from (3.96) we obtain the following which is (3.92). ▪

FINITE ELEMENT SCHEMES AND STABILITY ESTIMATES
We now provide DGFEMs for the approximation of solutions to (2.1) for k = 1, 2. For the case k = 1, we will not need to alter the bilinear form A 1 : V h × V h → R, given by (2.29). However, as mentioned in Section 2, obtaining an estimate in a H Δ -type norm (given by (2.31)) in the case k = 2 does not seem possible when considering curved finite elements, due to the form that the inverse inequality takes. This means that we must define a different bilinear form, which relies on a discrete analogue of the following identity Indeed, assuming that Ω is Lipschitz continuous, the above estimate follows from an application of integration by parts (twice), for i, j = 1, …, d, we see that for where the last equality is due to the fact that v| Ω = j v| Ω = 0 for j = 1, …, d. Summing (4.2) over all i, j = 1, …, d, we obtain (4.1). Furthermore (4.1) extends to u, v ∈ H 2 0 (Ω) by density, and, coupled with the Poincaré inequality, allows one to prove that the biharmonic problem (2.1) (for k = 2) is well posed in H 2 0 (Ω). Let us define the bilinear form C : V 2, h × V 2, h → R as follows: where n F is a fixed choice of unit normal to F,  F ≔ ∇ T ⋅ n F , and  ∶ R × R → R is defined by We now note the following consistency identity from [39] (c.f. [39], Lemma 4.1).

Lemma 4.1
Assume that Ω is piecewise C 2 and that ( h ) h>0 is a regular of order 1 meshes on Ω satisfying Assumption 3.21. Then, the bilinear form C : V 2, h × V 2, h → R satisfies the following consistency identity:

Numerical methods
We are now ready to provide the finite element schemes for the approximation of solutions to the Poisson, and biharmonic problem (2.1) for k = 1 and k = 2, respectively. Note that in the sequel we set V k, h ≔ V h, p , where we assume that p ≥ k.

Poisson problem
One seeks u 1, h ∈ V 1, h such that for all v h ∈ V 1, h , and we recall that B 1 , J 1 : V 1, h × V 1, h → R are defined as follows: where n F is a fixed choice of unit normal to F, and 1 F is a positive face dependent parameter.

Biharmonic problem
One seeks u 2, h ∈ V 2, h such that where B 2 , J 2 : V 2, h × V 2, h → R are defined as follows: F are positive face dependent terms to be provided, and C : V 2, h × V 2, h → R is defined by (4.3). Note that (4.6) is obtained by replacing u 2,h ), and including a tangential gradient penalty term in J 2 , in the definition of A 2 given by (2.29), and that the C(v h , u 2, h ) term results in A 2 defined by (4.6) being symmetric (while preserving the consistency of the scheme).
Remark 4.2 (Comparison to the method of [21]). In the method given by (4.6), we have taken V 2,h to be the space of discontinuous piecewise polynomials. However, if we instead use the space V h,c ≔ V 2,h ∩ H 1 0 (Ω), and assume that Ω is polygonal, then the scheme (4.6) coincides with the C 0 -interior penalty method proposed in [21] (without the second-order term, see (4.15) of [21]).

Stability estimates
We now provide stability estimates for the FEMs (4.5) and (4.6), which yield existence and uniqueness of a numerical solution. Let us recall the definitions of the norms that the stability estimates will hold in. We define the norms ‖⋅‖ h, k : V k, h → [0, ∞) as follows: for positive constants C *, k are to be determined.

Remark 4.3 (Trace estimates for jumps and averages). We note that for
where the multi-index satisfies | | ≤ m − 1, assuming that Ω is piecewise C m . Then, momentarily denoting {⋅} to be either the jump or average operator, it follows that for any In the sequel we shall utilize the above estimate several times, for various orders of , as it simplifies the exposition of the proofs.

Lemma 4.4
Assume that Ω ⊂ R d is piecewise C 2 , and that ( h ) h>0 is a regular of order 1 family of triangulations on Ω satisfying Assumption 3.21. Then, for any 1 ∈ (0, 1) there exists a constant 1 > 0 depending on 1 , C Tr , C I , and d, such that if Thus, there exists a unique u 1, h ∈ V 1, h that satisfies (4.5).

Lemma 4.5
Assume that Ω ⊂ R d is piecewise C 4 , and that ( h ) h>0 is a regular of order 3 family of triangulations on Ω satisfying Assumption 3.21. Then, for any 2 ∈ (0, 1) there exists a constant 2 > 0 depending on 2 , C Tr , C I , C P , and d, such that if

Proof.
Utilizing the trace and inverse estimates (3.41) and (3.47), and the Cauchy-Schwarz inequality with a parameter, we obtain for any 2 > 0, and Furthermore, a further application of the Cauchy-Schwarz inequality with a parameter yields .

ERROR ANALYSIS
We now use the consistency of the two schemes to prove optimal a priori error estimates for the numerical solutions u k, h , k = 1, 2, assuming that Ω is piecewise C m + 1 , m ∈ N, and that u k ∈ H k 0 (Ω)∩H 2k (Ω)∩ H s k (Ω;  h ), where s k = (s k K ) K∈ h and each s k K ≥ 2k, are the true solutions of (2.1) for k = 1, 2. That is, we prove the following estimate where t k K ≔ max{p + 1, m + 1, s k K }. Let us first recap on the approach we shall take, since this will in fact shorten the upcoming proofs. Let us take z k, h ∈ V k, h , k = 1, 2 to be arbitrary, denoting k, h ≔ u k − z k, h and k, h ≔ z k, h − u k, h , we see that Let us first estimate ‖ k, h ‖ h, k . Due to the interpolation estimate (3.65) (since the choice of z k, h ∈ V k, h is arbitrary) one can see that Since there are constants C k , k = 1, 2, satisfying 1 F ≤ C 1 and j F ≤ C 2 , j = 2, 3, 4, for all F ∈  i,b h , we have that Furthermore, Thus, from these two estimates, and an application of the interpolation estimate (3.65), we obtain .
(5.2) Applying (5.2) to (5.1) and taking square roots, we obtain and so, it remains to estimate ‖ k, h ‖ h, k , which relies upon the consistency of the schemes, and will be the objective of the next two proofs.
where t 1 K ≔ min{p + 1, m + 1, s 1 K }, u 1,h ∈ V 1,h is the unique solution of (4.5), and the constant C 1 depends upon the shape-regularity constant of the mesh, C I , C Tr , C  , d, s 1 K , and p, but not upon h K .
Proof. Let us take z h ∈ V 1, h , to be arbitrary, denoting h ≔ u 1 − z h and h ≔ z h − u 1, h , we see that Furthermore, from (5.3) for k = 1, we see that Then, applying the Cauchy-Schwarz inequality for vectors in R N , estimate (5.2), and the interpolation estimate (3.65), we obtain We again apply the Cauchy-Schwarz inequality for vectors in R N , estimate (5.2), the interpolation estimate (3.65), the trace estimate (3.41), and the inverse estimate (3.47) yielding

It then follows that
Applying the above estimate to (5.7), we obtain Dividing through by ‖ h ‖ h, 1 above, and applying the result in conjunction with (5.6) to (5.5) , (5.8) where t 2 K ≔ min{p + 1, m + 1, s 2 K }, u 2,h ∈ V 2,h is the unique solution of (4.6), and the constant C 2 depends upon the shape-regularity constant of the mesh, C P , C I , C Tr , C  , d, s 2 K , and p, but not upon h K .

NUMERICAL RESULTS
A discussion on the implementation of the numerical methods in Firedrake [46,47] is provided in [39].

Implementation
Software and code: The experiments in this section have been implemented in the most recent version of the Firedrake software [46,47] (as of July 3, 2018), which interfaces directly with PETSc [48,49] running through a Python interface [50,51]. There are two Firedrake scripts, Curved-Dirichlet-DGFEM.py (applicable to (4.5)), and DGFEM-curved-biharmonic.py (applicable to (4.6)) used to generate the experiments of this section is available in the Github repository: https:// github.com/ekawecki/Firedrake_Poisson_Biharmonic.
We took the penalty parameter 1 F = 10p 2 (obtained experimentally), where p is the polynomial degree of the space V comp h,p . For each polynomial degree p = 1, 2, 3, we successively refined the mesh quasiuniformly. We observe the predicted optimal convergence rate ‖u − u h ‖ h,1 = (h p ), as well as the optimal rate ‖u − u h ‖ L 2 (Ω) = (h p+1 ), with the true values and EOCs in brackets provided in Tables 1  and 2, respectively.

CONCLUSION
In the setting of curved finite elements, we have successfully reviewed several key estimates from theory of finite elements on polytopal domains, such as trace estimates, inverse estimates, discrete Poincaré-Friedrichs' inequalities, and optimal interpolation estimates in noninteger Hilbert-Sobolev norms, that are well known in the case of polytopal domains. Furthermore, we have proven curvature bounds for curved simplices, and utilized all of these estimates by providing stability, and a priori error analysis, of the IPDG method for the Poisson problem, originally introduced in [26], and for a variant of the h-version of the hp-DGFEM for the biharmonic problem introduced in [34]. In Section 6, we have provided numerical experiments for both the Poisson and biharmonic problem, where the domain is taken to be the unit disk. We implement a polynomial approximation of the domain, validating the a priori error estimates of Lemma 5.1. The estimates proven as part of this paper should serve useful for future applications to second-and fourth-order (as well as higher order) elliptic problems on curved domains, in particular, nondivergence form second-order elliptic equations.