Multiscale Finite Element Methods (MsFEMs) are now well-established finite element type approaches dedicated to multiscale problems. They first 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 differential equations. Non-intrusive MsFEM approaches are developed within the full generality of this framework, which may moreover be beneficial to steering software development and improving the theoretical understanding and analysis of MsFEMs.

§ INTRODUCTION

§ THE (INTRUSIVE) MULTISCALE FINITE ELEMENT METHOD

§.§ Discrete variational formulation

Let $d \geq 1$ denote the space dimension of interest and let $\Omega \subset \bbR^d$ be a bounded polytope.

Convexity can be assumed for elliptic regularity results to hold, for which we refer to <cit.>. 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. <ref> onwards, we will also consider more general problems, and we will mention other types of boundary conditions in Sec. <ref>.
More precisely, we focus here on the boundary value problem

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\varepsilon \nabla u^\varepsilon)}$
&in $\Omega$,
\\
&on $\partial \Omega$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-pde}
\end{equation}

where the diffusion tensor $A^\varepsilon \in L^\infty(\Omega,\ \bbR^{d \times d})$ satisfies the uniform bounds

\begin{gather}
\begin{aligned}
&\forall \, \xi\in\bbR^d,
\quad
m |\xi|^2 \leq \xi \cdot A^\varepsilon(x) \, \xi
\quad \text{a.e.~in } \Omega, \\
\text{and}
\qquad
&\forall \, \xi,\eta \in \bbR^d,
\quad
| \eta \cdot A^\varepsilon(x) \xi | \leq M \, |\xi| \, |\eta|
\quad \text{a.e.~in } \Omega,
\end{aligned}
\label{ass:bounds}
\end{gather}

for some $M \geq m > 0$ independent of $\varepsilon$. The right-hand side $f$ does not vary on the microscopic scale $\varepsilon$. We denote the diffusion tensor with a superscript $\varepsilon$ to keep in mind that $A^\varepsilon$ might be highly oscillatory on a typical length scale of size $\varepsilon$ much smaller than the diameter of $\Omega$ (assumed to be of order 1).
No further structural assumptions on $A^\varepsilon$ are made. In particular, $A^\varepsilon$ need not be the rescaling of a fixed periodic matrix of the form $A^\varepsilon(x) = A(x/\varepsilon)$. We will specialize to this periodic setting in Sec. <ref> 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\in L^2(\Omega)$ (rather than $f \in H^{-1}(\Omega)$, for which the problem (<ref>) 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 \in L^2(\Omega)$ cannot be relaxed. In this case, this is also explicitly stated.

Problem (<ref>) admits a unique solution in the space $H^1_0(\Omega)$. This solution is also characterized by the variational formulation

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{t?s}
\IEEEstrut
Find $u^\varepsilon \in H^1_0(\Omega)$ such that
\\
$a^{\varepsilon,\dif}(u^\varepsilon,\, v) = F(v)$
&for all $v \in H^1_0(\Omega)$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-vf}
\end{equation}

where the bilinear form $a^{\varepsilon,\dif}$ and the linear form $F$ are defined, for any $u,\, v \in H^1_0(\Omega)$, by

\begin{equation}
a^{\varepsilon, \dif}(u,\, v) =
\int_\Omega \nabla v \cdot A^\varepsilon \nabla u,
\quad
F(v) = \int_\Omega f v.
\label{eq:diffusion-bilin-form}
\end{equation}

The coercivity hypothesis in (<ref>) ensures that the bilinear form $a^{\varepsilon,\dif}$ is coercive on the space $H^1_0(\Omega)$. Then the Lax-Milgram Theorem <cit.> shows that (<ref>) is indeed well-posed.

The numerical approximation of (<ref>) with a finite element method starts by the introduction of a mesh $\mathcal{T}_H$ for $\Omega$. The subscript $H$ denotes the typical size of the mesh elements.
We assume $\mathcal{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 <cit.> 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 (<ref>) is obtained by restricting the equivalent formulation (<ref>) to a finite-dimensional subspace of $H^1_0(\Omega)$, typically consisting of functions that are piecewise polynomial on the mesh $\mathcal{T}_H$. We suppose that we are in the regime where $H$ is larger than or comparable to the microscale $\varepsilon$ and cannot be taken smaller due to computational limitations. In this case, it is well known that a Galerkin approximation of (<ref>) on, say, the standard (conforming) Lagrange $\Pone$ space on $\mathcal{T}_H$ provides only a poor, not to say an incorrect approximation of $u^\varepsilon$. See <cit.>, for instance, for an explicit example where the $\Pone$ approximation on a coarse mesh fails. At the same time, the use of a finite element method on a fine mesh of size $H \ll \varepsilon$ 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) <cit.>.

§.§ A simple multiscale finite element method

The MsFEM is a Galerkin approximation of (<ref>) 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^\varepsilon$ than a standard $\Pone$-approximation when $\varepsilon$ is smaller than $H$, and especially when $\varepsilon$ becomes asymptotically small. To begin with, we introduce here the simplest variant of the MsFEM, which originally appeared in <cit.>, before moving on to other MsFEM variants.

Let $x_1,\dots,\, x_{N_0}$ be an enumeration of the interior vertices of the mesh $\mathcal{T}_H$, i.e., the vertices that do not lie on $\partial \Omega$.
We denote by $\phiPone{i}$ the unique piecewise $\Pone$ function such that $\phiPone{i}(x_j) = \delta_{i,j}$ for all $1 \leq j \leq N_0$. (These are the basis functions for the standard $\Pone$ Lagrange finite element.) We define the multiscale basis functions $\phiEps{i}$ (for $1 \leq i \leq N_0$) by

\begin{equation}
\forall \ K \in \mathcal{T}_H,
\quad
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\varepsilon \nabla \phiEps{i})}$
&in $K$,
\\
&on $\partial K$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:MsFEM-basis}
\end{equation}

All these problems, on each mesh element $K$, are again well-posed by coercivity of $A^\varepsilon$ and the Lax-Milgram Theorem. The functions $\phiEps{i}$ so defined belong to the global space $H^1_0(\Omega)$ because the local boundary conditions on $\partial K$ imply continuity across all mesh elements $K$. It is also immediately seen that $\phiEps{i}$ is supported by exactly the same mesh elements as $\phiPone{i}$.

On each mesh element $K$, problem (<ref>) defines at most  non-trivial basis functions. Let $i_1,\dots,i_{d+1}$ be the indices of the vertices of $K$. It is easily inferred from (<ref>) that
\left. \phiEps{i_{d+1}} \right\vert_K
1 - \sum_{j=1}^d \left. \phiEps{i_j}\right\vert_K .
Thus, one only has to compute $d$ basis functions by the resolution of a PDE on $K$.

The multiscale approximation space is defined as $V_{H,0}^\varepsilon = \operatorname{span}\{ \phiEps{i} \mid 1 \leq i \leq N_0 \}$. This is a finite-dimensional space of the same dimension as the one used for a $\Pone$ Lagrange finite element approximation on the mesh $\mathcal{T}_H$. The MsFEM consists in computing the approximation $u^\varepsilon_H \in V_{H,0}^\varepsilon$ defined by the problem

\begin{equation}
\forall \, v_H^\varepsilon \in V_{H,0}^\varepsilon,
\qquad
a^{\varepsilon,\dif} \left(u^\varepsilon_H,\, v_H^\varepsilon \right) = F \left(v_H^\varepsilon \right).
\label{eq:diffusion-MsFEM}
\end{equation}

Since $V_{H,0}^\varepsilon$ is a subspace of $H^1_0(\Omega)$, the bilinear form $a^{\varepsilon,\dif}$ is coercive on $V_{H,0}^\varepsilon$ and the discrete problem (<ref>) is again well-posed by virtue of the Lax-Milgram Theorem.

The computation of the multiscale basis functions $\phiEps{i}$ is called the offline stage of the MsFEM and only has to be carried out once if (<ref>) has to be solved multiple times for different right-hand sides. Also note that all problems (<ref>) are independent of each other, and can thus be solved in parallel. In practice, the $\phiEps{i}$ are approximated numerically on a fine mesh of $K\in\mathcal{T}_H$ of mesh size $h \leq \varepsilon$ that resolves the oscillations of $A^\varepsilon$. 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 (<ref>) is called the online stage. The computational cost for this problem is the same as for a standard $\Pone$ approximation on the same mesh.
A further discussion of the practical implementation of the MsFEM is provided in Sec. <ref>. This discussion partially reproduces some elements of <cit.>. We include it here to clarify and motivate the developments in the sequel.

§.§ Intrusive workflow

The practical resolution of the global problem (<ref>) consists in the construction and resolution of the following linear system:

\begin{equation}
\mathds{A}^\varepsilon U^\varepsilon = \mathds{F}^\varepsilon,
\label{eq:linear-system-msfem}
\end{equation}

\begin{equation}
\forall \ 1 \leq i,\,j \leq N_0,
\quad
\mathds{A}^\varepsilon_{j,i} = a^{\varepsilon,\dif} \left( \phiEps{i},\, \phiEps{j} \right),
\quad
\mathds{F}^\varepsilon_j = F \left( \phiEps{j} \right),
\label{eq:linear-system-msbasis}
\end{equation}

where we recall that $N_0$ denotes the number of interior vertices of $\mathcal{T}_H$.
The MsFEM approximation $u^\varepsilon_H$ is given by
\sum_{i=1}^{N_0} U_i^\varepsilon \phiEps{i}.
The MsFEM can then be written (as it is traditionally presented) as in Algorithm <ref>. We use the notation $\displaystyle a^{\varepsilon,\dif}_K(u,v) = \int_K \nabla v \cdot A^\varepsilon \nabla u$ for all $u,v \in H^1(K)$ and we write $\displaystyle F_K(v) = \int_K fv$ for any $v \in L^2(K)$.

MsFEM approach for problem (<ref>) (see comments in the text)

[1]
Construct a mesh $\mathcal{T}_H$ of $\Omega$, denote $N_0$ the number of internal vertices and $\mathcal{N}(n,K)$ the global index of the vertex of $K \in \mathcal{T}_H$ that has local index $1 \leq n \leq d+1$ in $K$

Set $\mathds{A}^\varepsilon \coloneqq 0$ and $\mathds{F}^\varepsilon \coloneqq 0$

$K \in \mathcal{T}_H$
$1 \leq n \leq d+1$
Set $i \coloneqq \mathcal{N}(n,K)$
Solve for $\left. \phiEps{i} \right\vert_K$ in (<ref>)

$1 \leq l \leq d+1$
Set $j \coloneqq \mathcal{N}(l,K)$
$1 \leq n \leq d+1$
Set $i \coloneqq \mathcal{N}(n,K)$ and $\mathds{A}^\varepsilon_{j,i} \pluseq a_K^{\varepsilon,\dif}(\phiEps{i},\, \phiEps{j})$

Set $\mathds{F}^\varepsilon_{j} \pluseq F_K(\phiEps{j})$

Solve the linear system $\mathds{A}^\varepsilon U^\varepsilon = \mathds{F}^\varepsilon$

Obtain the MsFEM approximation $u^\varepsilon_H = \sum\limits_{i=1}^{N_0} U^\varepsilon_i \phiEps{i}$

Lines <ref>-<ref> of Algorithm <ref> (resp. <ref>-<ref>) constitute the offline (resp. online) stage of the MsFEM. Note that the computation of the stiffness matrix in line <ref> 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 $\phiEps{i}$ in line <ref>, one should apply Rem. <ref>. Only the online stage is to be repeated when problem (<ref>) is to be solved multiple times for various right-hand sides $f$.

Implementing Algorithm <ref> 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 <ref> of Algorithm <ref>), it is by no means obvious to adapt existing finite element codes based on generic approximation spaces (such as the spaces $V_H^L$ and $V_H^{CR}$ in Def. <ref>) to a different, problem-dependent choice of space such as $V_H^\varepsilon$. No analytic expressions for the basis functions $\phiEps{i}$ are available (and thus a fine mesh should be used to approximate them), the computation of $a_K^\varepsilon(\phiEps{i},\phiEps{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 ($\mathcal{N}(n,K)$ in Algorithm <ref>), the assembly of the global stiffness matrix $\{ \mathds{A}^\varepsilon_{j,i} \}_{1 \leq i,j \leq 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 (<ref>) by a $\Pone$ approximation on $\mathcal{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 (<ref>) requires the computation of the quantities

\begin{equation}
\mathds{A}^\varepsilon_{j,i}
\sum_{K\in\mathcal{T}_H} \int_K \nabla \phiEps{j} \cdot A^\varepsilon \nabla \phiEps{i},
\label{eq:stifness-msfem-msbasis}
\end{equation}

for all $1 \leq i,j \leq N_0$.

Following <cit.>, we rewrite the multiscale basis functions as

\begin{equation}
\tcboxmath{
\forall \, K \in \mathcal{T}_H,
\qquad
\phiEps{i}
\phiPone{i} + \sum_{\alpha=1}^d
\left.\left(\partial_\alpha \phiPone{i}\right) \right\vert_K \VK{\alpha}
\quad
\text{in } K,
\label{eq:diffusion-MsFEM-Vxy-grad}
\end{equation}

for all $1 \leq i \leq N_0$, where, for each mesh element $K$, we define the numerical corrector $\VK{\alpha} \in H^1_0(\Omega)$ ($1\leq\alpha\leq d$) as the function supported by $K$ that is the unique solution to the local problem

\begin{equation}
\tcboxmath{
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\varepsilon \nabla \VK{\alpha})}$
&$\operatorname{div}(A^\varepsilon e_ \alpha)$
&in $K$,
\\
&on $\partial K$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-MsFEM-correctors}
\end{equation}

Here, $e_\alpha$ denotes the $\alpha$-th canonical unit vector of $\bbR^d$.
The expansion (<ref>) is obtained upon rewriting (<ref>) as a PDE for $\phiEps{i}-\phiPone{i}$, and then using linearity of the PDE and the fact that $\nabla\phiPone{i}$ is constant in $K$ to show that
$ \displaystyle
\sum_{\alpha=1}^d \left.\left(\partial_\alpha \phiPone{i}\right) \right\vert_K \VK{\alpha}
$ is indeed the unique solution to this PDE.

Inserting (<ref>) for the trial and test functions in (<ref>) and again exploiting the fact that all $\phiPone{i}$ have piecewise constant gradients, we obtain

\begin{align*}
\mathds{A}^\varepsilon_{j,i}
\sum_{K \in \mathcal{T}_H}
\sum_{\alpha,\beta=1}^d
\left.\left(\partial_\beta \phiPone{j}\right)\right\vert_K \left(
\int_K \left( e_\beta + \nabla \VK{\beta} \right) \cdot A^\varepsilon \left( e_\alpha + \nabla \VK{\alpha} \right)
\right)
\left.\left(\partial_\alpha \phiPone{i}\right)\right\vert_K
\\
\sum_{K \in \mathcal{T}_H}
\sum_{\alpha,\beta=1}^d
\left.\left(\partial_\beta \phiPone{j}\right)\right\vert_K
\, a_K^{\varepsilon,\dif} \left(x^\alpha + \VK{\alpha}, x^\beta + \VK{\beta}\right) \,
\left.\left(\partial_\alpha \phiPone{i}\right)\right\vert_K.
\end{align*}

Next we define the piecewise constant effective diffusion tensor $\overline{A} \in \mathbb{P}_0(\mathcal{T}_H,\,\bbR^{d\times d})$ by

\begin{equation}
\left. \overline{A}_{\beta,\alpha} \right\vert_K
\frac{1}{|K|} a^{\varepsilon,\dif}_K \left(x^\alpha + \VK{\alpha},\, x^\beta + \VK{\beta} \right)
\quad \text{for each } K \in \mathcal{T}_H \text{ and each } 1 \leq \alpha,\beta \leq d,
\label{eq:diffusion-eff-msfem-gal}
\end{equation}

where $|K|$ denotes the measure of the mesh element $K$.
Then (<ref>) can be written as

\begin{equation}
\tcboxmath{
\mathds{A}^\varepsilon_{j,i} = \int_\Omega \nabla \phiPone{j} \cdot \overline{A} \, \nabla \phiPone{i}.
\label{eq:stiffness-msfem-p1basis}
\end{equation}

Motivated by (<ref>), we introduce the coarse-scale problem

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}\left(\,\overline{A}\, \nabla u\right)}$
&in $\Omega$,
\\
&on $\partial \Omega$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-pde-effective}
\end{equation}

and its Galerkin discretization with $\Pone$ Lagrange elements: with $V_{H,0} = \operatorname{span} \left\{\phiPone{i} \mid 1 \leq i \leq N_0 \right\}$ (note that the definition of $V_{H,0}$ will be generalized in Def. <ref>), find $u_{H,0} \in V_{H,0}$ such that

\begin{equation}
\forall \, v_H \in V_{H,0}, \quad \overline{a}^{\dif}(u_H,v_H) = F(v_H),
\label{eq:diffusion-FEM-effective}
\end{equation}

where the linear form $F$ is defined in (<ref>) and the bilinear form $\overline{a}^{\dif}$ is defined as

\begin{equation*}
\forall \, u, v \in H^1_0(\Omega), \quad \overline{a}^{\dif}(u,v) = \int_\Omega \nabla v \cdot \overline{A} \, \nabla u.
\end{equation*}

Problem (<ref>) equivalently writes

\begin{equation}
\mathds{A}^{\Pone} \, U^{\Pone} = \mathds{F}^{\Pone},
\label{eq:linear-system-p1}
\end{equation}

\begin{equation}
\forall \, 1 \leq i, j \leq N_0,
\quad
\mathds{A}^{\Pone}_{j,i} = \overline{a}^{\dif} \left(\phiPone{i},\phiPone{j} \right),
\quad
\mathds{F}^{\Pone}_j = F \left( \phiPone{j} \right).
\label{eq:linear-system-effective}
\end{equation}

Comparing the expressions (<ref>) and (<ref>), we deduce that $\mathds{A}^\varepsilon = \mathds{A}^\Pone$. In other words:

The stiffness matrix of the MsFEM problem (<ref>) is identical to the stiffness matrix of the $\Pone$ problem (<ref>).

This lemma immediately implies that the $\Pone$ problem (<ref>) is well-posed, since the MsFEM (<ref>) itself is well-posed.

Let us point out that problems (<ref>) and (<ref>) are defined entirely in terms of quantities that vary only on the macroscopic scale $H$. The finite element problem (<ref>) can thus be solved using a legacy code that is designed for standard FEMs.
Lemma <ref> then suggests including the $\Pone$ approximation (<ref>) of the effective, coarse-scale problem (<ref>) as an integral part of the MsFEM approach. We do so in Algorithm <ref> below.

The right-hand side vector $\mathds{F}^\varepsilon$ in (<ref>) is, in general, different from $\mathds{F}^\Pone$ in (<ref>). Indeed, we integrate the product of $f$ with highly oscillatory basis functions in the former problem and with $\Pone$ basis functions in the latter. The solutions $U^\varepsilon$ and $U^\Pone$ to (<ref>) and (<ref>), respectively, are thus different a priori.

§.§ Non-intrusive workflow

We propose the following non-intrusive MsFEM variant:

\begin{equation}
\text{Set } u_H^\varepsilon = u_H + \sum_{K\in\mathcal{T}_H} (\partial_\alpha u_H)\vert_K \, \VK{\alpha} \in V_{H,0}^\varepsilon \text{ where } u_H \in V_{H,0} \text{ is the unique solution to~\eqref{eq:diffusion-FEM-effective}}.
\label{eq:diffusion-MsFEM-noni}
\end{equation}

The MsFEM approximation $u_H^\varepsilon$ is well-defined, since we have seen above that problem (<ref>) is well-posed.

Note that the symbol $u^\varepsilon_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^\varepsilon_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 $\varphi$ on $K$, we have

\begin{equation}
\varphi(x) = \varphi \left(x_{c,K}\right) + \sum_{\alpha=1}^d \partial_\alpha \varphi \cdot \left( x^\alpha - x^\alpha_{c,K} \right),
\label{eq:P1-expansion}
\end{equation}

where $x^\alpha$ denotes the function that to a point $x\in\Omega$ associates its $\alpha$-th coordinate, and $x_{c,K}=(x^1_{c,K},\dots,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 $\partial_\alpha u_H$ element  by element at the end of the online stage, then $u^\varepsilon_H$ defined in (<ref>) can be computed element by element according to

\begin{equation}
\forall \, K \in \mathcal{T}_H,
\quad
u_H(x_{c,K}) + \sum\limits_{\alpha=1}^d \left.\left(\partial_\alpha u_H\right)\right\vert_K \left( x^\alpha - x^\alpha_{c,K} + \VK{\alpha}(x) \right)
\quad
\text{on } K,
\label{eq:msfem-diff-noni-post}
\end{equation}

without using the legacy code.

The above observations culminate in the computational approach presented in Algorithm <ref>. We can distinguish

* the offline stage consisting of lines <ref>-<ref>,

* the online stage being executed entirely in line <ref>,

* a post-processing step in line <ref>.

Non-intrusive MsFEM approach for problem (<ref>)

[1]
Let $\mathcal{T}_H$ be the mesh used by the legacy code

$K \in \mathcal{T}_H$

$1 \leq \alpha \leq d$
Solve for $\VK{\alpha}$ defined by (<ref>)

Compute $\overline{A}\vert_K$ defined by (<ref>)

Use the legacy code to solve for $u_H$ defined by (<ref>) and to save $\left\{u_H(x_{c,K})\right\}_{K\in\mathcal{T}_H}$ and $\left\{ (\partial_\alpha u_H) \vert_K \right\}_{K\in\mathcal{T}_H, \, 1 \leq \alpha \leq d}$

Obtain the MsFEM approximation $u^\varepsilon_H$ by (<ref>)

The superiority of Algorithm <ref> over the classical MsFEM Algorithm <ref> is that the global problem of the online stage can completely be constructed and solved by the use of a pre-existing $\Pone$ 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 $\Pone$ 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 $\mathds{F}^\Pone$ (see (<ref>)) for the global problem only requires a numerical quadrature on the coarse mesh and is therefore cheaper than the construction of $\mathds{F}^\varepsilon$ (see (<ref>)), involving the multiscale basis functions and requiring numerical quadratures at the microscale.

The part of the offline stage that manipulates fine meshes (lines <ref>-<ref>) and the post-processing step can be developed independently of the legacy code used in line <ref>.
The requirement for these fine-scale solvers is that they have access to the coarse mesh $\mathcal{T}_H$ used by the global solver.
Note also that the local problem (<ref>) is only indexed by the coarse mesh element $K$, in contrast to the local problem (<ref>) 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 (<ref>), this correspondence is not needed to compute $\overline{A}$, nor for the computation of the fine-scale solution $u^\varepsilon_H$ in (<ref>), 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 (<ref>) and the linear system solved for the non-intrusive MsFEM in (<ref>), are different in general. This motivates the comparison of the non-intrusive MsFEM approach (<ref>) to the following Petrov-Galerkin MsFEM:

\begin{equation}
\text{Find } u^\varepsilon_H \in V_{H,0}^\varepsilon \text{ such that }
a^{\varepsilon,\dif} \left(u_H^\varepsilon, \phiPone{j} \right)
F \left( \phiPone{j} \right) \quad \text{for all } 1 \leq j \leq N_{0},
\label{eq:diffusion-MsFEM-testP1}
\end{equation}

based on the trial space $V_{H,0}^\varepsilon$ and the test space $V_{H,0}$ for both the bilinear and the linear form. The following result was shown in <cit.>.

The non-intrusive MsFEM variant (<ref>) coincides with the Petrov-Galerkin MsFEM (<ref>).

The non-intrusive MsFEM approach is generalized in Sec. <ref> after the development of a general framework to define a wide variety of MsFEMs in Sec. <ref>. Lemma <ref> 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 <ref>.

§.§ 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^\varepsilon$ to a PDE as $\varepsilon$ tends to $0$. This relation to $H$-convergence provides an interesting interpretation of the effective tensor $\overline{A}$ introduced in (<ref>).

Let us suppose in this section (and in this section only, except for Sec. <ref>) that $A^\varepsilon(x) = A^\mathsf{per}(x/\varepsilon)$ for some bounded, $\bbZ^d$-periodic matrix $A^{\mathsf{per}}$ satisfying the coercivity property in (<ref>). In this case, the sequence of matrices $A^\varepsilon$ 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, <cit.> or <cit.> for details on periodic homogenization.

Due to $H$-convergence of $A^\varepsilon$, the functions $u^\varepsilon$, solution to (<ref>), converge to a limit function $u^\star$ (weakly in $H^1(\Omega)$, strongly in $L^2(\Omega)$) as $\varepsilon \to 0$. The homogenized limit $u^\star$ is the solution to the homogenized equation (<ref>) below.

Let $Q$ denote the unit cube of $\bbR^d$.
We introduce the corrector functions $w_1,\dots,\,w_d \in H^1_{per}(Q)$ solution to

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\mathsf{per} \nabla w_\alpha)}$
&in $\bbR^d$,
\\
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-correctors}
\end{equation}

which uniquely defines $w_\alpha$ up to an irrelevant additive constant. The entries of the (constant) homogenized diffusion tensor $A^\star \in \bbR^{d\times d}$ are given by

\begin{equation}
\int_Q (e_\beta + \nabla w_\beta) \cdot A^\mathsf{per}(e_\alpha + \nabla w_\alpha),
\qquad
1 \leq \alpha, \, \beta \leq d.
\label{eq:diffusion-hom-coef}
\end{equation}

The homogenized limit $u^\star$ of $u^\varepsilon$ is the unique solution in $H^1_0(\Omega)$ to the boundary value problem

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\star \nabla u^\star)}$
&in $\Omega$,
\\
&on $\partial \Omega$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:diffusion-hom-pde}
\end{equation}

The truncated reconstruction of $u^\star$ that is called the first-order two-scale expansion takes the form

\begin{equation}
u^{\varepsilon,1}(x) =
u^\star(x) +
\varepsilon \sum_{\alpha=1}^d \partial_\alpha u^\star(x) \, w_\alpha \left(\frac{x}{\varepsilon}\right).
\label{eq:expansion-1-hom}
\end{equation}

Under suitable regularity assumptions, the difference $u^\varepsilon - u^{\varepsilon,1}$ converges to $0$ strongly in $H^1(\Omega)$ as $\varepsilon\to0$. This property will be used for the convergence results in Sec. <ref>.

In the periodic setting, the expansion (<ref>) can be used to construct a numerical approximation of $u^\varepsilon$, without the need of any computations at the fine scale. This approximation is presumably valid only in the regime of very small parameters $\varepsilon$ and deteriorates if $\varepsilon$ 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 $\Omega$ and that depends on an effective tensor that is itself defined in terms of the corrector functions. Details can be found, e.g., in <cit.>. This prevents the $H$-convergence theory from being directly applicable for the numerical approximation of $u^\varepsilon$.

Numerical homogenization techniques, that draw their inspiration from the various elements above, offer an alternative for the approximation of $u^\varepsilon$ that can be applied in much more general contexts. We can see the similarities between the corrector functions $w_\alpha$ in (<ref>) and the numerical correctors $\VK{\alpha}$ in (<ref>). Note that the $\VK{\alpha}$ solve problems similar to (<ref>), but that they need to be solved at the microscale and on each mesh element $K$. Similarly, we note the resemblance between the reconstruction (<ref>) and the definition of $u_H^\varepsilon$ in (<ref>), and between the homogenized coefficient $A^\star$ defined in (<ref>) and the effective macroscopic coefficient $\overline{A}$ from (<ref>). However, contrary to $A^\star$, the MsFEM quantity $\overline{A}$ has to be computed on an element-by-element basis, and it is not necessarily constant throughout $\Omega$.
Finally, the MsFEM analogue of the homogenized problem (<ref>) is the resolution of the effective macroscale problem (<ref>).

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 (<ref>) in 2D posed on the unit square. Let us consider a mesh consisting of squares that are perfectly aligned with the periodicity of $A^\varepsilon$. We solve the corrector problems (<ref>) on all square mesh elements with periodic boundary conditions and subsequently compute the effective diffusion tensor $\overline{A}$ according to (<ref>).

In this case, the problems for the numerical correctors all reduce to (<ref>) and $\overline{A}$ is constant and equal to the homogenized coefficient $A^\star$ as defined by (<ref>). A $\mathbb{Q}_1$ discretization of the effective problem (<ref>) thus constitutes a non-intrusive MsFEM that is equivalent to the $\mathbb{Q}_1$ approximation of the homogenized equation (<ref>).

§ 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 (<ref>) (e.g. advection-diffusion equations) give rise to different choices of the local problem (<ref>), depending on which terms of the global PDE are included (see, for instance, <cit.>.)

The other reason is that, even for the pure diffusion problem (<ref>), the choice of the basis functions defined in (<ref>) has an important drawback. The definition of the multiscale basis functions requires a choice of arbitrary boundary conditions on the mesh element boundary $\partial K$, since the exact boundary condition satisfied by $u^\varepsilon$ is unknown. In (<ref>), 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^\varepsilon$ near $\partial K$ if $A^\varepsilon$ is highly oscillatory and the mesh $\mathcal{T}_H$ is coarse. Variations on the definition of the functions $\phiEps{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. <ref>.

The oversampling variant of the MsFEM was introduced along with the variant based on (<ref>) at the time of its first appearance in <cit.>. For this method, an oversampling domain $S_K$ is associated to each mesh element $K$ (details are provided in Sec. <ref>). The problems (<ref>) 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 $\partial 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 <cit.>.
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 $\partial 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. <ref>.

§.§ The non-intrusive approach

The intrusiveness of the specific MsFEM-lin variant introduced in Sec. <ref> is exemplary for all MsFEMs described in Sec. <ref>. It turns out that the non-intrusive MsFEM approach introduced in <cit.> and recalled in Sec. <ref> 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 <ref> (corresponding to the identities in boxes in Sec. <ref>).

The non-intrusive MsFEM follows from the expansion (<ref>), namely the expression of the multiscale basis function as a $\Pone$ 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 (<ref>) allows the preprocessing of the microstructure independently of the global approximation indices related to the finite element method;

* the expansion (<ref>) follows from the fact that $\nabla \phiPone{i}$ is piecewise constant combined with linearity of the local problems (<ref>);

* the stiffness matrix can be formulated in terms of a piecewise constant effective diffusion tensor in (<ref>) thanks to full localization of the corrector functions, the piecewise constant gradient of $\phiPone{i}$ in the expansion (<ref>) and bilinearity of the global problem (<ref>).

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. <ref>). 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. <ref> and Def. <ref>). This involves the definition of oversampling patches (for MsFEMs with oversampling, Def. <ref>), and an extension of the notion of degrees of freedom to define the boundary conditions for the numerical correctors (Def. <ref>, <ref> and <ref>) on oversampling patches. It is then possible to define the multiscale basis functions as a generalization of (<ref>) (see Def. <ref>) and finally to define the MsFEM for our general framework in Def. <ref>.

We note that our development of non-intrusive MsFEM approaches relies to a great extent on the fact that (<ref>), and its generalization (<ref>) in the general framework developed below, provide a description of the multiscale basis functions in terms of $\Pone$ basis functions, without the need of higher-order functions. Higher-order MsFEMs can be found in <cit.>. Possible analogues of (<ref>) 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 <cit.>.

§.§ 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 <cit.> within the Exa-Dune project <cit.>. For example, the distribution of local problems over multiple processors and subsequent coupling in a global problem are handled by designated software components <cit.>. 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 $\partial \Omega$ 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 <ref>). 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. <ref>. This is done in Def. <ref> and <ref> for the numerical correctors introduced in (<ref>), and in Def. <ref> for the expansion (<ref>) 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 (<ref>) (for a Petrov-Galerkin MsFEM) and (<ref>) (for a Galerkin MsFEM) in Sec. <ref>. 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^\varepsilon$ be a continuous bilinear form on $H^1(\Omega) \times H^1(\Omega)$.
We are interested in the solution to the problem

\begin{equation}
\text{Find } u^\varepsilon \in H^1_0(\Omega) \text{ such that }
a^\varepsilon(u^\varepsilon, v) = F(v)
\text{ for any } v \in H^1_0(\Omega),
\label{eq:gen-pb}
\end{equation}

where $F$ is defined as in (<ref>) for any $f\in L^2(\Omega)$. To ensure well-posedness of (<ref>), we suppose that the bilinear form $a^\varepsilon$ is coercive on $H^1_0(\Omega)$.
The bilinear form $a^\varepsilon$ may contain coefficients that oscillate on a microscopically small scale.

The oversampling and Crouzeix-Raviart variants of the MsFEM introduced in Sec. <ref> 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^\varepsilon$ is in fact defined on the broken Sobolev space $H^1(\mathcal{T}_H) \times H^1(\mathcal{T}_H)$. More precisely, we assume that we can represent it as $a^\varepsilon = \sum\limits_{K\in\mathcal{T}_H} a^\varepsilon_K$, where, for each $K \in \mathcal{T}_H$, $a^\varepsilon_K$ is a continuous bilinear form defined on $H^1(K) \times H^1(K)$.

To ensure well-posedness of MsFEMs, which may use nonconforming approximation spaces, coercivity on $H^1_0(\Omega)$ may be insufficient. Therefore, we add the following coercivity hypothesis for the bilinear forms $a_K^\varepsilon$:

\begin{equation}
\begin{IEEEeqnarraybox}[][c]{s}
\IEEEstrut
for all $K \in \mathcal{T}_H$, there exists $\alpha_K > 0$ such that
\\
\quad
\forall \, u \in H^1(K),
\quad
\geq
\alpha_K \, \lVert \nabla u \rVert_{L^2(K)}^2.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\label{eq:gen-pb-coer}
\end{equation}

In order to perform a convergence analysis, one also has to assume that the $\alpha_K$ are bounded from below by some $\tilde{\alpha}>0$ that does not depend on $H$. We provide convergence results in Sec. <ref> for the pure diffusion problem (<ref>), in which case we have $\alpha_K=m$ from (<ref>).

As an example, the introductory problem (<ref>) with the associated bilinear form $a^{\varepsilon,\dif}$ is covered by this framework as is made explicit in Example <ref> below. Other second-order PDEs that fit in our abstract variational formulation are given in Example <ref>.

The diffusion problem (<ref>) is covered by the abstract variational formulation above. Indeed, we can set

a^\varepsilon = a^{\varepsilon,\dif},

where $a^{\varepsilon,\dif}$ is the bilinear form defined in (<ref>). Further, we define

\coloneqq
\int_K \nabla v \cdot A^\varepsilon \nabla u,

for all $u,v \in H^1(\mathcal{T}_H)$, so that we have indeed $a^{\varepsilon,\dif} = \sum\limits_{K\in\mathcal{T}_H} a^{\varepsilon,\dif}_K$. Clearly, each $a^{\varepsilon,\dif}_K$ satisfies (<ref>) with $\alpha_K=m$ the coercivity constant from (<ref>).

The reaction-advection-diffusion equation,

\begin{equation*}
-{\operatorname{div}(A^\varepsilon \nabla u^\varepsilon)} + b\cdot\nabla u^\varepsilon + \sigma u^\varepsilon= f,
\end{equation*}

with a divergence-free advection field $b:\Omega \mapsto \bbR^d$ and a positive reaction coefficient $\sigma : \Omega \mapsto \bbR$,
can be modelled (under some regularity hypotheses that we do not state here) with the bilinear forms

\begin{equation*}
\int_K \nabla v \cdot A^\varepsilon \nabla u
+ v \, b \cdot \nabla u + \sigma u v.
\end{equation*}

However, these bilinear forms $a_K^\varepsilon$ do not satisfy (<ref>) even though the bilinear form $a^\varepsilon$ is coercive on $H^1(\Omega)$. To this end, a skew-symmetrized formulation of the transport term can be used.
The skew-symmetrized formulation uses the bilinear form

\begin{equation*}
\int_K \nabla v \cdot A^\varepsilon \nabla u
+ \frac{1}{2} v \, b \cdot \nabla u
- \frac{1}{2} u \, b \cdot \nabla v
+ \sigma u v,
\end{equation*}

which does satisfy (<ref>). However, the assumption (<ref>) is used for proving well-posedness of the MsFEM in Lemma <ref>, but both choices for $a_K^{\varepsilon}$ mentioned here can be studied in practice.
We refer e.g. to <cit.> for more details.
Within the general MsFEM framework, $b$ and $\sigma$ are allowed to be highly oscillatory, and this may impact the specific MsFEM strategy to be preferred.

§.§ Piecewise affine structure

In Sec. <ref>, 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.

Let a mesh $\mathcal{T}_H$ be given. The underlying  $\Pone$ space for the MsFEM, denoted $V_H$, is one of the following two spaces: the Lagrange approximation space

\begin{equation*}
V_H^{L} =
\{ v \in \Pone(\mathcal{T}_H) \mid v \text{ is continuous on } \Omega \},
\end{equation*}

in which case we shall refer to the associated MsFEM as the MsFEM-lin, or the Crouzeix-Raviart approximation space

\begin{equation*}
V_H^{CR} =
\left\{ v \in \Pone(\mathcal{T}_H) \mid \forall\ K \in \mathcal{T}_H, \ \forall e \in \mathcal{F}(K) \text{ such that } e \subset \Omega : \ \int_e \llbracket v \rrbracket = 0 \right\},
\end{equation*}

in which case the associated MsFEM shall be called the MsFEM-CR. We use the notation $\mathcal{F}(K)$ for the set of faces of $K$ and $\llbracket v \rrbracket$ denotes the jump of $v$ over the face $e$. The space $V_H^L$ is a subspace of $H^1(\Omega)$, but $V_H^{CR}$ is not. Note that no restrictions apply on faces lying on $\partial\Omega$.

We note that the underlying $\Pone$ space has the following property: if $v \in V_H$ is piecewise constant on the mesh $\mathcal{T}_H$, then $v$ is constant in $\Omega$. Contrary to the space $V_H^L$, functions in the Crouzeix-Raviart space $V_H^{CR}$ 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 $\Pone$ spaces that we introduced in Def. <ref>. The definition involves an arbitrary simplex $K$, which is typically an element of the mesh $\mathcal{T}_H$, or an associated oversampling patch (for the oversampling technique of the MsFEM) that we shall define in Def. <ref>. The latter is not always a simplex, and we extend Def. <ref> to such oversampling patches in Def. <ref> and <ref>.

A degree of freedom operator ( operator) $\Gamma$ associates to any simplex $K \subset \bbR^d$ and $v \in \Pone(K)$ a vector $\Gamma(K,v) \in \bbR^{d+1}$, whose components are called the degrees of freedom of $v$ on $K$, in such a way that the application $v \mapsto \Gamma(K,v)$ is a linear bijection from $\Pone(K)$ to $\bbR^{d+1}$. More precisely, $\Gamma(K,\cdot)$ will denote in the sequel one of the following two operators:

* ( operator for the MsFEM-lin.) Let $x_{0},\dots,x_{d}$ denote the vertices of $K$. We set

\begin{equation*}
\forall \, v \in \Pone(K),
\qquad
\Gamma^L(K,v) = \left( v(x_0),\dots,v(x_d) \right).
\end{equation*}

For $K \in \mathcal{T}_H$, the degree of freedom $[\Gamma^{L}(K,\cdot)]_j$ is said to be associated to the boundary if, for all $v \in \Pone(K)$, $[\Gamma^{L}(K,v)]_j = v(x)$ for a vertex $x$ of the mesh that lies on $\partial \Omega$.

* ( operator for the MsFEM-CR.) Let $e_0,\dots,e_d$ denote the faces of $K$. We set

\begin{equation*}
\forall \, v \in \Pone(K),
\qquad
\Gamma^{CR}(K,v) = \left( \frac{1}{|e_0|} \int_{e_0} v,\dots,\frac{1}{|e_d|} \int_{e_d} v \right).
\end{equation*}

For $K \in \mathcal{T}_H$, the degree of freedom $[\Gamma^{CR}(K,\cdot)]_j$ is said to be associated to the boundary if, for all $v \in \Pone(K)$, $\displaystyle [\Gamma^{CR}(K,v)]_j = \frac{1}{|e|} \int_{e} v$ for a face $e$ of the mesh that lies on $\partial \Omega$.

The  $\Pone$ test space is defined as

\begin{equation*}
V_{H,0} =
\left\{
v \in V_H \ \left\vert \
\begin{aligned}
&\forall \, K \in \mathcal{T}_H, \ \forall \, 1 \leq j \leq d+1, [\Gamma(K,v)]_j=0 \text{ if the degree} \\
&\text{of freedom } [\Gamma(K,\cdot)]_j \text{ is associated to the boundary}, \
\end{aligned}
\right.
\right\}.
\end{equation*}

The $\Pone$ test space is used in practice to approximate the subspace $H^1_0(\Omega)$ of $H^1(\Omega)$.
The degrees of freedom are defined element per element and are thus local.
Global properties of the underlying $\Pone$ space $V_H$ are most easily made explicit through the identification of a basis for $V_H$.

Let $V_H$ be an underlying $\Pone$ space as in Def. <ref>, and let $\Gamma$ be the associated  operator. We shall denote by $N$ the dimension of $V_H$. The $\Pone$ basis functions $\phiPone{1},\dots,\phiPone{N}$ are defined as follows:

* For the MsFEM-lin, let $x_1,\dots,x_N$ be an enumeration of the (internal and boundary) vertices of $\mathcal{T}_H$. Then $\phiPone{i}$ is defined by $\phiPone{i}(x_j) = \delta_{i,j}$ for all $1 \leq i,j \leq N$.

* For the MsFEM-CR, let $e_1,\dots,e_N$ be an enumeration of the (internal and boundary) faces of $\mathcal{T}_H$. Then $\phiPone{i}$ is defined by $\displaystyle \frac{1}{|e_j|} \int_{e_j} \phiPone{i} = \delta_{i,j}$ for all $1 \leq i,j \leq N$.

In both cases, these functions form a basis of the corresponding space $V_H$ of Def. <ref>.

§.§ Local problems

§.§.§ Oversampling patches

To replace the (standard) underlying $\Pone$ space by a space of the same (low) dimension, adapted to the microstructure of $a^\varepsilon$, we associate to each mesh element $K \in \mathcal{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.

Let $K\in\mathcal{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 $\rho\geq1$. The oversampling patch $S_K$ is defined as $S_K = S_K' \cap \Omega$.

See Fig. <ref> for an illustration of the construction of oversampling patches in dimension 2. In this work, we allow for the trivial homothety ratio $\rho=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.

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 $\Omega$ 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. <ref> is a simplex. However, for mesh elements close to the boundary $\partial \Omega$, 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 <cit.>.

§.§.§ Degrees of freedom on oversampling patches

Definition <ref> provides the definition of  operators on any simplex. For the MsFEM, we wish to compute multiscale functions on oversampling patches $S_K$, in which case Def. <ref> may be insufficient. We illustrated this in Fig. <ref>. 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 $\Pone$ basis function, we will still need a notion of  operator such that $\Gamma(S_K,\cdot)$ is a linear bijection from $\Pone(S_K)$ to $\bbR^{d+1}$. Therefore, we extend the definition of the degrees of freedom operators $\Gamma^L$ and $\Gamma^{CR}$ in Def. <ref> and <ref>.

Let $K\in\mathcal{T}_H$ and let $S_K$ be its associated oversampling patch. Let $x_0,\dots,x_d$ be a selection of $d+1$ distinct vertices of $S_K$. We define the  operator $\Gamma^L$ by

\begin{equation*}
\forall \, v \in \Pone(S_K),
\qquad
\Gamma^L(S_K,v) = \left( v(x_0),\dots,v(x_d) \right).
\end{equation*}

We note that any choice of $d+1$ nodal values unequivocally characterizes an affine function on $S_K$. Hence, $\Gamma^L(S_K,\cdot)$ is indeed a bijection. Since $\Gamma(S_K,\cdot)$ will only be used to identify elements of $\Pone(S_K)$ on the boundary of $S_K$, the particular choice of the vertices $x_0,\dots,x_d$ is unimportant in the above definition. Finally, when $S_K$ is a simplex, it has only $d+1$ vertices and Def. <ref> reduces to Def. <ref>.

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 $\mathcal{F}_a(S_K)$. The other faces of $S_K$ are referred to as the dilated faces, collected in the set $\mathcal{F}_d(S_K)$.
When the patch $S_K$ does not touch $\partial \Omega$, we have $\mathcal{F}_d(S_K)=\mathcal{F}(S_K)$ and $\mathcal{F}_a(S_K) = \emptyset$.
In Fig. <ref>, for example, the additional faces are exactly those faces that lie on $\partial\Omega$. This is not always the case, as is illustrated by Fig. <ref>.

For the definition of $\Gamma^{CR}(S_K,\cdot)$, we shall rely on the existence of $d+1$ dilated faces, because we need $\Gamma^{CR}(S_K,\cdot)$ to be a bijection between $\Pone(S_K)$ and $\bbR^{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. <ref>, the lower right dilated face falls outside $\Omega$ if the homothety ratio is too large, and the oversampling patch $S_K$ only has two dilated faces (edges here) and two additional faces.

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.

Let $K\in\mathcal{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,\dots,e_d$. We define the  operator $\Gamma^{CR}$ by

\begin{equation*}
\forall \, v \in \Pone(S_K),
\qquad
\Gamma^{CR}(S_K,v) = \left( \frac{1}{|e_0|} \int_{e_0} v,\dots,\frac{1}{|e_d|} \int_{e_d} v \right).
\end{equation*}

When $S_K$ is a simplex, we have $\mathcal{F}_d(S_K) = \mathcal{F}(S_K)$, and Def. <ref> coincides with the respective elements of Def. <ref>.

§.§.§ Numerical correctors: first oversampling strategy

We now provide the precise assumptions under which we will consider local problems, i.e., the analogues of (<ref>) defining the MsFEM-lin basis functions and the definition of the numerical correctors in (<ref>). In fact, since the numerical correctors play an essential role in the construction of non-intrusive MsFEM approaches, we define the numerical correctors first and use them to define the multiscale basis functions in Def. <ref>.

The two possible functional settings for all these constructions are provided by Def. <ref> and <ref>. 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 <ref> and <ref>).

Let $K\in\mathcal{T}_H$, let $S_K$ be its associated oversampling patch and let $\Gamma$ be a  operator from Def. <ref>. A subspace $V_K$ of $H^1(S_K)$ and bilinear form $s^\varepsilon_K : V_K \times V_K \mapsto \bbR$ are called sampling space and sampling form, respectively, if they satisfy the following:

* the space $V_K$ contains the space of affine functions $\Pone(S_K)$;

* the operator $\Gamma(S_K,\cdot)$ is well-defined on $V_K$;

* the -extended local problem: find $v \in V_K$ such that

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts}
\IEEEstrut
$s_K^\varepsilon(v, w)$
&$\langle g, w\rangle$
for all $w \in V_{K,0}$,
\\
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:sampling-space-pde-OS-dofe}
\end{equation}

has a unique solution for any $ g \in (H^1(S_K))'$. Here, $
V_{K,0} = \{ w \in V_K \mid \Gamma(S_K,w) = 0\}
$ is the sampling test space.

Problem (<ref>) is called `-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^\varepsilon$ shall be used to encode the oscillations of the bilinear form $a^\varepsilon$ 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^\varepsilon_K$ of the original problem (<ref>), or only some of them.
When the MsFEM was first proposed in <cit.>, it was suggested that $s_K^\varepsilon$ 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 <cit.>.

In the functional setting of Def. <ref>, the generalization of (<ref>) to define the numerical correctors for the general MsFEM framework is as follows.

We introduce for all $K\in\mathcal{T}_H$, for any $0\leq\alpha\leq d$, the function $\VSK{\alpha}{1} \in V_{K,0}$ as the unique solution to the corrector problem

\begin{equation}
\tcboxmath[colback=white, colframe=black]{
\forall \, w \in V_{K,0},
\quad
s_K^\varepsilon \left( \VSK{\alpha}{1}, w \right) =
\begin{cases}
\IEEEstrut
{-s_K^\varepsilon} \left( 1, w \right)
& \text{if } \alpha=0,
\\
{-s_K^\varepsilon} \left( x^\alpha - x^\alpha_{c,K}, w \right)
& \text{if } 1\leq\alpha\leq d,
\end{cases}
\label{eq:MsFEM-gen-correctors-dofe}
\end{equation}

The -extended numerical corrector $\VKOS{\alpha}{1}$ is defined as the restriction of $\VSK{\alpha}{1}$ to $K$, extended to all of $\Omega$ by $0$.

Note that the above definition introduces one more numerical corrector than introduced in (<ref>). The precise definition of the numerical correctors is chosen such that the analogous expansion of (<ref>) for the general framework (see (<ref>)) leads to a PDE for the multiscale basis functions analogous to (<ref>); we show this in Lemma <ref> and (for a second oversampling strategy introduced below) in Lemma <ref>. In the following example, we see that Def. <ref> is indeed a generalization of the numerical correctors defined by (<ref>) in Sec. <ref>.

[MsFEM-lin for diffusion problems]
We consider $V_H = V_H^{L}$ and $\Gamma = \Gamma^L$ from Def. <ref>.
For the diffusion problem (<ref>), we have $a^\varepsilon = a^{\varepsilon,\dif}$ and we set $s_K^\varepsilon = a_K^{\varepsilon,\dif}$ (see Example <ref>).
The sampling space for the MsFEM-lin is defined as

\begin{equation*}
\coloneqq
\left\{
v \in H^1(S_K) \ \mid \ \exists \, w \in \Pone(S_K) \text{ such that } v\vert_{\partial S_K} = w\vert_{\partial S_K}
\right\}.
\end{equation*}

Then the sampling test space $V^{L}_{K,0}$ is the space $H^1_0(S_K)$.
In this case, it holds $a_K^{\varepsilon,\dif}(1,w)=0$ for all $w \in V^{L}_{K,0}$.
Consequently, the -extended numerical corrector $\VKOS{0}{1}$ is identically equal to $0$; we obtain indeed exactly $d$ numerical correctors as in Sec. <ref>.
For the non-trivial numerical correctors, Def. <ref> corresponds to the weak formulation of the following boundary value problem:

\begin{equation*}
-{\operatorname{div}(A^\varepsilon \nabla \VSK{\alpha}{1})} = \operatorname{div}(A^\varepsilon e_\alpha) \
\text{ in } S_K,
\quad
\VSK{\alpha}{1} = 0 \
\text{ on } \partial S_K,
\end{equation*}

which is clearly well-posed.

[MsFEM-CR for diffusion problems]
Taking $a^\varepsilon$, $a^\varepsilon_K$ and $s_K^\varepsilon$ as in the previous example, we construct the MsFEM-CR with the sampling space $V_K^{CR} \coloneqq H^1(S_K)$.
With $V_H = V_H^{CR}$ and $\Gamma = \Gamma^{CR}$ from Def. <ref>, the -extended numerical corrector $\VKOS{\alpha}{1}$ is obtained from the boundary value problem:

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\varepsilon \nabla \VSK{\alpha}{1})}$
&$\operatorname{div}(A^\varepsilon e_\alpha)$
&in $S_K$,
\\
$\vec{n} \cdot A^\varepsilon \nabla \VSK{\alpha}{1}$
&$-{\vec{n} \cdot A^\varepsilon e_\alpha}$
&on each $h \in \mathcal{F}_a(S_K)$, \\
$\vec{n} \cdot A^\varepsilon \nabla \VSK{\alpha}{1}$
&$c_h - \vec{n} \cdot A^\varepsilon e_\alpha $
&on each $h \in \mathcal{F}_d(S_K)$, \\
$\displaystyle \frac{1}{|h|} \int_h \VSK{\alpha}{1}$
&for each $h \in \mathcal{F}_d(S_K)$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:msfem-cr-dofe}
\end{equation}

where $\vec{n}$ denotes the outward unit vector on $\partial 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 (<ref>), whereas the flux on the dilated faces of $S_K$ involves an additional constant, due to the fact that the test functions in $V_{K,0}^{CR}$ cannot take arbitrary values on the dilated faces. Indeed, their mean vanishes on these faces according to Def. <ref>.

When $S_K=K$ and when the faces of $K$ do not lie on $\partial \Omega$, this corresponds to the setting of the original MsFEM-CR defined in <cit.>. The latter work also provides an alternative characterization of the multiscale Crouzeix-Raviart space.

When a face $e$ of $K$ lies on $\Omega$, the basis functions that we will obtain do not satisfy $\phiEps{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 <cit.>. The MsFEM-CR with local boundary conditions as defined here was studied in <cit.>.

When $S_k=K$ (i.e., in the absence of oversampling), the  operator allows us to prescribe certain continuity properties on the faces of the mesh elements $K$. More precisely, when the MsFEM-lin with  operator $\Gamma^L$ is employed, the numerical correctors $\VKOS{\alpha}{1}$ vanish at the vertices of the mesh, and, with the correct choice of sampling space (see Example <ref>), they vanish on all faces of $K$ and are thus continuous on $\Omega$. When the MsFEM-lin with  operator $\Gamma^{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. <ref>, in the vein of the expansion (<ref>)) shows that this means that the continuity properties of the $\Pone$ basis functions of the underlying $\Pone$ 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 -extended local problems for our local computations, since the values on $\partial K$ are not controlled by the degrees of freedom $\Gamma(S_K,\cdot)$ on $\partial S_K$. Therefore, we introduce another variant of the local problems and the -extended numerical correctors in the next section.

§.§.§ Numerical correctors: second oversampling strategy

Let $K \in \mathcal{T}_H$ and let $V_K$ and $s_K^\varepsilon$ be a sampling space and sampling form, respectively, according to Def. <ref>. Additionally, suppose that the operator $\Gamma(K,\cdot)$ is well-defined on $V_K$. Then a -continuous local problem is to find $v\in V_K$ such that

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts}
\IEEEstrut
$s_K^\varepsilon(v, w)$
&$\langle g, w\rangle$
for all $w \in V_{K,0}$,
\\
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:sampling-space-pde-OS-dofc}
\end{equation}

for some $ g \in (H^1(S_K))'$.

Suppose any -continuous local problem in Def. <ref> is well-posed. Then we introduce, for all $K\in\mathcal{T}_H$ and all $0 \leq \alpha \leq d$, the functions $\VSK{\alpha}{2}$ as the unique functions in $V_K$ with $\Gamma\left(K, \VSK{\alpha}{2}\right)=0$ satisfying the corrector problem (<ref>). We define the -continuous numerical correctors $\VKOS{\alpha}{2}$ as the restriction of $\VSK{\alpha}{2}$ to $K$, extended to all of $\Omega$ by $0$.

We emphasize that the local problems of Def. <ref> and <ref> use test functions $w$ in the same space $V_{K,0}$. This means that the test functions satisfy $\Gamma(S_K,w)=0$, rather than $\Gamma(K,w)=0$.

Clearly, when $S_K=K$, there is no difference between the -extended and -continuous problems (<ref>) and (<ref>).
We shall in this case simply refer (<ref>) (or (<ref>)) as local problems, and we write $\VKOS{\alpha}{0} = \VK{\alpha}$ for the numerical correctors of MsFEMs without oversampling.

[MsFEM-lin for diffusion problems]
Continuing Example <ref>, consider now the -continuous numerical corrector $\VKOS{\alpha}{2}$. Equation (<ref>) solves the following problem: there exists $w \in \Pone(S_K)$ such that

\begin{equation*}
-{\operatorname{div}(A^\varepsilon \nabla \VKOS{\alpha}{2})} = \operatorname{div}(A^\varepsilon e_\alpha) \
\text{ in } S_K,
\quad
\VKOS{\alpha}{2} = w \
\text{ on } \partial S_K,
\quad
\VKOS{\alpha}{2} = 0 \
\text{ at the vertices of } K.
\end{equation*}

The boundary condition on $\partial S_K$ is `completed' by a condition at the vertices of $K$. Except when $A^\varepsilon$ is constant (and, consequently, $\phiEps{i}$ is affine on $S_K$), it is not evident whether a solution to this problem exists.

[MsFEM-CR for diffusion problems]
For the MsFEM-CR considered in Example <ref>, the -continuous numerical correctors satisfy the same problem (<ref>) as the -extended numerical correctors, but with the average condition replaced by $\displaystyle \frac{1}{|h|} \int_h \VKOS{\alpha}{2} = 0$ for each $h \in \mathcal{F}(K)$.
As we saw for the MsFEM-lin, this is not a standard boundary value problem on $S_K$.

Examples <ref> and <ref> show that a -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. <ref>. This is not the case for -continuous problems, for which one finds some boundary conditions on $ \partial 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 $\partial K$ that are explicitly imposed through the degrees of freedom on $K$ in (<ref>). Well-posedness is not obvious in general, and cannot always be deduced from well-posedness of the -extended counterpart (<ref>). See Examples <ref> and <ref> for concrete local problems. We now address the well-posedness of -continuous problems in more detail in Sec. <ref>.

§.§.§ Well-posedness of -continuous numerical correctors

We have seen in Examples <ref> and <ref> that -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 -continuous numerical correctors, and we use this strategy to discuss the well-posedness of the associated local problems.

In Def. <ref> we assume the well-posedness of -extended problems, and we have seen in Examples <ref> and <ref> that this is a natural assumption. It is also natural to assume that we can compute -extended numerical correctors numerically. We compute the -continuous numerical correctors from the -extended numerical correctors, by subtracting a linear combination of suitable functions $W^\beta$ from the -extended numerical correctors. The $W^\beta$ must all satisfy the homogeneous equation
s_K^\varepsilon \left( W^\beta , w \right) = 0
for all
w \in V_{K,0},
in order not to perturb the local problem (<ref>) that is already satisfied by both types of numerical correctors.
We shall use the functions $W^0 \coloneqq 1 + \VSK{0}{1}$ and $W^\beta \coloneqq x^\beta - x^\beta_{c,K} + \VSK{\beta}{1}$ for $1 \leq \beta \leq d$, where $\VSK{\beta}{1}$ is defined in Def. <ref>. The precise strategy is as follows.

Fix $0 \leq \alpha \leq d$. We look for coefficients $c_0^\alpha, \dots, c_d^\alpha$ such that
\displaystyle
\VSK{\alpha}{2}
\VSK{\alpha}{1}
- \sum_{\beta=0}^d c_\beta^\alpha \, W^\beta
on $K$, where we recall that $\VSK{\alpha}{2}$ is defined by Def. <ref>. Note that both sides of the equation clearly satisfy (<ref>). The desired equality thus holds if and only if
\displaystyle
\Gamma \left( K,
\VSK{\alpha}{1} - \sum_{\beta=0}^d c_\beta^\alpha \, W^\beta
\right)
Since the  operators are linear, this leads to the linear system

\begin{equation}
\underbrace{
\begin{bmatrix}
\vert & \vert & & \vert \\
\Gamma(K,W^0) & \Gamma(K,W^1) & \hdots & \Gamma(K,W^d) \\
\vert & \vert & & \vert \\
\end{bmatrix}
\displaystyle
\eqqcolon \mathds{M}
\begin{bmatrix}
c_0^\alpha \\ c_1^\alpha \\ \vdots \\ c_d^\alpha
\end{bmatrix}
\Gamma\left(K,\VSK{\alpha}{1}\right).
\label{eq:glue-system}
\end{equation}

Invertibility of the matrix $\mathds{M}$ is thus a sufficient condition for the existence of all -continuous numerical correctors, and the resolution of the linear system (<ref>) for each $\alpha$ (where all -extended numerical correctors are replaced by their numerical approximation) allows to compute the -continuous numerical correctors numerically.

Before studying the invertibility of the matrix $\mathds{M}$ in a few special cases, let us consider the matrix composed of the degrees of freedom on $S_K$, i.e., the matrix

\begin{equation*}
\widetilde{\mathds{M}}
\coloneqq
\begin{bmatrix}
\vert & \vert & & \vert \\
\Gamma(S_K,W^0) & \Gamma(S_K,W^1) & \hdots & \Gamma(S_K,W^d) \\
\vert & \vert & & \vert \\
\end{bmatrix}.
\end{equation*}

By definition of the functions $\VSK{\beta}{}$, we have $\Gamma(S_K,W^\beta) = \Gamma(S_K,x^\beta - x_{c,K}^\beta)$ for $1\leq \beta \leq d$, and $\Gamma(S_K, W^{0}) = \Gamma(S_K,1)$.
Note that the constant function together with the coordinate functions $x^\beta-x^\beta_{c,K}$ ($1 \leq \beta \leq d$) span $\Pone(S_K)$. Since $\Gamma(S_K,\cdot)$ is a bijection, the vectors $\Gamma(S_K, W^0),\dots,\Gamma(S_K,W^{d})$ are linearly independent. Hence the matrix $\widetilde{\mathds{M}}$ is invertible.
One may hope that the linear independence of the vectors $\Gamma(S_K,W^\beta)$ is preserved for the degrees of freedom on the interior boundary $\partial K$ instead of $\partial S_K$, yielding invertibility of $\mathds{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 $\mathds{M}$ in a few special cases. When $s_K^\varepsilon$ is the sampling form that was used in Example <ref> (corresponding to a diffusion problem; we will consider this case until the end of this section) and if $A^\varepsilon$ is constant, all numerical correctors vanish on $S_K$ and the foregoing argument for the matrix $\widetilde{\mathds{M}}$ shows invertibility of $\mathds{M}$.

In the periodic setting (see Sec. <ref>), even though $A^\varepsilon$ itself is not constant, its homogenized limit $A^*$ is. In this case, the $\VSK{\beta}{1}$ converge to zero weakly in $H^1(S_K)$. (We show this in Lemma <ref> in the absence of oversampling, but the argument can be generalized to -extended oversampling.)
Now consider the MsFEM-CR. The weak convergence of the $\VSK{\beta}{1}$ 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}(\partial K)$ in $L^2(\partial K)$ is compact, the $\VSK{\beta}{1}$ converge to $0$ strongly in $L^2$ on each face of $K$. Consequently, the degrees of freedom $\Gamma(K, \VSK{\beta}{1})$ (the averages on the faces of $K$) converge to zero as $\varepsilon\to0$. Thus, $\Gamma(K,W^0) \to \Gamma(K,1)$ and $\Gamma(K,W^\beta) \to \Gamma(K,x^\beta - x_{c,K}^\beta)$ as $\varepsilon\to0$ for all $1 \leq \beta \leq d$ and, by the above argument for the matrix $\widetilde{\mathds{M}}$, the matrix $\mathds{M}$ is invertible in this limit. By continuity of the determinant function, the matrix $\mathds{M}$ is invertible when $\varepsilon$ is small enough, and the -continuous basis functions exist in this regime.

The study of the -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 <cit.>, to see that the multiscale basis functions, obtained from the numerical correctors in Def. <ref> below, are in fact continuous for any bounded diffusion tensor. (See Example <ref> 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 -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 (<ref>) in terms of the numerical correctors. We recall that in Sec. <ref>, 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 <ref> and <ref>. Recall that $\phiPone{1},\dots,\,\phiPone{N}$ is a basis of the space $V_{H}$ (see Def. <ref>). We can suppose that the first $N_{0}$ basis functions form a basis of $V_{H,0}$. The following definition is the generalization of (<ref>) to the general MsFEM framework.

For each $i=1,\dots,N$, the multiscale basis function $\phiEps{i}$ is defined by

\begin{equation}
\tcboxmath[colback=white, colframe=black]{
\forall \, K \in \mathcal{T}_H,
\qquad
\left. \phiEps{i} \right\vert_K
\left. \phiPone{i} \right\vert_K + \phiPone{i}(x_{c,K}) \, \VKOS{0}{0} + \sum_{\alpha=1}^d \partial_\alpha \left(\left. \phiPone{i} \right\vert_K \right) \VKOS{\alpha}{0},
\label{eq:MsFEM-Vxy}
\end{equation}

where $\bullet = \mathsf{e}$ corresponds to -extended multiscale basis functions and $\bullet = \mathsf{c}$ corresponds to -continuous multiscale basis functions.

The -extended multiscale basis functions satisfy the variational formulation in the following lemma on the oversampling patches $S_K$.

Let $K$ be any mesh element and let $1 \leq i \leq N$. Consider an MsFEM with -extended basis functions. Define an extension of $\phiEps{i}$ from $K$ to $S_K$ by

\begin{equation}
\widehat{\phiEps{i}}
\widehat{\left.\phiPone{i}\right\vert_K} + \phiPone{i} \left(x_{c,K} \right) \, \VSK{0}{1} + \sum_{\alpha=1}^d \partial_\alpha \left(\left. \phiPone{i} \right\vert_K \right) \VSK{\alpha}{1},
\quad
\text{in } S_K,
\label{eq:MsFEM-VxyOS-dofe}
\end{equation}

where $\widehat{ \left. \phiPone{i} \right\vert_K }$ denotes the affine extension of $\left. \phiPone{i} \right\vert_K$ to $S_K$, and $\VKOS{\alpha}{1}$ is as in Def. <ref>.
Then $\widehat{\phiEps{i}}$ is the unique solution in $V_K$ to

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts}
\IEEEstrut
$s_K^\varepsilon \left(\widehat{\phiEps{i}}, w \right)$
for all $w \in V_{K,0}$,
\\
$\Gamma \left( S_K, \widehat{\phiEps{i}} \right)$
&$\Gamma \left( S_K , \widehat{\left.\phiPone{i}\right\vert_K} \right)$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:msbasis-OS-vf-dofe}
\end{equation}

In the case of the MsFEM-lin for the diffusion problem (<ref>), problem (<ref>) with $S_K=K$ coincides with the definition of the multiscale basis functions in (<ref>); see Example <ref>.

Problem (<ref>) has a unique solution in view of Def. <ref>. It thus suffices to show that $\widehat{\phiEps{i}}$ satisfies (<ref>).
Since the numerical correctors $\VKOS{\alpha}{1}$ belong to $V_{K,0}$ for all $1 \leq \alpha \leq N$, it is clear from (<ref>) that $\Gamma\left(S_K,\widehat{\phiEps{i}}\right)=\Gamma \left( S_K , \widehat{\left.\phiPone{i}\right\vert_K} \right)$.

Inserting (<ref>) into (<ref>) and applying (<ref>) to all $\VKOS{\alpha}{1}$, we find, for any test function $w\in V_{K,0}$,

\begin{align*}
s_K^\varepsilon \left(\widehat{\phiEps{i}}, w \right)
s_K^\varepsilon \left(\widehat{\phiEps{i}}, \widehat{\left.\phiPone{i}\right\vert_K} \right)
\phiPone{i} \left( x_{c,K} \right) s_K^\varepsilon \left( \widehat{\phiEps{i}}, \VSK{0}{1} \right)
\sum_{\alpha=1}^d
\left.\left(  \partial_\alpha \phiPone{i} \right) \right\vert_K s_K^\varepsilon \left( \widehat{\phiEps{i}}, \VSK{\alpha}{1}\right).\\
s_K^\varepsilon \left(\widehat{\phiEps{i}}, \widehat{\left.\phiPone{i}\right\vert_K} \right)
\phiPone{i} \left( x_{c,K} \right) s_K^\varepsilon \left( \widehat{\phiEps{i}}, 1 \right)
\sum_{\alpha=1}^d
\left.\left( \partial_\alpha  \phiPone{i} \right) \right\vert_K s_K^\varepsilon \left( \widehat{\phiEps{i}}, x^\alpha - x^\alpha_{c,K}\right)\\
s_K^\varepsilon \left(\widehat{\phiEps{i}}, \widehat{\left.\phiPone{i}\right\vert_K} \right)
s_K^\varepsilon \left(
\widehat{\phiEps{i}}, \,
\phiPone{i} \left( x_{c,K} \right)
+ \sum_{\alpha=1}^d \left.\left( \partial_\alpha \phiPone{i} \right) \right\vert_K \left( x^\alpha - x^\alpha_{c,K} \right)
\right)
\end{align*}

Here we use that $s_K^\varepsilon$ is a bilinear form on $V_K$, that all piecewise affine functions are contained in $V_K$ according to Def. <ref>, and the property that $\nabla \phiPone{i}$ is piecewise affine. Finally, we use (<ref>) for $\varphi=\widehat{\phiPone{i}}$ to conclude that

\begin{equation*}
s_K^\varepsilon \left(\widehat{\phiEps{i}}, w \right)
s_K^\varepsilon \left(\widehat{\phiEps{i}}, \widehat{\left.\phiPone{i}\right\vert_K} \right)
s_K^\varepsilon \left(\widehat{\phiEps{i}}, \widehat{\left.\phiPone{i}\right\vert_K} \right)
\end{equation*}

which establishes the desired variational formulation satisfied by $\widehat{\phiEps{i}}$.

If the -continuous problems (<ref>) are well-posed, we obtain by the same analysis the following result for -continuous multiscale basis functions.

Let $K$ be any mesh element and let $1 \leq i \leq N$. Assume that any -continuous local problem (<ref>) is well-posed. Consider an MsFEM with -continuous oversampling. Define an extension of $\phiEps{i}$ from $K$ to $S_K$ by

\begin{equation*}
\widehat{\phiEps{i}}
\widehat{\left.\phiPone{i}\right\vert_K} + \phiPone{i} \left(x_{c,K} \right) \, \VSK{0}{2} + \sum_{\alpha=1}^d \partial_\alpha \left(\left. \phiPone{i} \right\vert_K \right) \VSK{\alpha}{2},
\quad
\text{in } S_K,
\end{equation*}

where $\VKOS{\alpha}{2}$ is as in Def. <ref>.
Then $\widehat{\phiEps{i}}$ is the unique solution in $V_K$ to

\begin{equation}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts}
\IEEEstrut
$s_K^\varepsilon \left(\widehat{\phiEps{i}}, w \right)$
for all $w \in V_{K,0}$,
\\
$\Gamma \left( K,\widehat{\phiEps{i}} \right)$
&$\Gamma \left( K , \phiPone{i} \right)$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\label{eq:msbasis-OS-vf-dofc}
\end{equation}

[MsFEM-lin for diffusion problems]
In the setting of Example <ref>, any -extended multiscale basis function $\phiEps{i}$ for the MsFEM-lin constructed in (<ref>) is obtained, in each mesh element $K$, as the restriction of a function $\widehat{\phiEps{i}}$, which is the unique solution in $H^1(S_K)$ to

\begin{equation*}
-{\operatorname{div}(A^\varepsilon \nabla \widehat{\phiEps{i}})} = 0 \
\text{ in } S_K,
\quad
\widehat{\phiEps{i}} = \widehat{\phiPone{i}} \
\text{ on } \partial S_K.
\end{equation*}

For a -continuous basis function, $\widehat{\phiEps{i}}$ solves the same PDE in $S_K$, is affine on $\partial S_K$, and satisfies $\widehat{\phiEps{i}}(x_j) = \phiPone{i}(x_j)$ at all vertices $x_j$ of $K$.

[MsFEM-CR for diffusion problems]
In the continuation of Example <ref>, the -extended multiscale basis function $\phiEps{i}$ for the MsFEM-CR is the restriction to $K$ of $\widehat{\phiEps{i}}$, the unique solution in $H^1(S_K)$ to

\begin{equation*}
\left\{
\begin{IEEEeqnarraybox}[][c]{uts?s}
\IEEEstrut
$-{\operatorname{div}(A^\varepsilon \nabla \widehat{\phiEps{i}})}$
&in $S_K$,
\\
$\vec{n} \cdot A^\varepsilon \nabla \widehat{\phiEps{i}}$
&on each $h \in \mathcal{F}_a(S_K)$, \\
$\vec{n} \cdot A^\varepsilon \nabla \widehat{\phiEps{i}}$
&on each $h \in \mathcal{F}_d(S_K)$, \\
$\displaystyle \frac{1}{|h|} \int_h \widehat{\phiEps{i}}$
&$\displaystyle \frac{1}{|h|} \int_h \widehat{\phiPone{i}}$
&for each $h \in \mathcal{F}_d(S_K)$,
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\right.
\end{equation*}

where the constants $c_h$ are uniquely determined by the problem. For -continuous basis functions, the last condition is applied to the faces $h \in \mathcal{F}(K)$ (and all other conditions remain unchanged).

Our general framework allows two characterizations  of the multiscale basis functions, namely (<ref>) and (<ref>) or (<ref>), as was the case for the MsFEM studied in Sec. <ref> (where $\phiEps{i}$ is given by (<ref>) or (<ref>)). The essential advantage of (<ref>) is that the microscale is fully encoded in the numerical correctors $\VKOS{\alpha}{0}$, that can be computed element per element without any global information. In particular, the global index $i$ of the multiscale basis function $\phiEps{i}$ is irrelevant for the computation of the numerical correctors. The expression in (<ref>) is therefore the crucial relationship that we will employ to develop non-intrusive MsFEMs within the general framework in Sec. <ref>.

The second formulation of the multiscale basis functions, as solutions to the local problems (<ref>) or (<ref>), 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 $\Pone$ basis functions and the associated multiscale basis function. This is useful in particular for the well-posedness of the MsFEM, that we study in Lemma <ref>.

Our definition of the multiscale basis functions in (<ref>) is reminiscent of the Variational Multiscale Method, a framework developed in <cit.> to adapt Galerkin approximations on low-dimensional 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. <cit.>.

The first introduction of the MsFEM in <cit.> corresponds to the idea of oversampling with -continuous basis functions. Although their existence cannot be established in general, they are computed numerically by taking linear combinations of -extended basis functions (following an analogous strategy to the one we discussed in Sec. <ref>). The MsFEM with -extended basis functions is studied in the works <cit.> 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 <ref>, 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^\varepsilon$ and $V_{H,0}^\varepsilon$, as follows:

\begin{equation*}
V_H^\varepsilon = \left\{
\phiEps{i} \ \mid \ 1 \leq i \leq N
\right\},
\qquad
V_{H,0}^\varepsilon = \left\{
\phiEps{i} \ \mid \ 1 \leq i \leq N_0
\right\}.
\end{equation*}

Note that we only use $V_{H,0}^\varepsilon$ in the present section, because (<ref>) is posed with homogeneous Dirichlet conditions, but that the larger space $V_H^\varepsilon$ is used for more general boundary conditions (see Sec. <ref>). Applying (<ref>), we have the equivalent characterization in terms of the $\Pone$ space $V_H$,

\begin{equation*}
\left\{
\left.
v_H +
\sum_{K\in\mathcal{T}_H} \left(
v_H(x_{c,K}) \, \VKOS{0}{0}
+ \sum_{\alpha=1}^d
\partial_\alpha \left( \left. v_H \right\vert_K \right)
\VKOS{\alpha}{0}
\right)
\,\right\vert\,
v_H \in V_H
\right\}.
\end{equation*}

Let $V_H$ be an underlying $\Pone$ space defined in Def. <ref> with the associated  operator $\Gamma$ from Def. <ref> and Def. <ref>-<ref>. Define for each mesh element $K \in \mathcal{T}_H$ an oversampling patch (Def. <ref>), a sampling space and sampling form in accordance with Def. <ref>. Let the multiscale basis functions $\phiEps{i}$ be given as in Def. <ref>. Then a Multiscale Finite Element Method (MsFEM) for problem (<ref>) is: find $u^\varepsilon_H \in V_{H,0}^\varepsilon$ such that

\begin{equation}
\forall \, v_H^\varepsilon \in V_{H,0}^\varepsilon,
\qquad
\sum_{K\in\mathcal{T}_H} a_K^{\varepsilon} \left( u^\varepsilon_H,\, v_H^\varepsilon \right) = F \left( v_H^\varepsilon \right).
\label{eq:gen-MsFEM-OS}
\end{equation}

In the following lemma, we investigate the well-posedness of the MsFEM.

Consider an MsFEM without oversampling, or an MsFEM with oversampling using -continuous basis functions (assuming the associated basis functions are well-defined). When $a^\varepsilon$ satisfies (<ref>), the MsFEM (<ref>) has a unique solution.

Note that, with -continuous oversampling, but also without oversampling, the multiscale basis functions satisfy (<ref>). In particular, all degrees of freedom of $u_H^\varepsilon$ related to the boundary vanish. Also note that, the dimension of $V_{H,0}^\varepsilon$ being finite, it suffices to show that $u^\varepsilon_H=0$ is the unique solution to problem (<ref>) with $F=0$.

If $0 = F(u^\varepsilon_H) = a^\varepsilon(u_H^\varepsilon, u_H^\varepsilon)$, it follows from (<ref>) that $u^\varepsilon_H$ is piecewise constant.
Let us write $\displaystyle u_H^\varepsilon = \sum_{i=1}^{N_{0}} \alpha_i \, \phiEps{i}$ for some coefficients $\alpha_i\in\bbR$ and introduce the function $ \displaystyle u_H = \sum_{i=1}^{N_{0}} \alpha_i \, \phiPone{i} \in V_H$. Because of (<ref>), we have $\Gamma(K,u_H^\varepsilon) = \Gamma(K,u_H)$ for all mesh elements $K$. Since $u_H^\varepsilon$ is piecewise constant and $\Gamma(K,\cdot)$ is a bijection from $\Pone(K)$ to $\bbR^{d+1}$ (recall Def. <ref>), it follows that $u_H^\varepsilon=u_H$. In particular, the multiscale function $u_H^\varepsilon$ in fact belongs to the underlying $\Pone$ space $V_H$.

We remarked immediately below Def. <ref> that, for either of the two spaces $V_H=V_H^{L}$ or $V_H^{CR}$, the above implies that $u^\varepsilon_H$ is constant throughout $\Omega$. Since the degrees of freedom of $u^\varepsilon_H$ associated to the boundary vanish, we readily deduce that $u^\varepsilon_H=0$.

We do not know of the existence of a result on the well-posedness of MsFEMs with oversampling using -extended multiscale basis functions. In <cit.>, the authors establish an inf-sup result for a variant of the MsFEM-lin-OS with $\Pone$ test functions (see also Def. <ref>). This result is obtained for a periodic diffusion coefficient in the limit of sufficiently small $\varepsilon$.

§ 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. <ref>. We have seen in Lemma <ref> 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. <ref> and <ref> to obtain a non-intrusive MsFEM approach:

* the expansion (<ref>) allows to recast the matrix $\mathds{A}^\varepsilon$ of the linear system for the MsFEM as the matrix $\mathds{A}^\Pone$ associated to the $\Pone$ discretization of an effective problem;

* we approximate the right-hand side $\mathds{F}^\varepsilon$ of the MsFEM problem by the right-hand side $\mathds{F}^\Pone$ of this $\Pone$ discretization;

* the post-processing step (<ref>) applied to the $\Pone$ 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. <ref> is given by (<ref>) and that it can be rewritten in terms of the bilinear forms $a_K^\varepsilon$ satisfying (<ref>). Petrov-Galerkin variants of the multiscale finite element method with $\Pone$ test functions were previously studied in <cit.>. In our general MsFEM framework, the adaptation of Def. <ref> to a Petrov-Galerkin MsFEM is the following.

Let $V_H$ be an underlying $\Pone$ space defined in Def. <ref> with the associated  operator $\Gamma$ from Def. <ref> and Def. <ref>-<ref>. Define for each mesh element $K \in \mathcal{T}_H$ an oversampling patch (Def. <ref>), a sampling space and sampling form in accordance with Def. <ref>. Let the multiscale basis functions $\phiEps{i}$ be given as in Def. <ref>. Then a Petrov-Galerkin Multiscale Finite Element Method (PG-MsFEM) for problem (<ref>) is: find $u^\varepsilon_H \in V_{H,0}^\varepsilon$ such that

\begin{equation}
\forall \, v_H \in V_{H,0},
\qquad
\sum_{K\in\mathcal{T}_H} a_K^{\varepsilon} (u^\varepsilon_H,\, v_H) = F(v_H).
\label{eq:gen-MsFEM-OS-testP1}
\end{equation}

When confusion may arise, we shall refer to the MsFEM defined in Def. <ref> 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 (<ref>) to Lemma <ref>.

We now execute step <ref> of the summary of the non-intrusive MsFEM approach at the beginning of this section. The matrix $\mathds{A}^\varepsilon$ of the linear system associated to (<ref>) is defined by

\begin{equation}
\mathds{A}^\varepsilon_{j,i}
\sum_{K\in\mathcal{T}_H} a^\varepsilon_K \left(\phiEps{i}, \, \phiPone{j} \right),
\quad
1 \leq i,j \leq N_{0}.
\label{eq:system-msfem-gen-testP1}
\end{equation}

To find an effective $\Pone$ formulation with the same linear system, we will use the definition (<ref>) of the multiscale basis functions in the general framework, but first we combine it with (<ref>) applied to $\varphi=\phiPone{i}$ to rewrite (<ref>) as

\begin{align}
\left. \phiEps{i} \right\vert_K
\phiPone{i}(x_{c,K}) + \phiPone{i} \left( x_{c,K} \right) \VKOS{0}{0} + \sum_{\alpha=1}^d \left. \partial_\alpha \phiPone{i} \right\vert_K \left( x^\alpha - x^\alpha_{c,K} + \VKOS{\alpha}{0} \right)
\nonumber \\
\phiPone{i} \left( x_{c,K} \right) \, \VVK{0} + \sum_{\alpha=1}^d \left. \partial_\alpha \phiPone{i} \right\vert_K \VVK{\alpha},
\label{eq:MsFEM-VVxy}
\end{align}

\[
\VVK{0} \coloneqq 1+\VKOS{0}{0},
\quad
\VVK{\alpha} \coloneqq x^\alpha-x^\alpha_{c,K} + \VKOS{\alpha}{0},
\]

for all $1 \leq \alpha \leq d$ and each $K \in \mathcal{T}_H$. We recall that $\bullet \in \{\mathsf{e},\mathsf{c}\}$ indicates the choice of -extended or -continuous basis functions.
Inserting (<ref>) into (<ref>) for $\phiEps{i}$ and (<ref>) for $\varphi=\phiPone{j}$ yields

\begin{align*}
\mathds{A}^\varepsilon_{j,i}
\sum_{K\in\mathcal{T}_H} \left(
\phiPone{i}\left(x_{c,K}\right) \, a_K^\varepsilon\left(\VVK{0},1\right) \, \phiPone{j}\left(x_{c,K}\right)
+ \sum_{\alpha=1}^d \left.\left(\partial_\alpha\phiPone{i}\right)\right\vert_K \, a_K^\varepsilon \left(\VVK{\alpha},1\right) \, \phiPone{j}\left(x_{c,K}\right)
\right.
\nonumber\\
\sum_{\beta=1}^d \phiPone{i}\left(x_{c,K}\right) \, a_K^\varepsilon\left(\VVK{0},x^\beta-x_{c,K}^\beta\right) \, \left.\left(\partial_\beta\phiPone{j}\right)\right\vert_K
\nonumber\\
\left. \sum_{\alpha,\beta=1}^d
\left.\left(\partial_\alpha\phiPone{i}\right)\right\vert_K \, a^\varepsilon_K\left(\VVK{\alpha},x^\beta-x_{c,K}^\beta\right) \, \left.\left(\partial_\beta\phiPone{j}\right)\right\vert_K
\right),
\end{align*}

[ams equation]
|K|   ( M i   j )(x_c,K)
+ ∫_K j   B^1 ·∇i
+ i   B^2 ·∇j
+ ∇j ·A   ∇i,

where we have defined the effective mass $\overline{M}$, (adjoint) advection vector $\overline{B}^1$ and $\overline{B}^2$, and the effective diffusion tensor $\overline{A}$, for all $1 \leq \alpha, \beta \leq d$ and for each $K\in\mathcal{T}_H$, as

\begin{equation}
\begin{IEEEeqnarraybox}[][c]{us?us}
\IEEEstrut
$\left. \overline{M} \right\vert_K $
&$\displaystyle = \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{0},1\right)$,
&$\left. \overline{B}^1_\alpha \right\vert_K$
&$\displaystyle = \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{\alpha},1\right)$,
\\
$\left. \overline{B}^2_\beta \right\vert_K$
&$\displaystyle = \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{0},x^\beta-x_{c,K}^\beta\right)$,
&$\left. \overline{A}_{\beta,\alpha} \right\vert_K $
&$\displaystyle = \frac{1}{|K|} \, a^\varepsilon_K \left(\VVK{\alpha}, x^\beta-x_{c,K}^\beta\right)$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\label{eq:gen-MsFEM-eff-coef}
\end{equation}

Note that all integrals in (<ref>) can be computed exactly by evaluating the integrand at the centroid. With this quadrature rule, we observe that the term $|K| \left( \overline{M} \phiPone{i} \, \phiPone{j} \right)(x_{c,K})$ also equals the numerical approximation of the integral $\displaystyle \int_K \overline{M} \phiPone{i} \phiPone{j}$.

The new expression (<ref>) for the matrix of the linear system motivates us to introduce the effective bilinear forms $\overline{a}_K$ defined on $H^1(K) \times H^1(K)$ by

\begin{equation}
\label{eq:gen-MsFEM-eff-form}
\overline{a}_K(u,v)
\int_K \nabla v \cdot \overline{A} \, \nabla u
+ v \left(\overline{B}^1 \cdot \nabla u\right)
+ u \left(\overline{B}^2 \cdot \nabla v\right)
+ \overline{M} \, u \, v,
\quad
\text{for all } u,\,v \in H^1(K),
\end{equation}

and the associated $\Pone$ Galerkin approximation on the space $V_{H,0}$:

\begin{equation}
\text{Find } u_H\in V_{H,0} \text{ such that}
\sum_{K\in\mathcal{T}_H}
\overline{a}_K(u_H,v_H)
\quad \text{for all } v_H \in V_{H,0}.
\label{eq:gen-FEM-effective}
\end{equation}

This discrete problem leads to a linear system with the matrix

\begin{equation*}
\mathds{A}^\Pone_{j,i}
\overline{a} \left( \phiPone{i},\, \phiPone{j} \right)
\sum_{K\in\mathcal{T}_H} \overline{a}_K(\phiPone{i}, \, \phiPone{j}),
\quad
1 \leq i,j \leq N_0.
\end{equation*}

The identity (<ref>) thus implies the following result, which generalizes Lemma <ref> to the PG-MsFEM in the general framework.

The matrices $\mathds{A}^\varepsilon$ and $\mathds{A}^\Pone$ are identical if the integrals in (<ref>) are evaluated at the centroid of each $K$ for the computation of $\mathds{A}^\Pone$. Then the PG-MsFEM (<ref>) coincides with the resolution of the effective problem (<ref>) combined with the post-processing step

\begin{equation}
\label{eq:gen-MsFEM-post}
\left. u_H^\varepsilon \right\vert_K
u_H(x_{c,K}) \, \VVK{0} + \sum_{\alpha=1}^d \left. \partial_\alpha u_H \right\vert_K \VVK{\alpha}.
\end{equation}

Note that step <ref> of the summary at the beginning of this section is irrelevant for the PG-MsFEM. The computation of the right-hand side in (<ref>) is clearly part of any standard FEM software.

The computational approach described by Lemma <ref> naturally fits within the non-intrusive workflow of Algorithm <ref>. The numerical correctors on line <ref> are, of course, replaced by those of Def. <ref> or Def. <ref>. Line <ref> is replaced by the computation of all effective quantities in (<ref>). The online phase in line <ref> amounts to solving the $\Pone$ problem (<ref>), 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^\varepsilon_H$ from $u_H$ by virtue of (<ref>).

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. <ref>), we need to replace the $\Pone$ test space of the PG-MsFEM by the multiscale test space $V_{H,0}^\varepsilon$.
The matrix of the linear system associated to (<ref>) is given by

\begin{equation*}
\mathds{A}_{j,i}^{\varepsilon,\mathsf{G}}
\sum_{K\in\mathcal{T}_H} a_K^\varepsilon \left( \phiEps{i},\phiEps{j} \right),
\quad
1 \leq i,j \leq N_{0}.
\end{equation*}

Upon inserting (<ref>) for the test function $\phiEps{j}$, we find, for all $1 \leq i,j \leq N_0$,

\begin{equation*}
\mathds{A}_{j,i}^{\varepsilon,\mathsf{G}}
\mathds{A}_{j,i}^\varepsilon
+ \sum_{K\in\mathcal{T}_H} \left(
\phiPone{j}\left(x_{c,K}\right) \, a^\varepsilon_K\left(\phiEps{i}, \VKOS{0}{0}\right)
\sum_{\beta=1}^d \left.\left(\partial_\beta \phiPone{j}\right) \right\vert_K \, a^\varepsilon_K\left(\phiEps{i},\VKOS{\beta}{0}\right)
\right),
\end{equation*}

where $\mathds{A}^\varepsilon$ is the matrix of the Petrov-Galerkin MsFEM, see (<ref>) and (<ref>).

An effective formulation can again be derived by inserting (<ref>) for the $\phiEps{i}$. We obtain
[ams equation]
|K|   ( M^𝖦 i   j )(x_c,K)
+ ∫_K j   B^1,𝖦 ·∇i
+ i   B^2,𝖦 ·∇j
+ ∇j ·A^𝖦   ∇i,

where the effective mass, (adjoint) advection vectors and diffusion tensor are given by (using those defined in (<ref>))

\begin{equation}
\begin{IEEEeqnarraybox}[][c]{us?us}
\IEEEstrut
$\left. \overline{M}^\mathsf{G} \right\vert_K $
&$\displaystyle = \left. \overline{M} \right\vert_K + \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{0},\VKOS{0}{0}\right)$,
&$\left. \overline{B}^{1,\mathsf{G}}_\alpha \right\vert_K$
&$\displaystyle = \left. \overline{B}^{1}_\alpha \right\vert_K + \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{\alpha},\VKOS{0}{0}\right)$,
\\
$\left. \overline{B}^{2,\mathsf{G}}_\beta \right\vert_K$
&$\displaystyle = \left. \overline{B}^{2}_\beta \right\vert_K + \frac{1}{|K|} \, a_K^\varepsilon\left(\VVK{0},\VKOS{\beta}{0}\right)$,
&$\left. \overline{A}^\mathsf{G}_{\beta,\alpha} \right\vert_K$
&$\displaystyle = \left. \overline{A}_{\beta,\alpha} \right\vert_K + \frac{1}{|K|} \, a^\varepsilon_K \left(\VVK{\alpha}, \VKOS{\beta}{0}\right)$.
\IEEEstrut
\IEEEeqnarraynumspace %preserve the same alignment as obtained in an equation environment
\end{IEEEeqnarraybox}
\label{eq:gen-MsFEM-eff-coef-Gal}
\end{equation}

The above computations lead to the introduction of the effective bilinear form $\overline{a}^{\mathsf{G}} = \sum\limits_{K\in\mathcal{T}_H} \overline{a}^{\mathsf{G}}_K$ with

\begin{equation}
\label{eq:gen-MsFEM-eff-form-Gal}
\overline{a}_K^\mathsf{G}(u,v)
\int_K \nabla v \cdot \overline{A}^\mathsf{G} \, \nabla u
+ v \left(\overline{B}^{1,\mathsf{G}} \cdot \nabla u\right)
+ u \left(\overline{B}^{2,\mathsf{G}} \cdot \nabla v\right)
+ \overline{M}^\mathsf{G} \, u \, v.
\end{equation}

We formulate the following effective variational problem:

\begin{equation}
\text{Find } u_H\in V_{H,0} \text{ such that}
\sum_{K\in\mathcal{T}_H}
\overline{a}^{\mathsf{G}}_K(u_H,v_H)
\quad \text{for all } v_H \in V_{H,0}.
\label{eq:gen-FEM-effective-gal}
\end{equation}

The associated linear system has coefficients $\mathds{A}_{j,i}^{\Pone,\mathsf{G}} =  \overline{a}^{\mathsf{G}} \left( \phiPone{i},\phiPone{j} \right)$. We have the following analogue of Lemma <ref>, which generalizes Lemma <ref> to the G-MsFEM in the general framework.

The matrices $\mathds{A}^{\varepsilon,\mathsf{G}}$ and $\mathds{A}^{\Pone,\mathsf{G}}$ are identical if the integrals in (<ref>) are evaluated at the centroid of each $K$ in the computation of $\mathds{A}^{\Pone,\mathsf{G}}$.

Contrary to the matrices, the right-hand sides of the effective problem (<ref>) and the Galerkin MsFEM (<ref>) are not equal in general. We apply step <ref> 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:

\begin{equation}
\text{Find } u^\varepsilon_H \in V_{H,0}^\varepsilon \text{ such that }
\sum_{K \in \mathcal{T}_H}
a_K^{\varepsilon}\left(u_H^\varepsilon, \phiEps{j} \right)
F\left(\phiPone{j}\right) \quad \text{for all } 1 \leq j \leq N_{0}.
\label{eq:gen-MsFEM-noni}
\end{equation}

This problem is no longer a Galerkin approximation of (<ref>), because different test spaces are used for the bilinear and for the linear form. In view of Lemma <ref>, the non-intrusive MsFEM can equivalently be formulated as

\begin{equation*}
\text{compute } u_H \in V_{H,0} \text{ solution to~\eqref{eq:gen-FEM-effective-gal} and compute } u^\varepsilon_H \text{ from } u_H \text{ by~\eqref{eq:gen-MsFEM-post}},
\end{equation*}

provided all integrals in (<ref>) are evaluated at the centroid for the construction of the matrix of the linear system in (<ref>).

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 <ref>. For completeness, we provide the algorithm for the non-intrusive G-MsFEM in Algorithm <ref>.

Non-intrusive G-MsFEM for the general framework

[1]
Let $\mathcal{T}_H$ be the mesh used by the legacy code, let $\bullet \in \{\mathsf{e},\mathsf{c}\}$ be the chosen oversampling variant

$K \in \mathcal{T}_H$
$0 \leq \alpha \leq d$
Solve for the applicable $\VKOS{\alpha}{0}$ from Def. <ref> or <ref>
Compute the effective tensors defined in (<ref>)
Use the legacy code to construct the matrix $\mathds{A}^\Pone$ by evaluating (<ref>) at the centroid of each mesh element and to solve for $u_H$ defined by (<ref>)
Save $\left\{u_H(x_{c,K})\right\}_{K\in\mathcal{T}_H}$ and $\left\{ (\partial_\alpha u_H) \vert_K \right\}_{K\in\mathcal{T}_H, \, 1 \leq \alpha \leq d}$

Obtain the MsFEM approximation $u^\varepsilon_H$ from (<ref>)

The discussion surrounding Algorithm <ref> 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 -continuous oversampling, provide well-posedness of the G-MsFEM (<ref>) by Lemma <ref>. In this case, the non-intrusive approximation (<ref>) is also well-posed, for the matrices associated to both MsFEM variants are the same.
Regarding the PG-MsFEM (<ref>), 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 <ref> to the general framework.

Consider a G-MsFEM as defined by Def. <ref> without oversampling and suppose that the sampling form $s^\varepsilon_K$ equals the local bilinear form $a_K^\varepsilon$. Then the matrix associated to this G-MsFEM coincides with the matrix associated to the corresponding PG-MsFEM of Def. <ref>. Consequently, the non-intrusive Galerkin MsFEM (<ref>) coincides with the Petrov-Galerkin MsFEM of Def. <ref> and in particular, The Petrov-Galerkin MsFEM is well-posed.

To prove the lemma, we show that the matrices corresponding to the linear problems defined in (<ref>) and (<ref>) are equal. Using $s_K^\varepsilon = a_K^{\varepsilon,\dif}$, we have for all $1 \leq i, \, j \leq N_{0}$,

\begin{equation}
\mathds{A}^{\varepsilon,\mathsf{G}}_{j,i} - \mathds{A}^\varepsilon_{j,i}
\sum_{K\in\mathcal{T}_H} a_K^{\varepsilon,\dif}\left(\phiEps{i},\phiEps{j} - \phiPone{j}\right)
\sum_{K\in\mathcal{T}_H} s_K^\varepsilon \left(\phiEps{i}, \phiEps{j}-\phiPone{j} \right)
\label{eq:IBP-bubbles}
\end{equation}

Indeed, the multiscale basis functions satisfy $\Gamma(K,\phiEps{i}) = \Gamma(K,\phiPone{i})$ for all $K$, so that $\phiEps{j} - \phiPone{j} \in V_{K,0}$ (see (<ref>) with $S_K=K$ and recall Def. <ref> for the sampling test space $V_{K,0}^ \varepsilon$), and the variational problem in (<ref>) (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 non-intrusive MsFEM approaches. For more details, we refer to <cit.>.

Stabilized finite element formulations.
Stabilized finite element formulations add mesh-dependent terms to a discrete variational formulation (such as (<ref>)) to remove numerical instabilities, for example caused by sharp boundary layers of the exact solution. See <cit.> for such a variant of the MsFEM and see <cit.> for the stabilization of single-scale problems. The expansion (<ref>) 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 $\Pone$ 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 (<ref>). See e.g. <cit.>. An expansion of the kind (<ref>) 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 (<ref>) posed on the space $V_{H,0}$ and complemented with non-homogeneous Dirichlet conditions on $\partial\Omega$. This solution can directly be used to construct a multiscale approximation $u^\varepsilon_H \in V_H^\varepsilon$ from (<ref>). The translation of the Dirichlet condition to the MsFEM approximation is as follows: if -continuous oversampling is applied, the function $u^\varepsilon_H$ satisfies $[\Gamma(K,u_H^\varepsilon)]_j = [\Gamma(K,u_H)]_j$ for all degrees of freedom associated to the boundary. Here, $[\Gamma(K,u_H)]_j$ is determined by the legacy code. When -extended oversampling is used, the degrees of freedom associated to the boundary are equal to the sum of $[\Gamma(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 $\partial \Omega$, one solves a Galerkin approximation of the variational formulation in the space $V_H^\varepsilon$. The suitable adaptation of (<ref>) can be approximated by a non-intrusive Galerkin MsFEM following the same methodology as above. The effective $\Pone$ 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 (<ref>), the Neumann boundary condition in the effective problem is imposed on the effective flux $\vec{n} \cdot \overline{A} \nabla u_H$, where $\overline{A}$ is defined in (<ref>).

Parabolic equations.
When a parabolic equation is discretized in time, problems of the form (<ref>) 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^\varepsilon$, 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 $\Pone$ in (<ref>). 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 $\Pone$ 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^\varepsilon$ on the coarse scale. The HMM directly proposes to solve a $\Pone$ problem on the coarse scale, where effective coefficients of the $\Pone$ 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^\varepsilon$ at the coarse and the microscale by the use of multiscale basis functions, like the MsFEM. It is shown in <cit.> that a Petrov-Galerkin LOD method (also see <cit.>) can, with some additional approximations, be recast as the $\Pone$ 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 $\Pone$ 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 $\Pone$ 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^\varepsilon_K = a^{\varepsilon,\dif}_K$ defined in Example <ref>, and we choose the sampling form $s_K^\varepsilon = a^{\varepsilon,\dif}_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. <ref> and <ref> for the numerical correctors, Equation (<ref>) reduces to

\begin{equation}
a_K^{\varepsilon,\dif} \left( \VSK{\alpha}{0}, w \right)
-a_K^{\varepsilon,\dif} \left( x^\alpha, w \right),
\label{eq:diffusion-gen-correctors}
\end{equation}

for all $w\in V_{K,0}$ (where $V_{K,0}$ is the sampling test space for either the MsFEM-lin or the MsFEM-CR; see Examples <ref> and <ref>)
when $1 \leq \alpha \leq d$, whereas $\VSK{0}{0} = 0$. (The notation $\VK{\alpha}$ will be used in the absence of oversampling, see Rem. <ref>.) This means that $\VVK{0}=1$ in (<ref>). Consequently, regarding the formulation of the effective $\Pone$ problem, only the effective diffusion coefficient does not vanish in (<ref>) and (<ref>). Its definition in (<ref>) is identical to the formula in (<ref>) for the applicable choice of the numerical correctors.

The definition of the multiscale basis functions by (<ref>) reduces to (<ref>) (again upon replacing the numerical correctors $\VK{\alpha}$ by the relevant ones for the MsFEM under consideration). Hence, we can associate a multiscale counterpart in $V_H^\varepsilon$ to any $v_H \in V_H$, given by

\begin{equation}
v_H^\varepsilon = v_H
+ \sum_{K\in\mathcal{T}_H} \sum_{\alpha=1}^d
(\partial_\alpha v_H)\vert_K \, \VKOS{\alpha}{0}.
\label{eq:diffusion-MsFEM-Vxy}
\end{equation}

The non-intrusive MsFEM (<ref>) becomes

\begin{equation}
\text{Find } u^\varepsilon_H \in V_{H,0}^\varepsilon \text{ such that }
\sum_{K \in \mathcal{T}_H}
a_K^{\varepsilon,\dif}\left(u_H^\varepsilon, \phiEps{j} \right)
F\left(\phiPone{j}\right) \quad \text{for all } 1 \leq j \leq N_{0}.
\label{eq:diffusion-gen-MsFEM-noni}
\end{equation}

Lemma <ref> now amounts to the following.

Let `MsFEM' refer to the MsFEM-lin or the MsFEM-CR without oversampling. The non-intrusive Galerkin MsFEM (<ref>) coincides with the following Petrov-Galerkin MsFEM:

\begin{equation}
\text{Find } u^\varepsilon_H \in V_{H,0}^\varepsilon \text{ such that }
\sum_{K \in \mathcal{T}_H}
a_K^{\varepsilon,\dif} \left(u_H^\varepsilon, \phiPone{j} \right)
F \left( \phiPone{j} \right) \quad \text{for all } 1 \leq j \leq N_{0},
\label{eq:diffusion-gen-MsFEM-testP1}
\end{equation}

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 <ref>, <ref> and <ref> are generalizations of results of <cit.>, where the MsFEM-lin without oversampling is treated.

§.§ Convergence results

We estimate here the difference between the solutions to the (intrusive) Galerkin approximation (<ref>) and the non-intrusive MsFEM (<ref>), which coincides with the Petrov-Galerkin MsFEM (<ref>). We first show coercivity of the effective diffusion tensor $\overline{A}$.

Consider the MsFEM-lin or the MsFEM-CR without oversampling, or the MsFEM-CR with -continuous oversampling. The effective tensor $\overline{A}$ defined by (<ref>) with the appropriate numerical correctors satisfies

\begin{equation*}
\forall \, \xi \in \bbR^d,
\quad
m |\xi|^2 \leq \xi \cdot \overline{A} \, \xi.
\end{equation*}

Here, $m$ is the same coercivity constant as in (<ref>).

Let $\xi = (\xi_1,\dots,\,\xi_d) \in \bbR^d$, and let $K$ be any simplex of the mesh $\mathcal{T}_H$. We have

\begin{equation*}
|K| \, \xi \cdot \left. \overline{A} \right\vert_K \xi
\sum_{\alpha,\beta=1}^d a^{\varepsilon,\dif}_K \left(
\xi_\alpha \left(x^\alpha + \VK{\alpha}\right),\,
\xi_\beta \left(x^\beta + \VK{\beta} \right)
\right)
\int_K (\xi + \nabla \chi^\xi) \cdot A^\varepsilon (\xi + \nabla \chi^\xi),
\end{equation*}

denoting by $\chi^\xi$ the function $\displaystyle \chi^\xi = \sum_{\alpha=1}^d \xi_\alpha \VK{\alpha}$. Using (<ref>), we obtain

\begin{equation*}
|K| \, \xi \cdot \left. \overline{A} \right\vert_K \xi
\geq
m \int_K \left\vert \xi + \nabla \chi^\xi \right\rvert^2
\geq
m\,|K| \, |\xi|^2 +
2 m\,\int_K \xi \cdot \nabla \chi^\xi.
\end{equation*}

Using an integration by parts, we see that $\displaystyle \int_K \xi \cdot \nabla \chi^\xi = \int_{\partial K} \chi^\xi \, n\cdot\xi$, where $n$ is the unit outward normal vector on $\partial K$.
In the case of the MsFEM-lin, the function $\chi^\xi$ vanishes on $\partial K$. In the case of the MsFEM-CR with -continuous oversampling, or without oversampling, the function $\chi^\xi$ has average zero on each face of $K$. Since the factor $n\cdot\xi$ is constant on each face, the integral again vanishes. In conclusion, we have $\displaystyle \int_K \xi \cdot \nabla \chi^\xi = 0$.

We thus obtain the inequality
$    \xi \cdot \left. \overline{A} \right\vert_K \xi
\geq m |\xi|^2.
Since $K \in \mathcal{T}_H$ is arbitrary here, this shows coercivity of $\overline{A}$ and completes the proof.

Coercivity of the effective tensor $\overline{A}$ implies coercivity of the bilinear form $\overline{a}^\dif$ on $H^1_0(\Omega)$. By an application of the Lax-Milgram Theorem, we conclude that the (continuous) effective problem (<ref>) is well-posed for the MsFEM-lin and the MsFEM-CR without oversampling, and for the MsFEM-CR with -continuous oversampling.

The proof of the above lemma does not extend to the MsFEM-lin with oversampling, because there is no global information about $\VK{\alpha}$ on the faces of $K$.

The following lemma provides a variational characterization of the bijection (<ref>).

Consider the MsFEM-lin or the MsFEM-CR without oversampling. Let $v_H^\varepsilon\in V_H^\varepsilon$. The unique $v_H \in V_H$ for which (<ref>) holds, is the unique solution in $V_H$ to the problem

\begin{equation}
\forall \, w_H \in V_H,
\quad
\overline{a}^\dif(v_H,w_H)
\label{eq:diffusion-MsFEM-Vxy-var-char}
\end{equation}

In addition, we have, with the constants $m$ and $M$ from (<ref>), the estimate

\begin{equation*}
\lVert \nabla v_H \rVert_{L^2(\mathcal{T}_H)}
\leq
\frac{M}{m} \lVert \nabla v_H^\varepsilon \lVert_{L^2(\mathcal{T}_H)}.
\end{equation*}

Let $v_H \in V_H$ be the unique element of $V_H$ such that $v_H^\varepsilon$ and $v_H$ satisfy (<ref>). Take any $w_H \in V_H$. Using that $\nabla v_H$ and $\nabla w_H$ are piecewise constant, we compute

\begin{equation*}
\sum_{K\in\mathcal{T}_H} \sum_{\alpha,\beta=1}^d
(\partial_\beta w_H)\vert_K \,
a^{\varepsilon,\dif}_K \left( x^\alpha + \VK{\alpha}, x^\beta \right)
(\partial_\alpha v_H)\vert_K.
\end{equation*}

For the MsFEM without oversampling, the numerical correctors belong to the sampling test space $V_{K,0}$. We can thus use (<ref>) to obtain

\begin{equation*}
\forall \, 1 \leq \alpha, \, \beta \leq d,
\quad
a^{\varepsilon,\dif}_K \left( x^\alpha + \VK{\alpha}, x^\beta \right)
a^{\varepsilon,\dif}_K \left( x^\alpha + \VK{\alpha}, x^\beta + \VK{\beta} \right).
\end{equation*}

Using the definitions of $\overline{A}$ in (<ref>) and the of $\overline{a}^{\dif}$ in (<ref>) (we recall that these expressions hold true here upon replacing the numerical correctors by those under consideration), we conclude that

\begin{equation*}
\sum_{K\in\mathcal{T}_H} \sum_{\alpha,\beta=1}^d \int_K
\partial_\beta w_H \,
\overline{A}_{\beta,\alpha} \,
\partial_\alpha v_H
\overline{a}^{\dif}(w_H,v_H).
\end{equation*}

It follows that $v_H$ satisfies (<ref>). In addition, in view of the coercivity of $\overline{A}$ established in Lemma <ref> and by the Lax-Milgram Theorem, problem (<ref>) uniquely characterizes $v_H$.

The estimate on $v_H$ follows by testing the characterization (<ref>) against $w_H=v_H$. This yields, for any $K\in\mathcal{T}_H$,

\begin{align*}
m \left\lVert \nabla v_H \right\rVert_{L^2(K)}^2
\overline{a}^{\dif}(v_H,v_H)
a^{\varepsilon,\dif} \left( v_H^\varepsilon,v_H \right)
\int_K \nabla v_H \cdot A^\varepsilon \, \nabla v_H^\varepsilon
\leq
M \, \lVert \nabla v_H \rVert_{L^2(K)} \, \lVert \nabla v_H^\varepsilon \rVert_{L^2(K)}.
\end{align*}

The first inequality follows from coercivity of $\overline{A}$ and the second inequality from the upper bound on $A^\varepsilon$ in (<ref>) and the Cauchy-Schwarz inequality. We thus find that $\displaystyle \lVert \nabla v_H \rVert_{L^2(K)} \leq \frac{M}{m} \lVert \nabla v_H^\varepsilon \rVert_{L^2(K)}$.
The desired result is obtained upon squaring the inequality and summing over all $K\in\mathcal{T}_H$.

For the remainder of this section, we will consider MsFEMs without oversampling. Let $u^{\varepsilon,\mathsf{G}}_H$ denote the solution to the MsFEM approximation (<ref>) (we use the superscript $\mathsf{G}$ to stress that this is a Galerkin approximation) and let $u^{\varepsilon,\mathsf{PG}}_H$ denote the solution to the non-intrusive MsFEM (<ref>) (which is equivalent to the Petrov-Galerkin MsFEM (<ref>), since we do not apply the oversampling technique).

We first study the error $u^{\varepsilon,\mathsf{G}}_H - u^{\varepsilon,\mathsf{PG}}_H$ when $\varepsilon\to0$. In this case, we do not need a rate of convergence in $H$ and we shall relax the condition $f \in L^2(\Omega)$ to the condition $f \in H^{-1}(\Omega)$. Then the linear form $F$ in (<ref>) has to be redefined. Given $f \in H^{-1}(\Omega)$, there exist $f_0,f_1,\dots,f_d \in L^2(\Omega)$ such that
\displaystyle
\sum_{K\in\mathcal{T}_H} \left( \int_K f_0 \, v +
\sum_{\beta=1}^d \int_K f_\beta \, \partial_\beta v
\right),
which is in fact well-defined for any $v \in H^1(\mathcal{T}_H)$ and thus in particular on $V_H$, the underlying affine space for the MsFEM, and the multiscale space $V_H^\varepsilon$.

We consider in Lemma <ref> a sequence of diffusion tensors $A^\varepsilon$ that $H$-converges to a constant diffusion tensor. This means that $u^\varepsilon$ converges weakly in $H^1(\Omega)$ as $\varepsilon\to0$ towards a function $u^\star \in H^1_0(\Omega)$, solution to the homogenized problem (<ref>), and $A^\varepsilon\nabla u^\star \rightharpoonup A^\star \nabla u^\star$ weakly in $L^2(\Omega)$.

Consider the MsFEM-lin or the MsFEM-CR without oversampling. Suppose that $(A^\varepsilon)_{\varepsilon > 0}$ is a sequence of matrices satisfying (<ref>) that $H$-converges to a constant matrix. Let $f \in H^{-1}(\Omega)$. Then $\left\lVert u^{\varepsilon,\mathsf{G}}_H - u^{\varepsilon,\mathsf{PG}}_H \right\rVert_{H^1(\mathcal{T}_H)} \to 0$ as $\varepsilon\to0$.

A rate of convergence can be obtained under some additional structural assumptions on $A^\varepsilon$; see Lemma <ref>.

We need a few auxiliary results to establish Lemma <ref>. The first result below concerns the convergence of the numerical correctors as $\varepsilon\to0$.

Suppose that $A^\varepsilon$ $H$-converges to a constant homogenized tensor $A^\star$. Consider the MsFEM-lin or the MsFEM-CR without oversampling. Then, for all $K\in\mathcal{T}_H$ and all $1\leq\alpha\leq d$, we have $\VK{\alpha} \rightharpoonup 0$ weakly in $H^1(K)$ as $\varepsilon\to0$.

We introduce for each $\alpha=1,\dots,d$ the function $\tau^{\varepsilon,\alpha} = x^\alpha + \VK{\alpha}$. Then (<ref>) implies the PDE
-{\operatorname{div}(A^\varepsilon \nabla \tau^{\varepsilon,\alpha})}
$ in $K$.
The homogenized limit of this problem is
-{\operatorname{div}(A^\star \nabla \tau^{\star,\alpha})}
$ in $K$.
For the MsFEM-lin, the boundary conditions of the local problems impose $\tau^{\star,\alpha}=x^\alpha$ on $\partial K$. The boundary conditions associated to the MsFEM-CR are a constant flux $\vec{n}\cdot A^\star \nabla \tau^{\alpha,\star}$ on each face of $K$ and $\displaystyle \int_h \tau^{\star,\alpha} = \int_h x^\alpha$ for all faces $h$ of $K$. In both cases, the homogenized equation has a unique solution, which is easily seen to be $\tau^{\star,\alpha} = x^\alpha$, because $A^\star$ is constant. Therefore, $\tau^{\star,\alpha} \rightharpoonup x^\alpha$ weakly in $H^1(K)$. Subtracting again the function $x^\alpha$, we deduce the desired convergence.

We will also use the following result, which is a straightforward generalization of the extended Poincaré inequality in <cit.>.

Let $W$ be the subspace of $H^1(\mathcal{T}_H)$ defined by

\begin{equation*}
W = \left\{
v \in H^1(\mathcal{T}_H) \, \left\vert \,
\int_h \llbracket v \rrbracket = 0 \right. \text{ for each face } h \text{ of } \mathcal{T}_h, \,
\int_h v = 0 \text{ for each face } h \subset \partial \Omega
\right\}.
\end{equation*}

There exists a constant $C>0$ depending only on $\Omega$ but not on $H$ such that

\begin{equation*}
\forall \, v \in W,
\qquad
\lVert v \rVert_{L^2(\Omega)} \leq C \, \lVert \nabla v \rVert_{L^2(\mathcal{T}_H)}.
\end{equation*}

Note that the multiscale space $V_{H,0}^\varepsilon$ 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^{\varepsilon,\mathsf{G}}_H$ and $u^{\varepsilon,\mathsf{PG}}_H$.

Let $f \in H^{-1}(\Omega)$ and consider the MsFEM-lin or the MsFEM-CR without oversampling. Let $e^\varepsilon_H = u^{\varepsilon,\mathsf{G}}_H - u^{\varepsilon,\mathsf{PG}}_H$. There exists a unique $e_H^\Pone\in V_H$ and a linear combination of the numerical correctors, that we denote by $e_H^{\mathsf{osc}}$, such that $e^\varepsilon_H = e^\Pone_H + e^\mathsf{osc}_H$, and it holds, with the constants $m,M$ from (<ref>) and the constant $C$ from Lemma <ref>,

\begin{align}
a^{\varepsilon,\dif} \left(e^\varepsilon_H, \, e^\varepsilon_H \right)
\label{eq:estim-Gal-PG-Cea}\\
\left\lVert \nabla e^\mathsf{osc}_H \right\rVert_{L^2(K)}
\frac{M}{m} \left\lVert \nabla e^\Pone_H \right\rVert_{L^2(K)}
\quad
\text{for all } K \in \mathcal{T}_H,
\label{eq:estim-gradient-error3}\\
\left\lVert \nabla e^\mathsf{\Pone}_H \right\rVert_{L^2(\mathcal{T}_H)}
\frac{M}{m} \left\lVert \nabla e^\varepsilon_H \right\rVert_{L^2(\mathcal{T}_H)},
\label{eq:estim-gradient-error1}\\
\left\lVert \nabla e^\varepsilon_H \right\rVert_{L^2(\mathcal{T}_H)}
\sqrt{1+C^2} \frac{M^2}{m^3} \lVert F \rVert_{\mathcal{L}(H^1(\mathcal{T}_H))},
\label{eq:estim-gradient-error2}
\end{align}

where $\lVert \cdot \rVert_{\mathcal{L}(H^1(\mathcal{T}_H))}$ is the operator norm on $\mathcal{L}(H^1(\mathcal{T}_H))$.

Since the numerical approximations $u^{\varepsilon,\mathsf{G}}_H$ and $u^{\varepsilon,\mathsf{PG}}_H$ both belong to the multiscale approximation space $V_H^\varepsilon$, it follows that $e^\varepsilon_H \in V_H^\varepsilon$, and we are in a position to use (<ref>): there exists a unique $e^\Pone_H \in V_H$ such that

\begin{equation}