A Hybrid High-Order Method for Highly Oscillatory Elliptic Problems

1 Introduction

Over the last few years, many advances have been accomplished in the design of arbitrary-order polytopal discretization methods. Such methods are capable of handling meshes with polytopal cells, and possibly including hanging nodes. The use of polytopal meshes can be motivated by the increased flexibility, when meshing complex geometries, or when using agglomeration techniques for mesh coarsening (see, e.g., [7]). Classical examples of polytopal methods are the (polytopal) Finite Element Method (FEM) [46, 44], which typically uses non-polynomial basis functions to enforce continuity, and non-conforming methods such as the Discontinuous Galerkin (DG) [5, 16, 10] and the Hybridizable Discontinuous Galerkin (HDG) [15] methods. We also mention the Weak Galerkin (WG) [47] method (see [13] for its links to HDG).

More recently, new paradigms have emerged. One salient example is the Virtual Element Method (VEM) [9], which is formulated in terms of virtual (i.e., non-computed) conforming functions. The key idea is that the virtual space contains those polynomial functions leading to optimal approximation properties, whereas the remaining functions need not be computed (only their degrees of freedom need to be) provided some suitable local stabilization is introduced. The degrees of freedom in the VEM are attached to the mesh vertices, and, as the order of the approximation is increased, also to the mesh edges, faces, and cells. Another recent polytopal method is the Hybrid High-Order (HHO) method, which has been introduced for locking-free linear elasticity in [17], and for diffusion in [19]. The HHO method has been formulated originally as a non-conforming method, using polynomial unknowns attached to the mesh faces and cells. The HHO method has been bridged in [14] both to HDG (by identifying a suitable numerical flux trace), and to the non-conforming VEM considered in [6] (by identifying an isomorphism between the HHO degrees of freedom and a local virtual finite-dimensional space, which again contains those polynomial functions leading to optimal approximation properties). The focus here is on HHO methods. HHO methods offer several assets, including a dimension-independent construction, local conservativity, and attractive computational costs, especially in 3D. Indeed, the HHO stencil is more compact than for methods involving degrees of freedom attached to the mesh vertices, and static condensation allows one to eliminate cell degrees of freedom, leading to a global problem expressed in terms of face degrees of freedom only, whose number grows quadratically with the polynomial order, whereas the growth of globally coupled degrees of freedom is typically cubic for DG methods.

In this work, we are interested in elliptic problems featuring heterogeneous/anisotropic coefficients that are highly oscillatory. The case of slowly varying coefficients has already been treated in [18, 20], where error estimates tracking the dependency of the approximation with respect to the local heterogeneity/anisotropy ratios have been derived. Let Ω be an open, bounded, connected polytopal subset of ℝd, d∈{2,3}, and ε>0, supposedly much smaller than the diameter of the domain Ω, encode the highly oscillatory nature of the coefficients. We consider the model problem

where f∈L2⁢(Ω) is non-oscillatory, and 𝔸ε is an oscillatory, uniformly elliptic and bounded matrix-valued field on Ω. It is well known that the Hk+2-norm of the solution uε to problem (1.1) scales as ε-(k+1), meaning that monoscale methods (including the monoscale HHO method of order k≥0 of [18, 20]) provide an energy-norm decay of the error of order (hε)k+1. To be accurate, such methods must hence rely on a mesh resolving the fine scale, i.e., with size h≪ε. Since ε is supposedly much smaller than the diameter of Ω, an accurate approximation necessarily implies an overwhelming number of degrees of freedom. In a multi-query context, where the solution is needed for a large number of right-hand sides (e.g., an optimization loop, with f as a control and (1.1) as a distributed constraint), a monoscale solve is hence unaffordable. In that context, multiscale methods may be preferred. Multiscale methods aim at resolving the fine scale in an offline step, reducing the online step to the solution of a system of small size, based on oscillatory basis functions computed in the offline step, on a coarse mesh with size H≫ε. In a single-query context, multiscale methods are also interesting since they allow one to organize computations in a more efficient way.

Multiscale approximation methods on classical element shapes (such as simplices or quadrangles/hexahedra) have been analyzed extensively in the literature. Examples include, e.g., the multiscale Finite Element Method (msFEM) [34, 35, 23] (with energy-error bound of the form (ε12+H+(εH)12) in the periodic case), its variant using oversampling [34, 24] (with improved error bound of the form (ε12+H+εH) in the periodic case), or the Petrov–Galerkin variant of the msFEM using oversampling [36]. Let us also mention [3] (see also [33]), which is an extension to arbitrary orders of approximation of the classical msFEM (with error bound of the form (ε12+Hk+(εH)12) in the periodic case using H1-conforming finite elements of degree k≥1). These methods all rely on the assumption that a conforming finite element basis is available for the (coarse) mesh under consideration. Recent research directions essentially focus on the approximation of problems that do not assume scale separation, and on reducing and possibly eliminating the cell resonance error. One can cite, e.g., the Generalized msFEM (GmsFEM) [22], or the Local Orthogonal Decomposition (LOD) approach [32, 41]. We also mention that other paradigms exist to approximate oscillatory problems, like the Heterogeneous Multiscale Method (HMM) [21, 1].

On general polytopal meshes, the literature on multiscale methods is more scarce. For constructions in the spirit of the msFEM, one can cite the msFEM à la Crouzeix–Raviart of [39, 40], the so-called Multiscale Hybrid-Mixed (MHM) [4, 43] approach, and the (polynomial-based) method of [26] in the HDG context. Each one of these methods has its proper design, but they all share the same construction principles: they are based, more or less directly, on oscillatory basis functions that solve local Neumann problems with polynomial boundary data, and result in global systems (posed on the coarse mesh) that can be expressed in terms of face unknowns only. In the following, we will thus refer to those methods as skeletal-msFEM. The MHM approach actually presents a small difference with respect to the two other approaches since it is based on a hybridized primal formulation, which leads to consider flux-type unknowns at interfaces instead of potential-type unknowns; as a consequence, and in order to impose the compatibility constraint, one needs to solve a saddle-point global problem, whereas for the two other approaches, one ends up with a coercive problem. For the msFEM à la Crouzeix–Raviart, an error bound of the form (ε12+H+(εH)12) is proved in [39] in the periodic case. However, the analysis is led under the assumption that there exists a finite number of reference elements in the mesh sequence. For the MHM approach, which is designed in the same spirit, the same type of upper bound for the error is expected. Yet, in [43], the authors claim that their method is able to get rid of the resonance error (without oversampling); we clarify this issue in Remark A.3 below. For the HDG-like method, the analysis that is provided in [26] is sharp only in the regime H≪ε. As a consequence, there is, to date, no complete polytopal analysis available in the literature for skeletal-msFEM. Moreover, we observe that in the three methods, the discretization of the (non-oscillatory) right-hand side is realized in a somewhat suboptimal way, which can become a limiting issue in a multi-query context. In the msFEM à la Crouzeix–Raviart, the discretization is realized through a projection of the loading term onto the space spanned by the oscillatory basis functions. In the MHM and HDG-like approaches, the whole (local) space H1 is considered. In all cases, the approximation of the right-hand side does not take advantage of the fact that the latter is non-oscillatory. Let us mention, as another construction in the spirit of the msFEM, the work [38], which exploits in the DG context the ideas introduced in [3]. The drawback, which is inherent to DG methods, is the large size of the online systems. For constructions in the spirit of the GmsFEM, let us mention in the HDG context the contributions [25, 11] (that are based on [26]), and the work [42] in the WG context.

In this work, we devise a multiscale HHO (msHHO) method, which can be seen as a generalization (in particular to arbitrary orders of approximation) of the msFEM à la Crouzeix–Raviart of [39, 40]. Our contribution is twofold. First, we provide an analysis (in the periodic setting) of the method that is valid on general polytopal mesh sequences (in particular, we do not postulate the existence of reference elements); in that respect, this work presents the first complete polytopal analysis of a skeletal-msFEM. Note that considering general element shapes in the periodic setting is clearly not a good strategy (cf., e.g., [31]); however, this setting is not our final target. Second, taking advantage of the flexibility offered by the HHO framework, we improve on the existing methods. We introduce (polynomial) cell unknowns, that we use for the integration of the right-hand side (cf. Remarks 5.8 and 5.16 below). The non-oscillatory loading is hence discretized through a coarse-scale polynomial projection, while the size of the online system remains unchanged since the cell unknowns are locally eliminated in the offline step. Two versions of the msHHO method are proposed herein, both employing polynomials of arbitrary order k≥0 for the face unknowns. For the mixed-order msHHO method, the cell unknowns are polynomials of order (k-1) (if k≥1), whereas they are polynomials of order k≥0 for the equal-order msHHO method. The mixed-order msHHO method does not require stabilization, whereas a simple stabilization (which avoids computing additional oscillatory basis functions) is introduced in the equal-order case. We prove for both methods an energy-error estimate of the form (ε12+Hk+1+(εH)12)=:gk(H) in the periodic case. The analysis of the msHHO method differs from that of the monoscale HHO method since the local fine-scale space does not contain polynomial functions up to order (k+1); in this respect, our key approximation result is Lemma 4.5 below. With respect to [39], we also simplify the analysis and weaken the regularity assumptions (cf. Remark 4.6 below). Our analysis finally sheds new light on the relationship between the non-computed functions of the local virtual space and the associated local stabilization. To motivate the design and use of a high-order method, we note, as it was already pointed out in [3], that the upper bound gk⁢(H) is minimal for Hk=(ε12/2⁢(k+1))2/(2⁢k+3), hence as k≥0 increases, Hk increases whereas gk⁢(Hk) decreases. The msHHO method we devise is meant to be a first step in the design of an accurate and computationally effective multiscale approach on general meshes. The next step will be to address the resonance phenomenon and the more realistic setting of no scale separation.

The article is organized as follows. In Sections 2 and 3 we introduce, respectively, the continuous and discrete settings. In Section 4, we introduce the fine-scale approximation space, exhibiting its (oscillatory) basis functions and studying, locally, its approximation properties. In Section 5, we introduce the two versions of the msHHO method, analyze their stability, and derive energy-error estimates. In Section 6, we present some numerical illustrations in the periodic and locally periodic settings. Finally, in Appendix A we collect some useful estimates on the first-order two-scale expansion.

2 Continuous Setting

From now on, and in order to lead the analysis, we assume that the diffusion matrix 𝔸ε satisfies 𝔸ε(⋅)=𝔸(⋅/ε) in Ω, where 𝔸 is a symmetric and ℤd-periodic matrix field on ℝd. Letting Q:=(0,1)d, we define, for 1≤p≤+∞ and m∈ℕ⋆, the following periodic spaces:

with the classical conventions that Wperm,2⁢(Q) is denoted Hperm⁢(Q) and that the subscript “loc” can be omitted for p=+∞. Letting 𝒮d⁢(ℝ) denote the set of real-valued d×d symmetric matrices, we also define, for real numbers 0<a≤b,

We assume that there exist real numbers 0<α≤β such that

Assumption (2.1) ensures that 𝔸ε∈L∞⁢(Ω;ℝd×d) is such that 𝔸ε⁢(⋅)∈𝒮αβ a.e. in Ω for any ε>0, and hence guarantees the existence and uniqueness of the solution to (1.1) in H01⁢(Ω) for any ε>0. More importantly, assumption (2.1) ensures that the (whole) family (𝔸ε)ε>0 G-converges [2, Section 1.3.2] to some constant symmetric matrix 𝔸0∈𝒮αβ. Henceforth, we denote ρ:=βα≥1 the (global) heterogeneity/anisotropy ratio of both (𝔸ε)ε>0 and 𝔸0. Letting (𝒆1,…,𝒆d) denote the canonical basis of ℝd, the expression of 𝔸0 is known to read, for integers 1≤i,j≤d,

where, for any integer 1≤l≤d, the so-called corrector μl∈Hper1⁢(Q) is the solution with zero mean-value on Q to the problem

For further use, we also define the linear operator ℛε:Lperp⁢(Q)→Lp⁢(Ω), 1≤p≤+∞, such that, for any function χ∈Lperp⁢(Q), ℛε⁢(χ)∈Lp⁢(Ω) satisfies ℛε(χ)(⋅)=χ(⋅/ε) in Ω. In particular, for any integers 1≤i,j≤d, we have [𝔸ε]i⁢j=ℛε⁢(𝔸i⁢j). A useful property of ℛε is the relation ∂l⁡(ℛε⁢(χ))=1ε⁢ℛε⁢(∂l⁡χ), valid for any function χ∈Wper1,p⁢(Q) and any integer 1≤l≤d.

The homogenized problem reads

We introduce the so-called first-order two-scale expansion

Note that (uε-ℒε1⁢(u0)) does not a priori vanish on the boundary of Ω.

3 Discrete Setting

We denote by ℋ⊂ℝ+⋆ a countable set of meshsizes having 0 as its unique accumulation point, and we consider mesh sequences of the form (𝒯H)H∈ℋ. For any meshsize H∈ℋ, a mesh𝒯H is a finite collection of non-empty disjoint open polytopes (polygons/polyhedra) T, called elements or cells, such that Ω¯=⋃T∈𝒯HT¯ and H=maxT∈𝒯H⁡HT, HT standing for the diameter of the cell T. The mesh cells being polytopal, their boundary is composed of a finite union of portions of affine hyperplanes in ℝd called facets (each facet has positive (d-1)-dimensional measure). A closed subset F of Ω¯ is called a face if either

1. there exist T1,T2∈𝒯H such that F=∂⁡T1∩∂⁡T2∩Z, where Z is an affine hyperplane supporting a facet of both T1 and T2 (and F is termed interface), or

2. there exists T∈𝒯H such that F=∂⁡T∩∂⁡Ω∩Z, where Z is an affine hyperplane supporting a facet of both T and Ω (and F is termed boundary face).

Interfaces are collected in the set ℱHi, boundary faces in ℱHb, and we let ℱH:=ℱHi∪ℱHb. The diameter of a face F∈ℱH is denoted HF. For all T∈𝒯H, we define

the set of faces lying on the boundary of T; note that the faces in ℱT compose the boundary of T. For any T∈𝒯H, we denote by 𝒏∂⁡T the unit normal vector to ∂⁡T pointing outward T, and for any F∈ℱT, we let 𝒏T,F:=𝒏∂T∣F (by definition, 𝒏T,F is a constant vector on F).

We adopt the following notion of admissible mesh sequence; cf. [16, Section 1.4] and [20, Definition 2.1].

Two classical consequences of Definition 3.1 are that, for any mesh 𝒯H belonging to an admissible mesh sequence,

1. the quantity card⁢(ℱT) is bounded independently of the diameter HT for all T∈𝒯H (see [16, Lemma 1.41]),

2. mesh faces have a comparable diameter to the diameter of the cells to which they belong (see [16, Lemma 1.42]).

For any q∈ℕ, and any integer 1≤l≤d, we denote by ℙlq the linear space spanned by l-variate polynomial functions of total degree less than or equal to q. We let

Let a mesh 𝒯H be given. For any T∈𝒯H, ℙdq⁢(T) is composed of the restriction to T of polynomials in ℙdq, and for any F∈ℱH, ℙd-1q⁢(F) is composed of the restriction to F of polynomials in ℙdq (this space can also be described as the restriction to F of polynomials in ℙd-1q∘Θ-1, where Θ is any affine bijective mapping from ℝd-1 to the affine hyperplane supporting F). We also introduce, for any T∈𝒯H, the following broken polynomial space:

The term “broken” refers to the fact that no continuity is required between adjacent faces for functions in ℙd-1q⁢(ℱT). For any T∈𝒯H, we denote by (ΦTq,i)1≤i≤Ndq a set of basis functions of the space ℙdq⁢(T), and for any F∈ℱH, we denote by (ΦFq,j)1≤j≤Nd-1q a set of basis functions of the space ℙd-1q⁢(F). We define, for any T∈𝒯H and F∈ℱH, ΠTq and ΠFq as the L2-orthogonal projectors onto the spaces ℙdq⁢(T) and ℙd-1q⁢(F), respectively. Whenever no confusion can arise, we write, for all T∈𝒯H, all F∈ℱT, and all v∈H1⁢(T), ΠFq⁢(v) instead of ΠFq⁢(v∣F).

We conclude this section by recalling some classical results, that are valid for any mesh 𝒯H belonging to an admissible mesh sequence in the sense of Definition 3.1. For any T∈𝒯H and F∈ℱT, the trace inequalities

hold [16, Lemmas 1.46 and 1.49], as well as the local Poincaré inequality

where cP=π-1 for convex elements [8]; estimates in the non-convex case can be found, e.g., in [45]. Finally, proceeding as in [27, Lemma 5.6], one can prove using the above trace and Poincaré inequalities that

for integers 1≤s≤q+1 and 0≤m≤(s-1). All of the above constants are independent of the meshsize and can only depend on the underlying polynomial degree q, the space dimension d, and the mesh regularity parameter γ.

Henceforth, we use the symbol c to denote a generic positive constant, whose value can change at each occurrence, provided it is independent of the micro-scale ε, any meshsize HT or H, and the homogenized solution u0. We also track the direct dependency of the error bounds on the parameters α,β characterizing the spectrum of the diffusion matrix. The value of the generic constant c can depend on the space dimension d, the underlying polynomial degree, the mesh regularity parameter γ, and on some higher-order norms of the rescaling 𝔸β of the diffusion matrix or the correctors μl that will be made clear from the context.

4 Fine-Scale Approximation Space

Let k∈ℕ and let 𝒯H be a member of an admissible mesh sequence in the sense of Definition 3.1. In this section, we introduce the fine-scale approximation space on which we will base our multiscale HHO method. We first construct in Section 4.1 a set of cell-based and face-based basis functions, then we provide in Section 4.2 a local characterization of the underlying space, finally we study its approximation properties in Section 4.3.