A new approach to handle curved meshes in the hybrid high-order method

The hybrid high-order method is a modern numerical framework for the approximation of elliptic PDEs. We present here an extension of the hybrid high-order method to meshes possessing curved edges/faces. Such an extension allows us to enforce boundary conditions exactly on curved domains, and capture curved geometries that appear internally in the domain e.g. discontinuities in a diffusion coefficient. The method makes use of non-polynomial functions on the curved faces and does not require any mappings between reference elements/faces. Such an approach does not require the faces to be polynomial, and has a strict upper bound on the number of degrees of freedom on a curved face for a given polynomial degree. Moreover, this approach of enriching the space of unknowns on the curved faces with non-polynomial functions should extend naturally to other polytopal methods. We show the method to be stable and consistent on curved meshes and derive optimal error estimates in $L^2$ and energy norms. We present numerical examples of the method on a domain with curved boundary, and for a diffusion problem such that the diffusion tensor is discontinuous along a curved arc.


Introduction
In recent years, there has been a trend in the computational literature towards arbitrary order polytopal methods for the approximation of partial differential equations.Such methods have a greater flexibility in the mesh requirements and can capture more intricate geometric and physical details in the domain.Being of arbitrary order, they also benefit from better convergence rates with respect to the global degrees of freedom.A short list of such methods includes discontinuous Galerkin and hybridizable discontinuous Galerkin methods [13,18,24], virtual element methods [1,6,9,15], weak Galerkin methods [29], and polytopal finite elements [36].However, it is well known that any approximation method on a polytopal mesh of a smooth domain (i.e. with a first order representation of the boundary) will yield at best an order two convergence rate [35,37].Thus, any high-order method on curved domains requires a high-order (or exact) representation of the boundary for optimal convergence.
Developed in [23,25], hybrid high-order (HHO) schemes are modern polytopal methods for the approximation of elliptic PDEs.A key aspect of HHO is its applicability to generic meshes with arbitrarily shaped polytopal elements.This article focuses on the extension of HHO methods to allow for curved meshes, with unknowns that capture the geometry exactly, yet still achieve optimal convergence.While the approach is presented within an HHO framework for a diffusion problem, the key ideas are more general and can be extended to related polytopal methods and to other models such as linear elasticity, or the Stokes and Navier-Stokes equations.
There has been much work on the development of discontinuous Galerkin (DG) methods on curved meshes [12,14,28].We also make note of the article [31] which analyses several approaches to high-order finite element methods on curved meshes.However, for the aforementioned methods the problem is much simpler than for hybrid high-order methods due to the lack of unknowns on the mesh faces.The addition of unknowns on the mesh faces is one of the key benefits hybrid methods have over DG and other non-hybrid methods due to the strong enforcement of boundary conditions and the reduction of the global degrees of freedom via static condensation [22,Appendix B.3.2].
The article [5] proposes a virtual element method (VEM) in two dimensions for meshes possessing curved edges.For each curved edge the authors consider the space of polynomials on a linear reference segment in R and map this space onto the curved edge via a sufficiently smooth parameterisation.A similar approach is taken in the articles [4,21].A typical approach for hybridizable discontinuous Galerkin (HDG) methods on curved domains is to map the boundary data onto a polytopal sub-domain [cf.19,20].We make note of the articles [32,33] which also use this approach to curved boundaries.
While there has been some work on the development of hybrid high-order methods on curved meshes [8,10,11], the approach we take in this paper is quite different.Indeed, the article [8] approaches the issue of defining unknowns on curved faces by considering a polynomial mapping from a planar reference face onto the curved face.While this naturally requires the mesh faces to be polynomial, it also reduces the approximation order [7].Indeed, if the mapping onto the face has effective mapping order m (see [8,Equation (5) & Remark 1]), then defining face unknowns of degree l in the reference frame will yield approximation properties of at best order l m [8,Equation (8)].To recover the optimal approximation order observed for straight meshes, the degree of the face polynomials in the reference frame is increased by a factor of m, yielding a very large global stencil for high-order mappings.Moreover, approximation properties in the curved faces are unknown, and the authors assume them to be true [cf.8, Equation (9)] in order to obtain optimal error estimates.We also make note of the conference proceedings [34] which follows the same approach using reference frame polynomials to define unknowns on curved faces for a HDG method.
An alternative approach, first considered in [10,11], is to increase the polynomial degree of the element unknowns and weakly enforce the boundary or interface conditions without defining any unknowns on curved faces.This procedure has also been implemented for a fourth order bi-harmonic problem in a curved domain [26].Such an approach ensures stability of the system and that optimal convergence rates are achieved.However, this method does not capture the geometry exactly and requires a finely tuned Nitsche parameter to achieve stability and consistency [cf.30].Moreover, without unknowns defined on curved faces, it is not clear how to design an enriched method such as that proposed in [38], whereas the method devised in this paper works seamlessly with enrichment.
In this paper we take inspiration from the article [38] and consider unknowns on the faces to include the Neumann traces of higher order polynomials.We note that this approach does not consider reference elements or faces but rather directly defines non-polynomial spaces on curved faces.Such an approach is therefore more closely analogous to an enriched or extended method than it is to any of the previously mentioned methods of defining unknowns on curved faces.Using this approach, we are not restricted to consider polynomial faces, but can rather take any C 1 manifold.Moreover, the number of degrees of freedom on curved faces is strictly bounded above and does not grow arbitrarily large for high-order mappings.We are able to prove consistency of the scheme, and, by including the space of constant functions on the faces the method is shown to be stable.In Section 3 we prove optimal error estimates in energy and L 2 -norm, and in Section 5 we present a method for the design of quadrature rules on curved elements.The paper is concluded with some numerical tests in two dimensions in Section 6.

Model and Assumptions on the Mesh
We take a domain Ω ⊂ R d , d ≥ 2, and consider the Dirichlet-diffusion problem: where a(u, v) := (K∇u, ∇v) Ω and (v) := (f, v) Ω for some source term f ∈ L 2 (Ω) and diffusion tensor K assumed to be a symmetric, piecewise constant matrix-valued function satisfying, for all for two fixed real numbers 0 < K ≤ K.Here and in the following, (•, •) X is the L 2 -inner product of scalar-or vector-valued functions on a set X for its natural measure.We shall also denote by Let H ⊂ (0, ∞) be a countable set of mesh sizes with a unique cluster point at 0. For each h ∈ H, we partition the domain Ω into a mesh M h = (T h , F h ), where T h denotes the mesh elements and F h the mesh faces.
We suppose that the mesh elements T h are a disjoint set of bounded simply connected domains in R d with piece-wise C 1 boundary ∂T .We further suppose that Ω = T ∈T h T .
We suppose that the mesh faces F h are a disjoint set of non-intersecting, finite, (d − 1)dimensional C 1 manifolds which partition the mesh skeleton: T ∈T h ∂T = F ∈F h F , and that for each F ∈ F h there either exists two distinct elements T 1 , T 2 ∈ T h such that F ⊂ ∂T 1 ∩ ∂T 2 and F is called an internal face, or there exists one element T ∈ T h such that F ⊂ ∂T ∩ ∂Ω and F is called a boundary face.Interior faces are collected in the set F i h and boundary faces in the set F b h .The parameter h is given by h := max T ∈T h h T where, for X = T ∈ T h or X = F ∈ F h , h X denotes the diameter of X.We shall also collect the faces attached to an element T ∈ T h in the set F T := {F ∈ F h : F ⊂ T }.The unit normal to F ∈ F T pointing outside T is denoted by n T F , and n T : ∂T → R d is the unit normal defined by (n T )| F = n T F for all F ∈ F T .We note that as each F is C 1 , the normal n T F is well defined.It is also worth noting that the normal vector n T F will not be constant on curved faces.
We consider the following regularity assumption on the mesh elements.
Assumption 1 (Regular mesh sequence).There exists a constant > 0 such that, for each h ∈ H, every T ∈ T h is connected by star-shaped sets with parameter (see [22,Definition 1.41]).
Remark 1 (Assumptions on the mesh).Assumption 1 is taken from [27, Assumption 1] however we have removed the assumption that the faces are connected by star shaped sets.We shall also note that there is no requirement that the mesh elements be polytopal or for the mesh faces to be planar.
We further require that the elements of the mesh align with the discontinuities of the diffusion tensor, i.e., for each T ∈ T h , K| T := K T is a constant matrix.In an analogous manner to (1.2) we define quantities 0 < K T ≤ K T to satisfy and we define the local diffusion anisotropy ratio α T := K T K T . From hereon we shall write f g if there exists some constant C which is independent of the quantities f and g, the mesh diameter h, and of the diffusion tensor K, such that f ≤ Cg.
Under Assumption 1 the following continuous trace inequality holds: for all v ∈ H 1 (T ), A proof of (1.4) is provided in [27].We note that no assumption on T being polytopal is required.We also note the following inverse Sobolev inequality, a proof of which is provided for highly generic and potentially curved elements in [12,Lemma 4.23]: where we denote by P (T ) the space of polynomials on T of degree ≤ , ∈ N. Combining (1.4) and (1.5) yields the following discrete trace inequality:

Discrete Model
A standard hybrid high-order method on polytopal meshes defines the local discrete space as This makes sense on polytopal meshes where F is a (d−1)-dimensional hyperplane as there is no ambiguity in what is meant by P k (F ).Indeed, on such meshes it holds that P k (F ) = P k (Ω)| F = P k (Ω) d • n F .On curved meshes, it is not so obvious what the discrete space should be.We find that the appropriate local discrete space is that of where we define and n F is an arbitrary unit normal to the face F .The choice of unit normal n F does not affect the definition of P k (F ).We note that, even for a curved face, there is no ambiguity in the term P 0 (F ) as it represents the space of functions which are constant on the face F .We emphasise that as the unit normal n F is not constant, the space P k (F ) will be non-polynomial on curved faces.
Remark 2. If F is planar (that is, a (d−1)-dimensional hyperplane) then it holds that P k (F ) = P k (F ) and thus the discrete space in (2.1) coincides with the usual HHO space.
Remark 3. It suffices to take the space of unknowns on the faces as for stability and consistency to hold.However, we define the space as P k (F ) = P 0 (F ) + P k (Ω) d • n F for simpler implementation and robustness of more general models.
We shall denote by P k (F T ) the space and for a given We denote by π 0,k T : L 1 (T ) → P k (T ) and π 0,k F : L 1 (F ) → P k (F ) the L 2 -orthogonal projectors onto the spaces P k (T ) and P k (F ) respectively.We denote by π 1,k+1 K,T : L 1 (T ) → P k+1 (T ) the oblique elliptic projector onto the space P k+1 (T ) satisfying The following weighted inner-products and norms are taken from [27].The weighted innerproduct (•, (2.5) T ∇v| H r−1 (T ) d . (2.6) Lemma 1 (Approximation properties of π 1,k+1 K,T ).For all s = 1, . . ., k + 1 and v ∈ H k+2 (T ), Proof.A proof is provided by [27,Lemma 9].While that particular proof assumes the elements are polytopal, the proof only relies on [22, Theorem 1.50] which is provided for generic elements connected by star-shaped sets.
The interpolant T is defined by where π 0,k Lemma 2. The following commutation property holds: (2.9) Proof.It follows from the definitions of p k+1 K,T and where in the last two equalities we have integrated by parts and introduced the oblique elliptic projector using equation (2.4a).Taking Remark 4. We note that the commutation property (2.9) is the key result required to prove consistency of the scheme and relies on the fact that The additional condition that P 0 (F ) ⊂ P k (F ) is required for coercivity to hold.
We endow the discrete space U k T with the seminorm The local bilinear form a K,T : where s K,T : U k T × U k T → R is a local stabilisation term such that the following assumptions hold.
Assumption 2 (Local stabilisation term).The stabilisation term s K,T is a symmetric, positive semi-definite bilinear form that satisfies: 1. Stability and boundedness.For all v T ∈ U k T , An example of a stabilisation bilinear form satisfying Assumption 2 is provided in Section 4.
Lemma 3 (Consistency of s K,T ).Suppose s K,T : (2.14) Therefore, applying the upper bound in (2.12) and the definition (2.10) of Thus, we infer from Lemma 7 below that The proof follows from the approximation properties (2.7) of π 1,k+1 K,T .

Global Space and HHO Scheme
The global space of unknowns is defined as To account for the homogeneous boundary conditions, the following subspace is also introduced, For any v h ∈ U k h we denote its restriction to an element The HHO scheme reads: find where h : U k h,0 → R is a linear form defined as We define the discrete energy norm Thus, if v h a,K,h = 0 then it must hold that v T = v F = const for every T ∈ T h , F ∈ T h .However, we infer from the homogeneous boundary conditions that those constants must all be zero.

Error estimates
Theorem 5 (Consistency error).The consistency error E h (w; If such a w additionally satisfies w| T ∈ H k+2 (T ) for all T ∈ T h , the consistency error satisfies The global operators p k+1 K,h : such that their actions restricted to an element T ∈ T h are that of p k+1 K,T and π 1,k+1 K,T .The global interpolator Theorem 6 (Energy and L 2 error estimates).Let u ∈ H 1 0 (Ω) be the exact solution to equation (1.1) and suppose the additional regularity u ∈ H k+2 (T h ).Let u h be the exact solution to the discrete problem (2.18).Then the following error estimates hold: • L 2 estimate.Suppose additionally that the domain Ω is convex and K = I is the identity matrix, then optimal convergence in L 2 -norm holds: Proof of Theorems 5 and 6.The estimates (3.1) and (3.2) are provided in [27] and rely only on the design conditions stated in Assumption 2, the commutation property (2.9), the approximation properties of the elliptic projector (2.7), the consistency of s K,T (2.14), Lemma 4, as well as standard trace and inverse estimates provided in Section 1.1.
To prove (3.3) we require a slightly different approach to that of [22,Theorem 2.32].In particular, as π 0,k F is not a polynomial projector [22, Equation (2.78)] does not hold in our case.However, the remainder of the proof is the same so we only have to show that sup where z g is the solution to the dual problem a(v, z g ) = (g, v) Ω ∀v ∈ H 1 0 (Ω).As we have assumed Ω to be convex, the following elliptic regularity holds: Moreover, as K = I, the following equality established in the proof of [22,Lemma 2.18] holds true: The sum over the boundary term in (3.6) can be written as follows, As ∇π 1,k+1 K,T u • n T F ∈ P k (F ) we may drop the projector π 0,k F to write As ∇u ∈ H(div; Ω), the fluxes of u are continuous across every internal face F ∈ F i h .Therefore, as π 0,k F z g = 0 for all F ∈ F b h (due to z g = 0 on ∂Ω), it holds that Substituting back into (3.6)yields It follows from a Cauchy-Schwarz inequality and the consistency (2.14) that . It also follows from a Cauchy-Schwarz inequality, the continuous trace inequality (1.4) and the approximation properties (2.7) that and the proof follows from the elliptic regularity (3.5) and the bound g Ω ≤ 1.By a continuous trace inequality and a Poincaré-Wirtinger inequality The result holds due to the H 1 -approximation properties of the L 2 -projector [22, Lemma 1.43] which remain valid in curved domains.

Analysis of the stabilisation
We consider here the stabilisation bilinear form defined by however, the arguments we use to show robustness on curved meshes extend seamlessly to more general choices of stability such as those considered in [27, Section 4].It is clear that s K,T satisfies (2.13) so it remains to prove that (2.12) holds.
Proof.We first note the bound ∂T which follows from the ellipticity (1.3) of K T .Consider, by a triangle inequality First, we wish to bound the term v − π 0,k F T v ∂T .As π 0,k F T is the L 2 -orthogonal projector on P k (F T ), it minimises its respective norm.Therefore, we may replace π 0,k F T v with any element of P k (F T ).In particular, as It follows from the continuous trace inequality (1.4) and a Poincaré-Wirtinger inequality that Similarly, we apply the continuous trace inequality and a Poincaré-Wirtinger inequality on the term h where we have applied a triangle inequality to reach the conclusion.It follows from [22, Equation (1.77)] (which invokes [22, Equation (1.74)] which does not rely on the elements being polytopal) T .Thus, we can conclude that The proof follows by applying the ellipticity (1.3) of K T to yield Remark 6.We note that the inclusion P 0 (F ) ⊂ P k (F ) is crucial for the bound to hold, and without this inclusion, coercivity cannot hold.

Lemma 8 (Coercivity). It holds for all
Proof.It follows from the definition (2.10) of where we have added and subtracted π 0,k T p k+1 K,T v T to the volumetric term, and π 0,k T p k+1 K,T v T and π 0,k F T p k+1 K,T v T to the boundary term, and invoked triangle inequalities to reach the conclusion.Similar to the proof of Lemma 7, we apply the continuous trace inequality (1.4), a Poincaré-Wirtinger inequality (due to the zero mean value of π 0,k We can conclude from Lemma 7 that which combined with the definition of a K,T yields the result.

Lemma 9 (Boundedness). It holds for all
Proof.Consider by a triangle inequality and Lemma 7 Thus, we need to prove that It follows from the definition (2.3) of p k+1 K,T and an integration by parts that where we have applied Cauchy-Schwarz inequalities on both inner-products and the discrete trace inequality (1.6).The proof follows by simplifying (4.5) by |p k+1 K,T v T | K,H 1 (T ) and squaring.

Integration on curved domains
The design of integration methods on curved domains is an active area of research.In the recent article [2] a quadrature rule for curved domains is developed by considering a decomposition into triangular or rectangular pyramids T and a mapping T : [0, 1] d → T for each decomposition.
With knowledge of the Jacobian of such a mapping, integration can be performed on the preimage of each T .The article [17] develops an extension of the homogeneous integration rule developed in [16] by considering a curved triangulation of the domain and constructing a scaled boundary parameterisation on each curved triangle.Here, we also consider an extension of the homogeneous integration rule, but the approach we take is quite different.We avoid the need to split the curved domain into sub-regions and directly map the integral onto the boundary by constructing a Poincaré-type operator which inverts the divergence operator.Indeed, this operator was briefly mentioned in the appendix of [17], however, we develop the ideas here without a sub-triangulation, and independent of dimension.We begin with the formula developed in [16] to rewrite the integral onto the boundary of the element.This rule works by identifying a vector field such that ∇ • F = v for homogeneous functions v of degree q.Therefore, We would like to extend this rule to non-homogeneous functions.We begin by searching for a vector field of the form F = g r, such that ∇ • F = v where r denotes the unit vector in the radial direction.We find that the unknown function g must satisfy 1 where we denote by r = |x|.A solution is given by Thus, we have found an inverse divergence Therefore, an integral over the element T can be rewritten to its boundary as follows: We note that if v is a homogeneous function of degree q (that is, v(tx) = t q v(x)), then the inverse divergence formulae (5.1) and ( 5.3) coincide and thus so do the rules (5.2) and (5.4).
In this sense, the method can be considered an extension of the homogeneous integration rule developed in [16].
If we instead consider a vector field of the form F = g r0 where r0 is the unit radial direction from a shifted origin x 0 , we arrive at the more general formula ( Therefore, we may write This is very useful if the element contains one or more planar faces.For a vertex with coordinates ν, we can set x 0 = ν and it holds that (x − ν) • n T = 0 on any planar faces connected to the vertex ν.We note that if T is not star-shaped with respect to ν, then the integral 1 0 t d−1 v(tx + (1 − t)x 0 ) dt will pass through points outside of T .Thus, one would require a sufficiently smooth extension of v outside of T .However, for polynomials or functions analytic over Ω (such as an analytic source term), such an extension is trivial.

A quadrature rule for curved edges in two dimensions
For a given edge E, consider a parameterisation γ E : [t 0 , t 1 ] → E, t 0 < t 1 .Therefore, integration on curved edges is trivial: The above integral can easily be approximated with a one-dimensional Gaussian quadrature rule.In particular, let w i , x i , i = 1, . . ., N be the weights and abscissae associated with a quadrature rule on [0, 1].Then we can generate weights w E i and abscissae x E i on the edge E as follows: In practise, we generally store an arc length parameterisation for each edge and thus the term |γ E | is not required.

A quadrature rule for elements in two dimensions
In two dimensions the faces are edges and thus the boundary integral in (5.6) can be evaluated on each edge F ∈ F T using the rule described in (5.7).We let w F i and x F i , i = 1, . . ., N be the quadrature weights and abscissae associated with an edge F ∈ F T and w j , x j , j = 1, . . ., M be the weights and abscissae associated with a quadrature rule on [0, 1].We set ν to be the coordinate of a vertex of T connected to the highest number of straight edges in T .We then consider the quadrature rule That is, we store weights and abscissae x j x F i + (1 − x j )ν, for each i = 1, . . ., N , j = 1, . . ., M and on each edge F ∈ F T that is not a straight edge connected to the vertex ν.
If T is polygonal, then there always exists two straight edges connected to a vertex ν.Thus, the rule described by (5.8) consists of (|F T | − 2)N M quadrature points.If we consider a Gauss-Legendre rule on each edge which is exact for polynomials of degree k, then we require to take N = k+1 2 .However, for the inverse divergence formula (5.5) to reproduce polynomials of degree k exactly, we require to take M = k+2 2 due to the presence of the multiplier t. quadrature points.We note this is a slightly larger number of quadrature points than the usual (|F T | − 2) k+1 2 2 required by splitting the polygon T into (|F T | − 2) sub-triangles (an optimal sub-triangulation) and considering a Gauss-Legendre rule on each sub-triangle.However, (5.8) avoids the complex process of generating such a sub-triangulation.To avoid these additional quadrature points, one would need to consider a Gauss-Legendre rule on each edge, but a weighted Gaussian rule with the weight function w(t) = t for the integral (5.5).This is not explored further here.

A quadrature rule for elements in three dimensions
In three dimensions, a volumetric integral can be mapped onto the faces as follows, (5.9) Thus, given a quadrature rule for each face F ∈ F T , a quadrature rule for the element T can be developed analogously to the two-dimensional case.On each face x 0 ) dt.Take the planar region F ⊂ R 2 with (potentially curved) edges Ê F and a parameterisation γ F : F → F .It holds that where J( x) = det J t ( x)J ( x) and J is the Jacobian matrix of the map γ F .It then follows from (5.6) that (5.10) where n F Ê denotes the unit normal directed out of F and towards Ê.Therefore, given the parameterisation γ F and a parameterisation of each mapped edge Ê ∈ Ê F , the integral (5.10) can be evaluated analogously to the 2D case (5.8).

A note on planar faces
If the face F is planar, one can follow a procedure similar to that in [3] to rewrite the integrals on each face onto the edges E ∈ E F .We take γ F ( x) = x F + E x where x F is a point in the face F and E is an orthonormal matrix.Then it holds that J( x) ≡ 1 and Thus, we can map the integral (5.10) back to the edges of the face F as follows, However, as E is orthonormal it preserves distance and therefore it holds that (γ −1 where n F E denotes the unit normal directed out of F and towards E.Moreover, the mapping γ F is onto, so we can choose x0 such that γ F ( x0 ) = x F,0 for an arbitrary point x F,0 ∈ F .Therefore Again, we may choose x F,0 to be the vertex of the face F connected to the largest number of straight edges.The integral (5.11) is then evaluated in an identical manner as two-dimensional elements.

Implementation
The HHO method for curved edges is implemented using the open source C++ library PolyMesh [39].We generate curved meshes by first considering uniform Cartesian meshes and 'cutting' along a curve.The integrals are computed using the quadrature rule described by (5.8) where we take the one-dimensional integration rules to be Gauss-Legendre rules of degree 30.
A basis is formed for the space P k (F ) by first generating a spanning set by considering a canonical basis of P k (Ω) d and taking P 0 (F ) + P k (Ω) d • n F .The linearly dependent basis functions are removed algebraically using the FullPivLU class found in the Eigen library, with documentation available at https://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html.This requires a threshold to be set which determines the point at which pivots are considered to be numerically zero.We set this value to 10 −15 .We note that for sufficiently small h and large k this can result in certain linearly independent functions being removed from P k (F ).However, as these functions are 'close' to being linearly dependent, the method seems unaffected by their removal.The bases of both P k (F ) and P k (T ) are orthonormalised via a Gram-Schmidt process.
The relative error of the scheme is measured through the following three quantities: , where the norm • L 2 (T h ) is defined as the square-root of the sum of squares of • L 2 (T ) .We note that if the mesh conforms to the domain Ω then v L 2 (T h ) = v Ω for all v ∈ L 2 (Ω).
We consider here two sequence of meshes of the domain Ω.The curved meshes use an exact representation of the boundary, whereas the straight meshes take a piece-wise linear approximation of the boundary.The parameters of the mesh sequences are displayed in Table 1.Both sequences of meshes have the same parameters.Example curved meshes are plotted in Figure 1 and straight meshes are plotted in Figure 2  In Figure 3 we test both a curved HHO scheme and a classical HHO scheme (on straight meshes) with polynomial degrees given by k = 1 and k = 3.In both cases the curved HHO scheme on the fitted mesh observes significantly better convergence rates than the classical scheme on the straight mesh.While the scheme appears to converge optimally on curved meshes, it converges at most order 2 on straight meshes.
We consider Ω = {(x, y) : x 2 + y 2 < 1} to be the unit disc and K a piece-wise constant diffusion tensor given by .
We take R = 0.8, β 1 = 10 −6 and β 2 = 1 which corresponds to anisotropic diffusion in the region r < R, and a Poisson problem in r > R. We take the source term to be f ≡ 1.Again, we consider two sequences of meshes of the domain Ω.We take both sequences to fit the domain Ω exactly, however, the curved mesh we take to fit the discontinuity in K exactly and the straight mesh takes a piece-wise linear approximation of K.The mesh data is presented in Table 2.We note that both sequences of meshes have the same parameters.An example curved mesh and an example straight mesh is plotted in Figure 5.As we do not know the exact solution to this problem, we run the scheme on the finest curved mesh with k = 7.We denote by the discrete solution to this problem p k+1 K,h u h = u * h , which will play the role of the 'exact' solution.We measure the quantities Ω u * h ≈ 0.46006947 ; |u * h | H 1 (T h ) ≈ 0.80699766.

Mesh
We would then like to test the performance of the scheme on coarser meshes (both curved and straight) with smaller k by investigating the behaviour of We are less interested in the rate of convergence of these measures, but rather want to observe steady convergence, and investigate the difference between the two schemes.In Figure 6 we plot the quantities E 1 and E 2 against increasing polynomial degree k where we fix the mesh to be Mesh 1.It is clear that the scheme on the straight mesh, where we consider an approximate diffusion tensor, stops converging for k > 1 whereas the scheme on the curved mesh converges smoothly.In Figure 7 we test against decreasing mesh size h for polynomial degrees k.While for k = 1 the order of E 1 and E 2 are similar for both schemes (and at times smaller on the straight mesh), the convergence is much smoother on the curved mesh.For k = 3, the values are significantly smaller for the curved mesh, and the behaviour for the straight mesh does not differ much from the k = 1 case.This coincides with the previous observations that increasing k past 1 has little effect on the scheme when considering a piece-wise linear approximation of the discontinuity in the diffusion.
Straight mesh Curved mesh   Finally, in Figure 8 we show contour plots of the potential reconstructions of the discrete solutions on Mesh 1 with k = 7.We observe that the plot on the straight mesh seems to be distorted along the eigen vectors of K (that is, (1, 1) t and (1, −1) t ) when compared to the plot on the curved mesh.We also plot the absolute value of the difference between the two schemes and observe that this value seems to be of greatest magnitude around the discontinuity in the diffusion tensor.

)Remark 5 .
where the seminorm | • | H s (T h ) is defined as the square-root of the sum of squares of | • | H s (T ) for any s ∈ N. The L 2 -error estimate is stated with identity diffusion, corresponding to a Poisson problem.However, the result follows trivially (with a hidden constant depending additionally on the anisotropy of K) for any constant diffusion tensor K [cf.22, Remark 3.21].

Figure 1 :Figure 2 :
Figure 1: Example curved meshes used for the curved boundary test

Figure 5 :
Figure 5: Example meshes used for the heterogeneous diffusion test

Table 1 :
. Parameters of the mesh sequences used for the curved boundary test

Table 2 :
Parameters of the mesh sequences used for the heterogeneous diffusion test