2022 Non-intrusive implementation of a wide variety of Multiscale Finite Element Methods

Multiscale Finite Element Methods (MsFEMs) are now well-established ﬁnite element type approaches dedicated to multiscale problems. They ﬁrst compute local, oscillatory, problem-dependent basis functions that generate a suitable discretization space, and next perform a Galerkin approximation of the problem on that space. We investigate here how these approaches can be implemented in a non-intrusive way, in order to facilitate their dissemination within industrial codes or non-academic environments. We develop an abstract framework that covers a wide variety of MsFEMs for linear second-order partial diﬀerential equations. Non-intrusive MsFEM approaches are developed within the full generality of this framework, which may more-over be beneﬁcial to steering software development and improving the theoretical understanding and analysis of MsFEMs


Introduction
In this article, we consider highly oscillatory partial differential equations (PDEs) of the form L ε u ε = f, posed in a bounded domain Ω ⊂ R d , where L ε is a second-order linear differential operator with (possibly rough) coefficients that oscillate on a microscopic length scale of size ε much smaller than the diameter of Ω. (See Sec. 2 and 4 for a complete description of the problems that we study.)We seek a numerical approximation of u ε by applying a Galerkin approach.
It is well-known that standard (say P 1 ) finite element methods (FEMs) yield a poor approximation as a consequence of the highly oscillatory nature of the problem, unless a prohibitively expensive fine mesh is employed.An explicit example of this phenomenon is given, e.g., in [4,Example 1.1.].Dedicated multiscale approaches have thus been introduced, which provide a reasonably accurate approximation of u ε at a limited computational cost.Among the many multiscale approaches that have been proposed in the literature, we mention the Heterogeneous Multiscale Method (abbreviated HMM) [18], the Local Orthogonal Decomposition (LOD) [44] method, and the Multiscale Finite Element Method (MsFEM) [29], on which we focus here.We refer the reader to [1,19,36,4] for more comprehensive expositions of these multiscale methods.
The MsFEM is a finite element type approach that relies on the adaptation of the finite element space to the highly oscillatory differential operator, an idea that was first introduced in [5].The MsFEM was introduced in [29].It consists of two steps: 1.An "offline" stage, where highly oscillatory, problem-dependent basis functions are computed numerically as solutions to local problems (that mimic the reference problem on a subdomain).The local problems serve as a preprocessing step of the microstructure.
2. An "online" stage, where a Galerkin approximation of the reference problem, performed in the finite-dimensional space generated by the basis functions that are computed in the offline stage, is solved.This constitutes the global coupling between the local computations performed in the offline stage.
The MsFEM approach is particularly interesting for multi-query contexts, where the PDE of interest is to be solved repeatedly for multiple right-hand sides f (think e.g. of optimization problems or time-dependent problems where the time discretization results in a PDE in space to advance from one time step to the next).In this case, the basis functions, which depend on L ε , remain unchanged, so the offline stage is performed only once.The online stage solves a problem on a space of much lower dimension than a high fidelity space that fully resolves the microscale, resulting in a significant computational gain.
The specific choice of problem-dependent basis functions has led to various MsFEM variants in the literature (see, for instance, [29,19,38,39]).Although their implementation is rightfully considered relatively easy [19,45], these methods dictate an intrusive workflow, i.e., substantial modifications have to be made to the code of some existing finite element software to implement the MsFEM.For instance, all automatic computations involving integrals of standard (e.g.polynomial) basis functions have to be replaced by the corresponding computations for specialized basis functions for each specific problem.We explore this further in Sec.2.3.
The intrusive character of the MsFEM hinders the use of the method in many industrial contexts, because the time and tools required to adapt a legacy code that is currently in use may not be available.In this work, we propose a 'non-intrusive' MsFEM strategy that uses an existing legacy code for single-scale problems (based on standard finite elements), without any modifications, to obtain an accurate resolution of oscillatory PDEs.
The intrusiveness of the MsFEM is due to the following fact: the microstructure is preprocessed through the computation of specialized basis functions, thereby coupling the microstructure explicitly to the global numerical model.We circumvent this coupling with a novel formulation of the MsFEM basis functions.The resulting MsFEM strategy can be summarized as follows: 1.In the offline stage, so-called numerical correctors are computed as the solution to problems mimicking the reference problem on a subdomain.The numerical correctors are then used to average the microstructure in the form of effective, piecewise constant coefficients, leading to an effective PDE.
2. In the online stage, the effective PDE is solved by a standard FEM with the legacy code.
3. A post-processing stage is introduced to restore microscopic features in the macroscopic FEM result (obtained in the online stage) with the help of the numerical correctors computed in the offline stage.
Our non-intrusive MsFEM approach was introduced on a prototypical example in [12].In this article, we extend the findings of [12] to a large class of MsFEM variants.This requires the formulation of a general MsFEM framework, covering a generic definition of the MsFEM for the approximation of an abstract variational formulation of second-order linear PDEs.To the best of our knowledge, the question of how to make MsFEM approaches less intrusive has not been studied in the literature, except for the preliminary study in [12], and this work is a first step in that direction.We will comment on the intrusiveness of other multiscale methods in Sec.5.4.An overview of the contents of this article is as follows.In Sec. 2, we recall the basic principles of the FEM and the MsFEM and we explain the intrusive character of the MsFEM.We also review the non-intrusive MsFEM approach that was proposed in [12] for the simplest MsFEM variant on the example of a diffusion problem.We also highlight a link between the non-intrusive MsFEM approach and classical homogenization here.Then we summarize in Sec. 3 which properties of the MsFEM are essential for the non-intrusive workflow, and motivate the development of a general framework covering a wide variety of MsFEMs that we present in Sec. 4. We extend the nonintrusive MsFEM approach of [12] to our general MsFEM framework in Sec. 5.The non-intrusive MsFEM approach of [12] was found to be equivalent to a Petrov-Galerkin MsFEM (with P 1 test functions).This is no longer true for all MsFEMs covered by our general framework, and we obtain two non-intrusive MsFEMs: the Petrov-Galerkin MsFEM, which is completely equivalent to its non-intrusive implementation, and an approximate version of the Galerkin MsFEM that can be implemented in a non-intrusive way.The three essential formulas for the formulation of the nonintrusive MsFEM are highlighted in special boxes, both for the diffusion problem in Sec. 2 and for the general framework in Sec. 4 and 5.We then study the general MsFEM framework applied to diffusion problems in Sec.6, where we obtain a number of convergence results for the difference between the intrusive and non-intrusive MsFEM approaches.We conclude the article in Sec.7 by a numerical comparison of the intrusive and non-intrusive MsFEM approaches for diffusion problems, in order to assess the efficiency of our approaches for cases that are not covered by the convergence results of Sec. 6.Our results show that the Petrov-Galerkin MsFEM as well as the non-intrusive approximation of the Galerkin MsFEM are close to the original Galerkin MsFEM.Any possible additional error introduced by making the MsFEM non-intrusive is thus negligible.

Notation
In this article, we shall adopt standard notation for Sobolev spaces.In particular, the dual space of H 1 0 (Ω) is denoted H −1 (Ω).Further, for a given simplicial mesh T H of Ω, we use the notation H 1 (T H ) to denote the broken Sobolev space H 1 (T H ) = u ∈ L 2 (Ω) u| K ∈ H 1 (K) for all mesh elements K ∈ T H .
The standard norm for the space H 1 (Ω) is u H 1 (Ω) = u 2 L 2 (Ω) + ∇u 2 L 2 (Ω) and the correspond- . The space of functions whose restriction to each element of T H is a polynomial of degree k is denoted P k (T H ).
2 The (intrusive) multiscale finite element method

Discrete variational formulation
Let d ≥ 1 denote the space dimension of interest and let Ω ⊂ R d be a bounded polytope.Convexity can be assumed for elliptic regularity results to hold, for which we refer to [26].This technical assumption is not necessary for the algorithmic aspects of the MsFEM that are the main focus of this article.By way of example, we consider first the diffusion equation with homogeneous Dirichlet boundary conditions.In a second step, from Sec. 4 onwards, we will also consider more general problems, and we will mention other types of boundary conditions in Sec.5.3.More precisely, we focus here on the boundary value problem where the diffusion tensor A ε ∈ L ∞ (Ω, R d×d ) satisfies the uniform bounds for some M ≥ m > 0 independent of ε.The right-hand side f does not vary on the microscopic scale ε.We denote the diffusion tensor with a superscript ε to keep in mind that A ε might be highly oscillatory on a typical length scale of size ε much smaller than the diameter of Ω (assumed to be of order 1).No further structural assumptions on A ε are made.In particular, A ε need not be the rescaling of a fixed periodic matrix of the form A ε (x) = A(x/ε).We will specialize to this periodic setting in Sec.6.3 only to obtain convergence results, but this assumption is of no relevance for the practical implementation of the MsFEM.Let us also mention that none of the considerations in this article require symmetry of the diffusion tensor.Our development of non-intrusive MsFEMs also generalizes to linear systems of PDEs.The analysis we provide is also expected to extend to e.g. the system of linear elasticity up to some technicalities that we do not consider here.
For simplicity of exposition, we assume that f ∈ L 2 (Ω) (rather than f ∈ H −1 (Ω), for which the problem (1) is in fact well-posed).We do so to avoid unnecessary technicalities.Our proposed non-intrusive MsFEM carries over to the more general case.For some convergence results, the condition f ∈ L 2 (Ω) cannot be relaxed.In this case, this is also explicitly stated.
Problem (1) admits a unique solution in the space H 1 0 (Ω).This solution is also characterized by the variational formulation where the bilinear form a ε,diff and the linear form F are defined, for any u, v ∈ H 1 0 (Ω), by The coercivity hypothesis in (2) ensures that the bilinear form a ε,diff is coercive on the space H 1 0 (Ω).Then the Lax-Milgram Theorem [25,Theorem 5.8] shows that (3) is indeed well-posed.
The numerical approximation of (3) with a finite element method starts by the introduction of a mesh T H for Ω.The subscript H denotes the typical size of the mesh elements.We assume T H to be a simplicial, conformal mesh.For some convergence results, we shall assume quasi-uniformity.These assumptions are standard in finite element analysis.We refer, e.g., to [46,16,22] for a general exposition and various examples.Again, these regularity properties of the mesh do not have any impact on the implementation of the MsFEM on a given mesh.The regularity plays a role only to obtain convergence results.
A finite element method for (1) is obtained by restricting the equivalent formulation (3) to a finite-dimensional subspace of H 1 0 (Ω), typically consisting of functions that are piecewise polynomial on the mesh T H .We suppose that we are in the regime where H is larger than or comparable to the microscale ε and cannot be taken smaller due to computational limitations.In this case, it is well known that a Galerkin approximation of (3) on, say, the standard (conforming) Lagrange P 1 space on T H provides only a poor, not to say an incorrect approximation of u ε .See [4, Example 1.1], for instance, for an explicit example where the P 1 approximation on a coarse mesh fails.At the same time, the use of a finite element method on a fine mesh of size H ≪ ε might be unfeasible from a computational point of view because of its prohibitive computational cost.To remedy this issue, we shall next introduce the multiscale finite element method (MsFEM) [29,19].

A simple multiscale finite element method
The MsFEM is a Galerkin approximation of (3) for which the approximation space is adapted in order to achieve satisfactory accuracy even on a coarse mesh.The correct choice of approximation space yields a numerical approximation that is much closer to u ε than a standard P 1 -approximation when ε is smaller than H, and especially when ε becomes asymptotically small.To begin with, we introduce here the simplest variant of the MsFEM, which originally appeared in [29], before moving on to other MsFEM variants.
Let x 1 , . . ., x N 0 be an enumeration of the interior vertices of the mesh T H , i.e., the vertices that do not lie on ∂Ω.We denote by φ P 1 i the unique piecewise P 1 function such that φ P 1 i (x j ) = δ i,j for all 1 ≤ j ≤ N 0 .(These are the basis functions for the standard P 1 Lagrange finite element.)We define the multiscale basis functions All these problems, on each mesh element K, are again well-posed by coercivity of A ε and the Lax-Milgram Theorem.The functions φ ε i so defined belong to the global space H 1 0 (Ω) because the local boundary conditions on ∂K imply continuity across all mesh elements K.It is also immediately seen that φ ε i is supported by exactly the same mesh elements as φ P 1 i .Remark 1.On each mesh element K, problem (5) defines at most d + 1 non-trivial basis functions.Let i 1 , . . ., i d+1 be the indices of the vertices of K.It is easily inferred from (5) that . Thus, one only has to compute d basis functions by the resolution of a PDE on K.
The multiscale approximation space is defined as . This is a finite-dimensional space of the same dimension as the one used for a P 1 Lagrange finite element approximation on the mesh T H .The MsFEM consists in computing the approximation Since V ε H,0 is a subspace of H 1 0 (Ω), the bilinear form a ε,diff is coercive on V ε H,0 and the discrete problem ( 6) is again well-posed by virtue of the Lax-Milgram Theorem.
The computation of the multiscale basis functions φ ε i is called the offline stage of the MsFEM and only has to be carried out once if (1) has to be solved multiple times for different right-hand sides.Also note that all problems (5) are independent of each other, and can thus be solved in parallel.In practice, the φ ε i are approximated numerically on a fine mesh of K ∈ T H of mesh size h ≤ ε that resolves the oscillations of A ε .We omit these details here because they have no importance for the non-intrusive strategy that we shall propose later in this article.
The resolution of the global problem ( 6) is called the online stage.The computational cost for this problem is the same as for a standard P 1 approximation on the same mesh.A further discussion of the practical implementation of the MsFEM is provided in Sec.2.3.This discussion partially reproduces some elements of [12].We include it here to clarify and motivate the developments in the sequel.

Intrusive workflow
The practical resolution of the global problem (6) consists in the construction and resolution of the following linear system: with where we recall that N 0 denotes the number of interior vertices of T H .The MsFEM approxi- The MsFEM can then be written (as it is traditionally presented) as in Algorithm 1.We use the notation and we write Algorithm 1 MsFEM approach for problem (1) (see comments in the text) for 1 ≤ n ≤ d + 1 do

9:
Set j := N (l, K) end for 13: end for 15: end for 16: Solve the linear system Lines 1-12 of Algorithm 1 (resp.[13][14][15][16][17] constitute the offline (resp.online) stage of the MsFEM.Note that the computation of the stiffness matrix in line 11 only depends on the multiscale basis functions (and not on the right-hand side f ) and can therefore be carried out once and for all in the offline stage.Also note that, for an efficient computation of the φ ε i in line 6, one should apply Rem. 1.Only the online stage is to be repeated when problem (1) is to be solved multiple times for various right-hand sides f .Implementing Algorithm 1 in an industrial code is challenging.Indeed, the practical implementation of any finite element method relies on (i) the construction of a mesh, (ii) the construction of the linear system associated to the discrete variational formulation and (iii) the resolution of the linear system.An efficient implementation of the second step heavily relies on the choice of the discretization space.
Regarding the construction of the linear system (performed in line 11 of Algorithm 1), it is by no means obvious to adapt existing finite element codes based on generic approximation spaces (such as the spaces V L H and V CR H in Def.4.1) to a different, problem-dependent choice of space such as V ε H .No analytic expressions for the basis functions φ ε i are available (and thus a fine mesh should be used to approximate them), the computation of a ε K (φ ε i , φ ε j ) should be performed by quadrature rules on the fine mesh because the integrands are highly oscillatory, one should have at hand the correspondence between element and vertex indices of the coarse mesh (N (n, K) in Algorithm 1), the assembly of the global stiffness matrix {A ε j,i } 1≤i,j≤N 0 should be executed by a dedicated new piece of software, etc.To alleviate these obstacles, we introduce below a way of implementing the MsFEM that capitalizes on an existing code for solving (1) by a P 1 approximation on T H in the case of slowly varying diffusion coefficients.The three central identities for our approach that we aim to generalize to other MsFEMs in this article are framed in distinctive boxes.

Effective problem on the macroscopic scale
Let us consider the construction of the stiffness matrix of the MsFEM in more detail.The stiffness matrix defined in (8) requires the computation of the quantities for all 1 ≤ i, j ≤ N 0 .Following [12], we rewrite the multiscale basis functions as for all 1 ≤ i ≤ N 0 , where, for each mesh element K, we define the numerical corrector χ ε,α K ∈ H 1 0 (Ω) (1 ≤ α ≤ d) as the function supported by K that is the unique solution to the local problem Here, e α denotes the α-th canonical unit vector of R d .The expansion (10) is obtained upon rewriting (5) as a PDE for φ ε i − φ P 1 i , and then using linearity of the PDE and the fact that ∇φ K is indeed the unique solution to this PDE.
Inserting (10) for the trial and test functions in (9) and again exploiting the fact that all φ P 1 i have piecewise constant gradients, we obtain Next we define the piecewise constant effective diffusion tensor A ∈ P 0 (T H , R d×d ) by where |K| denotes the measure of the mesh element K. Then (9) can be written as Motivated by (13), we introduce the coarse-scale problem and its Galerkin discretization with P 1 Lagrange elements: with V H,0 = span φ P 1 i | 1 ≤ i ≤ N 0 (note that the definition of V H,0 will be generalized in Def.4.1), find u H,0 ∈ V H,0 such that where the linear form F is defined in (3) and the bilinear form a diff is defined as Problem (15) equivalently writes with Comparing the expressions ( 13) and ( 17), we deduce that A ε = A P 1 .In other words: Lemma 1.The stiffness matrix of the MsFEM problem (6) is identical to the stiffness matrix of the P 1 problem (15).
This lemma immediately implies that the P 1 problem ( 15) is well-posed, since the MsFEM (6) itself is well-posed.
Let us point out that problems (14) and (15) are defined entirely in terms of quantities that vary only on the macroscopic scale H.The finite element problem (15) can thus be solved using a legacy code that is designed for standard FEMs.Lemma 1 then suggests including the P 1 approximation (15) of the effective, coarse-scale problem (14) as an integral part of the MsFEM approach.We do so in Algorithm 2 below.
The right-hand side vector F ε in ( 8) is, in general, different from F P 1 in (17).Indeed, we integrate the product of f with highly oscillatory basis functions in the former problem and with P 1 basis functions in the latter.The solutions U ε and U P 1 to (7) and (16), respectively, are thus different a priori.

Non-intrusive workflow
We propose the following non-intrusive MsFEM variant: where u H ∈ V H,0 is the unique solution to (15).(18) The MsFEM approximation u ε H is well-defined, since we have seen above that problem (15) is well-posed.
Note that the symbol u ε H shall be used here and in the sequel for the solution to various MsFEMs variants to alleviate the notation.The exact MsFEM will be specified by the context.We will use distinct notation for different MsFEM variants when required for clarity.
The most efficient way to compute u ε H from u H is not as stated here, however.The evaluation of u H (x) may require the determination of the degrees of freedom associated to the simplex K to which x belongs.This demands the use of the internal mechanisms of the legacy code that is used to compute u H .The use of the legacy code can be avoided by expanding u H as follows.For any affine function ϕ on K, we have where x α denotes the function that to a point x ∈ Ω associates its α-th coordinate, and x c,K = (x 1 c,K , . . ., x d c,K ) is the centroid of K.If one uses the legacy code to store the values of u H (x c,K ) and ∂ α u H element by element at the end of the online stage, then u ε H defined in (18) can be computed element by element according to without using the legacy code.
The above observations culminate in the computational approach presented in Algorithm 2. We can distinguish (1) the offline stage consisting of lines 1-7, (2) the online stage being executed entirely in line 8, (3) a post-processing step in line 9.
The superiority of Algorithm 2 over the classical MsFEM Algorithm 1 is that the global problem of the online stage can completely be constructed and solved by the use of a pre-existing P 1 PDE solver.The only requirements for the legacy code are the functionality to provide piecewise constant diffusion coefficients to the solver and the existence of a procedure to store the value of the P 1 solution and its gradient at the centroids of the mesh.An additional advantage in the online stage is that the construction of the right-hand side F P 1 (see (17)) for the global problem only requires a numerical quadrature on the coarse mesh and is therefore cheaper than the construction of F ε (see (8)), involving the multiscale basis functions and requiring numerical quadratures at the microscale.
The part of the offline stage that manipulates fine meshes (lines 2-7) and the post-processing step can be developed independently of the legacy code used in line 8.The requirement for these fine-scale solvers is that they have access to the coarse mesh T H used by the global solver.Note Algorithm 2 Non-intrusive MsFEM approach for problem (1) 1: Let T H be the mesh used by the legacy code 2: for all K ∈ T H do 3: Solve for χ ε,α K defined by (11) 5: end for 6: Compute A| K defined by (12) 7: end for 8: Use the legacy code to solve for u H defined by (15) and to save {u 9: Obtain the MsFEM approximation u ε H by (20) also that the local problem ( 11) is only indexed by the coarse mesh element K, in contrast to the local problem (5) that is indexed both by the coarse mesh element K and the vertex index i.For the latter problems, one has to know, for each element K, the global index that corresponds to the vertices of K, a piece of information that may be difficult to access in a legacy code.For the problems (11), this correspondence is not needed to compute A, nor for the computation of the fine-scale solution u ε H in (20), both of which are entirely defined element-wise.

Interpretation of the non-intrusive MsFEM
We emphasized above that the right-hand sides of the linear system for the MsFEM in (8) and the linear system solved for the non-intrusive MsFEM in (17), are different in general.This motivates the comparison of the non-intrusive MsFEM approach (18) to the following Petrov-Galerkin MsFEM: based on the trial space V ε H,0 and the test space V H,0 for both the bilinear and the linear form.The following result was shown in [12].
The non-intrusive MsFEM approach is generalized in Sec. 5 after the development of a general framework to define a wide variety of MsFEMs in Sec. 4. Lemma 2 does not generalize to the full framework.We will see the conditions under which the non-intrusive approach leads to a Petrov-Galerkin MsFEM in Lemma 8.

Relation to homogenization theory
We highlight in this section the fact that many ingredients of our non-intrusive MsFEM approach are reminiscent of standard quantities of homogenization theory, or the theory of H-convergence, which studies the limit of a sequence of solutions u ε to a PDE as ε tends to 0. This relation to H-convergence provides an interesting interpretation of the effective tensor A introduced in (12).
Let us suppose in this section (and in this section only, except for Sec.6.3) that A ε (x) = A per (x/ε) for some bounded, Z d -periodic matrix A per satisfying the coercivity property in (2).In this case, the sequence of matrices A ε has a homogenized limit that is explicitly known.(An explicit characterization of the limit is not available for H-convergence in general.)We summarize the main results below.See, for instance, [10,47] or [2, Chapter 1] for details on periodic homogenization.
Due to H-convergence of A ε , the functions u ε , solution to (1), converge to a limit function u ⋆ (weakly in H 1 (Ω), strongly in L 2 (Ω)) as ε → 0. The homogenized limit u ⋆ is the solution to the homogenized equation (24) below.
Let Q denote the unit cube of R d .We introduce the corrector functions w 1 , . . ., which uniquely defines w α up to an irrelevant additive constant.The entries of the (constant) homogenized diffusion tensor A ⋆ ∈ R d×d are given by The homogenized limit u ⋆ of u ε is the unique solution in H 1 0 (Ω) to the boundary value problem The truncated reconstruction of u ⋆ that is called the first-order two-scale expansion takes the form Under suitable regularity assumptions, the difference u ε − u ε,1 converges to 0 strongly in H 1 (Ω) as ε → 0. This property will be used for the convergence results in Sec.6.3.
In the periodic setting, the expansion (25) can be used to construct a numerical approximation of u ε , without the need of any computations at the fine scale.This approximation is presumably valid only in the regime of very small parameters ε and deteriorates if ε grows.Moreover, in more general settings, the corrector functions are not local nor explicit, for their definition involves a PDE posed on the whole domain Ω and that depends on an effective tensor that is itself defined in terms of the corrector functions.Details can be found, e.g., in [3,2,43,47].This prevents the H-convergence theory from being directly applicable for the numerical approximation of u ε .
Numerical homogenization techniques, that draw their inspiration from the various elements above, offer an alternative for the approximation of u ε that can be applied in much more general contexts.We can see the similarities between the corrector functions w α in (22) and the numerical correctors χ ε,α K in (11).Note that the χ ε,α K solve problems similar to (22), but that they need to be solved at the microscale and on each mesh element K. Similarly, we note the resemblance between the reconstruction (25) and the definition of u ε H in (18), and between the homogenized coefficient A ⋆ defined in (23) and the effective macroscopic coefficient A from (12).However, contrary to A ⋆ , the MsFEM quantity A has to be computed on an element-by-element basis, and it is not necessarily constant throughout Ω.Finally, the MsFEM analogue of the homogenized problem (24) is the resolution of the effective macroscale problem (14).
Example 1.A very particular setting, although academic in nature and only useful for pedagogical purposes, actually leads to an MsFEM approximation that is exactly equivalent to a discretization of the periodic homogenization setting.Consider (1) in 2D posed on the unit square.Let us consider a mesh consisting of squares that are perfectly aligned with the periodicity of A ε .We solve the corrector problems (11) on all square mesh elements with periodic boundary conditions and subsequently compute the effective diffusion tensor A according to (12).In this case, the problems for the numerical correctors all reduce to (22) and A is constant and equal to the homogenized coefficient A ⋆ as defined by (23).A Q 1 discretization of the effective problem ( 14) thus constitutes a non-intrusive MsFEM that is equivalent to the Q 1 approximation of the homogenized equation (24).
3 Why develop a general framework?
In the sequel we develop a general framework for a wide variety of MsFEMs in an abstract setting.We motivate here why this general framework for MsFEMs is useful.

Local boundary conditions
First, let us explain why various different MsFEMs have been proposed in the literature.One reason is that different equations than (1) (e.g.advection-diffusion equations) give rise to different choices of the local problem (5), depending on which terms of the global PDE are included (see, for instance, [39].) The other reason is that, even for the pure diffusion problem (1), the choice of the basis functions defined in (5) has an important drawback.The definition of the multiscale basis functions requires a choice of arbitrary boundary conditions on the mesh element boundary ∂K, since the exact boundary condition satisfied by u ε is unknown.In (5), affine boundary conditions are imposed.In view of this choice, we shall refer to the MsFEM defined above as the 'MsFEM-lin'.
The MsFEM-lin cannot yield an accurate representation of u ε near ∂K if A ε is highly oscillatory and the mesh T H is coarse.Variations on the definition of the functions φ ε i have been proposed to improve the MsFEM.Here we summarize the ideas of oversampling and of MsFEM à la Crouzeix-Raviart, which together inspire the formulation of a general MsFEM framework in Sec. 4.
The oversampling variant of the MsFEM was introduced along with the variant based on (5) at the time of its first appearance in [29].For this method, an oversampling domain S K is associated to each mesh element K (details are provided in Sec.4.3.1).The problems (5) are solved on the larger domain S K rather than K, so the inadequate boundary conditions are pushed away from the actual mesh elements.To construct the multiscale basis functions, the resulting functions on S K are restricted to the actual mesh elements K and suitably combined around each vertex x i .The new multiscale basis functions oscillate on ∂K if the oversampling patch is taken large enough.We note that, in general, this strategy leads to discontinuous basis functions.Hence, the finite element space obtained is no longer conforming.
The MsFEM with Crouzeix-Raviart type boundary conditions for the local problems (which we shall abbreviate as 'MsFEM-CR') was introduced in [38].It uses basis functions associated to the edges of the mesh (in contrast to the MsFEM-lin presented above, and its oversampling variant, where basis functions are associated to the vertices of the mesh).A typical basis function satisfies the following on ∂K: the flux through each face of K is constant, and the constants are determined by the condition that the average of the basis function be 1 over one particular face and 0 over all other faces.Again, this is a way to avoid imposing any conditions on the trace of the basis function directly.The multiscale functions can thus be oscillatory on the faces of the mesh.As is the case for oversampling methods, the resulting finite element space is nonconforming.
All of these variations, applied to any MsFEM for linear second-order PDEs, are covered by the general MsFEM framework that we develop in Sec. 4.

The non-intrusive approach
The intrusiveness of the specific MsFEM-lin variant introduced in Sec.2.2 is exemplary for all Ms-FEMs described in Sec.3.1.It turns out that the non-intrusive MsFEM approach introduced in [12] and recalled in Sec.2.5 can also be generalized to all these MsFEM variants.We summarize the key ingredients that allow for the formulation of the non-intrusive MsFEM approach of Algorithm 2 (corresponding to the identities in boxes in Sec.2.4).
The non-intrusive MsFEM follows from the expansion (10), namely the expression of the multiscale basis function as a P 1 basis function and a linear combination of numerical correctors that are fully localized.We note that • the full localization of the numerical correctors defined in (11) allows the preprocessing of the microstructure independently of the global approximation indices related to the finite element method; • the expansion (10) follows from the fact that ∇φ P 1 i is piecewise constant combined with linearity of the local problems (5); • the stiffness matrix can be formulated in terms of a piecewise constant effective diffusion tensor in ( 13) thanks to full localization of the corrector functions, the piecewise constant gradient of φ P 1 i in the expansion (10) and bilinearity of the global problem (3).
These observations provide the main structure of the general framework.First, we choose an underlying, low-dimensional space of piecewise affine functions to which the MsFEM is associated (Def.4.1).This will be the standard conforming Lagrange space of order 1 (for the MsFEM-lin), or the Crouzeix-Raviart space of order 1 (for the MsFEM-CR).Second we need to formulate the local problems for the numerical correctors (Def.4.8 and Def.4.10).This involves the definition of oversampling patches (for MsFEMs with oversampling, Def.4.4), and an extension of the notion of degrees of freedom to define the boundary conditions for the numerical correctors (Def.4.2, 4.5 and 4.6) on oversampling patches.It is then possible to define the multiscale basis functions as a generalization of (10) (see Def. 4.11) and finally to define the MsFEM for our general framework in Def.4.12.
Remark 2. We note that our development of non-intrusive MsFEM approaches relies to a great extent on the fact that (10), and its generalization (33) in the general framework developed below, provide a description of the multiscale basis functions in terms of P 1 basis functions, without the need of higher-order functions.Higher-order MsFEMs can be found in [3,28].Possible analogues of (33) for such MsFEMs and the subsequent techniques to design a non-intrusive MsFEM variant are more involved and may be the topic of future work.See [11].

Other motivations for the general framework
Besides a unified formulation of our non-intrusive MsFEM approach, our general framework can also be beneficial to concrete code development for the MsFEM.Common features among various multiscale methods have previously been used to design flexible and efficient software for the implementation of such methods on the DUNE platform [7,6] within the Exa-Dune project [9].For example, the distribution of local problems over multiple processors and subsequent coupling in a global problem are handled by designated software components [8].Our work may contribute to the efficient implementation of all MsFEMs covered by our general framework in such a project and similar endeavours yet to come.
When formulating the general framework, we also clarify a few practical matters that are often left pending in the various research articles we are aware of.In particular, we give a rigorous definition of the oversampling procedure near the boundary ∂Ω of the global domain.
As we explore the general framework, we will also propose an MsFEM variant that has not yet appeared in the literature: the MsFEM-CR combined with the oversampling technique (see Example 9).We hope that our framework may also further the development of new MsFEM variants in an attempt to improve on the shortcomings of the methods known today.
Finally, the present study may also uncover a deeper understanding of MsFEMs by paving the way to a unified convergence analysis of different variants.This work is currently in preparation.

Abstract definition of the MsFEM
We develop here a general framework for multiscale finite element methods.The ultimate aim is to generalize the key identities of Sec.2.4.This is done in Def.4.8 and 4.10 for the numerical correctors introduced in (11), and in Def.4.11 for the expansion (10) of the multiscale basis functions.This allows the reformulation of the linear system of the MsFEM as the linear system of an effective problem in (41) (for a Petrov-Galerkin MsFEM) and ( 46) (for a Galerkin MsFEM) in Sec. 5.The other notions introduced in this section, although rather technical and abstract, are necessary tools to capture a wide variety of MsFEMs in our general framework.

The continuous problem
The abstract variational problem for our general MsFEM framework is as follows.Let a ε be a continuous bilinear form on H 1 (Ω) × H 1 (Ω).We are interested in the solution to the problem where F is defined as in (4) for any f ∈ L 2 (Ω).To ensure well-posedness of (26), we suppose that the bilinear form a ε is coercive on H 1 0 (Ω).The bilinear form a ε may contain coefficients that oscillate on a microscopically small scale.
The oversampling and Crouzeix-Raviart variants of the MsFEM introduced in Sec.3.1 show that we need to accommodate for approximation spaces with discontinuities at the interfaces.This requires some additional assumptions on the formulation of the abstract problem.We suppose that the bilinear form a ε is in fact defined on the broken Sobolev space H 1 (T H ) × H 1 (T H ).More precisely, we assume that we can represent it as To ensure well-posedness of MsFEMs, which may use nonconforming approximation spaces, coercivity on H 1 0 (Ω) may be insufficient.Therefore, we add the following coercivity hypothesis for the bilinear forms a ε K : for all K ∈ T H , there exists In order to perform a convergence analysis, one also has to assume that the α K are bounded from below by some α > 0 that does not depend on H.We provide convergence results in Sec.6 for the pure diffusion problem (1), in which case we have α K = m from (2).
As an example, the introductory problem (1) with the associated bilinear form a ε,diff is covered by this framework as is made explicit in Example 2 below.Other second-order PDEs that fit in our abstract variational formulation are given in Example 3.
Example 2. The diffusion problem (1) is covered by the abstract variational formulation above.Indeed, we can set a ε = a ε,diff , where a ε,diff is the bilinear form defined in (4).Further, we define with α K = m the coercivity constant from (2).
Example 3. The reaction-advection-diffusion equation, with a divergence-free advection field b : Ω → R d and a positive reaction coefficient σ : Ω → R, can be modelled (under some regularity hypotheses that we do not state here) with the bilinear forms However, these bilinear forms a ε K do not satisfy (27) even though the bilinear form a ε is coercive on H 1 (Ω).To this end, a skew-symmetrized formulation of the transport term can be used.The skew-symmetrized formulation uses the bilinear form which does satisfy (27).However, the assumption ( 27) is used for proving well-posedness of the MsFEM in Lemma 5, but both choices for a ε K mentioned here can be studied in practice.We refer e.g. to [35,40,11] for more details.Within the general MsFEM framework, b and σ are allowed to be highly oscillatory, and this may impact the specific MsFEM strategy to be preferred.

Piecewise affine structure
In Sec.2.2, we have seen that the relation between multiscale basis functions and piecewise affine functions is essential for the development of our non-intrusive MsFEM.For the MsFEM definition in the general framework, we start by choosing such a structure in the following definition.Definition 4.1.Let a mesh T H be given.The underlying P 1 space for the MsFEM, denoted V H , is one of the following two spaces: the Lagrange approximation space in which case we shall refer to the associated MsFEM as the MsFEM-lin, or the Crouzeix-Raviart approximation space in which case the associated MsFEM shall be called the MsFEM-CR.We use the notation F (K) for the set of faces of K and v denotes the jump of v over the face e.The space V L H is a subspace of H 1 (Ω), but V CR H is not.Note that no restrictions apply on faces lying on ∂Ω.
We note that the underlying P 1 space has the following property: if v ∈ V H is piecewise constant on the mesh T H , then v is constant in Ω. Contrary to the space V L H , functions in the Crouzeix-Raviart space V CR H are discontinuous in general.They are continuous, however, at the centroids of all faces of the mesh.
For standard finite elements, the notion of degrees of freedom allows to characterize any finite element function.The idea of the MsFEM is to preserve this notion of degrees of freedom (in a suitable way made precise below) in the definition of a multiscale approximation space, while adapting the piecewise affine structure to the microstructure of the PDE.We formalize this notion for the two underlying P 1 spaces that we introduced in Def.4.1.The definition involves an arbitrary simplex K, which is typically an element of the mesh T H , or an associated oversampling patch (for the oversampling technique of the MsFEM) that we shall define in Def.4.4.The latter is not always a simplex, and we extend Def.4.2 to such oversampling patches in Def.4.5 and 4.6.Definition 4.2.A degree of freedom operator (DOF operator) Γ associates to any simplex K ⊂ R d and v ∈ P 1 (K) a vector Γ(K, v) ∈ R d+1 , whose components are called the degrees of freedom of v on K, in such a way that the application v → Γ(K, v) is a linear bijection from P 1 (K) to R d+1 .More precisely, Γ(K, •) will denote in the sequel one of the following two operators: 1. (DOF operator for the MsFEM-lin.)Let x 0 , . . ., x d denote the vertices of K. We set for a vertex x of the mesh that lies on ∂Ω.
2. (DOF operator for the MsFEM-CR.)Let e 0 , . . ., e d denote the faces of K. We set For K ∈ T H , the degree of freedom [Γ CR (K, •)] j is said to be associated to the boundary for a face e of the mesh that lies on ∂Ω.
The P 1 test space is defined as •)] j is associated to the boundary, .
The P 1 test space is used in practice to approximate the subspace H 1 0 (Ω) of H 1 (Ω).The degrees of freedom are defined element per element and are thus local.Global properties of the underlying P 1 space V H are most easily made explicit through the identification of a basis for V H . Definition 4.3.Let V H be an underlying P 1 space as in Def.4.1, and let Γ be the associated DOF operator.We shall denote by N the dimension of V H .The P 1 basis functions φ P 1 1 , . . ., φ P 1 N are defined as follows: • For the MsFEM-lin, let x 1 , . . ., x N be an enumeration of the (internal and boundary) vertices of T H . Then φ P 1 i is defined by φ P 1 i (x j ) = δ i,j for all 1 ≤ i, j ≤ N .
• For the MsFEM-CR, let e 1 , . . ., e N be an enumeration of the (internal and boundary) faces of T H . Then φ P 1 i is defined by In both cases, these functions form a basis of the corresponding space V H of Def.4.1.

Oversampling patches
To replace the (standard) underlying P 1 space by a space of the same (low) dimension, adapted to the microstructure of a ε , we associate to each mesh element K ∈ T H an oversampling patch.It serves to avoid imposing artificial, non-oscillatory boundary conditions on K directly when computing numerical correctors to process the microstructure.
Definition 4.4.Let K ∈ T H be any mesh element and let S ′ K be a simplex obtained from K by homothety around the centroid of K with homothety ratio ρ ≥ 1.The oversampling patch S K is defined as S K = S ′ K ∩ Ω. See Fig. 1 for an illustration of the construction of oversampling patches in dimension 2. In this work, we allow for the trivial homothety ratio ρ = 1.In this case, the patch S K coincides with K.
We will call an MsFEM without oversampling an MsFEM for which all oversampling patches satisfy S K = K.Otherwise, the MsFEM is called an MsFEM with oversampling.We speak simply of an MsFEM when there are no assumptions on the oversampling patches.
Figure 1: Oversampling patches for MsFEM in 2D.Left: the patch for the mesh element K is obtained from K by homothety.Right: The triangle S ′ K partially lies outside the domain Ω and the oversampling patch S K is not homothetic to K. It is not even a triangle.
For most mesh elements K, the patch S K in Def.4.4 is a simplex.However, for mesh elements close to the boundary ∂Ω, alternative constructions should be considered.We have not found any explicit description of such a construction in the literature.This complicates the reproducibility of the method as well as a rigorous convergence analysis.The precise definitions of this section provide a first step to address these issues.A fully rigorous convergence analysis of the MsFEM with oversampling as described here is the subject of ongoing investigations [11].

Degrees of freedom on oversampling patches
Definition 4.2 provides the definition of DOF operators on any simplex.For the MsFEM, we wish to compute multiscale functions on oversampling patches S K , in which case Def.4.2 may be insufficient.We illustrated this in Fig. 1b.Indeed, the number of vertices/faces of the oversampling patch may be larger than d + 1.In order to associate a multiscale basis function to every P 1 basis function, we will still need a notion of DOF operator such that Γ(S K , •) is a linear bijection from P 1 (S K ) to R d+1 .Therefore, we extend the definition of the degrees of freedom operators Γ L and Γ CR in Def.4.5 and 4.6.Definition 4.5.Let K ∈ T H and let S K be its associated oversampling patch.Let x 0 , . . ., x d be a selection of d + 1 distinct vertices of S K .We define the DOF operator Γ L by We note that any choice of d + 1 nodal values unequivocally characterizes an affine function on S K .Hence, Γ L (S K , •) is indeed a bijection.Since Γ(S K , •) will only be used to identify elements of P 1 (S K ) on the boundary of S K , the particular choice of the vertices x 0 , . . ., x d is unimportant in the above definition.Finally, when S K is a simplex, it has only d + 1 vertices and Def.4.5 reduces to Def. 4.2.
To generalize the notion of degrees of freedom for the Crouzeix-Raviart space to non-simplicial patches, we need to introduce some additional notation.On the boundary of a non-simplicial oversampling patch, we can identify some faces that collapse to a single vertex if we shrink S K to K. We call these faces the additional faces and denote the set containing them by F a (S K ).The other faces of S K are referred to as the dilated faces, collected in the set F d (S K ).When the patch S K does not touch ∂Ω, we have F d (S K ) = F (S K ) and F a (S K ) = ∅.In Fig. 2a, for example, the additional faces are exactly those faces that lie on ∂Ω.This is not always the case, as is illustrated by Fig. 2b.
For the definition of Γ CR (S K , •), we shall rely on the existence of d + 1 dilated faces, because we need Γ CR (S K , •) to be a bijection between P 1 (S K ) and R d+1 .This imposes a constraint on the choice of the homothety ratio used to construct S K .For example, in the case of Fig. 2a, the lower right dilated face falls outside Ω if the homothety ratio is too large, and the oversampling patch S K only has two dilated faces (edges here) and two additional faces.
Figure 2: Non-simplicial oversampling patches in 2D.The dilated edges of the patch S K , those that 'correspond' to the edges of the original triangle K, are dashed.Definition 4.6.Let K ∈ T H and let S K be its associated oversampling patch.We assume that S K has d + 1 dilated faces, and we denote hem by e 0 , . . ., e d .We define the DOF operator Γ CR by When S K is a simplex, we have F d (S K ) = F (S K ), and Def.4.6 coincides with the respective elements of Def.4.2.

Numerical correctors: first oversampling strategy
We now provide the precise assumptions under which we will consider local problems, i.e., the analogues of ( 5) defining the MsFEM-lin basis functions and the definition of the numerical correctors in (11).In fact, since the numerical correctors play an essential role in the construction of nonintrusive MsFEM approaches, we define the numerical correctors first and use them to define the multiscale basis functions in Def.4.11.
The two possible functional settings for all these constructions are provided by Def.4.7 and 4.9.Both definitions present a different implementation of the oversampling technique.These definitions involve a 'sampling space', whose name is inspired by the idea that only a limited number of local problems will be solved to encode the microstructure of the PDE in the numerical model.The choice of sampling space has to accommodate for the boundary conditions that one wishes to impose on the numerical correctors and basis functions (e.g.essential or natural; see Examples 4 and 5).Definition 4.7.Let K ∈ T H , let S K be its associated oversampling patch and let Γ be a DOF operator from Def. 4.2.A subspace V K of H 1 (S K ) and bilinear form s ε K : V K × V K → R are called sampling space and sampling form, respectively, if they satisfy the following: 1. the space V K contains the space of affine functions P 1 (S K ); 2. the operator Γ(S K , •) is well-defined on V K ; 3. the DOF -extended local problem: find v ∈ V K such that has a unique solution for any g ∈ (H 1 (S K )) ′ .Here, Problem ( 28) is called 'DOF -extended' because the degrees of freedom, controlling the boundary conditions associated to the local problem, are imposed on the oversampling patch S K rather than the (generally smaller) mesh element K.
The sampling form s ε K shall be used to encode the oscillations of the bilinear form a ε and thus the microstructure of the problem in the multiscale finite element functions.There is some flexibility in choosing the sampling form; one may to choose to include all the same terms as those in the bilinear form a ε K of the original problem (26), or only some of them.When the MsFEM was first proposed in [29], it was suggested that s ε K should include those terms that correspond to the highest-order terms of the PDE that is to be solved.In the context of the advection-diffusion equation, one may thus choose to include only the diffusion terms, but also the advection term can be included in our MsFEM framework.Both options have been studied e.g. in [39,40].
In the functional setting of Def.4.7, the generalization of (11) to define the numerical correctors for the general MsFEM framework is as follows.
Definition 4.8.We introduce for all K ∈ T H , for any 0 ≤ α ≤ d, the function χ ε,α,e S K ∈ V K,0 as the unique solution to the corrector problem The DOF -extended numerical corrector χ ε,α,e K is defined as the restriction of χ ε,α,e S K to K, extended to all of Ω by 0. Note that the above definition introduces one more numerical corrector than introduced in (11).The precise definition of the numerical correctors is chosen such that the analogous expansion of (11) for the general framework (see (33)) leads to a PDE for the multiscale basis functions analogous to (5); we show this in Lemma 3 and (for a second oversampling strategy introduced below) in Lemma 4. In the following example, we see that Def.4.8 is indeed a generalization of the numerical correctors defined by (11) in Sec.2.4.
Example 4 (MsFEM-lin for diffusion problems).We consider V H = V L H and Γ = Γ L from Def. 4.2.For the diffusion problem (1), we have a ε = a ε,diff and we set . The sampling space for the MsFEM-lin is defined as Then the sampling test space V L K,0 is the space H 1 0 (S K ).In this case, it holds a ε,diff K (1, w) = 0 for all w ∈ V L K,0 .Consequently, the DOF -extended numerical corrector χ ε,0,e K is identically equal to 0; we obtain indeed exactly d numerical correctors as in Sec.2.4.For the non-trivial numerical correctors, Def.4.8 corresponds to the weak formulation of the following boundary value problem: which is clearly well-posed.
Example 5 (MsFEM-CR for diffusion problems).Taking a ε , a ε K and s ε K as in the previous example, we construct the MsFEM-CR with the sampling space where n denotes the outward unit vector on ∂S K and c h is a constant whose value is uniquely determined by the above problem.We note that the condition for the flux on the additional faces of S K is entirely determined by the right-hand side in (29), whereas the flux on the dilated faces of S K involves an additional constant, due to the fact that the test functions in V CR K,0 cannot take arbitrary values on the dilated faces.Indeed, their mean vanishes on these faces according to Def. 4.7.
When S K = K and when the faces of K do not lie on ∂Ω, this corresponds to the setting of the original MsFEM-CR defined in [38].The latter work also provides an alternative characterization of the multiscale Crouzeix-Raviart space.When a face e of K lies on Ω, the basis functions that we will obtain do not satisfy φ ε e = 0 on e, but only satisfy a weak boundary condition in the average sense on e.This does not correspond to the original definition of the MsFEM-CR in [38,37].The MsFEM-CR with local boundary conditions as defined here was studied in [17,42,34].
When S k = K (i.e., in the absence of oversampling), the DOF operator allows us to prescribe certain continuity properties on the faces of the mesh elements K.More precisely, when the MsFEMlin with DOF operator Γ L is employed, the numerical correctors χ ε,α,e K vanish at the vertices of the mesh, and, with the correct choice of sampling space (see Example 4), they vanish on all faces of K and are thus continuous on Ω.When the MsFEM-lin with DOF operator Γ CR is considered, we obtain weak continuity of the numerical correctors over all faces of the mesh.The definition of the multiscale basis functions that we give below (see Def. 4.11, in the vein of the expansion ( 10)) shows that this means that the continuity properties of the P 1 basis functions of the underlying P 1 space are not perturbed when building the multiscale basis functions.
In the general case, when the oversampling patch S K is larger than K, we cannot preserve any of these continuity properties if we use DOF -extended local problems for our local computations, since the values on ∂K are not controlled by the degrees of freedom Γ(S K , •) on ∂S K .Therefore, we introduce another variant of the local problems and the DOF -extended numerical correctors in the next section.

Numerical correctors: second oversampling strategy
Definition 4.9.Let K ∈ T H and let V K and s ε K be a sampling space and sampling form, respectively, according to Def. 4.7.Additionally, suppose that the operator Γ(K, •) is well-defined on V K .Then a DOF -continuous local problem is to find v ∈ V K such that for some g ∈ (H 1 (S K )) ′ .
Definition 4.10.Suppose any DOF -continuous local problem in Def.4.9 is well-posed.Then we introduce, for all K ∈ T H and all 0 ≤ α ≤ d, the functions χ ε,α,c S K as the unique functions in V K with Γ K, χ ε,α,c S K = 0 satisfying the corrector problem (29).We define the DOF -continuous numerical correctors χ ε,α,c K as the restriction of χ ε,α,c S K to K, extended to all of Ω by 0. We emphasize that the local problems of Def.4.7 and 4.9 use test functions w in the same space V K,0 .This means that the test functions satisfy Γ(S K , w) = 0, rather than Γ(K, w) = 0. Remark 3. Clearly, when S K = K, there is no difference between the DOF -extended and DOF -continuous problems (28) and (31).We shall in this case simply refer (28) (or (31)) as local problems, and we write χ ε,α,• K = χ ε,α K for the numerical correctors of MsFEMs without oversampling.
Example 6 (MsFEM-lin for diffusion problems).Continuing Example 4, consider now the DOFcontinuous numerical corrector χ ε,α,c K .Equation (36) solves the following problem: there exists w ∈ P 1 (S K ) such that The boundary condition on ∂S K is 'completed' by a condition at the vertices of K. Except when A ε is constant (and, consequently, φ ε i is affine on S K ), it is not evident whether a solution to this problem exists.
Example 7 (MsFEM-CR for diffusion problems).For the MsFEM-CR considered in Example 5, the DOF -continuous numerical correctors satisfy the same problem (30) as the DOF -extended numerical correctors, but with the average condition replaced by As we saw for the MsFEM-lin, this is not a standard boundary value problem on S K .
Examples 4 and 5 show that a DOF -extended local problem is typically equivalent to a PDE with boundary conditions on S K .Under reasonable assumptions, these problems have a unique solution as required by Def.4.7.This is not the case for DOF -continuous problems, for which one finds some boundary conditions on ∂S K (because the degrees of freedom of test functions in V K,0 are prescribed on S K ) and another set of conditions on ∂K that are explicitly imposed through the degrees of freedom on K in (31).Well-posedness is not obvious in general, and cannot always be deduced from well-posedness of the DOF -extended counterpart (28).See Examples 6 and 7 for concrete local problems.We now address the well-posedness of DOF -continuous problems in more detail in Sec.4.3.5.

Well-posedness of DOF-continuous numerical correctors
We have seen in Examples 6 and 7 that DOF -continuous local problems lead to non-standard boundary conditions.This poses not only a theoretical issue, but also a computational challenge.To complete our study of the general MsFEM framework, we now present a computational strategy to obtain the DOF -continuous numerical correctors, and we use this strategy to discuss the wellposedness of the associated local problems.
In Def.4.7 we assume the well-posedness of DOF -extended problems, and we have seen in Examples 4 and 5 that this is a natural assumption.It is also natural to assume that we can compute DOF -extended numerical correctors numerically.We compute the DOF -continuous numerical correctors from the DOF -extended numerical correctors, by subtracting a linear combination of suitable functions W β from the DOF -extended numerical correctors.The W β must all satisfy the homogeneous equation s ε K W β , w = 0 for all w ∈ V K,0 , in order not to perturb the local problem (29) that is already satisfied by both types of numerical correctors.We shall use the functions W 0 := 1 + χ ε,0,e S K and W β := x β − x β c,K + χ ε,β,e S K for 1 ≤ β ≤ d, where χ ε,β,e S K is defined in Def.4.8.The precise strategy is as follows.
Fix 0 ≤ α ≤ d.We look for coefficients c α 0 , . . ., c α d such that χ ε,α,c where we recall that χ ε,α,c S K is defined by Def.4.10.Note that both sides of the equation clearly satisfy (29).The desired equality thus holds if and only if Γ Since the DOF operators are linear, this leads to the linear system Invertibility of the matrix M is thus a sufficient condition for the existence of all DOF -continuous numerical correctors, and the resolution of the linear system (32) for each α (where all DOFextended numerical correctors are replaced by their numerical approximation) allows to compute the DOF -continuous numerical correctors numerically.Before studying the invertibility of the matrix M in a few special cases, let us consider the matrix composed of the degrees of freedom on S K , i.e., the matrix By definition of the functions χ ε,β,c S K , we have Γ(S K , W β ) = Γ(S K , x β − x β c,K ) for 1 ≤ β ≤ d, and Γ(S K , W 0 ) = Γ(S K , 1).Note that the constant function together with the coordinate functions ) is a bijection, the vectors Γ(S K , W 0 ), . . ., Γ(S K , W d ) are linearly independent.Hence the matrix M is invertible.One may hope that the linear independence of the vectors Γ(S K , W β ) is preserved for the degrees of freedom on the interior boundary ∂K instead of ∂S K , yielding invertibility of M. We found this to hold for all numerical tests that we performed, involving both the MsFEM-lin and the MsFEM-CR.
We can prove invertibility of M in a few special cases.When s ε K is the sampling form that was used in Example 4 (corresponding to a diffusion problem; we will consider this case until the end of this section) and if A ε is constant, all numerical correctors vanish on S K and the foregoing argument for the matrix M shows invertibility of M.
In the periodic setting (see Sec. 2.7), even though A ε itself is not constant, its homogenized limit A * is.In this case, the χ ε,β,e S K converge to zero weakly in H 1 (S K ).(We show this in Lemma 13 in the absence of oversampling, but the argument can be generalized to DOF -extended oversampling.)Now consider the MsFEM-CR.The weak convergence of the χ ε,β,e S K in H 1 (S K ) ensures weak convergence on each face of K in the H 1/2 -norm by continuity of the trace operator.Since the embedding of H 1/2 (∂K) in L 2 (∂K) is compact, the χ ε,β,e S K converge to 0 strongly in L 2 on each face of K. Consequently, the degrees of freedom Γ(K, χ ε,β,e S K ) (the averages on the faces of K) converge to zero as ε → 0. Thus, Γ(K, W 0 ) → Γ(K, 1) and Γ(K, W β ) → Γ(K, x β − x β c,K ) as ε → 0 for all 1 ≤ β ≤ d and, by the above argument for the matrix M, the matrix M is invertible in this limit.By continuity of the determinant function, the matrix M is invertible when ε is small enough, and the DOF -continuous basis functions exist in this regime.
The study of the DOF -continuous numerical correctors for the MsFEM-lin is more delicate, since pointwise operations are involved in evaluating the degrees of freedom, which are ill-defined on H 1 (S K ).One can invoke the De Giorgi-Nash result, which can be found e.g. in [25,Theorem 8.22], to see that the multiscale basis functions, obtained from the numerical correctors in Def.4.11 below, are in fact continuous for any bounded diffusion tensor.(See Example 8 for a definition of the multiscale basis functions for the MsFEM-lin independent of the numerical correctors.)Pointwise evaluation is then justified.It would therefore be convenient to study the DOF -continuous basis functions directly, without the intermediate step of the numerical correctors.We omit the details here.

The multiscale basis functions
We can now define the multiscale basis functions for the approximation of the abstract problem (26) in terms of the numerical correctors.We recall that in Sec.2.2, the numerical correctors were derived from the definition of the basis functions.We give an equivalent definition of the multiscale basis functions, independent of the numerical correctors, in Lemmas 3 and 4. Recall that φ P 1 1 , . . ., φ P 1 N is a basis of the space V H (see Def. 4.3).We can suppose that the first N 0 basis functions form a basis of V H,0 .The following definition is the generalization of (10) to the general MsFEM framework.Definition 4.11.For each i = 1, . . ., N , the multiscale basis function φ ε i is defined by where • = e corresponds to DOF -extended multiscale basis functions and • = c corresponds to DOF -continuous multiscale basis functions.
The DOF -extended multiscale basis functions satisfy the variational formulation in the following lemma on the oversampling patches S K .Lemma 3. Let K be any mesh element and let 1 ≤ i ≤ N .Consider an MsFEM with DOF-extended basis functions.Define an extension of φ ε i from K to S K by where φ P 1 i K denotes the affine extension of φ P 1 i K to S K , and χ ε,α,e K is as in Def.4.8.Then In the case of the MsFEM-lin for the diffusion problem (1), problem (35) with S K = K coincides with the definition of the multiscale basis functions in (5); see Example 8.
Proof.Problem (35) has a unique solution in view of Def.4.7.It thus suffices to show that φ ε i satisfies (35).Since the numerical correctors χ ε,α,e K belong to V K,0 for all 1 ≤ α ≤ N , it is clear from (35) Inserting (34) into (35) and applying ( 29) to all χ ε,α,e K , we find, for any test function w ∈ V K,0 , Here we use that s ε K is a bilinear form on V K , that all piecewise affine functions are contained in V K according to Def. 4.7, and the property that ∇φ P 1 i is piecewise affine.Finally, we use ( 19) for ϕ = φ P 1 i to conclude that which establishes the desired variational formulation satisfied by φ ε i .
If the DOF -continuous problems (31) are well-posed, we obtain by the same analysis the following result for DOF -continuous multiscale basis functions.
where χ ε,α,c K is as in Def.4.10.Then Example 8 (MsFEM-lin for diffusion problems).In the setting of Example 4, any DOF -extended multiscale basis function φ ε i for the MsFEM-lin constructed in ( 35) is obtained, in each mesh element K, as the restriction of a function φ ε i , which is the unique solution in H 1 (S K ) to For a DOF -continuous basis function, φ ε i solves the same PDE in S K , is affine on ∂S K , and satisfies φ ε i (x j ) = φ P 1 i (x j ) at all vertices x j of K.
Example 9 (MsFEM-CR for diffusion problems).In the continuation of Example 5, the DOFextended multiscale basis function φ ε i for the MsFEM-CR is the restriction to K of φ ε i , the unique solution in H 1 (S K ) to where the constants c h are uniquely determined by the problem.For DOF -continuous basis functions, the last condition is applied to the faces h ∈ F (K) (and all other conditions remain unchanged).
Our general framework allows two characterizations of the multiscale basis functions, namely ( 33) and ( 35) or (36), as was the case for the MsFEM studied in Sec.2.2 (where φ ε i is given by ( 5) or ( 10)).The essential advantage of ( 33) is that the microscale is fully encoded in the numerical correctors χ ε,α,• K , that can be computed element per element without any global information.In particular, the global index i of the multiscale basis function φ ε i is irrelevant for the computation of the numerical correctors.The expression in ( 33) is therefore the crucial relationship that we will employ to develop non-intrusive MsFEMs within the general framework in Sec. 5.
The second formulation of the multiscale basis functions, as solutions to the local problems ( 35) or (36), provides a more direct interpretation of the multiscale basis functions in terms of the sampling form chosen.It also gives a relation between the degrees of freedom of the P 1 basis functions and the associated multiscale basis function.This is useful in particular for the wellposedness of the MsFEM, that we study in Lemma 5.
Remark 4. Our definition of the multiscale basis functions in (33) is reminiscent of the Variational Multiscale Method, a framework developed in [32,33] to adapt Galerkin approximations on lowdimensional spaces to the presence of multiscale features.In this context, our formulation of the MsFEM also exhibits a link with residual-free bubbles, see e.g.[14,13,33].
Remark 5.The first introduction of the MsFEM in [29] corresponds to the idea of oversampling with DOF -continuous basis functions.Although their existence cannot be established in general, they are computed numerically by taking linear combinations of DOF -extended basis functions (following an analogous strategy to the one we discussed in Sec.4.3.5).The MsFEM with DOFextended basis functions is studied in the works [20,30] dealing with the convergence analysis of the MsFEM-lin with oversampling.Let us also note that the combination of Crouzeix-Raviart MsFEM and oversampling has, to the best of our knowledge, not yet been proposed in the literature.This method, for which the basis functions are given explicitly in Example 9, is a natural by-product of the identification of the abstract MsFEM framework.

The global problem
We can now define the multiscale trial and test spaces, respectively V ε H and V ε H,0 , as follows: Note that we only use V ε H,0 in the present section, because ( 26) is posed with homogeneous Dirichlet conditions, but that the larger space V ε H is used for more general boundary conditions (see Sec. 5.3).
Applying (33), we have the equivalent characterization in terms of the P 1 space V H , Definition 4.12.Let V H be an underlying P 1 space defined in Def.4.1 with the associated DOF operator Γ from Def. 4.2 and Def.4.5-4.6.Define for each mesh element K ∈ T H an oversampling patch (Def.4.4), a sampling space and sampling form in accordance with Def.4.7.
Let the multiscale basis functions φ ε i be given as in Def.4.11.Then a Multiscale Finite Element Method (MsFEM) for problem (26) is: find In the following lemma, we investigate the well-posedness of the MsFEM.
Lemma 5. Consider an MsFEM without oversampling, or an MsFEM with oversampling using DOF-continuous basis functions (assuming the associated basis functions are well-defined).When a ε satisfies (27), the MsFEM (37) has a unique solution.
Proof.Note that, with DOF -continuous oversampling, but also without oversampling, the multiscale basis functions satisfy (36).In particular, all degrees of freedom of u ε H related to the boundary vanish.Also note that, the dimension of V ε H,0 being finite, it suffices to show that u ε H = 0 is the unique solution to problem (37) with Because of (36), we have Γ(K, u ε H ) = Γ(K, u H ) for all mesh elements K. Since u ε H is piecewise constant and Γ(K, •) is a bijection from P 1 (K) to R d+1 (recall Def.4.2), it follows that u ε H = u H .In particular, the multiscale function u ε H in fact belongs to the underlying P 1 space V H .We remarked immediately below Def.4.1 that, for either of the two spaces V H = V L H or V CR H , the above implies that u ε H is constant throughout Ω.Since the degrees of freedom of u ε H associated to the boundary vanish, we readily deduce that u ε H = 0.
We do not know of the existence of a result on the well-posedness of MsFEMs with oversampling using DOF -extended multiscale basis functions.In [30], the authors establish an inf-sup result for a variant of the MsFEM-lin-OS with P 1 test functions (see also Def. 5.1).This result is obtained for a periodic diffusion coefficient in the limit of sufficiently small ε.

Non-intrusive MsFEM for the general framework
We show in this section how to develop a non-intrusive approach for the general MsFEM framework of Sec. 4. We have seen in Lemma 2 that for a particular MsFEM variant, the non-intrusive Galerkin MsFEM approach coincides with a Petrov-Galerkin MsFEM.This does not hold for all MsFEMs in the general framework.We first develop a non-intrusive MsFEM approach for a Petrov-Galerkin MsFEM in the general framework.The non-intrusive approach for the Petrov-Galerkin MsFEM is equivalent to the Petrov-Galerkin MsFEM itself.In a second step, we introduce a non-intrusive approximation of the Galerkin MsFEM.Before doing so, let us summarize the main steps of Sec.2.4 and 2.5 to obtain a non-intrusive MsFEM approach: (1) the expansion (10) allows to recast the matrix A ε of the linear system for the MsFEM as the matrix A P 1 associated to the P 1 discretization of an effective problem; (2) we approximate the right-hand side F ε of the MsFEM problem by the right-hand side F P 1 of this P 1 discretization; (3) the post-processing step (20) applied to the P 1 approximation of the effective problem yields the MsFEM approximation.

The Petrov-Galerkin MsFEM
We recall that the abstract continuous problem for which we developed the MsFEM in Sec. 4 is given by (26) and that it can be rewritten in terms of the bilinear forms a ε K satisfying (27).Petrov-Galerkin variants of the multiscale finite element method with P 1 test functions were previously studied in [30,28].In our general MsFEM framework, the adaptation of Def.4.12 to a Petrov-Galerkin MsFEM is the following.Definition 5.1.Let V H be an underlying P 1 space defined in Def.4.1 with the associated DOF operator Γ from Def. 4.2 and Def.4.5-4.6.Define for each mesh element K ∈ T H an oversampling patch (Def.4.4), a sampling space and sampling form in accordance with Def.4.7.Let the multiscale basis functions φ ε i be given as in Def.4.11.Then a Petrov-Galerkin Multiscale Finite Element Method (PG-MsFEM) for problem (26) When confusion may arise, we shall refer to the MsFEM defined in Def.4.12 as the Galerkin MsFEM (G-MsFEM).To study well-posedness of the PG-MsFEM, it is most convenient to relate this method to the G-MsFEM.Therefore, we postpone well-posedness of (38) to Lemma 8.
We now execute step (1) of the summary of the non-intrusive MsFEM approach at the beginning of this section.The matrix A ε of the linear system associated to (38) is defined by To find an effective P 1 formulation with the same linear system, we will use the definition (33) of the multiscale basis functions in the general framework, but first we combine it with (19) applied to ϕ = φ P 1 i to rewrite (33) as where Λ ε,0 K , for all 1 ≤ α ≤ d and each K ∈ T H .We recall that • ∈ {e, c} indicates the choice of DOF -extended or DOF -continuous basis functions.Inserting (40) into (39) for φ ε i and ( 19) for ϕ = φ P 1 j yields where we have defined the effective mass M , (adjoint) advection vector B 1 and B 2 , and the effective diffusion tensor A, for all 1 ≤ α, β ≤ d and for each K ∈ T H , as Note that all integrals in (41) can be computed exactly by evaluating the integrand at the centroid.With this quadrature rule, we observe that the term |K| M φ P 1 i φ P 1 j (x c,K ) also equals the numerical approximation of the integral K M φ P 1 i φ P 1 j .The new expression (41) for the matrix of the linear system motivates us to introduce the effective bilinear forms a K defined on and the associated P 1 Galerkin approximation on the space V H,0 : This discrete problem leads to a linear system with the matrix The identity (41) thus implies the following result, which generalizes Lemma 1 to the PG-MsFEM in the general framework.
Lemma 6.The matrices A ε and A P 1 are identical if the integrals in (43) are evaluated at the centroid of each K for the computation of A P 1 .Then the PG-MsFEM (38) coincides with the resolution of the effective problem (44) combined with the post-processing step Note that step (2) of the summary at the beginning of this section is irrelevant for the PG-MsFEM.The computation of the right-hand side in (38) is clearly part of any standard FEM software.
The computational approach described by Lemma 6 naturally fits within the non-intrusive workflow of Algorithm 2. The numerical correctors on line 4 are, of course, replaced by those of Def.4.8 or Def.4.10.Line 6 is replaced by the computation of all effective quantities in (42).The online phase in line 8 amounts to solving the P 1 problem (44), where all integrations to construct the matrix of the linear system are to be performed by evaluation at the centroid.(This is not the case for the construction of the right-hand side, however.)Finally, in the post-processing phase we construct u ε H from u H by virtue of ( 45).Next we generalize the above expansions to design a non-intrusive approximation of the G-MsFEM.

The non-intrusive Galerkin MsFEM
For the G-MsFEM (introduced in Def.4.12), we need to replace the P 1 test space of the PG-MsFEM by the multiscale test space V ε H,0 .The matrix of the linear system associated to (37) is given by Upon inserting (33) for the test function φ ε j , we find, for all 1 ≤ i, j ≤ N 0 , where A ε is the matrix of the Petrov-Galerkin MsFEM, see (39) and ( 41).An effective formulation can again be derived by inserting (40) for the φ ε i .We obtain where the effective mass, (adjoint) advection vectors and diffusion tensor are given by (using those defined in (42)) The above computations lead to the introduction of the effective bilinear form We formulate the following effective variational problem: The associated linear system has coefficients A P 1 ,G j,i = a G φ P 1 i , φ P 1 j .We have the following analogue of Lemma 6, which generalizes Lemma 1 to the G-MsFEM in the general framework.
Lemma 7. The matrices A ε,G and A P 1 ,G are identical if the integrals in (48) are evaluated at the centroid of each K in the computation of A P 1 ,G .
Contrary to the matrices, the right-hand sides of the effective problem (49) and the Galerkin MsFEM (37) are not equal in general.We apply step (2) formulated at the beginning of this section: the right-hand side of the G-MsFEM is approximated by the right-hand side of the effective problem to obtain an approximate, but non-intrusive, MsFEM.The non-intrusive G-MsFEM becomes: This problem is no longer a Galerkin approximation of (1), because different test spaces are used for the bilinear and for the linear form.In view of Lemma 7, the non-intrusive MsFEM can equivalently be formulated as compute u H ∈ V H,0 solution to (49) and compute u ε H from u H by ( 45), provided all integrals in (48) are evaluated at the centroid for the construction of the matrix of the linear system in (49).The latter formulation of the non-intrusive MsFEM immediately suggests how to effectively implement the non-intrusive MsFEM in a non-intrusive way similar to Algorithm 2. For completeness, we provide the algorithm for the non-intrusive G-MsFEM in Algorithm 3.

Algorithm 3 Non-intrusive G-MsFEM for the general framework
1: Let T H be the mesh used by the legacy code, let • ∈ {e, c} be the chosen oversampling variant Compute the effective tensors defined in (47) 7: end for 8: Use the legacy code to construct the matrix A P 1 by evaluating (48) at the centroid of each mesh element and to solve for u H defined by (49) The discussion surrounding Algorithm 2 regarding the advantages for the implementation of this non-intrusive MsFEM approach also applies here.
Let us now comment on the well-posedness of the MsFEMs for the general framework introduced above.We recall that the hypotheses of the general framework without oversampling, or with DOFcontinuous oversampling, provide well-posedness of the G-MsFEM (37) by Lemma 5.In this case, the non-intrusive approximation (50) is also well-posed, for the matrices associated to both MsFEM variants are the same.Regarding the PG-MsFEM (38), we can only establish well-posedness if the associated matrix coincides with the matrix of the corresponding Galerkin MsFEM.This is stated in the following lemma, which generalizes Lemma 2 to the general framework.Lemma 8. Consider a G-MsFEM as defined by Def.4.12 without oversampling and suppose that the sampling form s ε K equals the local bilinear form a ε K .Then the matrix associated to this G-MsFEM coincides with the matrix associated to the corresponding PG-MsFEM of Def.5.1.Consequently, the non-intrusive Galerkin MsFEM (50) coincides with the Petrov-Galerkin MsFEM of Def.5.1 and in particular, The Petrov-Galerkin MsFEM is well-posed.Proof.To prove the lemma, we show that the matrices corresponding to the linear problems defined in (50) and (38) are equal.Using s ε K = a ε,diff K , we have for all 1 ≤ i, j ≤ N 0 , Indeed, the multiscale basis functions satisfy Γ(K, φ ε i ) = Γ(K, φ P 1 i ) for all K, so that φ ε j − φ P 1 j ∈ V K,0 (see (35) with S K = K and recall Def.4.7 for the sampling test space V ε K,0 ), and the variational problem in (35) (with S K = K) shows that the above quantity vanishes.

Further extensions of the non-intrusive MsFEM
We sketch some other FEM settings to which we have applied the above strategy to develop nonintrusive MsFEM approaches.For more details, we refer to [11].
Stabilized finite element formulations.Stabilized finite element formulations add meshdependent terms to a discrete variational formulation (such as (37)) to remove numerical instabilities, for example caused by sharp boundary layers of the exact solution.See [39] for such a variant of the MsFEM and see [15,31,46] for the stabilization of single-scale problems.The expansion (40) can also be inserted in these additional terms to find a non-intrusive implementation of the associated MsFEM.
Petrov-Galerkin formulations.Other test spaces than the P 1 space V H,0 can be considered in Petrov-Galerkin formulations.An example would be to use multiscale test functions that locally solve the adjoint problem rather than the direct problem that is used for the multiscale functions in (5).See e.g.[23].An expansion of the kind (40) can still be used to find a non-intrusive formulation, with a suitably adapted definition of the numerical correctors.
Non-homogeneous Dirichlet conditions.Suppose that a legacy FEM code can provide a solution to an effective problem such as (15) posed on the space V H,0 and complemented with non-homogeneous Dirichlet conditions on ∂Ω.This solution can directly be used to construct a multiscale approximation u ε H ∈ V ε H from (45).The translation of the Dirichlet condition to the MsFEM approximation is as follows: if DOF -continuous oversampling is applied, the function ] j for all degrees of freedom associated to the boundary.Here, [Γ(K, u H )] j is determined by the legacy code.When DOF -extended oversampling is used, the degrees of freedom associated to the boundary are equal to the sum of [Γ(K, u H )] j and a perturbation due to the fact that the degrees of freedom of the numerical correctors do not vanish.
Neumann conditions.To apply Neumann conditions on ∂Ω, one solves a Galerkin approximation of the variational formulation in the space V ε H .The suitable adaptation of (37) can be approximated by a non-intrusive Galerkin MsFEM following the same methodology as above.The effective P 1 approximation that is obtained corresponds to the resolution of an effective PDE with Neumann conditions, for which a legacy code can be used.In the case of the diffusion problem (1), the Neumann boundary condition in the effective problem is imposed on the effective flux n • A∇u H , where A is defined in (12).
Parabolic equations.When a parabolic equation is discretized in time, problems of the form (26) are typically obtained for each time step, but with a right-hand side that depends on the solution of the previous time step.This term belongs to the space V ε H , so it varies on the microscale and cannot be integrated numerically by the legacy code that operates on the coarse mesh.The non-intrusive strategy of the foregoing sections cannot be applied directly to find a non-intrusive MsFEM.In the vein of our non-intrusive approach, one could introduce an additional approximation upon replacing the multiscale solution of the previous time step by its underlying P 1 in (40).The effect of this approximation is beyond the scope of the present work.

Intrusiveness of other multiscale methods
Some work on the formulation of effective P 1 problems in multiscale methods, and the related question of non-intrusive approaches, can be found in the literature.We discuss here the case of the HMM and the LOD method.
First, the HMM is less intrusive than the original MsFEM, because its main objective is to approximate u ε on the coarse scale.The HMM directly proposes to solve a P 1 problem on the coarse scale, where effective coefficients of the P 1 problem are defined in terms of the solutions to local problems.This workflow corresponds to our non-intrusive MsFEM approach, and when the local problems of the HMM coincide with the computation of the numerical correctors introduced in this work, the HMM and the MsFEM for the pure diffusion problem are identical.For more general problems, there is an important difference between the two methods.In the MsFEM, the form of the effective equation and the definition of the effective coefficients follows directly from the choice of basis functions, and thus from the choice of local problems.For the HMM, the local problems and the effective equation are formulated independently, and the link between the two is only justified heuristically, drawing inspiration from homogenization theory.
The LOD method aims at approximating u ε at the coarse and the microscale by the use of multiscale basis functions, like the MsFEM.It is shown in [24] that a Petrov-Galerkin LOD method (also see [21]) can, with some additional approximations, be recast as the P 1 discretization of an appropriate coarse-scale problem.This opens the way to non-intrusive implementations in the spirit of the present article.The LOD method and the MsFEM notably differ in the fact that the LOD basis functions are defined on a patch around the vertices of the mesh that should generally be taken larger than the support of the associated P 1 functions.In contrast, the MsFEM uses fully localized basis functions (even though they may have been computed using oversampling patches), each of which has the same support as the corresponding P 1 basis functions.

Comparison of the classical and non-intrusive MsFEM for diffusion problems
We study in this section a particular setting within the general MsFEM framework, namely that of MsFEMs for diffusion problems.We set in this section a ε K = a ε,diff K defined in Example 2, and we choose the sampling form s ε K = a ε,diff K .

The general framework for diffusion problems
For the convenience of the reader, we first give an explicit description of the simplifications of the general framework in the diffusion setting.In Def.4.8 and 4.10 for the numerical correctors, Equation ( 29) reduces to for all w ∈ V K,0 (where V K,0 is the sampling test space for either the MsFEM-lin or the MsFEM-CR; see Examples 4 and 5) when 1 ≤ α ≤ d, whereas χ ε,0,• S K = 0. (The notation χ ε,α K will be used in the absence of oversampling, see Rem. 3.) This means that Λ ε,0 K = 1 in (40).Consequently, regarding the formulation of the effective P 1 problem, only the effective diffusion coefficient does not vanish in (42) and (47).Its definition in (47) is identical to the formula in (12) for the applicable choice of the numerical correctors.
The definition of the multiscale basis functions by (33) reduces to (10) (again upon replacing the numerical correctors χ ε,α K by the relevant ones for the MsFEM under consideration).Hence, we can associate a multiscale counterpart in V ε H to any v H ∈ V H , given by The non-intrusive MsFEM (50) becomes Lemma 8 now amounts to the following.
Lemma 9. Let 'MsFEM' refer to the MsFEM-lin or the MsFEM-CR without oversampling.The non-intrusive Galerkin MsFEM (54) coincides with the following Petrov-Galerkin MsFEM: We will specify for all results in this section to which specific MsFEMs they apply among the MsFEM-lin and the MsFEM-CR, with or without oversampling.Lemmas 10, 11 and 16 are generalizations of results of [12], where the MsFEM-lin without oversampling is treated.

Convergence results
We estimate here the difference between the solutions to the (intrusive) Galerkin approximation (37) and the non-intrusive MsFEM (54), which coincides with the Petrov-Galerkin MsFEM (55).We first show coercivity of the effective diffusion tensor A.
Lemma 10.Consider the MsFEM-lin or the MsFEM-CR without oversampling, or the MsFEM-CR with DOF-continuous oversampling.The effective tensor A defined by (12) with the appropriate numerical correctors satisfies Here, m is the same coercivity constant as in (2).
Proof.Let ξ = (ξ 1 , . . ., ξ d ) ∈ R d , and let K be any simplex of the mesh T H .We have Using an integration by parts, we see that where n is the unit outward normal vector on ∂K.In the case of the MsFEM-lin, the function χ ξ vanishes on ∂K.In the case of the MsFEM-CR with DOF -continuous oversampling, or without oversampling, the function χ ξ has average zero on each face of K. Since the factor n • ξ is constant on each face, the integral again vanishes.In conclusion, we have We thus obtain the inequality ξ • A K ξ ≥ m|ξ| 2 .Since K ∈ T H is arbitrary here, this shows coercivity of A and completes the proof.
Coercivity of the effective tensor A implies coercivity of the bilinear form a diff on H 1 0 (Ω).By an application of the Lax-Milgram Theorem, we conclude that the (continuous) effective problem ( 14) is well-posed for the MsFEM-lin and the MsFEM-CR without oversampling, and for the MsFEM-CR with DOF -continuous oversampling.Remark 6.The proof of the above lemma does not extend to the MsFEM-lin with oversampling, because there is no global information about χ ε,α K on the faces of K.The following lemma provides a variational characterization of the bijection (53).
In addition, we have, with the constants m and M from (2), the estimate Proof.Let v H ∈ V H be the unique element of V H such that v ε H and v H satisfy (53).Take any w H ∈ V H . Using that ∇v H and ∇w H are piecewise constant, we compute For the MsFEM without oversampling, the numerical correctors belong to the sampling test space V K,0 .We can thus use (52) to obtain Using the definitions of A in (12) and the of a diff in (15) (we recall that these expressions hold true here upon replacing the numerical correctors by those under consideration), we conclude that It follows that v H satisfies (56).In addition, in view of the coercivity of A established in Lemma 10 and by the Lax-Milgram Theorem, problem (56) uniquely characterizes v H .
The estimate on v H follows by testing the characterization (56) against w H = v H .This yields, for any The first inequality follows from coercivity of A and the second inequality from the upper bound on A ε in (2) and the Cauchy-Schwarz inequality.We thus find that ∇v . The desired result is obtained upon squaring the inequality and summing over all K ∈ T H .
For the remainder of this section, we will consider MsFEMs without oversampling.Let u ε,G H denote the solution to the MsFEM approximation (37) (we use the superscript G to stress that this is a Galerkin approximation) and let u ε,PG H denote the solution to the non-intrusive MsFEM (54) (which is equivalent to the Petrov-Galerkin MsFEM (55), since we do not apply the oversampling technique).
We first study the error u ε,G H −u ε,PG H when ε → 0. In this case, we do not need a rate of convergence in H and we shall relax the condition f ∈ L 2 (Ω) to the condition f ∈ H −1 (Ω).Then the linear form F in (3) has to be redefined.Given f ∈ H −1 (Ω), there exist f 0 , f 1 , . . ., f d ∈ L 2 (Ω) such that which is in fact well-defined for any v ∈ H 1 (T H ) and thus in particular on V H , the underlying affine space for the MsFEM, and the multiscale space V ε H .We consider in Lemma 12 a sequence of diffusion tensors A ε that H-converges to a constant diffusion tensor.This means that u ε converges weakly in H 1 (Ω) as ε → 0 towards a function u ⋆ ∈ H 1 0 (Ω), solution to the homogenized problem (24), and A ε ∇u ⋆ ⇀ A ⋆ ∇u ⋆ weakly in L 2 (Ω).
Lemma 12. Consider the MsFEM-lin or the MsFEM-CR without oversampling.Suppose that (A ε ) ε>0 is a sequence of matrices satisfying (2) that H-converges to a constant matrix.Let f ∈ H → 0 as ε → 0.
Remark 7. A rate of convergence can be obtained under some additional structural assumptions on A ε ; see Lemma 18.
We need a few auxiliary results to establish Lemma 12.The first result below concerns the convergence of the numerical correctors as ε → 0.
Lemma 13.Suppose that A ε H-converges to a constant homogenized tensor A ⋆ .Consider the MsFEM-lin or the MsFEM-CR without oversampling.Then, for all K ∈ T H and all 1 ≤ α ≤ d, we have χ ε,α K ⇀ 0 weakly in H 1 (K) as ε → 0. Proof.We introduce for each α = 1, . . ., d the function τ ε,α = x α + χ ε,α K .Then (52) implies the PDE −div(A ε ∇τ ε,α ) = 0 in K.The homogenized limit of this problem is −div(A ⋆ ∇τ ⋆,α ) = 0 in K.For the MsFEM-lin, the boundary conditions of the local problems impose τ ⋆,α = x α on ∂K.The boundary conditions associated to the MsFEM-CR are a constant flux n • A ⋆ ∇τ α,⋆ on each face of K and h τ ⋆,α = h x α for all faces h of K.In both cases, the homogenized equation has a unique solution, which is easily seen to be τ ⋆,α = x α , because A ⋆ is constant.Therefore, τ ⋆,α ⇀ x α weakly in H 1 (K).Subtracting again the function x α , we deduce the desired convergence.
We will also use the following result, which is a straightforward generalization of the extended Poincaré inequality in [22,Lemma 3.31].Lemma 14.Let W be the subspace of H 1 (T H ) defined by There exists a constant C > 0 depending only on Ω but not on H such that Note that the multiscale space V ε H,0 is contained in W for both the MsFEM-lin and the MsFEM-CR without oversampling.Finally, we provide a number of useful bounds for the difference between u ε,G H and u where Proof.Since the numerical approximations u ε,G H and u Applying Lemma 11 to v ε H = e ε H ,we immediately obtain (59).Now recall that the numerical correctors are defined by (52).Using the fact that ∇e P 1 H is piecewise constant, this implies that e osc H satisfies the following variational problem in each K ∈ T H : Without oversampling, it holds χ ε,α K ∈ V K,0 for each 1 ≤ α ≤ d, so e osc H can be used as a test function here.With the bounds in (2), implying continuity and coercivity of a ε,diff K , we obtain (58).Next using (61), we can write We deduce from (51) that a ε,diff u ε,PG H , e osc H = 0. Since e ε H can be used as a test function in the discrete problem (37) where C is the Poincaré constant from Lemma 14.Now applying (58) and (59) on the right, and using coercivity of a ε,diff on the left, we find from which we deduce (60).
We next study the convergence of u ε,G H − u ε,PG H as H → 0. To this end, we return to the original hypotheses of Sec. 2, i.e., f ∈ L 2 (Ω).Note that for the next result, the additional convergence hypothesis of Lemma 12 for A ε is not needed.Lemma 16.Consider the MsFEM-lin or the MsFEM-CR without oversampling.Assume that f ∈ L 2 (Ω).Then there exists a constant C independent of ε, H and f such that To prove this lemma, we will use some Poincaré-Friedrichs inequalities, for which we refer e.g. to [38,Lemma 4.3], [22,Lemma B.66 For the MsFEM-lin, it holds that χ ε,α K = 0 on ∂K for all mesh elements K and all 1 ≤ α ≤ d, and it follows that e osc H = 0 on the boundaries of all mesh elements.In the case of the MsFEM-CR, it holds that h χ ε,α K = 0 for all faces h of the mesh and all α.Since ∂ α e P 1 H is constant on each mesh element K, we also have One more time using the lower bound in (2), we find The proof is concluded by application of Lemma 14 to e ε H .

Convergence results in the periodic setting
We now study the MsFEM-lin applied to the periodic setting introduced in Sec.2.7 in some more detail.To the best of our knowledge, all convergence results known for the MsFEM are obtained in this periodic setting (see e.g.[20,30,19,3,28,38,37,41,40]).The analysis in these works relies on the explicit description of the microstructure that we summarized in Sec.2.7.In particular, recall the existence of a homogenized diffusion coefficient given by ( 23) and the first-order twoscale expansion (25).We emphasize, however, that the application of the MsFEM does not require the periodic setting, nor does it even suppose the PDE under consideration to be embedded in a sequence of PDEs for a family of parameters ε that tend 0. Applying the MsFEM to a sequence of matrices A ε = A per (•/ε), we obtain a sequence of effective tensors A(ε).Each A(ε) is defined by (12) for a fixed value of ε.We have the following convergence result.
Proof.We fix a mesh element K ∈ T H . First observe that A(ε) and A ⋆ satisfy for each 1 ≤ α, β ≤ d, in view of the variational formulations satisfied by χ ε,α K (solution to the PDE (11)) and w α (solution to the PDE (22)).Now let τ ε,α = x α + χ ε,α K .In view of Lemma 13, τ ε,α ⇀ τ ⋆,α as ε → 0 weakly in H 1 (K).Writing the two-scale expansion (25) of τ ε,α , we thus have, when ε is small, and the difference tends to zero in H 1 (K) as ε → 0. Inserting this convergence in (64), we deduce that lim The convergence to the mean on the unit cube in the last equality follows from the Q-periodicity of the function e β • A per (e α + ∇w α ).
The following lemma studies the convergence of u ε,G H − u ε,PG H towards 0 as ε → 0 for the MsFEMlin.As was stated in Rem. 7, thanks to the periodic setting, we now obtain a rate for the convergence stated in Lemma 12.
Lemma 18.Let f ∈ L 2 (Ω).Suppose that the family of meshes (T H ) H>0 is quasi-uniform.Consider the MsFEM-lin without oversampling.For A ε = A per (•/ε) sufficiently regular, we have where the constant C depends on the dimension d and the constants m, M in (2), but not on ε, H or f .
Proof.Let e ε H = u ε,G H −u ε,PG H . Lemma 15 applies, so we can use (57) and a Cauchy-Schwarz inequality to find a ε,diff (e ε H , e ε H ) ≤ f L 2 (Ω) e osc H L 2 (Ω) ≤ f L 2 (Ω) Next we seek a bound on χ ε,α K in L 2 (K).Using (11) and ( 22), we have div Since χ ε,α K vanishes on ∂K (recall that we consider the MsFEM-lin without oversampling), the maximum principle [25,Theorem 8.1] yields When A per is sufficiently regular, the corrector functions w α are uniformly bounded.Then the mesh regularity provides a constant C such that for each K ∈ T H and each 1 ≤ α ≤ d, we have Since all χ ε,α K have disjoint supports, we can use the latter estimate to bound . (66) The last inequality relies on the quasi-uniformity of the mesh.We insert (66) combined with (59) into (65) to find a ε,diff (e ε H , e ε H ) ≤ Cε f L 2 (Ω) ∇e ε H L 2 (Ω) .
Applying the coercivity property in (2) on the left-hand side, we obtain the desired result.
The classical error estimate for the Galerkin MsFEM approach (6) is obtained in the periodic setting and under some regularity assumption on A per and on the homogenized limit u ⋆ .The bound obtained in [19,Theorem 6.5] reads ≤ C H + ε + ε/H , for some C independent of ε and H. Lemma 18 shows that the same estimate holds true for u ε,PG H , the Petrov-Galerkin MsFEM approximation, under the correct regularity assumptions.We note that the bound for u ε,PG H can also be inferred from Lemma 16.However, since the MsFEM is applied in the regime where ε < H, the result of Lemma 18 is more precise, thanks to the extra structural assumptions made on the diffusion tensor A ε .
A reference solution u ε h is computed on a uniform 1024×1024 mesh T h by means of a standard P 1 finite element method using FreeFEM++ [27].The FreeFEM++ scripts to perform all different MsFEMs can be found in the GitHub repository at https://github.com/RBiezemans/MsFEM-in-FreeFEM.
We compare the reference solution u ε h to MsFEM solutions obtained on a coarse mesh T H for varying H.The mesh T H is a uniform 1/H × 1/H triangulation of Ω.We test the MsFEM-lin and the MsFEM-CR using the sampling operator s ε K = a ε,diff K .All oversampling methods in this section use a homothety ratio of 3 for the construction of the oversampling patches in Def.4.4, and DOF -continuous basis functions, which ensure certain continuity properties on the boundary of the mesh elements.A precise definition of the associated basis functions can be found in Examples 8 and 9.The mesh T h is a refinement of T H for all values of H. Therefore, for each K ∈ T H , we use the corresponding submesh of T h (consisting of all triangles included in K) for the numerical approximation of the numerical correctors in (29) by P 1 Lagrange finite elements.
We first compare the approximations u ε,G H and u ε,G-ni H for varying H in Fig. 3. Without oversampling (OS), the approximation u ε,G-ni H equals u ε,PG H due to Lemma 8. We also report the error committed by the G-MsFEM.We observe that, without oversampling, the difference u ε,G H − u ε,G-ni H is much smaller than this error.As a result, the errors obtained with the G-MsFEM and its nonintrusive approximation are of the same size.Indeed, the error of the non-intrusive G-MsFEM-lin deviates from the error of the G-MsFEM-lin by at most 0.05% for all tests that we report here.For the MsFEM-CR, this is at most 1.2%.In both cases, the two MsFEM variants thus have practically the same accuracy.This is in agreement with the theoretical result of Lemma 16.The estimates obtained in Sec.6 do not apply to MsFEMs with oversampling.From Fig. 3, we can see that the difference u ε,G H − u ε,G-ni H is still small with respect to the error committed by the G-MsFEM when oversampling is applied.The approximation errors for the non-intrusive G-MsFEMs with oversampling differ by at most 1.3% from the error of the G-MsFEM.Let us also point out the qualitative and quantitative similarities between the performance of the MsFEM for the periodic and the nonperiodic diffusion coefficient.Although the homogenized coefficient of A ε,np has a more complicated e. in Ω, and ∀ ξ, η ∈ R d , |η • A ε (x)ξ| ≤ M |ξ| |η| a.e. in Ω,

Lemma 4 .
Let K be any mesh element and let 1 ≤ i ≤ N .Assume that any DOF-continuous local problem (31) is well-posed.Consider an MsFEM with DOF-continuous oversampling.Define an extension of φ ε i from K to S K by

H d/ 2 ∂ α e P 1 H K 2 ≤ 2 K∈T H d α=1 ∂ α e P 1 H 2 L 2 = Cε 2 ∇e P 1 H 2 L 2
Cε (K) (Ω) ε,PG H . Lemma 15.Let f ∈ H −1 (Ω) and consider the MsFEM-lin or the MsFEM-CR without oversampling.Let e ε H = u ε,G H − u ε,PG H .There exists a unique e P 1 H ∈ V H and a linear combination of the numerical correctors, that we denote by e osc H , such that e ε H = e P 1 H + e osc H , and it holds, with the constants m, M from (2) and the constant C from Lemma 14, and e P 1 H in (55), we have a ε,diff u ε,G H , e ε H − a ε,diff u ε,PG H , e P 1 H = F (e osc H ), which shows (57).It follows that a ]. and recall the results of Lemma 15.We have e ε H = e P 1 H + e osc H (see(61)), and (57) provides, for f ∈ L 2 (Ω), the equality a ε,diff (e ε H , e ε H ) = (f, e osc H ) L 2 (Ω) .Hence, by the Cauchy-Schwarz inequality, a ε,diff (e ε H , e ε H ) ≤ f L 2 (Ω) e osc H L 2 (Ω) .