{'fraction_non_alphanumeric': 0.04043345426771052, 'fraction_numerical': 0.012793584113041817, 'mean_word_length': 4.306231003039514, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-intrusive implementation of Multiscale Finite Element Methods: an illustrative example

Multiscale Finite Element Methods (MsFEM) are finite element type approaches dedicated to multiscale problems. They first compute local, oscillatory, problem-dependent basis functions which generate a specific 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.


Introduction
We consider the highly oscillatory diffusion problem − div (A ε ∇u ε ) = f in Ω, u ε = 0 on ∂Ω, in a bounded domain Ω ⊂ R d , where the diffusion coefficient A ε is assumed to oscillate on a typical length scale of size ε much smaller than the diameter of Ω, and to satisfy the usual ellipticity assumptions (see Section 2 below for details). We seek a numerical approximation of (1) by applying a Galerkin approach. It is well-known that standard (say P 1 ) finite element methods yield an approximation of poor accuracy, as a consequence of the highly oscillatory nature of the problem, unless a prohibitively expensive fine mesh is employed. Dedicated multiscale approaches have thus been introduced, which provide a reasonably accurate approximation of (1) for a limited computational cost. Among the many multiscale approaches that have been proposed in the literature, we mention the Heterogeneous Multiscale Method (henceforth abbreviated as HMM) [1], the Local Orthogonal Decomposition method (LOD) [2], and the Multiscale Finite Element Method (MsFEM) on which we focus here (see [7,5,8] for a comprehensive exposition).
Specifically, the MsFEM method is a finite element type method, which consists of two steps: • an "offline" phase, where highly oscillatory, problem-dependent (but right-hand side independent) basis functions are numerically computed as solutions to local problems (that mimick the reference problem on a subdomain); • an "online" phase, where a Galerkin approximation of (1), performed in the finitedimensional space generated by the basis functions computed in the offline phase, is solved.
The MsFEM approach is particularly interesting in multi-query contexts, where (1) is to be repeatedly solved for multiple right-hand sides f (think e.g. of optimization problems, or of time-dependent problems where (1) would be the typical equation to solve to advance from one time step to the next). In this case, the offline phase is only performed once and a significant computational gain is thus achieved. Several MsFEM variants exist, depending on the specific definition of the basis functions. Although their implementation is rightfully considered to be relatively easy (see e.g. [5,10]), they are all definitely intrusive. They indeed require to change the finite element basis set and adjust it to the problem at hand. Our aim in this article is to investigate how these approaches can be adapted (possibly at the price of a marginal loss in their efficiency) so that they become as little intrusive as possible, thereby allowing to use only a legacy, single-scale software (based on standard finite elements) to recover an accurate approximation of (1). We believe that such an endeavor will eventually facilitate the dissemination of MsFEM methods within industrial, non academic codes.
The article is organized as follows. In Section 2, on the example of the diffusion problem (1), we introduce the simplest MsFEM approach (namely the so-called linear MsFEM) in a Galerkin setting (see Algorithm 1), and present its non-intrusive variant, namely Algorithm 2. We also outline the relation between this variant and a Petrov-Galerkin variant of the MsFEM. Some elements of numerical analysis (along with an illustrative numerical result) are provided in Section 3, to estimate the additional error introduced by the non-intrusive implementation of the method. It is well-known that the linear MsFEM variant considered in Sections 2 and 3 is outperformed by several MsFEM variants. For pedagogical purposes, we have deliberately chosen to present our ideas on this simple variant and to collect in Section 4 some concluding remarks on the many possible extensions of our methodology to design non-intrusive implementations of existing approaches.
We postpone a comprehensive presentation of our procedure, with applications to several MsFEM variants, several equations and various boundary conditions, to an upcoming article [4] (see also [3]).

Non-intrusive implementation of MsFEM: a simple case
We provide problem (1) with all the usual assumptions that make it well-posed at the continuous level and amenable to a classical Galerkin approximation. In particular, we assume that A ε ∈ L ∞ (Ω, R d×d ) satisfies the bounds for some M ≥ m > 0 independent of ε, and that f does not oscillate on the microscopic scale ε. Note that no further structural assumptions on A ε are made (in particular, A ε need not be of the form A(·/ε) for a fixed rescaled function A). We respectively denote by the bilinear and linear forms associated to the variational formulation of (1). We seek a numerical approximation of (1) by applying an MsFEM type Galerkin approach. To this end, we introduce a conformal simplicial mesh T H of Ω (i.e., made of triangles if d = 2, tetrahedra if d = 3) and denote by V H the usual conformal P 1 approximation space on

The Multiscale Finite Element Method
The MsFEM is a Galerkin approximation that adapts the finite-dimensional approximation space in order to obtain a satisfactory accuracy even on a coarse mesh. Let x 1 , . . . , x Nv be the interior vertices (i.e., the N v vertices that do not lie on ∂Ω) of T H , and let φ P1 i be the unique element of V H such that φ P1 i (x j ) = δ i,j for all 1 ≤ j ≤ N v . For any 1 ≤ i ≤ N v , we define the multiscale basis function φ ε i ∈ H 1 0 (Ω), which is supported by the exact same mesh elements as φ P1 i , by The multiscale approximation space is defined as and it has the same dimension as V H . The MsFEM approach then consists in computing the approximation u ε H ∈ V ε H defined by the problem Since the space V ε H is problem-dependent, we can hope (and this is indeed the case) the approximation (5) to capture the exact solution much better than a P 1 approximation on the same mesh (even on a mesh of size H of the order of ε).
Remark 1 (Offline vs. online phase). In practice, the local problems (4) need to be approximated, for instance using a standard P 1 approximation on a fine mesh (of K ∈ T H ) of mesh size h ≤ ε which resolves the oscillations of A ε . We omit here this additional discretization. All what follows can readily be extended to the case when only some numerical approximation of φ ε i is available.
The computation of {φ ε i } 1≤i≤Nv is called the offline phase of the MsFEM. All the problems (4) on different mesh elements are independent of each other, and can thus be solved in parallel. On the other hand, the term online phase is used for the resolution of the global problem (5), where the number of degrees of freedom is the same as in a standard P 1 approximation on V H .
Algorithm 1 MsFEM approach for Problem (1) (see comments in the text) Solve for φ ε i in (4) 4: end for 5: Set A ε := 0 and F ε := 0 6: for all K ∈ T H do 7: end for 12: end for 14: end for Implementing Algorithm 1 in an industrial code is challenging. Indeed, the practical implementation of any finite element method relies on (i) the construction of a mesh, (ii) the construction of the linear system associated to the discrete variational formulation and (iii) the resolution of the linear system. An efficient implementation of the second step heavily relies on the choice of the discretization space. Regarding the construction of the linear system (performed in line 10 of Algorithm 1), it is by no means obvious to adapt existing finite element codes based on generic spaces like V H to a different, problem-dependent choice of space such as V ε H . No analytic expressions for the basis functions φ ε i are available (and thus a fine mesh should be used to approximate them), the computation of a ε K φ ε i , φ ε j and F K φ ε 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 in the coarse mesh, the assembly of the global stiffness matrix A ε j,i 1≤i,j≤Nv should be manually performed, . . . To alleviate these obstacles, we shall next introduce a way of implementing MsFEM that capitalizes on an existing code for solving (1) by a P 1 approximation on T H in the case of slowlyvarying diffusion coefficients.
To the best of our knowledge, the question of how to make MsFEM approaches less intrusive has not been studied in the literature, and this work (that will be complemented in the future works [4,3]) is a first step in that direction. On the other hand, for some other multiscale approaches (including HMM and LOD), this question has been (at least partially) addressed. By construction, HMM methods are less invasive, since they primarily aim at approximating u ε on the coarse scales. The first step of these methods somewhat consists in building an effective, slowly-varying diffusion coefficient (which plays the role of the matrix A introduced in (11) below), which can next be used in any single-scale solver. The LOD approach, which, similarly to MsFEM, aims at approximating u ε on the coarse and the fine scales, is also invasive in general, since it also introduces adapted basis functions. As shown in [6], the LOD can be recast as a P 1 discretization of an appropriate single-scale problem, an observation which opens the way to non-intrusive implementations. Note however that, despite its current intrusiveness, MsFEM has its own advantages over other multiscale approaches: it directly aims at approximating the oscillatory solution u ε (including its fine scale details), and makes use of fully localized basis functions to do so.

Equivalent problem on the macroscopic scale
Our starting point for reducing intrusiveness in the above MsFEM implementation is the following key observation. On any K ∈ T H , by linearity of the definition (4) of φ ε i in terms of φ P1 i , and because the finite element space V H consists of functions that are piecewise affine (that is to say, ∇φ P1 i is constant in each K), we have the expansion for any basis function φ ε i of V ε H . Here, χ ε,α K ∈ H 1 0 (K) is defined as the solution to the local problem where e α denotes the α-th canonical unit vector of R d . Considering indeed the right-hand side of (6), we compute that its gradient in K is The right-hand side of (6) therefore satisfies the PDE in (4). It also satisfies the boundary conditions in (4), in view of the boundary conditions in (7). Since the solution to (4) is unique, we deduce the identity (6).
For each K ∈ T H and each 1 ≤ α ≤ d, we extend χ ε,α K by 0 outside K, thereby obtaining a function χ ε,α K ∈ H 1 0 (Ω). We then deduce from (6) that there is a one-to-one correspondence between functions in V ε H and functions in V H . More precisely, for any and conversely, for We now consider the linear system associated to (5) and insert therein the expansion (8) for the multiscale basis functions. The solution to (5) with (8) and using that ∇φ P1 i is piecewise constant, we write Inserting this relation in the definition of the matrix elements A ε j,i , we obtain where x α = x · e α is the α-th coordinate. We now define the piecewise constant matrixvalued field A ∈ P 0 (T H , R d×d ) by where |K| denotes the area or volume of the mesh element K. Using (7) and the bounds (2) satisfied by A ε , it is easy to show (see [4]) that A satisfies the following uniform lower and upper bounds: for any ξ and η in R d , m |ξ| 2 ≤ ξ · A(x)ξ and η · A(x)ξ ≤ M 1 + M m |ξ| |η| a.e. in Ω.
Motivated by (10) and (11), we introduce the coarse-scale problem and its P 1 Galerkin discretization: find u H ∈ V H such that where the linear form F is defined by (3) and the bilinear form a A is defined by

Problem (14) equivalently writes
with We then deduce from (10) that where we recall that the piecewise constant matrix A given by (11), and therefore the stiffness matrix A P1 , depends on the fine scale oscillations of A ε . The above calculations yield the following result.
Lemma 2. The stiffness matrix A ε in (9) of the MsFEM problem (5) is identical to the stiffness matrix A P1 in (15) of the P 1 problem (14).
We note that the right-hand vector F ε in (9) is in general different from the right-hand side vector F P1 in (15), since we integrate f against highly oscillatory basis functions in the former problem and against P 1 basis functions in the latter. The solutions U ε and U P1 to (9) and (15), respectively, are thus a priori different.
The above observations suggest to use the identity (17) of the stiffness matrices to replace the MsFEM discrete problem (9) by the discrete problem (15) stemming from the P 1 Galerkin approximation of the single scale problem (13), which itself may be easily implemented in legacy codes. This results in our non-intrusive MsFEM strategy presented in Algorithm 2. We can distinguish there the offline phase (in lines 1-14) and the online phase (in lines 15-20). The few lines that differ between Algorithm 2 and a standard P 1 algorithm are highlighted in blue. The lines in black are already present in standard codes and are written in Algorithm 2 for the sake of completeness.
The superiority of Algorithm 2 over the classical MsFEM Algorithm 1 is that the global problem of the online phase (including its right-hand side) can be completely constructed and solved using a pre-existing P 1 PDE solver. The only requirements in the legacy code are the ability to provide piecewise constant diffusion coefficients to the solver  N (n, K) the global index of the vertex of K ∈ T H that has local index 1 ≤ n ≤ d + 1 in K 2: for all K ∈ T H do 3: Solve for χ ε,α K defined by (7) 5: end for 6: Compute A| K defined by (11) 7: end for end for 15: The requirement for these fine-scale solvers is that they have access to the coarse mesh T H of the global solver and that they can evaluate u H anywhere in Ω (which is useful in line 20). Note also that the fine-scale problem (7) is only indexed by the coarse mesh element K, in contrast to the fine scale problem (4), which is indexed both by the coarse mesh element K and the vertex index i. In the latter case, one has to know, for each element K, the global number of the element vertices, a piece of information which may be difficult to access to in a legacy code. In the former case, this correspondence is not needed to compute A and u ε H , which are entirely defined element-wise.
Remark 3 (Link with homogenization theory). The above non-intrusive implementation involves quantities which are of course reminiscent of standard quantities introduced in homogenization. The problem (7), the effective diffusion matrix (11), the single-scale problem (13) and the expansion (8) ressemble the corrector problem, the homogenized matrix, the homogenized problem and the two-scale expansion of the oscillatory solution, respectively. In the same spirit, in the case when A ε is the rescaling of some periodic matrix, it can be shown (see [4]) that the effective diffusion matrix A converges to the homogenized matrix when ε → 0.
As shown by (17), the matrix A P1 in (15) is identical to the matrix A ε in (9). However, in general, the right-hand sides F ε in (9) and F P1 in (15) are different. This motivates the introduction of the following Petrov-Galerkin variant of the MsFEM: find u ε,P G Note that, in contrast to (5), we take the test functions in the P 1 space V H rather than j the stiffness matrix of the resulting linear system. Of course, the right-hand side vector of this linear system is equal to F P1 defined by (16).
Lemma 4. The stiffness matrices A ε and A ε,P G of the problems (5) and (18) are identical (and thus identical to A P1 ). The right-hand side vector of (18) is identical to the right-hand side vector F P1 of the problem (14). The solution u ε,P G H to the Petrov-Galerkin MsFEM (18) thus coincides with the function u ε H computed by Algorithm 2. Proof. To prove this lemma, it is enough to show that the stiffness matrices of (5) and (18) are equal. Using an integration by parts, we compute where n is the unit outward normal vector on ∂K. The two terms above vanish in view of (4). This implies the identity of the stiffness matrices of (5) and (18) and the well-posedness of the Petrov-Galerkin approximation (18).
To summarize, the above procedure to go from the MsFEM problem (5) to its nonintrusive implementation described in Algorithm 2 is based on the following steps: 1. we use the linearity of the problem and the fact that gradients of P 1 basis functions are constant in each mesh element to establish the identity (6); 2. we can then recast the MsFEM stiffness matrix A ε as the stiffness matrix A P1 of the P 1 discretization of an appropriate problem; 3. we approximate the right-hand side F ε of the MsFEM problem by a right-hand side F P1 which can be computed in a manner consistent with a P 1 discretization; 4. we postprocess the P 1 solution to obtain an approximation of the reference solution.
We nowhere use in these steps that the discrete problem we are actually solving (here (15)) in fact corresponds to some discretization (here, a Petrov-Galerkin discretization) of the reference problem. This correspondence is useful to estimate the error introduced by the non-intrusive implementation (a task we perform in Section 3 below), but it does not necessarily hold for other variants of MsFEM (e.g. the oversampling variant briefly mentioned in Section 4 below) and it is not required to put the non-intrusive approach in action. We expect the computational cost of Algorithm 2 to be smaller than that of Algorithm 1, since F P1 is cheaper to compute than F ε (we do not need to use a quadrature rule on the fine mesh). On the other hand, Algorithm 2 (and non-intrusive implementations in general) may introduce additional numerical errors. We estimate these in Section 3.
(provided by Algorithm 2). Assume that f ∈ L 2 (Ω). There exists a constant C independent of ε, H and f such that The classical error estimate for the Galerkin MsFEM approach (5) is obtained in the literature under the assumption that A ε is actually the rescaling of a given periodic matrix. The bound in the estimate reads C (H + √ ε+ ε/H) for some C independent of ε and H (see e.g. [5]). Lemma 6 and a triangle inequality show that the same convergence rate holds for the Petrov-Galerkin MsFEM approach (18), of course under the same periodicity assumption.
Proof. Let e ε H = u ε,G H − u ε,P G H . Since the numerical approximations u ε,G H and u ε,P G Upon testing against v = e osc edges of the coarse mesh elements, since the exact solution oscillates along the edges while the numerical approximation does not. To overcome this drawback, several alternative definitions of the multiscale basis functions have been proposed, leading to different MsFEM variants, including the oversampling variant, introduced in [7] and nowadays considered to be a reference MsFEM variant. We will investigate in [4,3] how the nonintrusive implementation procedure presented here can be extended to that case. We note that oversampling multiscale basis functions are again defined in terms of φ P1 i in a linear manner, which should allow to obtain a formula analogous to (6). The stiffness matrices of the Galerkin and the Petrov-Galerkin approximations are in general different (in contrast to the case studied here, see Lemma 4), but formal arguments show that the two of them can be expressed as the stiffness matrix obtained using a standard P 1 approximation of the single-scale problem (13), for an appropriately defined piecewise constant effective matrix A. We refer to [4,3] for definite conclusions.
To outline the versality of the above procedure leading to a non-intrusive implementation, we eventually make a few remarks on the advection-diffusion problem − div(A ε ∇u ε ) + b · ∇u ε = f in Ω, u ε = 0 on ∂Ω.
The MsFEM basis functions may be defined by (4) or by a similar equation including the advection term (see e.g. [9,3]). Oversampling may also be used. In all cases, for the same reasons as above, an identity of the type (6) should hold, which allows to express the stiffness matrix of the approach as the stiffness matrix obtained using a standard P 1 finite element approximation of a single-scale problem containing diffusion and advection terms, with appropriate definitions of a diffusion matrix A and an advection field b that are both piecewise constant. We again refer to [4,3] for details.