{'fraction_non_alphanumeric': 0.04977264348039818, 'fraction_numerical': 0.015239031584121912, 'mean_word_length': 4.045886730053741, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interplay between diffusion anisotropy and mesh skewness in Hybrid High-Order schemes

We explore the effects of mesh skewness on the accuracy of standard Hybrid High-Order (HHO) schemes for anisotropic diffusion equations. After defining a notion of regular skewed mesh sequences, which allows, e.g., for elements that become more and more elongated during mesh refinement, we establish an error estimate in which we precisely track the dependency of the local multiplicative constants in terms of the diffusion tensor and mesh skewness. This dependency makes explicit an interplay between the local diffusion properties and the distortion of the elements. We then provide several numerical results to assess the practical convergence properties of HHO for highly anisotropic diffusion or highly distorted meshes. These tests indicate a more robust behaviour than the theoretical estimate indicates.


Introduction
The last few years have seen a increased interest in novel discretisation methods, for diffusion equations, that support polytopal meshes (made of general polygons/polyhedra) and allow for arbitrary approximation orders: Hybridisable Discontinuous Galerkin methods [4], Virtual Element Methods (VEM) [3], Weak Galerkin Methods [11], etc. The Hybrid High-Order (HHO) method [7,8] is one of these arbitrary-order polytopal methods, and shares with the aforementioned ones the hybrid structure of unknowns (contrary to Discontinuous Galerkin methods [6]), that is, unknowns located in the elements and on their faces. We refer to the introduction of [8] for a thorough review of the literature on polytopal methods. The HHO method can be seen as a high-order extension of the Hybrid Mimetic Mixed method [9] and, contrary to some other polytopal methods, it has a flux formulation that makes it a Finite Volume method. Additionally, the design of HHO schemes is dimension-independent and has an enhanced compliance with the physics due to the construction of local problem-dependent reconstruction operators.
HHO schemes have been applied to and analysed for a variety of models (see [8] and references therein), with error estimates that have an explicit dependency on the physical data. These estimates are however obtained for "regular" polytopal meshes, that is, meshes whose elements are "isotropic" (not elongated in any particular direction, and whose faces have a diameter comparable to their elements' diameters). In this work we analyse and numerically test the HHO scheme for highly anisotropic diffusion equations and families of distorted meshes, that no longer satisfy the usual regularity conditions. We consider the archetypal linear diffusion model where Ω is a polytopal domain of R d , K : Ω → R d×d sym is a symmetric bounded uniformly coercive diffusion tensor, and f ∈ L 2 (Ω). The solution to (1) is taken in the classical weak sense.
A review of historical or polytopal numerical methods on distorted meshes is out of this paper's scope. We however mention the recent works [2,12] about the VEM on anisotropic meshes, which present numerical results for a Poisson problem with internal layer, and derive approximation properties of the relevant interpolators. The novelty of our work, besides considering the HHO method instead of the VEM, is to establish complete error estimates (not just interpolator approximation properties) that take into account not only the distortion of the mesh, but also the high anisotropy of the diffusion tensor and the subtle interplay between these two features. The approach used here can be adapted to other methods, such as VEM, to yield error estimates that account for this interplay.
This paper is organised as follows. The concept of regular skewed mesh sequences, for which the error analysis will be carried out, is introduced in Section 2; these meshes can have very elongated elements, provided that some local linear map transforms them into isotropic elements. The oblique elliptic projector is at the core of HHO schemes; its approximation properties on skewed elements are presented in Section 3, and are used in Section 4 to perform the error analysis of HHO schemes on skewed meshes. This analysis is based on local transports of each skewed element T into an isotropic element T ; this transport identifies a new diffusion tensor on T , whose anisotropy properties dictate the contribution of T to the global error estimate. The error estimate stated in Theorem 10 therefore highlights how the diffusion anisotropy and the mesh skewness are combined in the multiplicative constants. This approach has the added advantage of leading to an error estimate that is as optimal as the standard error estimate for anisotropic diffusion models on regular (non-skewed) mesh sequences. In Section 5, we perform a series of tests to evaluate the practical impact of high diffusion anisotropy and mesh skewness on the accuracy of HHO schemes. Some of the conclusions drawn from these tests are predicted by the error estimate but, overall, the HHO scheme is found to be more robust with respect to the diffusion anisotropy and mesh skewness than what the theoretical analysis seems to indicate. A conclusion is provided in Section 6.
Notations. The Euclidean norm of a vector ξ ∈ R d is denoted by |ξ|. If L : (R d ) s → R is an s-linear map, we define the norm of L by For X an open subset of R n , n ∈ {d, d − 1}, (·, ·)X and · X denote respectively the L 2 (X)-or L 2 (X) n -inner product and norm. Letting D s v be the s-th order differential of v, the H s (X)-seminorm of a function v ∈ H s (X) is |v| H s (X) := Ns(D s v) X .

Regular skewed mesh sequences
Let us first briefly recall the definition of polytopal mesh, referring to [8, Section 1.1] for details. A polytopal mesh of Ω is a couple M h = (T h , F h ) where T h is a collection of disjoint polytopes T -the elements-such that Ω = ∪T ∈T h T , and F h is the set of mesh faces whose closures form a partition of ∪T ∈T h ∂T , and such that each face is contained in one or two elements boundaries. Mesh faces can be different from the geometrical faces of the polytopes, the latter being possibly cut in two mesh faces in case of non-conforming mesh [8, Fig. 1.2]. The diameter of a subset X of R d is denoted by hX . The index h in M h is the meshsize h = maxT ∈T h hT . For T ∈ T h , we let FT be the set of faces F ∈ F h such that ∂T = ∪F ∈F T F . The outer normal to T on F ∈ FT is nT F . A matching simplicial mesh of T ∈ T h is a polytopal mesh of T made of simplices and whose faces correspond to the geometrical simplicial faces.
We now define the concept of regular skewed mesh sequence, which allows for elements that become more and more stretched along the sequence, provided that each element can be linearly mapped onto an "isotropic" element, that satisfies the regularity conditions of a standard regular mesh sequence [8, Definition 1.9].
Definition 1 (Regular skewed mesh sequence). Let H ⊂ (0, +∞) be a countable set with 0 as only accumulation point. For each h ∈ H, let M h be a polytopal mesh and φ h = (φT )T ∈T h be a family of isomorphisms of R d . The sequence (M h , φ h ) h∈H is a regular skewed mesh sequence if there exists ∈ (0, 1) such that, for all h ∈ H and all T ∈ T h , the following properties hold: 3. For all τ ∈ T T , it holds h T ≤ hτ and hτ ≤ rτ , where rτ is the inradius of τ .
Remark 2 (Comparison with [12]). The notion of regular skewed mesh sequence is close to the notion of regular anisotropic mesh of [12], in particular in the usage of maps from skewed elements to isotropic elements. A noticeable difference, however, is the requirement in [12] that two neighbouring elements T, T must have similar isotropy (that is, the corresponding mappings φT , φ T must be close in a proper measure); this is due to the type of interpolators considered in [12], which are adapted to VEM and therefore require to compute averaged values around each vertex. Such a requirement of similar isotropy for neighbouring elements is absent from Definition 1, which is geared towards methods -such as HHO-whose interpolators are L 2 -projections on cell and face polynomials; as a consequence, this definition allows for example for meshes with layers of very thin rectangles neighbouring layers of squares.
In the rest of the paper, we consider a regular skewed mesh sequence (M h , φ h ) h∈H with parameter , and we write a b if a ≤ Cb with C > 0 depending only on Ω and and, when the inequality involves H s seminorms, also on the exponent s. We write a ≈ b if a b and b a. We also make the following assumption.
Assumption 3 (Piecewise constant diffusion tensor). For all h ∈ H, the diffusion tensor K is piecewise constant on T h . For any T ∈ T h we set KT = K |T .
Let T be an element of one of the meshes M h . If x ∈ T we set x = φT (x) ∈ T . The gradient (resp. differential) with respect to x is denoted by ∇ (resp. D). For . We also set JφT = |det φT |, and define Jφ T |F as the absolute value of the determinant of the restriction φ T |F : HF → H F , where HX denotes the hyperplane generated by X = F or F ; Jφ T |F can be computed using any pairs of orthonormal bases in F and F . Letting φ t T be the transpose of φT , the relevant diffusion tensor on T is: The maximal and minimal eigenvalues of K φ, T are denoted by K φ, T and K φ, T .
2. Transport of L 2 -inner products and norms. For all w, z in L 2 (T ) or L 2 (T ) d , For all F ∈ FT and w, z ∈ L 2 (F ), For all w, z ∈ H 1 (T ), Since T contains a ball of radius 2 h T (Point 3 in Definition 1) and T has diameter hT ≤ −1 h T , we see that φ −1 T maps a ball of radius 2 h T into a ball of radius −1 h T . Hence, where A is the matrix of φ T |F in (BF , BF ). In particular, JφT = |det A|λ = Jφ T |F λ and thus λ = JφT /Jφ T |F . Transposing the matrix above gives the matrix of φ t T in the orthonormal bases ( B, B). Since the last vector of B (resp. B) is n T F (resp. nT F ), reading the last column of this transposed matrix gives φ t T n T F = λnT F and proves (3). 2. Simple changes of variables (in T or F ) establish (4) and (5).
, an induction on s shows that, for all ξ 1 , . . . , ξ s ∈ R d , We infer that . . , s, and the proof of (6) is complete The relation (7) is obtained transposing (9) for s = 1. The second relation in (8) follows from the first with z = w. To prove this first relation, apply (4) to KT ∇w and ∇z instead of w and z, use the fact that KT is constant and invoke (7) to write

Oblique elliptic projector on skewed elements
Here, T is a generic element of M h . Fix a polynomial degree ≥ 0 and recall the definition in [8, Section 3.2.1] of the oblique elliptic projector π 1, The approximation properties of the oblique elliptic projector form an essential component of the analysis of HHO schemes for (1). To establish these approximation properties, let us first describe how the elliptic projector is transported through φT .
Lemma 5 (Transport of the elliptic projector). Letting π 1, K,φ, T be the oblique elliptic projector on T for the tensor K φ, T defined by (2), it holds Proof. Take w ∈ P (T ) and write, using the definition (10) of π 1, K,T , the transport relation (8) applied to (v, w) instead of (w, z), and the definition of π 1, K,φ, T together with w ∈ P ( T ), Using Proposition 6 (Approximation properties of the elliptic projector on skewed elements). For all s ∈ {1, . . . , + 1} and all v ∈ H s (T ), and, if s ≥ 2, for all F ∈ FT , Proof. Since T satisfies Points 2 and 3 in Definition 1, [8, Theorem 3.3] yields The volumetric (16) and trace (17) estimates are obtained transporting these estimates with (12). We start with the volumetric estimate. Using (12) and (8) we have Hence, applying (18) and using the estimate h T hT (see Point 1 in Definition 1), By the definition of the H s -seminorm, the relation (6) and the transport (4) give Plugged into (20), this concludes the proof of (16). We now turn to (17). The transport relations (12), (7) and (5) together with the definition (2) of K φ, T yield Estimate (17) follows plugging this relation into (19), using (21) and recalling (15) and that h T hT .
Remark 7 (Optimality of the approximation properties). This proof shows that (16) and (17) come from the corresponding inequalities (18) and (19) for isotropic elements, and from (21), itself derived from (6). The latter inequality is optimal in the sense that, for any φT , there are functions w for which it is an equality. Hence, the approximation properties (16) and (17) for skewed elements are as optimal as the corresponding approximation properties for isotropic elements.

Analysis of HHO schemes on skewed meshes
We briefly recall the construction of HHO schemes for (1) (referring to [8, Chapter 3.1] for a comprehensive presentation), and establish key properties for proving error estimates on skewed meshes. In the following, k ≥ 0 is a fixed polynomial degree.

Local space and potential reconstruction
For T ∈ T h , the local space of unknowns is For isotropic elements, this norm is usually defined using hF instead of dT F , see [8, Section 3.1.3.2]. The choice made in (22) ensures, for skewed elements, optimal estimates in terms of φT . The potential reconstruction p k+1 Lemma 8 (Transport of potential reconstruction). It holds where v T = ( vT , ( vF )F ∈F T ) ∈ U k T is the transported v T , and p k+1 K,φ, T is the potential reconstruction on T for the diffusion tensor K φ, T .
Proof. For all w ∈ P k+1 (T ), where we have used in this order the transport relation (8), the definition (23) of p k+1 K,T , the transport relations (8) and (5), and the fact that KT and nT F are constant. Invoking (7), (3) and (2), we have and (26) gives the conclusion following from the definition of p k+1 K,φ, T . Since w is arbitrary in P k+1 ( T ), this proves that p k+1 K,T v T and p k+1 K,φ, T v T have the same gradient. Using (4) and (24) we also see that they have same average on T , which concludes the proof of (25).

Local bilinear form
The difference operators δ k K,T : U k T → P k (T ) and, for F ∈ FT , δ k K,T F : where, for X = T or F , π 0,k X : L 2 (X) → P k (X) is the L 2 (X)-orthogonal projection. We note that, for any w ∈ L 2 (X), The local stabilisation bilinear form is given by: The local HHO bilinear form a K,T : U k T × U k T → R is then defined by: In the right-hand side above, the first term is responsible for the consistency of the bilinear form, while the addition of the second term ensures the stability and boundedness property stated in the following proposition. Other choices of s K,T are possible [8, Assumption 3.9], and, on isotropic meshes, the factor dT F in this stabilisation bilinear form can be replaced by hF .
Proposition 9 (Stability and boundedness of a K,T ). It holds Proof.

HHO scheme and error estimate
The global discrete space of unknowns is obtained patching local spaces and enforcing homogeneous Dirichlet boundary conditions: This global bilinear form defines the energy norm such that, The HHO scheme for (1) is written: find u h ∈ U k h,0 such that This scheme is well-posed, and is a Finite Volume scheme in the sense that it has a flux formulation [8,Lemma 3.17]. Our main result is the following theorem.
Theorem 10 (Discrete energy error estimate for HHO schemes on skewed meshes). Assume that the weak solution u ∈ H 1 0 (Ω) to (1) is such that, for some r ∈ {0, . . . , k}, u |T ∈ H r+2 (T ) for all T ∈ T h . Let u h ∈ U k h,0 be the solution to the HHO scheme (38). Then, it holds where α K,φ, T is the anisotropy ratio of K φ, T , defined by α K,φ, T := Remark 11 (Optimality of the error estimate). Following Remark 7, Estimate (39) is as optimal with respect to the mesh skewness as the corresponding estimate [8, Theorem 3.18], for isotropic meshes, is optimal with respect to the diffusion tensor.
Remark 12 (Interplay between mesh skewness and diffusion anisotropy). Assume for simplicity that d = 2 and that, for any T ∈ T h , there is an orthonormal basis in which where |u| H r+2 (T h ) is the usual broken H r+2 -seminorm of u. The first term in the righthand side of (41) encodes the interaction between the skewness of the mesh elements and the local anisotropy of the diffusion tensor.
Theorem 10. Applying the 3rd Strang lemma [5], we have where The following relation is established in the proof of [8, Lemma 3.15]: Let T nT F , using Cauchy-Schwarz inequalities and |K where we have used (17) (with = k + 1 and s = r + 2) and the norm equivalence (31) Since T is isotropic and I k T u = I k T u (owing to (28)), the consistency properties [8, Lemma 3.10] of s K,φ, T and the relations h T hT and (21) yield Plug this estimate and (44) into (43), use α K,φ, T ≥ 1 and recall (42) to conclude.

Numerical evaluation of the effects of diffusion anisotropy and mesh skewness
We provide here a series of numerical results, on the domain Ω = (0, 1) 2 (and with nonhomogeneous Dirichlet boundary conditions -see [8,Section 2.4] for the adaptation of the scheme (38) to this case), to assess the practical optimality of the error estimate (39) and its consequence (41) in cases of highly anisotropic diffusion tensor and/or skewed mesh families. The accuracy of the HHO scheme is measured through the following two relative errors: where · a,K,h is defined by (37), and · 1,h is the diffusion-independent discrete H 1 -norm obtained adding together the local seminorms (22) with K = I d , that is: The numerical tests have been performed using the code "HHO-Diffusion" in the C++ open source library HArDCore [1]. This library provides generic tools for implementing 2D and 3D numerical methods with unknowns made of polynomials on the edges/faces and cells of the mesh; it also naturally handles generic polygonal and polyhedral meshes. The variety of possible tests to assess the practical efficiency of the scheme (38) with anisotropic diffusion/skewed meshes is infinite, given the numerous possible parameters (polynomial degrees k, diffusion tensors, exact solutions, type of meshes, etc). We only report a few relevant results here, but all the meshes and data used in the tests below are available in HArDCore for the interested reader to run additional tests.

Test A: anisotropic diffusion tensor
This test focuses on the effect of an anisotropic and heterogeneous diffusion tensor. For λ ∈ {10 −6 , 1, 10 6 }, we consider the tensor and fix the exact solution u(x, y) = cos(πx) cos(πy); the source term and boundary conditions are computed from this solution. Since (∂xu) |y=0.5 = 0, we still have ∇·(K∇u) ∈ L 2 (Ω) despite the discontinuity of K along y = 0.5. We consider a family of locally refined meshes from the FVCA5 benchmark [10] (see Fig. 1), for which the setting of Remark 12 holds with λT = λ and aT = bT = 1; the estimate (41) therefore predicts a dependency of the energy error on max(λ, λ − 1 2 ). The results for k = 1, 3 are presented in Fig. 2; tests with other polynomial degrees present the same trend. The energy error E a,K,h appears to depend much less on the anisotropy ratio than predicted; the error E 1,h shows a more pronounced dependency on the tensor anisotropy, especially for low degrees where a factor of about 30 is seen on the finest mesh between λ = 1 and λ = 10 −6 , 10 6 .

Test B: skewed mesh
In this test, we study the impact of the mesh skewness. We take K = I d and the exact solution u(x, y) = cos(πx) cos(πy). The meshes are (mostly) hexagonal, and more and more skewed as their size decreases (see Fig. 3). The results in Fig. 4 show a clear loss of rate of convergence, compared to the expected rate for isotropic meshes.
To estimate more precisely the effect of mesh skewness, we introduce the flatness factor defined as fl h := maxT ∈T h flT with flT := h T ρ T , where ρT is the radius of the largest ball centred at the centre of mass of T and contained in T . The skewness of the considered meshes comes from the large flatness factors flT of some elements T . It is easy to convince oneself that this setting is compatible with Remark 12 with λ = 1, aT = flT and bT = 1. As a consequence, (41) predicts an upper bound      Table 1 an evaluation of the rates of growth of E h /h k+1 with respect to fl h . Estimate (46) tells us that, at least for the energy error, this rate should be at a maximum of 2. As can be seen in Table 1, the actual rates are much smaller than 2, and both errors are less sensitive to the mesh flatness than (46) predicts; the diffusion-independent norm E 1,h is the least sensitive of both. Table 1 also reports the condition numbers (CN) in 1-norm for the statically condensed system. For regular mesh sequences, CNs of HHO systems grow as h −2 . Here, the growth is in h −4 (but the CNs do not depend much on k). The additional power of 2 could come from a factor fl 2 (since, here, fl ∼ 1/h). Further analysis and tests are however necessary to reach a definitive conclusion, and it should also be noted that the meshes considered here contain a large portion of skewed elements; the condition numbers could be reduced for meshes with a smaller portion of distorted cells.

Test C
We assess here the interplay between mesh skewness and diffusion anisotropy, taking K(x, y) = diag(10 6 , 1), and u(x, y) = cos(πx) cos(πy) as before. We consider two families of meshes: regular hexagonal, and skewed hexagonal with flatness factor multiplied by two from one mesh member to the next; see Fig. 5.
The setting of Remark 12 is valid with (aT , bT ) = (1, fl h ) with fl h ≤ 10 3 for the considered meshes; (41) thus predicts a bound E a,K,h 10 6 fl −1 h h k+1 |u| H k+2 (Ω) . For the skewed meshes, we have fl h ∼ 1/h and we therefore expect a better rate of convergence than the usual h k+1 rate for isotropic meshes. Fig. 6 confirms this improvement, albeit in a non-uniform way.
The improvement is clearer if we superimpose the errors for the families of regular and skewed meshes, see Figs. 7a and 7b. For a given meshsize, selecting a mesh that is stretched in the direction of strong diffusion improves the convergence in both norms; this gain is valid for all degrees, but especially prominent for the lowest-order case k = 0 (for which, at the considered meshsizes, there is no apparent convergence on non-stretched meshes). In Figs. 7c and 7d the same errors are plotted against the number of globally coupled degrees of freedom, which for HHO schemes correspond to the edge unknowns (the element unknowns can be eliminated by static condensation [8,Appendix B]). In terms      of errors vs. number of degrees of freedom, the gain in using skewed meshes is less clear, except for k = 0; the reason is that meshes entirely made of stretched elements usually have, for a given meshsize, more edges than regular meshes.

Conclusion
We presented a theoretical and numerical study of the accuracy and robustness of the classical HHO method, when applied to anisotropic diffusion equations on distorted meshes. We defined a notion of mesh sequences that accepts in particular elements that become more and more elongated as the mesh is refined, and we established an error estimate that tracks the dependency of the constants with respect to the local diffusion anisotropy and elements skewness. We then presented the results of several numerical tests designed to explore the optimality of the error estimate. These results indicate that some behaviours highlighted by the theoretical estimate (such as the interplay between diffusion anisotropy and mesh skewness) are perceptible in practical numerical results, but they also show that this estimate appears to be pessimistic in its prediction of the behaviour of the error in case of strong anisotropy or skewness. Further work remains to be done to obtain more optimal estimates in terms of dependency with respect to the tensor anisotropy (this only has to be done for non-skewed meshes, as our approach would then provide an optimal estimate for skewed meshes). An aspect that is not covered by our definition of regular skewed mesh sequences is the case of small edges/faces in otherwise isotropic elements; another approach has to be adopted to derive error estimates in such situations. Finally, even though the standard HHO scheme displays some level of robustness on distorted meshes, it would be interesting to develop a variant that is specifically adapted to such meshes, and leads to better condition numbers than the standard method.