Optimized Oversampling

## Abstract

Abstract In this work we address the Multiscale Spectral Generalized Finite Element Method (MS-GFEM) developed in [I. Babuška and R. Lipton, Multiscale Modeling and Simulation 9 (2011), pp. 373–406]. We outline the numerical implementation of this method and present simulations that demonstrate contrast independent exponential convergence of MS-GFEM solutions. We introduce strategies to reduce the computational cost of generating the optimal oversampled local approximating spaces used here. These strategies retain accuracy while reducing the computational work necessary to generate local bases. Motivated by oversampling we develop a nearly optimal local basis based on a partition of unity on the boundary and the associated A-harmonic extensions.

## 1 Introduction

Modern science and technology demand increasingly efficient methods for solution of multiscale heterogeneous problems.
These problems arise in structural analysis for aerospace and infrastructure or naturally in the study of biological or geological structures. The primary computational challenge in all of these problems is due to the extreme degrees of freedom associated with parameterizing material heterogeneity. Multiscale numerical methods provide one way to solve this problem by invoking local and independent computations that enable a global solve involving a drastically reduced number of degrees of freedom.

A typical and important example is local stress analysis inside an epoxy fiber reinforced composite structure. Here the fibers have a diameter on the order of 10 microns and are densely distributed within an epoxy filler. A critical numerical analysis of such a problem is only possible on massively parallel computers. Hence a special and very parallelizable numerical method is needed. This paper addresses such a method.
In this computational investigation we work with the generalized finite element method (GFEM) introduced in [] and expanded on in [], [], and []. This approach is a partition of unity method (PUM) [], which utilizes the results of many independent and local computations carried out across the computational domain. Essentially this is a domain decomposition and the computational domain is partitioned into an overlapping collection of preselected subsets $\omega_{i}$, $i=1,2,...m$. Finite-dimensional approximation spaces $V_{\omega_{i}}$ are constructed over each subset using local information. In this paper the finite dimensional space of local solutions to the problem is first solved over a slightly larger domain $\omega^{*}_{i}$ and generates the finite dimensional space $V_{\omega^{*}_{i}}$. In the second step we restrict the solutions to the smaller set (patch) $\omega_{i}\subset\subset\omega_{i}^{*}$. In this way the local space $V_{\omega_{i}}$ is now given by the restriction of elements in $V_{\omega^{*}_{i}}$ to $\omega_{i}$. This way of generating the local basis is called oversampling. As is the case in all GFEM schemes, each subspace is computed independently and the full global solution is obtained by solving a global (macro) system which can be orders of magnitude smaller than the system corresponding to a direct application of the finite element method (FEM) to the full structure.

In this paper we carry out numerical investigations using optimal local bases introduced in [] and []. These bases are proved to have the best approximation properties over all local bases constructed by oversampling A-harmonic functions [], []. Here an A-harmonic function $v$ satisfies

$$ ${\rm div}(A(x)\nabla v)=0\hbox{ on $\omega^{\ast}$},$ (1) $$

where $A(x)$ is a given $L^{\infty}$ coefficient satisfying the usual coercivity and boundedness conditions.
The N-dimensional optimal local basis is shown to be the span of the first $N$ eigenfunctions associated with the singular values of the restriction operator used in the oversampling. These local bases have been theoretically proved to have approximation error that asymptotically approaches an exponential decay with the dimension of the approximation space ([] and []). The convergence result applies to any elastic or conductivity tensor field with $L^{\infty}$ coefficients satisfying coercivity and boundedness conditions.
The associated scheme introduced in [] and [] is the Multiscale-Spectral GFEM (MS-GFEM) and has global approximation error that converges exponentially asymptotically in the energy norm.

As discussed above the ultimate motivation behind the MS-GFEM and its exponential convergence rate is its potential for use in large parallel implementations. The method has two essential parts. A. Parallel construction of the local non-polynomial shape functions leading to an exponential approximation of the actual solution in the energy norm. Here we note that each local basis is defined over separate patches and can be computed independently on separate processors using local memory. B. Efficient parallel solution of the global system of linear algebraic equations. This part is made efficient noting that
global basis vectors associated with each patch are only influenced by immediate neighboring patches and a Schwartz alternating method [], [] can be applied as a preconditioner for an exponentially convergent iterative scheme. The current paper concentrates on A while B is discussed in section [8](#S8) and its details will be presented in a forthcoming paper. Another very expensive aspect of large implementations is mesh generation. Here the GFEM allows each patch to be meshed independently.

In this work we provide a numerical implementation of MS-GFEM and our aim is threefold. We focus on a simple two patch domain decomposition to illustrate 1) contrast independent convergence of the local basis and global basis, 2) new methods for inexpensively generating local basis functions, and 3) identification of alternate local bases with good approximation properties. The computational domain is composed of heat conducting particles included within a connected second heat conducting phase. The connected phase is often referred to as the matrix phase. Computations are carried out here for a large range of contrasts between particle and matrix conductivities. For inclusions that do not touch and have smooth boundary our computations show that the MS-GFEM approximations converge exponentially and independently of the contrast between phases in the preasymptotic regime. This is consistent with the simulations given in [] for a medium with separated holes. This is also consistent with the more recent numerical investigation carried out for scalar and elastic problems using the optimal local basis functions in []. The theoretical proof of contrast independent exponential convergence for MS-GFEM applied to heterogeneous elastic and conduction problems is given in a forthcoming publication.

The primary numerical work in generating the best local basis is not in the solution of the eigenvalue problem but in the numerical generation of the local A-harmonic subspace used to construct the optimal local basis.
In section [6](#S6) we introduce a method to reduce the computational cost of generating optimal local basis functions. We investigate the inexpensive construction of A-harmonic subspaces for numerically generating the optimal local basis. The optimal $N$-dimensional approximation space is given by the span of the optimal local basis given by the first $N$ eigenfunctions associated with the largest $N$ singular values of the restriction operator acting on A-harmonic functions. These eigenfunctions are seen to define the Kolmogorov $n$-width [] and are referred to as the $n$-width eigenfunctions.

We explore the use of discrete A-harmonic extensions of hat functions defined on the boundary of $\omega^{\ast}$ necessary to generate these eigenfunctions numerically. Here we are motivated by the fact that longer wavelength boundary data applied to $\partial\omega^{*}$ penetrates further into the subdomain $\omega$ than short wavelength “oscillatory” boundary data [] and []. In light of the exponential decay of the $n$-width eigenfunctions it becomes clear that the space spanned by the lower $n$-width modes are generated by linear combinations of A-harmonic extensions of low frequency or slowly changing boundary data. A-harmonic extensions of more oscillatory boundary data decay faster into the domain and thus contribute only to the higher $n$-width modes. Thus we use the A-harmonic extension of a collection of hat functions with a relatively large support set that taken together form a partition of unity on the boundary. Here the support of the boundary data given by hat functions can extend over tens or hundreds of boundary nodes associated with the discrete A-harmonic extension of the boundary data. Our numerical experiments reveal the size of support of the boundary hat function that can be used for accurate yet inexpensive computation of the optimal local basis. It turns out that the number of independent boundary hat functions can be roughly the same as the number of discrete $n$-width basis functions needed for a good approximation. The numerical experiments presented in Section [6](#S6) show that this method retains accuracy while effectively reducing the computational work necessary to generate local bases.
Motivated by our results in section [6](#S6) we use the A-harmonic extensions of the independent boundary hat functions restricted to the subdomain $\omega$ as elements of a local basis. In Section [7](#S7) we find that the relative global error of solution using these as local shape functions in an oversampled-GFEM scheme developed here while not as small still remains comparable to that of our MS-GFEM. This can be deduced from the numerical experiments of Section [6](#S6) that show that for a given number of local basis functions the span of this type of local basis is “nearly the same” as the span of the $n$-widths used in MS-GFEM. In forthcoming work we illustrate this reason in rigorous mathematical terms for scalar problems in terms of Tchebycheff approximations.

Related recent work is motivated by the theory of low rank approximations based on randomized SVD [] or rSVD. The work [] outlines a reduction of the dimension of the space used to generate $n$-width eigenfunctions numerically by sampling randomly generated basis functions. Another strategy [] provides a method based on rSVD for choosing boundary data and using their A-harmonic extensions as a sub optimal basis and demonstrating a nearly exponential decay in the expected value of the error. In a recent independent development the numerical investigations [] show that local bases constructed by A-harmonically extending traces of harmonic polynomials deliver exponential decay of the error.

We conclude the introduction noting that there are several other strategies to reduce the computational work for multiscale numerical implementations.
One way to address multiscale computation is to exploit the local periodicity or stochasticity of a microstructure and make use of homogenization theory []. These methods include the original Multiscale FEM [] and the variational multiscale method []. In the absence of local structure one can consider approaches to numerical homogenization for rough coefficients (i.e., $L^{\infty}$ coefficients); however, the lack of local structure naturally degrades the efficiency and convergence rates. Nevertheless the size of the global solve can be reduced and several methods have been proposed for $L^{\infty}$ coefficients that offer significant dimension reduction. There is a huge literature and a variety of methods have been introduced. These methods include upscaling based on harmonic coordinates and elliptic inequalities [], [], elliptic solvers based on H-matrices [, ], explicit solution of local computations through Bayesian numerical homogenization [],
dimension reduction methods based on global changes of coordinates and MS-FEM for upscaling porous media flows [], [], the heterogeneous multiscale methods [], [], [], and an adaptive coarse scale–fine scale projection method []. Additional contemporary methods include numerical homogenization based on the flux norm for $L^{\infty}$ coefficients [], rough polyharmonic splines [], subgrid upscaling methods [] and global Galerkin projection schemes for problems with $L^{\infty}$ coefficients and homogeneous Dirichlet boundary data []. For a coarse mesh of diameter $H$ local bases that deliver order $H$ convergence with $O((log{(1/H)})^{d+1}$ approximation functions are developed in []. For comparison the method presented here and in [], [], show that the coarse mesh can be fixed arbitrarily and for a given relative error $\tau$ one needs $O((log{(1/\tau)})^{d+1}$ $n$-width approximation functions.

## 2 Problem Formulation

### 2.1 Variational Formulation of the Problem

To fix ideas we consider the scalar problem over a bounded domain $\Omega\subset\mathbb{R}^{2}$ with piecewise $C^{1}$-boundary given by

$$ $-\operatorname{div}(A(x)\nabla u(x))=f(x),\text{ }x\in\Omega$ (2) $$

with Neumann boundary conditions prescribed on the boundary $\partial\Omega_{N}\subset\partial\Omega$

$$ $n\cdot A(x)\nabla u(x)=g(x),\text{ }x\in\partial\Omega_{N},$ (3) $$

where $n$ is the unit outer normal vector, and Dirichlet boundary conditions on $\partial\Omega_{D}\subset\partial\Omega$

$$ $u(x)=h(x)\text{, }x\in\partial\Omega_{D},$ (4) $$

such that $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$ and $\overline{\partial\Omega_{D}}\cup\overline{\partial\Omega_{N}}=\partial\Omega$. Here $A(x)$ is the $2\times 2$ conductivity matrix with rough coefficients $a_{ij}(x)\in L^{\infty}(\Omega)$, and satisfies the standard ellipticity and boundedness conditions

$$ $0\leq\alpha v\cdot v\leq A(x)v\cdot v\leq\beta v\cdot v<\infty,\text{ }\forall x\in\Omega\text{ and }v\in\mathbb{R}^{2}.$ (5) $$

The unique weak solution of ([2](#S2.E2)) belongs to the convex space

$$ $H^{1}_{D}(\Omega)=\{u\in H^{1}(\Omega):u=h\text{ on }\partial\Omega_{D}\},$ (6) $$

and satisfies

$$ $B\left(u,v\right)=F(v)$ (7) $$

for all $v$ in the energy space

$$ $H^{1}_{0D}=\{v\in H^{1}(\Omega):v=0\text{ on }\partial\Omega_{D}\},$ (8) $$

where

$$ $B(u,v)=\int_{\Omega}A(x)\nabla u\cdot\nabla v\,dx\text{ and }F(v)=\int_{\Omega}fv\,dx+\int_{\partial\Omega_{N}}gv\,ds.$ $$

The energy norm is given by $\lVert u\rVert_{\mathcal{E}(\Omega)}=(B(u,u))^{1/2}$.

### 2.2 MS-GFEM

In this paper the numerical solution to ([7](#S2.E7)) is computed using MS-GFEM. We outline the a priori convergence estimates and motivation for the method. To illustrate the ideas we consider a computational domain $\Omega\subset\mathbb{R}^{2}$ containing smooth inclusions separated by a prescribed minimum distance. The MS-GFEM is a version of GFEM [] and is a domain decomposition given by a partition of unity of overlapping subdomains. Over each subdomain a local approximation space of shape functions is constructed that is adapted to the coefficient of the PDE restricted to the subdomain. In MS-GFEM the local shape functions are characterized by an optimal local approximation space obtained from oversampling. The local computations associated with construction of local approximation spaces are independent and can be performed in parallel. The resulting global stiffness matrix can be several orders of magnitude smaller than the stiffness matrix obtained by applying FEM directly ([]).

We begin by describing the partition of unity.
Let $\{\omega_{i}\}_{i=1}^{m}$ be a collection of open sets covering the domain $\Omega$ such that $\cup_{i=1}^{N}\omega_{i}=\Omega$ and let $\phi_{i}\in C^{1}(\omega_{i})$, $i=1,\ldots,m$, be a partition of unity subordinate to the open covering. Here the maximum number of sets $\omega_{i}$ that any point $x$ in $\Omega$ can belong to is at most $\kappa$. The partition of unity functions satisfy the following properties:

$$ $\displaystyle 0\leq\phi_{i}\leq 1,\,\,\,\,\,i=1,\ldots,N,$ $\displaystyle\phi_{i}(x)=0,\,x\in\Omega\setminus\omega_{i},\,\,\,\,\,i=1,\ldots,N,$ $\displaystyle\sum_{i}\phi_{i}(x)=1\,\,\,\,\,\forall x\in\Omega,$ $\displaystyle\max_{x\in\Omega}\lvert\phi_{i}(x)\rvert\leq C_{1},\,\,\,\,\,i=1,\ldots,N,$ $\displaystyle\max_{x\in\Omega}\lvert\nabla\phi_{i}(x)\rvert\leq\frac{C_{2}}{\operatorname{diam}(\omega_{i})},\,\,\,\,\,i=1,\ldots,N,$ $$

where $C_{1}$ and $C_{2}$ are bounded positive constants and $\operatorname{diam}(\omega_{i})$ is the diameter of the set $\omega_{i}$. The partition of unity functions are chosen to be flat-topped so that local approximation spaces are linearly independent and for ensuring good conditioning of the global stiffness matrix ([]).

Now we describe the construction of the optimal local approximation spaces that characterize the MS-GFEM method. These spaces are constructed through oversampling and are spanned by spectral bases given in terms of the eigenspaces of a restriction operator on A-harmonic functions described below, see [] and [].
To construct the local bases first let $\{\omega_{i}^{*}\}_{i=1}^{N}$ be a second collection of open sets such that each $\omega_{i}$ is contained within the larger open set $\omega_{i}^{*}$. For interior domains we require $\operatorname{dist}(\partial\omega^{\ast}_{i},\partial\omega_{i})>0$. For $\omega_{i}$ touching the boundary of $\Omega$, i.e., $\partial\omega_{i}\cap\partial\Omega\not=\emptyset$, we require $\operatorname{dist}(\partial\omega_{i}^{*}\cap\Omega,\partial\omega_{i}\cap\Omega)>0$. In what follows we will refer to the subdomains $\omega_{i}$ and $\omega_{i}^{*}$ as patches. The local approximations are constructed from a local affine space spanned by a particular solution (see equation ([11](#S2.E11)) below) and a local approximation space $V^{m_{i}}_{\omega_{i}}$. Here the dimension of the local approximation space is denoted by $m_{i}$.
On each subdomain $\omega_{i}$ the local approximation space is created by first constructing a finite dimensional space of functions on $\omega_{i}^{*}$ denoted by $V^{m_{i}}_{\omega_{i}^{*}}$ and then restricting them to $\omega_{i}$ to form the approximation space $V^{m_{i}}_{\omega_{i}}$. In MS-GFEM the local shape functions describing a basis for $V^{m_{i}}_{\omega_{i}}$ are the eigenfunctions associated with the singular values of the restriction operator (see section [2.3](#S2.SS3)) acting on A-harmonic functions, i.e., functions $\xi$ that satisfy

$$ $-\operatorname{div}(A(x)\nabla\xi(x))=0\text{ on }\omega_{i}^{*}.$ (9) $$

The space of A-harmonic functions defined on $\omega_{i}^{*}$ is written $H_{A}(\omega_{i}^{*})$. For patches that do not share a boundary with $\Omega$ functions are chosen such that they are equivalent up to a constant. The associated quotient space is written $H_{A}(\omega_{i}^{*})/\mathbb{R}$. For boundary patches sharing Neumann conditions on $\partial\omega_{i}^{*}\cap\partial\Omega_{N}$, the local functions are taken to satisfy homogeneous Neumann boundary conditions on $\partial\omega_{i}^{*}\cap\partial\Omega_{N}$. For boundary patches sharing non-homogeneous Dirichlet conditions on $\partial\omega_{i}^{*}\cap\partial\Omega_{D}$, the local functions are taken to satisfy homogeneous Dirichlet boundary conditions on $\partial\omega_{i}^{*}\cap\partial\Omega_{D}$. Here all boundary patches either share a Dirichlet boundary or Neumann boundary but not both. The construction of A-harmonic approximation spaces for boundary patches are described in detail in [] and []. For interior and boundary patches with Neumann data the local approximation space is augmented with the constant functions. In all cases the dimension of the local approximation space over $\omega_{i}$ is denoted by $m_{i}$.
The global approximation space is constructed from the local approximation spaces and is defined by

$$ $V^{N}=\left\{\sum_{i=i}^{N}\phi_{i}\xi_{i}:\xi_{i}\in V_{\omega_{i}}^{m_{i}}\right\},$ (10) $$

and one verifies as in [] that $V^{N}$ is a subspace of the energy space $H_{0D}^{1}(\Omega)$.
The boundary data ([3](#S2.E3)), ([4](#S2.E4)), and the right hand side of ([2](#S2.E2)) are satisfied by local particular solutions. These solutions $\chi_{i}\in H^{1}(\omega_{i}^{*})$ are defined by

$$ $-\operatorname{div}(A(x)\nabla\chi_{i}(x))=f(x),\text{ }x\in\omega_{i}^{*},$ (11) $$

with Dirichlet data $\chi_{i}=0$ on $\partial\omega_{i}^{*}$ for interior patches. For boundary patches sharing non-homogeneous Neumann data the particular solutions $\chi_{i}$ also satisfy the boundary data $\eqref{eq:globalneumann}$ and $\chi_{i}=0$ on $\partial\omega_{i}^{*}\cap\Omega$. For boundary patches sharing non-homogeneous Dirichlet data the particular solutions $\chi_{i}=\chi_{i}^{R}+\chi_{i}^{D}$ where $\chi_{i}^{R}$ solves ([11](#S2.E11)) with $\chi_{i}^{R}=0$ on $\partial\omega_{i}^{*}$ and $\chi_{i}^{D}$ satisfy the boundary data $\eqref{eq:globaldirichlet}$ on $\partial\omega_{i}^{*}\cap\partial\Omega$ and homogeneous Neumann data on $\partial\omega^{*}_{i}\cap\Omega$ and ([11](#S2.E11)) with $f=0$. The global particular solution $u^{F}$ is then defined by pasting together the local particular solutions, i.e.,

$$ $u^{F}=\sum_{i=1}^{N}\phi_{i}\chi_{i}.$ (12) $$

The finite dimensional approximate solution of ([7](#S2.E7)) is posed over the convex space $K^{N}=V^{N}+u^{F}$. The problem becomes a variational inequality over the convex space $K^{N}$. We seek a solution $u^{G}\in V^{N}$ to the following problem for all $v\in V^{N}$

$$ $B(u^{G},v)=F(v)-B(u^{F},v).$ (13) $$

The MS-GFEM approximate solution to ([7](#S2.E7)) is given by $u_{0}=u^{G}+u^{F}$ and

$$ $u_{0}=\operatorname*{arg\,min}\left\{\|u-w\|_{{\mathcal{E}}(\Omega)}\,:\,w\in V^{N}\right\}$ (14) $$

follows from Galerkin orthogonality.
The existence of a unique solution $u_{0}\in K^{N}$ follows from the standard theory of variational inequalities see, e.g., []. Given any function $u$ in $H^{1}_{D}(\Omega)$ and trial field $u_{T}=\sum_{i=1}^{N}\phi_{i}\xi_{i}+u^{F}$ a nontrivial calculation following the proofs of Theorems 3.2 and 3.3 of [] shows that if the local error on every $\omega_{i}$ satisfies

$$ $\|u-(\xi_{i}+\chi_{i})\|_{{\mathcal{E}}{(\omega_{i}})}\leq\epsilon_{i}\|u\|_{{\mathcal{E}}{(\omega_{i}^{*}})}\leq\epsilon_{i}\|u\|_{\mathcal{E}(\Omega)},$ (15) $$

for some $\xi_{i}$ in $V_{\omega_{i}}^{m_{i}}$ and local particular solution $\chi_{i}$,
then the global (i.e., total) approximation error of the trial field is bounded by

$$ $\|u-u_{T}\|_{{\mathcal{E}}(\Omega)}\leq(\kappa)^{1/2}C\left(\sum_{i}^{N}\epsilon_{i}^{2}\right)^{1/2}\|u\|_{{\mathcal{E}}{(\Omega})},$ (16) $$

where $C>0$ is independent of $u$ in $H^{1}_{D}(\Omega)$. Now from Galerkin minimality ([14](#S2.E14)) together with ([16](#S2.E16)) we see that the global error of the Galerkin solution $u_{0}$ is controlled by the local errors, this is the hallmark of the GFEM method. The calculations behind ([15](#S2.E15)) and ([16](#S2.E16)) for MS-GFEM are provided in the appendix for completeness. In the next section we show how to construct optimal finite dimensional local approximation spaces to produce local approximations $\xi_{i}$ that satisfy ([15](#S2.E15)) when the local particular solution $\chi_{i}$ is given.

### 2.3 Optimal Local Approximation Spaces and MS-GFEM

We now describe the optimal local approximation spaces developed in [] and []. We begin by restricting attention to a single subdomain $\omega$ and omit the subscript. The oversampling problem is defined over the larger subdomain $\omega^{*}\supset\omega$. The space of A-harmonic functions $H_{A}(\omega^{*})$ is the space defined as

$$ $H_{A}(\omega^{*})=\{v\in H^{1}(\omega^{*}):B(v,w)=0\text{ for all }w\in H^{1}_{0}(\omega^{*})\}.$ (17) $$

Note we can add a constant function to any function in $H_{A}(\omega^{*})$ and still satisfy ([17](#S2.E17)). Hence we consider the quotient space made up of all equivalence classes of functions in $H_{A}(\omega^{*})$ that are the same up to a constant and denote this as
$H_{A}(\omega^{*})/\mathbb{R}$. Elements $u$ of $H_{A}(\omega^{*})/\mathbb{R}$ restricted to $\omega$ can be approximated using an optimal local spectral basis [] and []. The restriction operator $P$ is defined by $Pu(x)=u(x)$ for $x$ in $\omega$. The restriction is a compact map from
$H_{A}(\omega^{*})/\mathbb{R}$ into $H_{A}(\omega)/\mathbb{R}$ and the optimal local basis on $\omega$ is given by the span of the eigenfunctions $\{\xi_{j}\}_{i=1}^{\infty}$ associated with the singular values $\{\lambda_{j}\}_{j=1}^{\infty}$ of the restriction operator

$$ $P^{*}P\xi_{j}=\lambda_{j}\xi_{j},$ (18) $$

see [].
This space provides local shape functions and we see that MS-GFEM is a spectral approximation method. The optimal local approximation space is given by the span of the spectral basis and denoted by $V^{m}_{\omega}={\rm span}\{\xi_{1},\xi_{2},\ldots,\xi_{m}\}$.
The optimality of this oversampling space is seen from the theory of Kolmogorov $n$-widths, []. The best accuracy of approximation over all $n$-dimensional subspaces $S(n)\subset H_{A}(\omega)/\mathbb{R}$ of a function $u\in H_{A}(\omega^{*})/\mathbb{R}$ is measured by the Kolmogorov $n$-width as

$$ $d_{n}(\omega,\omega^{*})=\inf_{S(n)}\sup_{u\in H_{A}(\omega^{*})/\mathbb{R}}\inf_{v\in S(n)}\frac{\lVert Pu-v\rVert_{\mathcal{E}(\omega)}}{\lVert u\rVert_{\mathcal{E}(\omega^{*})}},$ (19) $$

and most importantly $d_{n}(\omega,\omega^{\ast})=\sqrt{\lambda_{n+1}}$. Thus the square root of the $n^{th}$ eigenvalue gives the approximation error incurred when approximating the A-harmonic part of the solution using $V_{\omega}^{m}$.
The optimal subspace achieving the infimum of ([19](#S2.E19)) for $n=m$ is written as $V_{\omega}^{m}$ and is defined as the span of the first $m$ eigenfunctions satisfying ([18](#S2.E18)), see []. These eigenfunctions can be computed explicitly and are characterized by the following theorem.

###### Theorem 1 ( [ 5 ] , Theorem 3.2) .

The optimal approximation space is given by $V^{m}_{\omega}=\operatorname{span}\{\xi_{1},\ldots,$ $\xi_{m}\}$, where $\xi_{j}=P\phi_{j}$ and $\lambda_{j}$ and $\phi_{j}$ are the largest $m$ eigenvalues and corresponding eigenfunctions that satisfy

$$ $(\phi_{j},\delta)_{\mathcal{E}(\omega)}=\lambda_{j}(\phi_{j},\delta)_{\mathcal{E}(\omega^{*})}\text{ }\forall\delta\in H_{A}(\omega^{*}).$ (20) $$

The functions $\{\phi_{j}\}_{j=1}^{\infty}$ form a complete orthonormal set in $H_{A}(\omega^{*})/\mathbb{R}$.

###### Definition 1 .

The optimal local basis $\{\xi_{i}\}_{i=1}^{m}$ is referred to as the “$n$-width” functions or spectral basis and $V^{m}_{\omega}={\rm span}\{\xi_{1},\xi_{2},\ldots,\xi_{m}\}$ as the spectral subspace.

An asymptotically exponential decay of the $n$-widths $\sqrt{\lambda_{m+1}}$ for $n=m$ is proved in [].

###### Theorem 2 .

The accuracy has nearly exponential decay for $m$ sufficiently large, i.e., for any small $\epsilon$ one has

$$ $d_{m}(\omega,\omega^{*})=\sqrt{\lambda_{m+1}}\leq e^{-m^{\frac{1}{1+d}-\epsilon}}.$ (21) $$

It follows that for $u_{A}\in H_{A}(\omega^{*})$ there exists $\xi\in V_{\omega}^{m}$ such that

$$ $\lVert u_{A}-\xi\rVert_{\mathcal{E}(\omega)}\leq\sqrt{\lambda_{m+1}}\lVert u_{A}\rVert_{\mathcal{E}(\omega^{*})}.$ (22) $$

Now for a solution $u$ of ([2](#S2.E2)) in $H^{1}_{D}(\Omega)$ we construct a trial $u_{T}$ with a $\xi_{i}$ on each $\omega_{i}^{*}$ such that ([15](#S2.E15)) holds. Here we will write $\lambda_{m+1}^{i}$ to correspond to domain $\omega_{i}^{*}$. Note first that for any $\chi_{i}$ and element $\xi_{i}\in V^{m_{i}}_{\omega_{i}}$ we have

$$ $u-(\xi_{i}+\chi_{i})=(u-\chi_{i})-\xi_{i},$ (23) $$

where $u-\chi_{i}$ is A-harmonic. Thus from Theorem [2](#Thmtheorem2) we can choose a $\xi_{i}$ such that

$$ $\lVert u-\chi_{i}-\xi_{i}\rVert_{\mathcal{E}(\omega_{i})}\leq\sqrt{\lambda^{i}_{m+1}}\lVert(u-\chi_{i})\rVert_{\mathcal{E}(\omega_{i}^{*})}.$ (24) $$

Now note from ([11](#S2.E11)) that for interior patches

$$ $\displaystyle\|\chi_{i}\|^{2}_{\mathcal{E}(\omega_{i}^{*})}$ $\displaystyle=\int_{\omega_{i}^{*}}A(x)\nabla\chi_{i}\cdot\nabla\chi_{i}\,dx$ (25) $\displaystyle=\int_{\omega_{i}^{*}}f\chi_{i}\,dx=\int_{\omega_{i}^{*}}A(x)\nabla u\cdot\nabla\chi_{i}\,dx\leq\|u\|_{\mathcal{E}(\omega_{i}^{*})}\|\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})},$ $$

so

$$ $\displaystyle\|\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}\leq\|u\|_{\mathcal{E}(\omega_{i}^{*})},$ (26) $$

and a simple calculation shows

$$ $\displaystyle\|u-\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}\leq 2\|u\|_{\mathcal{E}(\omega_{i}^{*})}.$ (27) $$

From this we get

$$ $\lVert u-\chi_{i}-\xi_{i}\rVert_{\mathcal{E}(\omega_{i})}\leq\epsilon_{i}\lVert u\rVert_{\mathcal{E}(\omega_{i}^{*})}.$ (28) $$

with $\epsilon_{i}=2\sqrt{\lambda^{i}_{m+1}}$.
Proceeding in the same way for boundary patches sharing Neumann data we note $u-\chi_{i}$ is A-harmonic and choose $\xi_{i}$ to obtain

$$ $\lVert u-\chi_{i}-\xi_{i}\rVert_{\mathcal{E}(\omega_{i})}\leq\epsilon_{i}\lVert u\rVert_{\mathcal{E}(\omega_{i}^{*})},$ (29) $$

with $\epsilon_{i}=2\sqrt{\lambda^{i}_{m+1}}$.
For Dirichlet data similarly $u-\chi_{i}$ is A-harmonic and we may choose $\xi_{i}$ to obtain

$$ $\lVert u-\chi_{i}-\xi_{i}\rVert_{\mathcal{E}(\omega_{i})}\leq\epsilon_{i}\lVert u\rVert_{\mathcal{E}(\omega^{*}_{i})},$ (30) $$

with $\epsilon_{i}=2\sqrt{5\lambda^{i}_{m+1}}$.
These calculations are provided in the Appendix for completeness.
The global error estimate ([16](#S2.E16)) now follows from ([28](#S2.E28)), ([29](#S2.E29)), and ([30](#S2.E30)).
Therefore it is evident that the a priori estimates showing asymptotic exponential accuracy are obtained when using the optimal local bases in the construction of the global approximation space $V^{N}$. The global error on $\Omega$ is controlled by the local errors and we have an a priori asymptotic exponential decay of the MS-GFEM approximation.
The numerical experiments given here show an immediate exponential decay of the $n$-widths for this class of coefficients, see section [4.6](#S4.SS6). The numerical experiments also show that the global approximation error is exponentially decaying independent of contrast, see sections [4](#S4), [5](#S5).

We summarize noting that the $n$-width eigenfunctions are the optimal oversampled approximation functions and are computed for the particular geometry and microstructure $\Omega$. They carry information about the local fields and allow for a drastic reduction in the dimension of the global stiffness matrix. The total computational complexity of this method is the same as a direct application of the FEM method and is of the same order as the number of inclusions in the domain $\Omega$, see []. But now most all of the work is done in parallel and offline. This results in a global solve given by the inversion of a matrix that is orders of magnitude smaller than the total number of inclusions.

###### Theorem 1 ( [ 5 ] , Theorem 3.2) .

###### Definition 1 .

###### Theorem 2 .

## 3 Computational Method

We outline the implementation of the MS-GFEM algorithm and describe the construction of the global stiffness matrix, particular solutions, and right hand side. Recall that we wish to solve ([7](#S2.E7)) using the optimal local bases pasted together over the computational domain $\Omega\subset\mathbb{R}^{2}$. As data we are given a matrix of material properties $A(x)$ together with Neumann and Dirichlet boundary conditions specified over the boundary of $\Omega$. We begin by first creating a standard finite element mesh on the entire domain $\Omega$ which will be used in all computations. Future implementations will use independent meshes over each patch $\omega_{i}$ enhancing the parallel nature of the local computations.

The steps in the MS-GFEM algorithm are:

- 1.
Define a covering of $\Omega$ by square patches $\{\omega_{i}\}_{i=1}^{N}$ and expand each patch $\{\omega_{i}^{*}\}_{i=1}^{N}$ such that $\omega_{i}\subset\omega_{i}^{*}$ are concentric and the ratio of side edges of $\omega_{i}$ to $\omega_{i}^{*}$ is less than 1.
- 2.
Construct partition of unity functions $\{\phi_{i}\}_{i=1}^{N}$ subordinate to $\{\omega_{i}\}_{i=1}^{N}$. In this work we use the standard bilinear and quadratic shape functions over triangles and squares to construct the partition of unity functions.
- 3.
Solve local particular solutions $\{\chi_{i}\}_{i=1}^{N}$ satisfying ([11](#S2.E11)). The boundary data on interior patches is given by $\chi_{i}=0$ on $\partial\omega_{i}^{*}$. On boundary patches $\chi_{i}$ satisfies Neumann boundary data ([3](#S2.E3)) or Dirichlet data ([4](#S2.E4)) on $\partial\omega_{i}^{*}\cap\partial\Omega$ with $\chi_{i}=0$ on $\partial\omega_{i}^{*}\cap\Omega$.
- 4.
Construct right hand side of ([13](#S2.E13)), denoted as $\mathbf{r}$.
- 5.
Construct finite dimensional local spaces $S^{n_{i}}_{\omega_{i}^{*}}\subset H_{A}(\omega_{i}^{*})$. Here $S^{n_{i}}_{\omega_{i}^{*}}$ is thought of a discrete approximation to $H_{A}(\omega^{*}_{i})$. The dimension of these spaces are $n_{i}$ and the elements of this space are A-harmonic extensions of boundary data given by the piece wise linear hat functions defined on $\partial\omega_{i}^{*}$. Together these piece wise linear functions form a partition of unity on $\partial\omega_{i}^{*}$.
- 6.
Solve the local spectral problems over $S^{n_{i}}_{\omega_{i}^{*}}$ to build the local spectral bases that generate the spaces $V_{\omega_{i}}^{m_{i}}$ with $n_{i}>m_{i}$.
- 7.
Construct the global stiffness matrix $\mathrm{G}$ in ([13](#S2.E13)) with entries given by $B(u,v)$ for $u,v\in V^{N}$.
- 8.
Solve the Galerkin formulation given by $\mathrm{G}\mathbf{x}=\mathbf{r}$.
- 9.
The approximate solution to ([6](#S2.E6)) is given by $u_{0}=u^{G}+u^{F}$.

We elaborate on the elements of this algorithm below and illustrate the computations for a simple partition of unity comprised of two subdomains in section [4](#S4). For details in computing finite dimensional subspaces of $H^{0}_{A}(\omega_{i}^{*})$ and the generation local spectral bases see section [6](#S6).

### 3.1 Construction of Local Solution Spaces S ω i ∗ n i subscript superscript 𝑆 subscript 𝑛 𝑖 superscript subscript 𝜔 𝑖 S^{n_{i}}_{\omega_{i}^{*}} .

The material domain $\Omega$ is covered by square or hexahedral patches $\omega_{1},\ldots,\omega_{N}$ with an expanded covering $\omega_{1}^{*},\ldots,\omega_{N}^{*}$ such that $\omega_{i}\subset\omega_{i}^{*}$.
We use the following notations. Let $\mathcal{H}=\{E_{e}\}_{e=1}^{n_{el}}$ be a finite element mesh of $\Omega$, with elements $E_{e}$, which conforms to the boundaries of each patch $\omega_{i}$ and $\omega_{i}^{*}$, $i=1,\dots,N$. The piecewise bilinear and quadratic shape functions are denoted $H_{n}(\alpha)$ where ${\alpha}=(\alpha_{1},\alpha_{2})$ are coordinates on the reference element and ${x}=(x_{1},x_{2})$ are the Cartesian coordinates on the domain $\Omega$. The map $\alpha_{e}:E_{e}\rightarrow\square$ is the change of coordinates map from the $e$-th element $E_{e}$ to the reference element $\square$. We begin by generating a discrete A-harmonic solution space $S^{n_{i}}_{\omega_{i}^{*}}\subset H_{A}(\omega_{i}^{*})$ of A-harmonic FEM solutions of

$$ $\left\{\begin{aligned} -\operatorname{div}(A(x)\nabla w_{i}^{k}(x))=0\text{ for }x\in\omega_{i}^{*}\\ w_{i}^{k}(x)=h_{i}^{k}(x)\text{ for }x\in\partial\omega_{i}^{*}.\end{aligned}\right.$ (31) $$

Here $h_{i}^{k}(x)$ is the $k$-th hat function defined by $h_{i}^{k}=1$ at the $k$-th boundary node on $\partial\omega_{i}^{*}\cap\Omega$ and $h_{i}^{k}=0$ at all other boundary nodes. Here $1\leq k\leq n_{i}$. For boundary patches, $\partial\omega_{i}^{*}\cap\partial\Omega\neq\emptyset$, only hat functions corresponding to interior nodes, $\partial\omega_{i}^{*}\cap\Omega$, are used as boundary data, and $h_{i}^{k}(x)$ satisfies homogeneous Neumann data on $\partial\omega_{i}^{*}\cap\partial\Omega_{N}$ and homogeneous Dirichlet data on $\partial\omega_{i}^{*}\cap\partial\Omega_{D}$.
We write the $k$-th A-harmonic function over $\omega_{i}^{*}$ as

$$ ${w}_{i}^{k}({x})=\sum_{e=1}^{n_{el}}\sum_{l=1}^{n_{en}}\mathbf{w}_{i,e}^{k,l}H_{l}({\alpha}_{e}\left({x})\right),$ (32) $$

where $\mathbf{w}_{i,e}^{k}$ is the vector of nodal values at the $e$-th element, $1\leq e\leq n_{el}$, contained in $\omega_{i}^{*}$, with $\mathbf{w}_{i,e}^{k,l}$ the component at the $l$-th node, $1\leq l\leq n_{en}$.

### 3.2 Construction of Local Spectral Bases, V ω i m i subscript superscript 𝑉 subscript 𝑚 𝑖 subscript 𝜔 𝑖 V^{m_{i}}_{\omega_{i}}

Now we construct the spectral basis that delivers the local approximation with optimal approximation properties. We introduce the span $V_{\omega_{i}^{*}}^{m_{i}}=\operatorname{span}\left\{\xi_{i}^{1},\dots\xi_{i}^{m_{i}}\right\}$, where the $\xi_{i}^{j}$ are the eigenfunctions corresponding to the $m_{i}$ largest eigenvalues of

$$ $\mathrm{Q}_{i}\mathbf{x}=\lambda\mathrm{P}_{i}\mathbf{x}.$ (33) $$

Here the entries of the matrices $\mathrm{Q}_{i}$ and $\mathrm{P}_{i}$ are given as

$$ $\displaystyle\mathrm{Q}_{i}^{jk}$ $\displaystyle=(w_{i}^{j},w_{i}^{k})_{\mathcal{E}(\omega_{i})}=\int_{\omega_{i}}A\nabla w_{i}^{j}\nabla w_{i}^{k}dx=\int_{\partial\omega_{i}}({n}\cdot A\nabla w_{i}^{j})w_{i}^{k}ds,$ (34) $\displaystyle\mathrm{P}_{i}^{jk}$ $\displaystyle=(w_{i}^{j},w_{i}^{k})_{\mathcal{E}(\omega_{i}^{*})}=\int_{\omega_{i}^{*}}A\nabla w_{i}^{j}\nabla w_{i}^{k}dx=\int_{\partial\omega_{i}^{*}}({n}\cdot A\nabla w_{i}^{j})w_{i}^{k}ds,$ (35) $$

where we have used $S^{n_{i}}_{\omega_{i}^{*}}$ as an approximate basis for $H_{A}(\omega_{i}^{*})$ in the local spectral problems defined by ([20](#S2.E20)). These functions give the discrete approximation to the $n$-width functions when $\xi^{j}_{i}$’s are restricted to $\omega_{i}$ and the space $V^{m_{i}}_{\omega_{i}}=\operatorname{span}\left\{\xi_{i}^{1},\dots\xi_{i}^{m_{i}}\right\}$ gives the discrete approximation to the spectral subspace on $\omega_{i}$. The $k$-th entry of the vector $\mathbf{x}$ in ([33](#S3.E33)) is defined to be the coefficient of the $k$-th basis function $w_{i}^{k}(x)\in S^{n_{i}}_{\omega_{i}^{*}}$. Thus, the $q$-th, $1\leq q\leq m_{i}$, eigenfunction over $\omega_{i}^{*}$ is written as

$$ ${\xi}_{i}^{q}({x})=\sum_{k=1}^{m_{i}}\xi_{i}^{q,k}{w}_{i}^{k}({x}),$ (36) $$

with coefficients $\mathbf{x}_{k}=\xi_{i}^{q,k}$.
Here the eigenfunctions are listed as $\{\xi_{i}^{1},\xi_{i}^{2},\ldots,\xi_{i}^{m_{i}}\}$, where the corresponding eigenvalues are listed in descending order $1>\lambda_{i}^{1}\geq\lambda_{i}^{2}\geq\cdots\lambda_{i}^{m_{i}}>0$ ($\mathrm{P}_{i}$ and $\mathrm{Q}_{i}$ are symmetric positive definite).

Combining ([32](#S3.E32)) and ([36](#S3.E36)) gives the $n$-width functions written in terms of the underlying finite element discretization

$$ ${\xi}_{i}^{q}({x})=\sum_{e=1}^{n_{el}}\sum_{k=1}^{m_{i}}\sum_{l=1}^{n_{en}}\xi_{i}^{q,k}\mathbf{w}_{i,e}^{k,l}H_{l}({\alpha}_{e}({x})).$ $$

The entries of the spectral matrices ([34](#S3.E34)) may be computed quickly by rewriting the middle integrals in ([34](#S3.E34)) as

$$ $\displaystyle Q_{i}^{jk}$ $\displaystyle=\int_{\omega}A\nabla w^{j}\cdot\nabla(v_{\omega}+\hat{w}^{k})\,ds$ $\displaystyle P_{i}^{jk}$ $\displaystyle=\int_{\omega^{*}}A\nabla w^{j}\cdot\nabla(v_{\omega^{*}}+\tilde{w}^{k})\,ds$ $$

where $w^{k}=v_{\omega}+\hat{w}^{k}$ on $\omega$ with $v_{\omega}\in H^{1}_{0}(\omega)$ vanishing on $\partial\omega$ and $w^{k}=v_{\omega^{*}}+\tilde{w}^{k}$ with $v_{\omega^{*}}\in H^{1}_{0}(\omega^{*})$ vanishing on $\partial\omega^{*}$. In the finite element discretization $\hat{w}^{k}$ is zero on all nodes in the interior of $\omega$ and $\hat{w}^{k}=w^{k}$ on nodes lying on $\partial\omega$. Similarly $\tilde{w}^{k}$ is zero on all nodes interior to $\omega^{*}$ and $\tilde{w}^{k}=w^{k}$ on nodes lying on $\partial\omega^{*}.$ Then the spectral matrices may be computed by

$$ $\displaystyle Q_{i}^{jk}$ $\displaystyle=\int_{\omega}A\nabla w^{j}\cdot\nabla(\hat{w}^{k})\,ds$ $\displaystyle P_{i}^{jk}$ $\displaystyle=\int_{\omega^{*}}A\nabla w^{j}\cdot\nabla(\tilde{w}^{k})\,ds$ $$

where the integrand is nonzero only on elements away from the boundary. This follows since $H_{A}^{0}\perp H^{1}_{0}$ in the $\lVert\cdot\rVert_{\mathcal{E}}$-norm by definition. This integral reduces the number of elements to be integrated over from quadratic to linear in two dimensions and cubic to quadratic in three dimensions. The reduction in memory usage is the same, allowing for a much larger FEM mesh to be used underneath these computations.

### 3.3 Global Stiffness Matrix and Right Hand Side

Using the bases for the local trial fields $V^{m_{i}}_{\omega_{i}}$ the global trial field $V^{N}$ is defined as in ([10](#S2.E10)) by

$$ $V^{N}=\operatorname{span}\{\phi_{i}\xi_{i}^{q}:1\leq i\leq N,\text{ }1\leq q\leq m_{i},\text{ and }\xi_{i}^{q}\in V^{m_{i}}_{\omega_{i}}\}.$ $$

The local spaces $V^{m_{i}}_{\omega_{i}}$ are augmented by the constant functions for interior patches or patches with Neumann boundary conditions.
For boundary patches on which Dirichlet boundary conditions are enforced, $\partial\omega_{i}\cap\partial\Omega_{D}\neq\emptyset$, constant functions need not be added to $V^{m_{i}}_{\omega_{i}}$.
The Galerkin discretization defined by ([13](#S2.E13)) gives the matrix equation $\mathrm{G}\mathbf{x}=\mathbf{r}$ where the stiffness matrix $\mathrm{G}$ is defined entrywise by

$$ $\mathrm{G}_{iqjr}=(\phi_{i}{\xi}_{i}^{q},\phi_{j}{\xi}_{j}^{r})_{\mathcal{E}(\Omega)}.$ (37) $$

The $ij$-th block of $G$ corresponds to cross products of functions from $\phi_{i}V^{m_{i}}_{\omega_{i}}$ and $\phi_{j}V^{m_{j}}_{\omega_{j}}$. The $iq$-th element of the unknown vector $\mathbf{x}$ is then the coefficient of the global trial function $\phi_{i}\xi_{i}^{q}$, so that we have the part of the solution corresponding to the local A-harmonic functions written

$$ $u^{G}=\sum_{i=1}^{N}\sum_{q=1}^{m_{i}}\mathbf{x}_{i}^{q}\phi_{i}\xi_{i}^{q}.$ (38) $$

Likewise, the right hand side of ([13](#S2.E13)) is written entrywise as

$$ $\mathbf{r}_{iq}=F(\phi_{i}\xi_{i}^{q})-B(u^{F},\phi_{i}\xi_{i}^{q})=\int_{\Omega}f\phi_{i}\xi_{i}^{q}\,dx+\int_{\partial\Omega_{N}}g\phi_{i}\xi_{i}^{q}\,dS-\int_{\Omega}u^{F}\phi_{i}\xi_{i}^{q}\,dx.$ (39) $$

With $\mathrm{G}\mathbf{x}=\mathbf{r}$ defined as above, solving for $\mathbf{x}$ gives the A-harmonic part of the solution ([38](#S3.E38)).
The approximate solution to ([2](#S2.E2)) is given by

$$ $u_{0}=u^{G}+u^{F}=\sum_{i=1}^{N}\sum_{q=1}^{m_{i}}\mathbf{x}_{i}^{q}\phi_{i}\xi_{i}^{q}+\sum_{i=1}^{N}\phi_{i}\chi_{i}$ $$

for local particular solutions $\chi_{i}$ as in ([11](#S2.E11)).

## 4 Numerical Implementation for a Benchmark Problem

For the implementation of the MS-GFEM algorithm multiple partial differential equations must be solved. It should be noted that any method to solve PDEs may be used, here we use the finite element method. The problem domain $\Omega$ has an underlying finite element mesh and overlapping patches $\omega_{i}$ to which the mesh conforms. All computations are performed over the finite elements. A concrete numerical example of the MS-GFEM algorithm is demonstrated below.

### 4.1 Problem Formulation

All computations in this work are done on the same rectangular domain with 100 separated circular inclusions with small variations in radii (figure [1](#S4.F1)). For demonstration of the numerical implementation we consider the following problem.

$$ $\left\{\begin{aligned} -\operatorname{div}(A(x_{1},x_{2})\nabla u(x_{1},x_{2}))&=0,\text{ }(x_{1},x_{2})\in\Omega\\ u(x_{1},x_{2})&=0\text{, }(x_{1},x_{2})\text{ on the left of }\partial\Omega\\ u(x_{1},x_{2})&=1\text{, }(x_{1},x_{2})\text{ on the right of }\partial\Omega\\ {n}\cdot A\nabla u(x_{1},x_{2})&=0,\text{ }(x_{1},x_{2})\text{ on the top and bottom of }\partial\Omega.\end{aligned}\right.$ (40) $$

Dirichlet boundary conditions are applied on the left and right of the domain while the top and bottom have homogeneous Neumann boundary conditions.
The matrix material has property $A_{1}$ and the inclusion material has property $A_{2}$. They are given as

$$ $A_{1}=\begin{pmatrix}\alpha&0\\ 0&\alpha\end{pmatrix},\,\,\,A_{2}=\begin{pmatrix}\beta&0\\ 0&\beta\end{pmatrix}.$ $$

In the section [5](#S5) we will study a large range of material contrasts given by the ratio $\alpha/\beta$. For the purpose of illustration we shall take $\alpha=1$ and $\beta=100$ in this section.

Figure: Figure 1: Rectangular domain of size $20\times 10$ with 100 circular inclusions for problem ([40](#S4.E40)).
Refer to caption: /html/1910.08856/assets/inclusions_01_100total.png

### 4.2 Mesh Selection

All computations in this work are performed on the same rectangular domain. Here we will compare the accuracy of the MS-GFEM implementation to the standard FEM method. In order to compare the methods we select the underlying mesh used for the MS-GFEM method by solving for the standard FEM solution for the same problem over progressively finer meshes for the highest contrast used in our study. The FEM error for a given mesh is measured using the a posteriori estimate

$$ $\epsilon_{h}^{2}=\frac{\lVert u_{h/2}\rVert^{2}-\lVert u_{h}\rVert^{2}}{\lVert u\rVert^{2}}\frac{1}{1-2^{-2\gamma}},$ (41) $$

where we use bilinear elements and assume $\epsilon_{h}^{2}=Ch^{2\gamma}$ for $1/2\leq\gamma\leq 1$ (see []). Here the a posteriori estimate ([41](#S4.E41)) is standard and we have chosen $\gamma=1/2$. The finest mesh which gives an acceptable FEM solution over all material contrasts will be the same one used for comparison of the MS-GFEM method with the FEM method. Here we will carry out the study for contrasts ranging from $1/1000$ to $1000$. From Table [1](#S4.T1) we choose the finest mesh with $1,560,058$ elements as the common mesh in which to compare MS-GFEM solutions. The a posteriori estimate for the FEM using finer and finer mesh seeds is shown in figure [2](#S4.F2). For the purposes of this study we denote the FEM solution with the finest mesh as the overkill solution for a given contrast and denote this solution by $u$. The MS-GFEM solution will now be computed over this same mesh. The following subsections describe the MS-GFEM implementation.

Figure: Figure 2: Log plot of a posteriori error estimate against the log of the number of mesh elements The green curve is the a posteriori contrast for the matrix conductivity 1000 and the inclusion 1, while the blue curve is for matrix conductivity 1 and inclusion conductivity 1000. For other contrasts the convergence graphs are essentially the same.
Refer to caption: /html/1910.08856/assets/x1.png

**Table 1: Log of a posteriori errors computed according to ([41](#S4.E41)).**
|  | A posteriori Error by Contrast (fiber:matrix) |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| # Elements | Log(# Elements) | 100:1 | 200:1 | 1000:1 | 1:100 | 1:200 | 1:1000 |
| 6686 | 3.825 | – | – | – | – | – | – |
| 24187 | 4.384 | -0.94312 | -0.94368 | -0.94418 | -0.64733 | -0.64384 | -0.64103 |
| 95017 | 4.978 | -1.20170 | -1.20145 | -1.20129 | -0.93778 | -0.93424 | -0.93140 |
| 388189 | 5.589 | -1.48816 | -1.48767 | -1.48730 | -1.23698 | -1.23345 | -1.23061 |
| 1560058 | 6.193 | -1.78071 | -1.78001 | -1.77948 | -1.54125 | -1.53772 | -1.53489 |

### 4.3 Construction of Partition of Unity

The domain is a $20\times 10$ rectangle with 100 circular inclusions with small variation in radii. To fix ideas we work with a two patch covering $\{\omega_{1},\omega_{2}\}$ of $\Omega$ by a rectangular patch $\omega_{1}=12\times 6$ and a rectangular annulus $\omega_{2}=20\times 10-8\times 4$. The expanded patches are defined as $\omega_{1}^{*}=16\times 8$ and $\omega_{2}^{*}=20\times 10-4\times 2$. The covering is shown in Figure [3](#S4.F3).

Figure: Figure 3: Top row: The patches $\omega_{1}$, $\omega_{1}^{*}$, respectively. Bottom row: $\omega_{2}$, $\omega_{2}^{*}$, respectively. Red indicates the subdomain, green is the area of $\Omega$ outside the subdomain.
Refer to caption: /html/1910.08856/assets/n_om_1_20004.png

To construct the partition of unity function $\phi_{1}$ over $\omega_{1}$, the overlap between $\omega_{1}$ and $\omega_{2}$ is determined and split into two sets as shown in Figure [4](#S4.F4). On $\omega_{1}-\omega_{2}$ (the $8\times 4$ green rectangle on the bottom left of Figure [3](#S4.F3)) the partition of unity function is 1, outside $\omega_{1}$ it is 0. On the overlap shown on the left of Figure [4](#S4.F4), $\phi_{1}$ is linear. On the corners of the overlap, shown on the right, it is bilinear. The values of $\phi_{1}$ are computed at the nodes and after that interpolated to the integration points using the finite element shape functions. The partition of unity function $\phi_{2}$ is then just $\phi_{2}=1-\phi_{1}$. The partition of unity functions and their derivatives are shown in Figure [5](#S4.F5).

Figure: Figure 4: Sets in the overlap of $\omega_{1}\cap\omega_{2}$ shown in red. On the left region $\phi_{1}$ is linear. On the right $\phi_{1}$ is bilinear.
Refer to caption: /html/1910.08856/assets/pou_edges_smoothed.png

Figure: Figure 5: The partition of unity functions $\phi_{1}$ (top left), $\phi_{2}$ (top right), and the derivatives with respect to $x_{1}$ (bottom left) and $x_{2}$ (bottom right) over $\Omega$.
Refer to caption: /html/1910.08856/assets/pou_omega_1.png

### 4.4 Local Particular Solutions

We find the local particular solutions $\chi_{1}$ and $\chi_{2}$ over $\omega_{1}^{*}$ and $\omega_{2}^{*}$, respectively. On boundaries of the patches $\partial\omega_{i}^{*}$ which do not coincide with the outer boundary $\partial\Omega$ we enforce $\chi_{i}=0$ and on boundaries which coincide with outer boundaries the given boundary conditions in ([40](#S4.E40)) are applied. Since the PDE is homogeneous and the boundary data is homogeneous Dirichlet data, $\chi_{1}\equiv 0$.
For $\omega_{2}^{*}$ the problem becomes

$$ $\left\{\begin{aligned} -\operatorname{div}(A\nabla\chi_{2}(x_{1},x_{2}))&=0\text{, }(x_{1},x_{2})\in\omega_{2}^{*}\\ -n\cdot A\nabla\chi_{2}(x_{1},x_{2})&=0\text{, }(x_{1},x_{2})\in\partial\omega_{2}^{*}\cap\Omega\\ \chi_{2}(x_{1},x_{2})&=0\text{, }(x_{1},x_{2})\text{ on the left side of }\partial\omega_{2}^{*}\cap\partial\Omega_{D}\\ \chi_{2}(x_{1},x_{2})&=1\text{, }(x_{1},x_{2})\text{ on the right side of }\partial\omega_{2}^{*}\cap\partial\Omega_{D}\\ -{n}\cdot A\nabla\chi_{2}(x_{1},x_{2})&=0\text{, }(x_{1},x_{2})\in\partial\omega_{2}^{*}\cap\partial\Omega_{N}.\end{aligned}\right.$ (42) $$

The particular solution $\chi_{2}$ is shown in Figure [6](#S4.F6). The finite element analysis for $\chi_{2}$ is used in the assembly of the right hand side, see ([39](#S3.E39)).

Figure: Figure 6: The particular solution $\chi_{2}$ displayed on $\omega_{2}^{*}$. Here $\chi_{1}\equiv 0$ since $\omega_{1}^{*}$ is an interior domain ($\partial\Omega\cap\partial\omega_{1}^{*}=\emptyset$) and the PDE has 0 right hand side.
Refer to caption: /html/1910.08856/assets/particular_solution_omega_2.png

### 4.5 Construction of Local Solution Spaces S ω i ∗ n i subscript superscript 𝑆 subscript 𝑛 𝑖 superscript subscript 𝜔 𝑖 S^{n_{i}}_{\omega_{i}^{*}}

Next we generate the local spaces $S^{n_{1}}_{\omega_{1}^{*}}$ and $S^{n_{2}}_{\omega_{2}^{*}}$. Here we use as boundary data the piece wise linear hat functions defined on $\partial\omega_{i}^{*}$. Taken together these piece wise linear functions form a partition of unity on $\partial\omega_{i}^{*}$. To make an element of $S^{n_{1}}_{\omega^{*}_{1}}$ we consider each of the $1\leq j\leq n_{1}$ nodes on $\partial\omega_{1}^{*}$. For each node we solve the discrete Dirichlet problem

$$ $\left\{\begin{aligned} -\operatorname{div}(A\nabla w_{1}^{j}(x,y))&=0,\text{ }(x,y)\in\omega_{1}^{*}\\ w_{1}^{j}(x,y)&=1\text{, }\text{on the }j\text{-th node of }\partial\omega_{1}^{*}\\ w_{1}^{j}(x,y)&=0\text{, }\text{every other node on }\partial\omega_{1}^{*}.\end{aligned}\right.$ (43) $$

the collection of A-harmonic functions constructed this way gives the basis for $S^{n_{1}}_{\omega^{*}_{1}}$.
On $\omega_{2}^{*}$, $S^{n_{2}}_{\omega_{2}^{*}}$ consists of functions corresponding to hat functions on $\partial\omega_{2}^{*}\cap\Omega$ as boundary data. Homogeneous Dirichlet and Neumann conditions are enforced on $\partial\omega_{2}^{*}\cap\partial\Omega_{D}$ and $\partial\omega_{2}^{*}\cap\partial\Omega_{N}$, respectively. For each of the $1\leq j\leq n_{2}$ nodes on $\partial\omega_{2}^{*}\cap\Omega$ the following problem is solved.

$$ $\left\{\begin{aligned} -\operatorname{div}(A\nabla w_{2}^{j}(x,y))&=0,\text{ }(x,y)\in\omega_{2}^{*}\\ w_{2}^{j}(x,y)&=1\text{, }\text{on the }j\text{-th node of }\partial\omega_{2}^{*}\cap\Omega\\ w_{2}^{j}(x,y)&=0\text{, }\text{every other node on }\partial\omega_{2}^{*}\cap\Omega\\ w_{2}^{j}(x,y)&=0\text{, }(x,y)\in\partial\omega_{2}^{*}\cap\partial\Omega_{D}\\ -{n}\cdot A\nabla w_{2}^{j}(x,y)&=0,\text{ }(x,y)\in\partial\omega_{2}^{*}\cap\partial\Omega_{N}.\end{aligned}\right.$ (44) $$

Figure [7](#S4.F7) shows two computed functions $w_{i}^{j}(x,y)$ over $\omega_{1}^{*}$ and $\omega_{2}^{*}$.

Figure: Figure 7: Two of the computed functions $w_{1}^{j}(x,y)$ (left) and $w_{2}^{j}(x,y)$ (right) shown over $\omega_{1}^{*}$ and $\omega_{2}^{*}$, respectively, using a hat function with a single node as support for boundary data. We have re-scaled the region near the boundary datum to make the field visible.
Refer to caption: /html/1910.08856/assets/hat_function_1_node_20006_1c.jpg

The figure shows that the functions $w_{i}^{j}(x,y)$ are localized and decay rapidly away from the boundary when the hat function used as boundary data has support of only one node. We note that if we consider broader hat functions with support over several boundary nodes instead of just one node then the associated solutions of ([43](#S4.E43)) and ([44](#S4.E44)) have energy densities that extend and are nonzero further inside $\omega_{1}^{*}$ and $\omega_{2}^{*}$. To see this we provide a computation using a hat function with support over 301 nodes on the boundary of $\omega_{1}^{*}$ and 101 nodes on the boundary of $\omega_{2}^{*}$ instead of one node, results are shown in Figure [8](#S4.F8). In section [6](#S6) we investigate the use of hat function boundary data with different support to generate the local solution spaces. In this approach we pick a collection of hat functions all of the same width and taken together form a partition of unity on the boundary. It is clear that using fewer boundary datum reduces the number of problems to solve in generating $S^{n_{i}}_{\omega^{*}_{i}}$. However the interesting observation found in section [6](#S6) is that the reduction of the number of boundary datum used in generating the space done in this way does not significantly effect the convergence rate of the spectral basis. This is because A-harmonic extensions of boundary data with large support on the boundary have energy that decays slowly into the domain relative to A-harmonic extensions of boundary data with small support sets, see Figures [7](#S4.F7) and [8](#S4.F8) and [] and []. We say that solutions with slow energy decay penetrate into the domain.

Figure: Figure 8: An example of A-harmonic functions $w_{i}^{j}$ on $\omega_{1}^{*}$ (left) and $\omega_{2}^{*}$ (right) using a hat function with 301 nodes (left) and 101 nodes (right) as support for boundary data.
Refer to caption: /html/1910.08856/assets/hat_function_301_node_20006_1.png

In conclusion the equations ([43](#S4.E43)) and ([44](#S4.E44)) will be solved using bilinear or quadratic FEM, the functions $w_{i}^{j}$ are saved at the nodes of the underlying finite element mesh and represented by finite element shape functions and the stiffness matrices $P_{i}$ and $Q_{i}$ are assembled.

### 4.6 Construction of Local Spectral Bases

Now we use $S_{\omega_{i}^{*}}^{n_{i}}$ to construct the local spectral basis with span given by the space $V_{\omega_{i}^{*}}^{m_{i}}={\rm span}\left\{\xi_{i}^{1},\dots\xi_{i}^{m_{i}}\right\}$, where the $\xi_{i}^{j}$ are the eigenfunctions corresponding to the $m_{i}$ largest eigenvalues of

$$ $\mathrm{Q}_{i}\mathbf{x}=\lambda\mathrm{P}_{i}\mathbf{x}$ (45) $$

with entries of $\mathrm{P}$ and $\mathrm{Q}$ defined as in ([34](#S3.E34)) and ([35](#S3.E35)). For this study the integration is performed over the domain rather than the boundary using the underlying finite elements, and, due to symmetry, only for the entries $j\geq k$. Equation ([34](#S3.E34)) is written as

$$ $\int_{\omega_{i}}A\nabla w_{i}^{j}\nabla w_{i}^{k}\,dx=\sum_{e=1}^{N^{e}_{\omega_{i}}}{(\mathbf{w}_{i,e}^{j})^{T}}\mathrm{K}_{e}\mathbf{w}_{i,e}^{k},$ (46) $$

where $T$ denotes the transpose, $N^{e}_{\omega_{i}}$ is the number of elements in patch $\omega_{i}$, $\mathbf{w}_{i,e}^{j}$ and $\mathbf{w}_{i,e}^{k}$ are the vectors of nodal values of the functions $w_{i}^{j}$ and $w_{i}^{k}$ at the element $e$, and $\mathrm{K}_{e}$ is the element stiffness matrix.
The integral ([35](#S3.E35)) over $\omega_{i}^{*}$ is computed by adding ([46](#S4.E46)) to the integral over the set $\omega_{i}^{*}-\omega_{i}$. The summation is done in parallel using OpenMP. To ensure that the matrix $\mathrm{P}_{i}$ is positive definite, each row $j$ of $\mathrm{P}_{i}$ and $\mathrm{Q}_{i}$ is normalized by $\mathrm{P}_{i}^{jj}$.

Figure: Figure 9: Log plot of eigenvalues for $\omega_{1}$ (left) and $\omega_{2}$ (right) for conductivity in the matrix 100 and in the fibers 1.
Refer to caption: /html/1910.08856/assets/x2.png

The generalized eigenvalue problems are then solved using routines from the Intel^® Math Kernel Library (MKL). The problems are reduced to standard symmetric eigenvalue problems. In the following description of the methodology, we neglect the index $i$ keeping in mind that the problem has to be solved for each patch. Considering that the matrix $\mathrm{P}$ is symmetric and positive definite, a Cholesky factorization is performed

$$ $\mathrm{P}=\mathrm{U}^{T}\mathrm{U}$ (47) $$

where $\mathrm{U}$ is an upper triangular matrix. Equation ([47](#S4.E47)) is then inserted in ([45](#S4.E45))

$$ $\mathrm{Q}\mathbf{x}=\lambda\mathrm{U}^{T}\mathrm{U}\mathbf{x}.$ $$

The inverse of $\mathrm{U}^{T}$ is computed and multiplied on both sides. On the left side there is also an identity matrix $\mathrm{I}=\mathrm{U}^{-1}\mathrm{U}$ inserted between $\mathrm{Q}$ and $\mathbf{x}$

$$ ${\mathrm{U}^{T}}^{-1}\mathrm{Q}\mathrm{U}^{-1}\mathrm{U}\mathbf{x}=\lambda\mathrm{U}\mathbf{x}.$ (48) $$

With $\mathbf{y}=\mathrm{U}\mathbf{x}$ and $\mathrm{C}={\mathrm{U}^{T}}^{-1}\mathrm{Q}\mathrm{U}^{-1}$ ([48](#S4.E48)) becomes a standard eigenvalue problem

$$ $\mathrm{C}\mathbf{y}=\lambda\mathbf{y}.$ (49) $$

The eigenvalues for problem ([49](#S4.E49)) are the same as for ([45](#S4.E45)), the eigenvectors $\mathbf{x}$ are obtained by solving $\mathrm{U}\mathbf{x}=\mathbf{y}$.

The eigenvalues for the two patches are shown in Figure [9](#S4.F9). Eigenfunctions corresponding to eigenvalues which are smaller than $10^{-12}$ are discarded. One can see in the graph on the right that the eigenvalues come in pairs, due to symmetry, except for the first one of $\omega_{2}$. Figure [10](#S4.F10) shows $\xi_{1}^{4}$, $\xi_{1}^{8}$, $\xi_{2}^{1}$, and $\xi_{2}^{6}$. In addition to the local spectral basis functions, $V^{m_{1}}_{\omega_{1}^{*}}$ is augmented by the constant functions.

Figure: Figure 10: Local spectral basis functions $\xi_{1}^{5}$ (top left) and $\xi_{1}^{7}$ (top right) over $\omega_{1}^{*}$ and $\xi_{2}^{1}$ (bottom left) and $\xi_{2}^{6}$ (bottom right) over $\omega_{2}^{*}$.
Refer to caption: /html/1910.08856/assets/ef_omega1_005.png

### 4.7 Solving the Global System, G ​ 𝐱 = 𝐫 G 𝐱 𝐫 \mathrm{G}\mathbf{x}=\mathbf{r}

The global stiffness matrix and right hand side entries, ([37](#S3.E37)) and ([39](#S3.E39)), are computed using integration on the underlying finite element mesh similarly to ([46](#S4.E46)) using the saved element stiffness matrices from the particular solution runs. The constructions of the global stiffness matrix and right hand side are performed in parallel using OpenMP.
The global system $\mathrm{G}\mathbf{x}=\mathbf{r}$ is solved using the Pardiso solver from the Intel^® MKL. The solution vector $\mathbf{x}$ is used to recover $u^{G}$ as in ([38](#S3.E38)). The solution $u^{G}$ is shown in Figure [11](#S4.F11) (left). The particular solution $u^{F}$ over the whole domain $\Omega$ consisting of the particular solutions over the patches glued together with the partition of unity functions is shown in Figure [11](#S4.F11) (right).

Figure: Figure 11: The solutions $u^{G}$ (left) and $u^{F}$ (right).
Refer to caption: /html/1910.08856/assets/u_G__20006.png

The final solution is then recovered from $u_{0}=u^{G}+u^{F}$ shown in Figure [12](#S4.F12) (left). The solution provided directly by a finite element method is displayed in Figure [12](#S4.F12) (right) for comparison with the MS-GFEM solution.

Figure: Figure 12: The final solution $u_{0}$ computed by the MS-GFEM (left) and an FEM (right).
Refer to caption: /html/1910.08856/assets/u_MSGFEM_20006.png

The MS-GFEM computed solution has energy $\lVert u_{0}\rVert_{\mathcal{E}(\Omega)}^{2}=0.799227$. The FEM, “overkill,”solution has energy given by $\lVert u\rVert_{\mathcal{E}(\Omega)}^{2}=0.799227.$ The relative error of the MS-GFEM solution compared to the FEM solution is given by

$$ $\frac{\lVert u-u_{0}\rVert_{\mathcal{E}(\Omega)}}{\lVert u\rVert_{\mathcal{E}(\Omega)}}=2.52\times 10^{-5}.$ $$

## 5 Contrast Independent Convergence Study

In this section we present numerical results indicating that the MS-GFEM converges to the “overkill” FEM solution at an exponential rate that is independent of contrast between material properties of the matrix and the inclusions.
To demonstrate the contrast independent convergence rate we consider again problem ([40](#S4.E40)). We compute the MS-GFEM solution $u_{0}$ for different material contrasts between the matrix and inclusion materials. We adopt the convention $A:B$ where the first number is the matrix conductivity and the number following the colon is the inclusion conductivity. We carry out the simulation for 1:1, 1:1.5, 1:2, 1:10, 1:100, 1:200, and 1:1000 and 1.5:1, 2:1, 10:1, 100:1, 200:1, and 1000:1.
All errors are computed relative to the FEM “overkill” solution $u$ for each choice of contrast. The relative error with respect to the energy norm is computed as

$$ $E=\frac{\lVert u_{0}-u\rVert_{\mathcal{E}(\Omega)}^{2}}{\lVert u\rVert_{\mathcal{E}(\Omega)}^{2}}.$ (50) $$

Tables [2](#S5.T2) and [3](#S5.T3) show the relative errors of the MS-GFEM solution $u_{0}$ computed across contrasts. The number of spectral basis functions from $\omega_{1}$ and $\omega_{2}$ used in computing $u_{0}$ are listed under their respective headings in the tables.

**Table 2: Relative error versus dimension of spectral basis for different contrasts computed over a mesh of 1,560,058 elements. Here the A-harmonic local space $S^{n_{i}}_{\omega_{i}}$ is generated by hat functions defined on $\partial\omega^{*}_{i}$ of width 1 node.**
| $\omega_{1}$ | $\omega_{2}$ | $1000:1$ | $200:1$ | $100:1$ | $10:1$ | $2:1$ | $1.5:1$ | $1:1$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 5 | 2 | 1.25E-01 | 1.25E-01 | 1.24E-01 | 1.23E-01 | 1.08E-01 | 1.01E-01 | 8.98E-02 |
| 10 | 4 | 2.64E-02 | 2.61E-02 | 2.58E-02 | 2.29E-02 | 1.81E-02 | 1.58E-02 | 2.52E-02 |
| 15 | 6 | 6.11E-03 | 6.04E-03 | 5.98E-03 | 6.60E-03 | 7.80E-03 | 6.86E-03 | 1.75E-03 |
| 20 | 8 | 1.86E-03 | 1.85E-03 | 1.86E-03 | 3.35E-03 | 3.73E-03 | 2.78E-03 | 1.65E-03 |
| 25 | 10 | 1.05E-03 | 1.06E-03 | 1.08E-03 | 1.72E-03 | 2.18E-03 | 1.74E-03 | 1.28E-03 |
| 30 | 12 | 2.23E-04 | 2.83E-04 | 3.60E-04 | 1.31E-03 | 1.56E-03 | 1.08E-03 | 3.58E-04 |
| 35 | 14 | 8.79E-05 | 1.19E-04 | 1.64E-04 | 8.61E-04 | 1.10E-03 | 7.30E-04 | 5.33E-05 |
| 40 | 16 | 3.20E-05 | 4.63E-05 | 7.06E-05 | 4.82E-04 | 6.61E-04 | 4.43E-04 | 4.12E-05 |
| 45 | 18 | 2.72E-05 | 3.69E-05 | 5.43E-05 | 3.35E-04 | 4.15E-04 | 2.63E-04 | 2.77E-05 |
| 50 | 20 | 2.42E-05 | 3.04E-05 | 4.30E-05 | 2.37E-04 | 1.94E-04 | 1.23E-04 | 1.22E-05 |
| 55 | 22 | 2.42E-05 | 2.75E-05 | 2.86E-05 | 1.19E-04 | 9.26E-05 | 5.48E-05 | 6.89E-06 |
| 60 | 24 | 2.51E-05 | 2.49E-05 | 2.55E-05 | 4.39E-05 | 4.18E-05 | 2.87E-05 | 6.70E-06 |
| 65 | 26 | 2.68E-05 | 2.57E-05 | 2.52E-05 | 3.52E-05 | 2.53E-05 | 1.69E-05 | 6.25E-06 |
| 70 | 28 | 2.81E-05 | 2.67E-05 | 2.59E-05 | 2.59E-05 | 2.04E-05 | 1.46E-05 | 5.91E-06 |

**Table 3: Relative error versus dimension of spectral basis for different contrasts computed over a mesh of 1,560,058 elements. Here the A-harmonic local space $S^{n_{i}}_{\omega_{i}}$ is generated by hat functions defined on $\partial\omega_{i}^{*}$ of width 1 node.**
| $\omega_{1}$ | $\omega_{2}$ | $1:1000$ | $1:200$ | $1:100$ | $1:10$ | $1:2$ | $1:1.5$ | $1:1$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 5 | 2 | 2.89E+00 | 1.23E+00 | 8.49E-01 | 2.16E-01 | 8.01E-02 | 8.11E-02 | 8.98E-02 |
| 10 | 4 | 1.39E+00 | 6.08E-01 | 4.34E-01 | 1.17E-01 | 2.94E-02 | 2.35E-02 | 2.52E-02 |
| 15 | 6 | 8.48E-01 | 4.40E-01 | 3.31E-01 | 9.30E-02 | 1.52E-02 | 9.46E-03 | 1.75E-03 |
| 20 | 8 | 2.09E-01 | 1.09E-01 | 6.83E-02 | 1.17E-02 | 5.28E-03 | 3.35E-03 | 1.65E-03 |
| 25 | 10 | 2.87E-02 | 1.48E-02 | 1.12E-02 | 8.01E-03 | 3.13E-03 | 1.98E-03 | 1.28E-03 |
| 30 | 12 | 4.05E-03 | 3.52E-03 | 3.42E-03 | 4.34E-03 | 1.85E-03 | 1.14E-03 | 3.58E-04 |
| 35 | 14 | 2.57E-03 | 2.52E-03 | 2.52E-03 | 2.57E-03 | 1.12E-03 | 6.99E-04 | 5.33E-05 |
| 40 | 16 | 9.46E-04 | 9.65E-04 | 1.01E-03 | 1.29E-03 | 7.21E-04 | 4.31E-04 | 4.12E-05 |
| 45 | 18 | 5.21E-04 | 4.72E-04 | 4.32E-04 | 4.81E-04 | 3.61E-04 | 2.43E-04 | 2.77E-05 |
| 50 | 20 | 4.55E-04 | 3.93E-04 | 3.66E-04 | 2.66E-04 | 1.31E-04 | 8.51E-05 | 1.22E-05 |
| 55 | 22 | 5.59E-05 | 1.35E-04 | 3.18E-04 | 2.04E-04 | 8.73E-05 | 5.36E-05 | 6.89E-06 |
| 60 | 24 | 3.40E-05 | 4.89E-05 | 6.34E-05 | 1.67E-04 | 5.33E-05 | 3.06E-05 | 6.70E-06 |
| 65 | 26 | 2.63E-05 | 3.33E-05 | 3.86E-05 | 4.86E-05 | 2.81E-05 | 1.78E-05 | 6.25E-06 |
| 70 | 28 | 2.50E-05 | 2.24E-05 | 2.33E-05 | 3.37E-05 | 1.83E-05 | 1.22E-05 | 5.91E-06 |

**Table 4: Relative errors versus dimension of spectral basis of the MS-GFEM solution using quadratic finite elements on a mesh with 388,189 elements and 1,156,882 nodes. Here the A-harmonic local space $S^{n_{i}}_{\omega_{i}}$ is generated by hat functions defined on $\partial\omega_{i}^{*}$ of width 3 nodes.**
| $\omega_{1}$ | $\omega_{2}$ | $1000:1$ | $100:1$ | $1:1$ | $1:100$ | $1:1000$ |
| --- | --- | --- | --- | --- | --- | --- |
| 5 | 2 | 1.25E-01 | 1.24E-01 | 9.03E-02 | 8.02E-01 | 2.75E+00 |
| 10 | 4 | 2.64E-02 | 2.59E-02 | 2.49E-02 | 2.45E-01 | 7.74E-01 |
| 15 | 6 | 6.08E-03 | 5.97E-03 | 1.74E-03 | 1.40E-01 | 3.80E-01 |
| 20 | 8 | 1.85E-03 | 1.83E-03 | 1.63E-03 | 5.84E-02 | 2.11E-01 |
| 25 | 10 | 1.03E-03 | 1.05E-03 | 1.48E-03 | 1.32E-02 | 2.07E-02 |
| 30 | 12 | 2.08E-04 | 2.87E-04 | 3.71E-04 | 4.74E-03 | 5.50E-03 |
| 35 | 14 | 8.32E-05 | 1.23E-04 | 3.92E-05 | 2.44E-03 | 3.19E-03 |
| 40 | 16 | 3.19E-05 | 5.88E-05 | 3.79E-05 | 7.73E-04 | 7.48E-04 |
| 45 | 18 | 1.97E-05 | 3.93E-05 | 3.02E-05 | 3.06E-04 | 3.89E-04 |
| 50 | 20 | 9.85E-06 | 2.56E-05 | 1.57E-05 | 2.87E-04 | 3.66E-04 |
| 55 | 22 | 2.98E-06 | 9.70E-06 | 3.23E-06 | 2.45E-04 | 3.45E-05 |
| 60 | 24 | 1.45E-06 | 6.43E-06 | 2.54E-06 | 4.27E-05 | 2.29E-05 |
| 65 | 26 | 7.42E-07 | 4.22E-06 | 1.81E-06 | 1.96E-05 | 1.03E-05 |
| 70 | 28 | 6.05E-07 | 2.37E-06 | 1.04E-06 | 1.04E-05 | 7.39E-06 |

Figure: Figure 13: Log plot of relative error of MS-GFEM solution (computed with bilinear FEM) versus number of local basis functions used. Each line is for a different matrix conductivity and inclusions with conductivity 1.
Refer to caption: /html/1910.08856/assets/1hat_mat_greater_fib-contrast_1000_100-50-20-10-2-1_v2.png

Figure: Figure 14: Log plot of relative error of MS-GFEM solution (computed with bilinear FEM) versus number of local basis functions used. Each line is for matrix of conductivity 1 and inclusions of different conductivity.
Refer to caption: /html/1910.08856/assets/1hat_mat_less_fib-contrast_1-1p1-1p3-1p5-2-10-100-200-1000_v2.png

Figure: Figure 15: Comparison of quadratic and linear finite elements for computing MS-GFEM global relative errors. Log plot of relative error versus number of local basis functions used. Each line is the relative error for quadratic or linear elements of contrast listed with matrix conductivity followed by inclusion conductivity.
Refer to caption: /html/1910.08856/assets/quad_vs_lin_mat_greater_fib-contrast_1000_100-1_v2.png

In Figure [13](#S5.F13) the log of the relative errors is plotted against the dimension of the spectral basis for inclusions with conductivity 1 and matrix with conductivity 1, 2, 10, 20, 50, 100, and 1000.
Figure [14](#S5.F14) shows the same for inclusions with conductivity 1, 1.1, 1.3, 1.5, 2, 10, 100, 200, and 1000, and matrix with conductivity 1. These figures show the exponential convergence of the error seen with respect to the dimension of the local approximation space. Here the contrast between matrix and inclusions is seen to not influence the convergence rate.

Next we increase the order of finite elements to quadratic to see what effect it has on the rate of convergence of the relative error versus dimension of the local spectral bases. The relative errors of the MS-GFEM solution using quadratic finite elements on a mesh with 388,189 elements and 1,156,882 nodes are now used. The spectral bases are now generated by an A-harmonic approximation space $S_{\omega_{i}^{*}}^{n_{i}}$ with boundary data given by hat functions defined on $\partial\omega^{*}_{i}$ that are supported over three boundary nodes. The convergence rates for the associated MS-GFEM for several differenct contrasts between matrix and particles are displayed in Figure [15](#S5.F15). The convergence rate for MS-GFEM using quadratic elements is also compared to the convergence rate for MS-GFEM using linear elements in Figure [15](#S5.F15). The simulations show that the change of element type has no effect of the convergence rates associated with the dimension of the local spectral basis up to a dimension of about fifty. Beyond fifty the convergence rate of MS-GFEM using linear elements the convergence rate goes to zero, whereas the convergence rate for MS-GFEM with quadratic elements continues to converge exponentially with dimension of up to eighty. The relative error decreases below $10^{-5}$ for all cases using spectral bases of dimension greater than seventy. Here the convergence for MS-GFEM using linear elements flattens and is zero after dimension fifty due to the linear FEM approximation error incurred in the generation of the discrete A-harmonic spaces $S_{\omega^{*}_{i}}^{m_{i}}$.

### 5.1 Eigenvalue Decay and Size of Oversampling Domain

The key feature of MS-GFEM is the construction of optimal local solution spaces through oversampling. The size of the oversampled patch $\omega^{*}$ as compared to $\omega$ has a direct impact on the decay rate of the eigenvalues which bound the local errors (cf. Theorem [2](#Thmtheorem2)). We investigate the size of the oversampling domain $\omega^{*}$ relative to the domain $\omega$ used in approximation. We first consider oversampeling domain $\omega^{*}$ given by a disk of radius $R$ and the smaller concentric disk $\omega$ of radius $\rho$, $R>\rho$. We look at the A-harmonic space of solutions of ([9](#S2.E9)) with a homogeneous coefficient $A(x)=I$ for $\omega_{i}^{*}=\omega^{*}$. A straight forward calculation shows that the $n$-width eigenvalues are given by

$$ $\lambda_{n}=\left(\frac{\rho}{R}\right)^{2n}.$ (51) $$

Now keeping in mind that the earlier numerical results given here for particle composites illustrate the contrast independent convergence of the MS-GFEM, and so the $n$-width eigenvalues should also be independent of contrast. Therefore we argue that a relationship like $\lambda_{n}=\left(\frac{\rho}{R}\right)^{2n}$ should also hold for the heterogeneous case. Our computations show that this is approximately true. To illustrate the point we consider rectangular domains with coefficient matrix $A(x)$ associated with particle composites as in Figure [1](#S4.F1) with contrast $1$ in the matrix and $1000$ in the inclusions. We denote by $\rho$ and $R$ the longest side lengths of $\omega$ and $\omega^{*}$, respectively. Figure [16](#S5.F16) shows the decay of the $n$-width eigenvalues with respect the ratio $\rho/R$ to for three values $16/18$, $12/18$, and $6/18$. We see from these simulations that there is good agreement between the formula ([51](#S5.E51)) given by the solid lines and the numerically computed convergence rate given by the data points.

Figure: Figure 16: Log plot of eigenvalues for different subdomains $\omega_{1}$ holding $\omega_{1}^{*}$ fixed. In this case the side lengths of $\omega_{1}$ are proportionately changed. Conductivity 1 in matrix and 1000 in inclusions. Eigenfunctions generated with 19-node width hat functions as boundary data for $H_{A}(\omega^{*})$.
Refer to caption: /html/1910.08856/assets/x4.png

## 6 Reduction of Computational Work by Penetration; an Effective Computation of Finite Dimensional Subspaces of H A ​ ( ω i ∗ ) subscript 𝐻 𝐴 superscript subscript 𝜔 𝑖 H_{A}(\omega_{i}^{*})

The primary numerical work in generating the best local basis is not in the solution of the eigenvalue problem but in the numerical generation of the local A-harmonic subspace over which the eigenvalue problem is solved.
Numerically generating the A-harmonic subspaces $S^{n_{i}}_{\omega_{i}^{*}}\subset H_{A}(\omega_{i}^{*})$ can be computationally expensive for fine meshes. Defining the boundary data as in ([31](#S3.E31)) with a single boundary element for the support of the hat function $h_{i}^{j}$ results in as many problems to solve as there are boundary nodes on $\partial\omega_{i}^{*}$. Most importantly, we see from Figure [7](#S4.F7) that the resulting function $w_{i}^{j}$ given by the A-harmonic extension of the boundary hat functions with support over one boundary element has very small support and gradient inside $\omega^{\ast}$. However for A-harmonic extensions of boundary hat functions having 301 nodes as support we see from Figure [8](#S4.F8) that the support of this function penetrates further into in $\omega^{\ast}$. This is the phenomenon of penetration where A-harmonic extensions of boundary data with less oscillation have a gradient that penetrates further into $\omega^{*}$, than A-harmonic extensions of more oscillatory boundary data, see [] and []. In light of the exponential decay of the $n$-width eigenfunctions in the energy norm over $\omega$ it becomes clear that the space spanned by the lower $n$-width modes are generated by linear combinations of A-harmonic extensions of coarse boundary data, ie., boundary hat functions of large support. This motivates the usage of significantly wider hat functions on the boundary such that their A-harmonic extension will have significant variation that penetrates into $\omega$. This means that less functions need to be used for the computation of the local space $S^{n_{i}}_{\omega^{*}_{i}}$ if we use boundary data given by hat functions having support over several boundary nodes.

We see from Theorem [2](#Thmtheorem2) that the local error controlled by the $n$-width eigenvalues, decay at worst exponentially so the number of local spectral basis functions $\xi_{i}^{j}$ can be reduced with little effect on the error. For the example of section [4](#S4), the eigenvalues decay rapidly, see Figure [9](#S4.F9), and the 100-th and 50-th eigenvalues dip below $10^{-15}$ for $\omega_{1}$ and $\omega_{2}$, respectively. Therefore only a handful of the associated $n$-width eigenfunctions need to be used in the local approximation spaces to achieve excellent accuracy; the rest can be discarded.

Thus if we want to approximate $M$ dimensional local subspaces using $n$ spectral basis functions ($n<M$) we are motivated to use $M$ wider boundary boundary hat functions and their A harmonic extensions. This results in fewer problems to solve and smaller spectral matrices in ([33](#S3.E33)) since the number of A-harmonic functions $w_{i}^{j}$ can be of the same order of wide hat functions used for boundary data. This gives a strategy for reducing the computational work in computing the local basis and for computing the spectral basis.

When the number of boundary hat functions is reduced from $N$ to $M$ by defining boundary hat functions with larger support then the computational time for constructing $S^{m_{i}}_{\omega_{i}^{*}}$ is reduced by a factor of $\frac{M}{N}$. The time for computing each entry of the matrices ([34](#S3.E34)) and ([35](#S3.E35)) is approximately the same. Due to symmetry of the matrices, only the upper triangular part and the diagonal need to be computed. If N hat functions were used then $\frac{N(N+1)}{2}$ entries have to be computed. If that number is reduced to $M$ then

$$ $\frac{M(M+1)}{2}=\frac{N\rho\left(N\rho+1\right)}{2}$ $$

entries have to be computed, where $\rho=M/N$. For large values of $N$

$$ $\frac{M(M+1)}{2}\sim\frac{N^{2}\rho^{2}}{2}$ $$

and the computational time for generating the matrices is reduced by a factor of $\left(\frac{M}{N}\right)^{2}$. The number of hat functions used for each successive widening progresses as $\{N,N/2,N/3,N/4,\ldots\}$. For example, if the support of the hat functions is increased from 1 to 3 nodes, then the total number of boundary value problems ([31](#S3.E31)) to solve is $M=N/2$ and the time to fill the spectral matrices is about $1/4$ the time as with single node hat functions. We illustrate these observations with the following example.

We solve the problem ([40](#S4.E40))
as before
on the domain shown in Figure [1](#S4.F1) using hat functions with varying widths to generate the local solution spaces $S^{n_{i}}_{\omega_{i}^{*}}$. We label solutions $u_{k},\ k=1,3,5,7,9,\dots$, where $k$ represents the number of nodes in the support of the hat function used as boundary data for problem ([31](#S3.E31)).

The matrix material has conductivity $A_{1}$ and the inclusion material has conductivity $A_{2}$. We compute solutions to ([40](#S4.E40)) using two different sets of material properties. For Case 1

$$ $A_{1}=\begin{pmatrix}1&0\\ 0&1\end{pmatrix}\ \ \ A_{2}=\begin{pmatrix}100&0\\ 0&100\end{pmatrix}$ $$

and for Case 2 we reverse the material properties to

$$ $A_{1}=\begin{pmatrix}100&0\\ 0&100\end{pmatrix}\ \ \ A_{2}=\begin{pmatrix}1&0\\ 0&1\end{pmatrix}.$ $$

The relative error of the final solutions $u_{k}$ as compared to the “overkill” FEM solution $u$ is computed with respect to the energy norm as

$$ $E_{k}=\frac{\lVert u_{k}-u\rVert_{\mathcal{E}(\Omega)}^{2}}{\lVert u\rVert_{\mathcal{E}(\Omega)}^{2}}.$ (52) $$

The relative global errors versus the support of the hat functions used to generate the local approximation spaces are reported in Table [5](#S6.T5) and are graphed in Figure [17](#S6.F17). The result is that the global error is insensitive to the support of the local hat functions up to a width of twenty five nodes. Here the “Error 1” column lists the relative error for a matrix of conductivity $1$ and inclusion of conductivity $500$. The “Error 2” column lists the relative errors for a matrix of conductivity $500$ and inclusion of conductivity $1$. The table also shows the size of the spectral matrices $\mathrm{P}$ and $\mathrm{Q}$ and the number of eigenfunctions used in the computation of $u_{k}$.

**Table 5: Relative error for global MS-GFEM solutions using hat functions with different numbers of boundary nodes for support. The first column indicates the number of nodes in the support of the hat functions. The first error is measured in the energy norm forwhen the matrix material has conductivity 1 and the inclusions 100, and the second errors are when the conductivities are switched. Also the size of the spectral matrices $\mathrm{P}$ and $\mathrm{Q}$ (number of hat functions used), and the number of local spectral basis functions used in the approximation (chosen such that $10^{-12}\leq\lambda\leq 1$).**
|  | Eigenvalue System Size | Number of Eigenfunctions |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| # | Error 1 | Error 2 | $\omega_{1}$ | $\omega_{2}$ | $\omega_{1}$ | $\omega_{2}$ |
| 1 | 4.185E-05 | 5.004E-05 | 3840 | 961 | 124 | 49 |
| 3 | 4.147E-05 | 5.007E-05 | 1920 | 481 | 124 | 49 |
| 5 | 4.135E-05 | 5.002E-05 | 1280 | 321 | 124 | 49 |
| 7 | 4.120E-05 | 5.003E-05 | 960 | 241 | 124 | 49 |
| 9 | 4.122E-05 | 5.008E-05 | 768 | 193 | 124 | 49 |
| 11 | 4.109E-05 | 4.998E-05 | 640 | 161 | 124 | 49 |
| 13 | 4.113E-05 | 5.001E-05 | 548 | 138 | 124 | 49 |
| 15 | 4.108E-05 | 4.984E-05 | 480 | 121 | 124 | 49 |
| 17 | 4.108E-05 | 4.984E-05 | 426 | 107 | 124 | 49 |
| 19 | 4.109E-05 | 4.969E-05 | 384 | 97 | 124 | 49 |
| 21 | 4.077E-05 | 4.951E-05 | 349 | 88 | 123 | 49 |
| 23 | 4.124E-05 | 4.950E-05 | 320 | 81 | 123 | 49 |
| 25 | 4.126E-05 | 4.936E-05 | 295 | 74 | 123 | 47 |

Figure: Figure 17: Plot of MS-GFEM relative global error versus the number of nodes in the support of the boundary hat functions used to generate the local spaces $S_{\omega_{i}^{*}}^{n_{i}}\subset H_{A}(\omega_{i}^{*})$. Here the two cases of matrix to inclusion conductivities of 100 to 1 and 1 to 100 are shown.

In this example the error ([52](#S6.E52)) is seen to be constant, see Figure [17](#S6.F17), while the size of the spectral matrices decays rapidly for the first few sets of wider hat functions. The number of eigenfunctions selected for the computation begins decreasing after widening the boundary hat functions beyond 9 nodes. Collecting these observations, we conclude that this method of generating the local A-harmonic space $S^{n_{i}}_{\omega_{i}^{*}}$ is an effective method to reduce computational costs while maintaining acceptable error tolerance.

## 7 Span of Hat Functions Used Directly as Local Approximation Space

In the last section we exploited the penetration of fields inside $\omega^{*}_{i}$ obtained by A-harmonically extending hat functions of large support on the boundary. Motivated by this we investigate the direct use of such functions as local approximations on $\omega_{i}\subset\omega^{*}_{i}$. We call the GFEM method using these “oversampled” local basis the oversampled-GFEM. We investigate the decay rate of the global error using oversampled-GFEM. Here we prescribe hat functions on the boundary $\partial\omega^{*}_{i}$ and compute their A-harmonic extensions into $\omega^{*}_{i}$ as in ([43](#S4.E43)) and ([44](#S4.E44)). Then the basis for the local approximation space is given by the restriction of these functions to $\omega_{i}\subset\omega_{i}^{*}$. As in the last section we take the hat functions on the boundary as elements of a partition of unity on the boundary of $\omega^{*}_{i}$. We try coarser to finer boundary partitions of unity to generate our local subspace for the oversampled-GFEM. We find as we increase the number of basis elements by refining the boundary partition of unity we generate global approximate solutions with accuracy that are nearly exponentially decreasing with the dimension of the local basis. Tables [6](#S7.T6), [7](#S7.T7), and [8](#S7.T8) provides a comparison of the relative global error obtained using the same dimension of local approximation space given by $m_{i}$ for the oversampled-GFEM and $V_{\omega_{i}}^{m_{i}}$ for the MS-GFEM. Each table addresses a different material contrast. Here we are comparing with MS-GFEM that use $n$-width eigenfunctions generated from a local high dimensional approximation space $S^{n_{i}}_{\omega^{*}_{i}}$ with $n_{i}>>m_{i}$ that is spanned by A-harmonic extensions of hat functions with support over one boundary element. The tables show that for all contrasts the MS-GFEM does better than oversampled-GFEM at reducing the relative error in the global approximation as the number of bases elements are decreased. However as the dimension of the local approximating space for oversampled-GFEM and MS-GFEM is reduced we see that the oversampled-GFEM delivers global solutions with relative error comparable to MS-GFEM. This numerical study is carried out for the matrix-inclusion composite with contrasts $1:100$ and $100:1$, $1:500$ and $500:1$, and $1:1000$ and $1000:1$. Here the coarseness of the hat function A-harmonic local basis used in the oversampled-GFEM is measured by the support of the boundary hat functions, see Tables [6](#S7.T6), [7](#S7.T7), and [8](#S7.T8). The relative error versus total degrees of freedom given in Tables [6](#S7.T6), [7](#S7.T7), and [8](#S7.T8) are plotted in Figure [18](#S7.F18). Here the top left chart in Figure [18](#S7.F18) shows the exponential decay of oversampled-GFEM with respect to global degrees of freedom for different material contrasts and the top right chart displays the exponential decay of MS-GFEM with respect to global degrees of freedon for different contrasts. The bottom left chart in Figure [18](#S7.F18) plots oversampled-GFEM together with MS-GFEM and their decay with global degrees of freedom with matrix conductivity equal to unity and different particle conductivity’s and the lower right hand chart plots the decay of the error of both oversampled-GFEM and MS-GFEM with respect to global degrees of freedom for different matrix conductivity’s and particle conductivity unity. Note that both MS-GFEM and oversampled-GFEM converge exponentially with respect to the global degrees of freedom for all contrasts but MS-GFEM gives the most rapid decrease in approximation error with respect to global degrees of freedom.

**Table 6: Relative global errors for oversampled-GFEM and MS-GFEM, both computed and compared to the FEM overkill solution (relative energy norm error). Contrasts 100:1 and 1:100. The columns $\omega_{1},\,\omega_{2}$ are the number of shape functions used in the corresponding subdomain. The local spaces $S_{\omega_{i}}^{n_{i}}$ are composed of elements obtained by A-harmonically extending boundary hat functions containing a single-node.**
|  | 100:1 | 1:100 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Hat Widths | $\omega_{1}$ | $\omega_{2}$ | GFEM | MS-GFEM | GFEM | MS-GFEM |
| 69 | 109 | 28 | 2.89e-05 | 2.91e-05 | 4.63e-05 | 2.25e-05 |
| 99 | 76 | 20 | 1.11e-04 | 2.75e-05 | 1.50e-04 | 4.34e-05 |
| 129 | 59 | 15 | 1.51e-04 | 3.14e-05 | 1.57e-04 | 1.61e-04 |
| 159 | 48 | 13 | 7.28e-04 | 5.03e-05 | 5.19e-04 | 4.90e-04 |
| 189 | 40 | 11 | 7.38e-04 | 1.11e-04 | 2.07e-03 | 1.21e-03 |
| 219 | 34 | 9 | 2.74e-03 | 3.73e-04 | 3.64e-03 | 7.74e-03 |
| 249 | 30 | 8 | 4.40e-03 | 1.15e-03 | 5.73e-03 | 9.77e-03 |
| 279 | 27 | 7 | 5.55e-03 | 4.52e-03 | 1.17e-02 | 1.36e-02 |
| 309 | 24 | 7 | 1.29e-02 | 4.72e-03 | 7.85e-03 | 1.63e-02 |
| 339 | 22 | 6 | 1.34e-02 | 4.88e-03 | 1.11e-02 | 3.80e-02 |
| 369 | 20 | 6 | 1.22e-02 | 5.02e-03 | 1.28e-02 | 6.86e-02 |
| 399 | 19 | 5 | 2.66e-02 | 1.57e-02 | 1.53e-02 | 1.09e-01 |
| 499 | 15 | 4 | 3.50e-02 | 1.64e-02 | 4.93e-02 | 3.53e-01 |
| 599 | 12 | 4 | 2.04e-02 | 1.86e-02 | 8.39e-02 | 3.94e-01 |

**Table 7: Relative global errors for oversampled-GFEM and MS-GFEM, both computed and compared to the FEM overkill solution (relative energy norm error). Contrasts 500:1 and 1:500. The columns $\omega_{1},\,\omega_{2}$ are the number of shape functions used in the corresponding subdomain. The local spaces $S_{\omega_{i}}^{n_{i}}$ are composed of elements obtained by A-harmonically extending boundary hat functions containing a single-node.**
|  | 500:1 | 1:500 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Hat Widths | $\omega_{1}$ | $\omega_{2}$ | GFEM | MS-GFEM | GFEM | MS-GFEM |
| 69 | 109 | 28 | 6.81E-03 | 3.06E-05 | 3.41E-05 | 4.61E-05 |
| 99 | 76 | 20 | 1.04E-04 | 2.81E-05 | 1.33E-04 | 4.10E-05 |
| 129 | 59 | 15 | 1.31E-04 | 2.75E-05 | 1.49E-04 | 2.08E-04 |
| 159 | 48 | 13 | 6.99E-04 | 3.51E-05 | 6.92E-04 | 5.46E-04 |
| 189 | 40 | 11 | 6.65E-04 | 1.04E-04 | 2.60E-03 | 1.28E-03 |
| 219 | 34 | 9 | 2.70E-03 | 3.40E-04 | 3.92E-03 | 1.16E-02 |
| 249 | 30 | 8 | 4.39E-03 | 1.07E-03 | 7.23E-03 | 1.70E-02 |
| 279 | 27 | 7 | 5.49E-03 | 4.50E-03 | 1.81E-02 | 2.10E-02 |
| 309 | 24 | 7 | 1.29E-02 | 4.74E-03 | 1.03E-02 | 2.80E-02 |
| 339 | 22 | 6 | 1.34E-02 | 4.91E-03 | 1.63E-02 | 7.18E-02 |
| 369 | 20 | 6 | 1.22E-02 | 5.00E-03 | 1.56E-02 | 1.74E-01 |
| 399 | 19 | 5 | 2.68E-02 | 1.57E-02 | 2.34E-02 | 3.39E-01 |
| 499 | 15 | 4 | 3.54E-02 | 1.65E-02 | 9.79E-02 | 7.59E-01 |
| 599 | 12 | 4 | 2.05E-02 | 1.87E-02 | 1.61E-01 | 8.58E-01 |

**Table 8: Relative global errors for oversampled-GFEM and MS-GFEM, both computed and compared to the FEM overkill solution (relative energy norm error). Contrasts 1000:1 and 1:1000. The columns $\omega_{1},\,\omega_{2}$ are the number of shape functions used in the corresponding subdomain. The local spaces $S_{\omega_{i}}^{n_{i}}$ are composed of elements obtained by A-harmonically extending boundary hat functions containing a single-node.**
|  | 1000:1 | 1:1000 |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Hat Widths | $\omega_{1}$ | $\omega_{2}$ | GFEM | MS-GFEM | GFEM | MS-GFEM |
| 69 | 109 | 28 | 2.99E-05 | 3.09E-05 | 9.71E-04 | 5.62E-05 |
| 99 | 76 | 20 | 1.17E-04 | 2.85E-05 | 5.87E-03 | 4.65E-05 |
| 129 | 59 | 15 | 1.30E-04 | 2.80E-05 | 1.48E-04 | 2.29E-04 |
| 159 | 48 | 13 | 6.95E-04 | 3.49E-05 | 7.71E-04 | 5.64E-04 |
| 189 | 40 | 11 | 6.55E-04 | 1.05E-04 | 2.87E-03 | 1.36E-03 |
| 219 | 34 | 9 | 2.70E-03 | 3.38E-04 | 4.09E-03 | 1.26E-02 |
| 249 | 30 | 8 | 4.39E-03 | 1.06E-03 | 7.91E-03 | 2.24E-02 |
| 279 | 27 | 7 | 5.48E-03 | 4.49E-03 | 2.14E-02 | 2.60E-02 |
| 309 | 24 | 7 | 1.29E-02 | 4.75E-03 | 1.17E-02 | 3.60E-02 |
| 339 | 22 | 6 | 1.34E-02 | 4.91E-03 | 1.94E-02 | 9.93E-02 |
| 369 | 20 | 6 | 1.22E-02 | 5.00E-03 | 1.79E-02 | 2.21E-01 |
| 399 | 19 | 5 | 2.68E-02 | 1.57E-02 | 3.03E-02 | 5.44E-01 |
| 499 | 15 | 4 | 3.55E-02 | 1.65E-02 | 1.36E-01 | 1.06E+00 |
| 599 | 12 | 4 | 2.05E-02 | 1.87E-02 | 2.19E-01 | 1.20E+00 |

Figure: Figure 18: Relative errors of oversampled-GFEM and MS-GFEM solutions for different contrasts versus total number of shape functions used. Errors computed relative to the FEM overkill solution.
Refer to caption: /html/1910.08856/assets/x6.png

## 8 Iterative solution implementation

The prime motivation behind the MS-GFEM and its contrast independent exponential convergence rate is its potential use in large parallel implementations. As expressed in the introduction the method has two essential parts. A. Parallel construction of local shape functions leading to the exponential approximation of the actual solution in the energy norm. Here each local basis is defined over separate patches and is computed independently on separate processors using local memory. B. Efficient parallel solution of the global system of linear algebraic equations. Recall that we have the partition of unity $\{\phi_{i}\}_{i=1}^{N}$ subordinate to the overlapping patches $\omega_{i}$, with $\Omega=\cup_{i=1}^{N}\omega_{i}$ and $\sum_{i=1}^{N}\phi_{i}=1$.
We consider the span $V^{m_{i}}_{\omega_{i}}$ of $n$-width functions $\xi_{j}^{i},\,j=1,\ldots,m_{i}$ defined on $\omega_{i}^{\ast}$ restricted to $\omega_{i}$. Set ${V^{m_{i}}_{\omega_{i}}}^{\ast}=\phi_{i}V^{m_{i}}_{\omega_{i}}$ and the global approximation space is

$$ $V^{N}=\sum_{i=1}^{N}{V^{m_{i}}_{\omega_{i}}}^{\ast}.$ (53) $$

We use the partition of unity structure of the global approximation $V^{N}$ space to iteratively solve the global problem applying the Schwarz method for overlapping domains. Because each patch overlaps only with its neighbours we can develop an iterative scheme that can be implemented in parallel. Proceeding in this way we arrive at the preconditioned Richardson iteration

$$ $\mathrm{B}_{N}(\mathbf{x}^{n+1}-\mathbf{x}^{n})+\mathrm{G}\mathbf{x}^{n}=\mathbf{r}_{N},\,n=0,1,\ldots,x^{0}\in\mathbb{R}^{N},$ (54) $$

where $B_{N}$ is taken to be an additive two level Schwartz preconditioner [], []. This scheme is implemented in parallel and the iteration process converges exponentially. Here we will assume a fixed number of cores per patch. For this implementation the time for convergence is anticipated to be constant for a fixed ratio between the size of the problem and the number of cores, so it weakly scales with the number of cores used in the computation []. This iterative method will be addressed in detail in a forthcoming paper.

## 9 Conclusions

The MS-GFEM numerical method is shown to be a viable method for the computation of local fields in heterogeneous materials with any contrast between material properties. A full implementation of the method along with the theoretical foundation of the method were discussed. A method for dramatically reducing computational costs in computing local trial fields given by the $n$-width functions was presented along with numerical results supporting the method. We remark that A-harmonic extensions with long wavelength boundary data when restricted to subdomains have energies on these subdomains that decay slower than short wavelength data as the subdomain boundary moves away from where the data is prescribed, see [] and []. Motivated by this fact we have proposed to use A-harmonic extensions of boundary partitions of unity where the boundary hat functions have large support as local bases. Our numerical investigations show that local basis generated this way provide nearly exponential decay in the relative global error as the dimension of the approximation space is increased. Use of these bases provides a further reduction in computational cost. In future work the implementation of MS-GFEM discussed here will be further parallelized for use on large scale computers. For much larger problems with many patches greater savings from computational costs are anticipated by using coarser hat functions on the boundary to generate local approximation spaces. We finish by noting that in [] related work has been recently carried out showing the contrast independence of the local spectral basis functions on their domains and that the global solutions converge exponentially in the $L^{2}$ norm. There it is found that the local $n$-width basis functions outperform all other contemporary choices of local trial fields constructed by oversampling.

## 10 Appendix

For completeness we show how to estimate the local error for patches that share a common boundary with the computational domain.
For patches with Neumann boundary data ([3](#S2.E3)) we introduce the space $\mathcal{U}=\{v\in H^{1}(\omega^{*}_{i}):\,v=0\hbox{ on }\partial{\omega^{*}_{i}}\cap\Omega\}$. The local particular solution $\chi_{i}$ of ([11](#S2.E11)) with ([3](#S2.E3)) solves the equivalent variational problem: $\chi_{i}\in\mathcal{U}$ and

$$ $\displaystyle\int_{\omega_{i}^{*}}A(x)\nabla\chi_{i}\cdot\nabla v\,dx=\int_{\omega^{*}_{i}}vf\,dx+\int_{\partial\omega^{*}_{i}\cap\partial\Omega}vg\,ds,$ (55) $$

for all $v$ in $\mathcal{U}$.
Next observe that the actual solution $u$ also solves

$$ $\displaystyle\int_{\omega_{i}^{*}}A(x)\nabla u\cdot\nabla v\,dx=\int_{\omega^{*}_{i}}vf\,dx+\int_{\partial\omega^{*}_{i}\cap\partial\Omega}vg\,ds,$ (56) $$

for all $v$ in $\mathcal{U}$ so

$$ $\displaystyle\int_{\omega_{i}^{*}}A(x)\nabla\chi_{i}\cdot\nabla v\,dx=\int_{\omega_{i}^{*}}A(x)\nabla u\cdot\nabla v\,dx$ (57) $$

for all $v$ in $\mathcal{U}$. On choosing $v=\chi_{i}$ in ([57](#S10.E57)) and applying Cauchy’s inequality we get

$$ $\displaystyle\|\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}\leq\|u\|_{\mathcal{E}(\omega_{i}^{*})},$ (58) $$

and as before a simple calculation shows

$$ $\displaystyle\|u-\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}\leq 2\|u\|_{\mathcal{E}(\omega_{i}^{*})}.$ (59) $$

and ([29](#S2.E29)) follows.

For Dirichlet boundary data on $\partial\omega_{i}^{*}\cap\partial\Omega$ we use linearity and write
$\chi_{i}=\chi_{i}^{R}+\chi_{i}^{D}$, where $\chi_{i}^{R}\in H^{1}_{0}(\omega_{i}^{*})$ solves

$$ $\displaystyle\int_{\omega_{i}^{*}}A(x)\nabla\chi_{i}^{R}\cdot\nabla v\,dx=\int_{\omega^{*}_{i}}vf\,dx,$ (60) $$

for all $v$ in $H_{0}^{1}(\omega_{i}^{*})$, and $\chi_{i}^{D}=u$ on $\partial\omega_{i}^{*}\cap\partial\Omega$, where $n\cdot A(x)\nabla\chi_{i}^{D}=0$ on $\partial\omega_{i}^{*}\cap\Omega$, and solves

$$ $\displaystyle\operatorname{div}(A(x)\nabla\chi_{i}^{D})=0,$ (61) $$

in $\omega_{i}^{*}$. To conclude we will show that

$$ $\displaystyle\|\chi_{i}^{D}\|_{\mathcal{E}(\omega_{i}^{*})}\leq\|u-\chi_{i}^{R}\|_{\mathcal{E}(\omega^{*}_{i})}\leq 2\|u\|_{\mathcal{E}(\omega^{*}_{i})}.$ (62) $$

This follows from

$$ $\displaystyle\|\chi_{i}^{D}\|_{\mathcal{E}(\omega_{i}^{*})}^{2}=\inf\left\{\|w\|_{\mathcal{E}(\omega_{i}^{*})}^{2}:\,w\in\mathcal{A}_{ad}\right\}$ (63) $\displaystyle\hbox{where }\mathcal{A}_{ad}=\left\{w\in H^{1}(\omega_{i}^{*}),\,\operatorname{div}(A(x)\nabla w)=0,\hbox{ in }\omega_{i}^{*},\,\,w=u\hbox{ on }\partial\omega_{i}^{*}\cap\partial\Omega\right\}$ $$

To see ([63](#S10.E63)) write any $w\in\mathcal{A}_{ad}$ as $\chi_{i}^{D}+\delta$, $\delta=w-\chi_{i}^{D}$ to observe $(\chi_{i}^{D},\delta)_{\mathcal{E}(\omega_{i}^{*})}=0$ so

$$ $\displaystyle\|\chi_{i}^{D}\|^{2}_{\mathcal{E}(\omega_{i}^{*})}\leq\|w\|^{2}_{\mathcal{E}(\omega_{i}^{*})}=\|\chi_{i}^{D}\|^{2}_{\mathcal{E}(\omega_{i}^{*})}+\|\delta\|^{2}_{\mathcal{E}(\omega_{i}^{*})}$ (64) $$

and since $u-\chi_{i}^{R}\in\mathcal{A}_{ad}$ the first inequality of ([62](#S10.E62)) follows. Note as before $\|\chi_{i}^{R}\|_{\mathcal{E}(\omega_{i}^{*})}\leq\|u\|_{\mathcal{E}(\omega_{i}^{*})}$ and together with the triangle inequality this gives the second inequality of ([62](#S10.E62)).
On collecting results we conclude using the triangle inequality that

$$ $\displaystyle\|\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}=\|\chi_{i}^{D}+\chi_{i}^{R}\|_{\mathcal{E}(\omega_{i}^{*})}\leq 3\|u\|_{\mathcal{E}_{(}\omega^{*}_{i})},$ (65) $$

and as before a simple calculation shows

$$ $\displaystyle\|u-\chi_{i}\|_{\mathcal{E}(\omega_{i}^{*})}\leq 2\sqrt{5}\|u\|_{\mathcal{E}(\omega_{i}^{*})}.$ (66) $$

and ([30](#S2.E30)) follows.

For completeness we finish by showing that ([15](#S2.E15)) implies ([16](#S2.E16)). Our calculation follows Theorem 3.3 and remarks 3.4 and 3.5 of []. We expand

$$ $\displaystyle\|u-u_{T}\|^{2}_{{\mathcal{E}}(\Omega)}$ $\displaystyle=\int_{\Omega}\,A\nabla(u-u^{T})\cdot\nabla(u-u^{T})\,dx$ (67) $\displaystyle=\int_{\Omega}\,A\nabla(\sum_{i=1}^{N}\phi_{i}(u-\xi_{i}-\chi_{i})\cdot\nabla(\sum_{i=1}^{N}\phi_{i}(u-\xi_{i}-\chi_{i})\,dx$ $\displaystyle=\int_{\Omega}\,A(\sum_{i=1}^{N}k_{i})\cdot(\sum_{i=1}^{N}k_{i})\,dx$ $$

with

$$ $\displaystyle k_{i}=[\nabla\phi_{i}(u-\xi_{i}-\chi_{i})+\phi_{i}\nabla(u-\xi_{i}-\chi_{i})].$ (68) $$

It follows that

$$ $\displaystyle\|u-u_{T}\|^{2}_{{\mathcal{E}}(\Omega)}$ $\displaystyle\leq 2\int_{\Omega}\,A(\sum_{i=1}^{N}\nabla\phi_{i}(u-\xi_{i}-\chi_{i})\cdot(\sum_{i=1}^{N}\nabla\phi_{i}(u-\xi_{i}-\chi_{i}))\,dx$ (69) $\displaystyle+2\int_{\Omega}\,A(\sum_{i=1}^{N}\phi_{i}\nabla(u-\xi_{i}-\chi_{i})\cdot(\sum_{i=1}^{N}\phi_{i}\nabla(u-\xi_{i}-\chi_{i}))\,dx$ $\displaystyle\leq 2\kappa\sum_{i=1}^{N}\int_{\omega_{i}}\,A\nabla\phi_{i}(u-\xi_{i}-\chi_{i})\cdot\nabla\phi_{i}(u-\xi_{i}-\chi_{i}))\,dx$ $\displaystyle+2\kappa\sum_{i=1}^{N}\int_{\omega_{i}}\,A\phi_{i}\nabla(u-\xi_{i}-\chi_{i})\cdot\phi_{i}\nabla(u-\xi_{i}-\chi_{i}))\,dx$ $\displaystyle\leq C(\kappa\sum_{i=1}^{N}\|u-\xi_{i}-\chi_{i}\|_{L_{A}^{2}({\omega_{i}})}+\kappa\sum_{i=1}^{N}\|u-\xi_{i}-\chi_{i}\|_{\mathcal{E}(\omega_{i})}),$ $$

where for isotropic material properties $A(x)=Ia(x)$ we have $\|w\|_{L^{2}_{A}(\omega_{i})}=(\int_{\omega_{i}}\,aw^{2}\,dx)^{1/2}$. For interior patches and patches on the boundary with Neumann data the local spaces $V^{m_{i}}_{\omega_{i}}$ contain constant functions and we can choose them so that the weighted Poincaré inequality holds

$$ $\displaystyle\|u-\xi_{i}-\chi_{i}\|^{2}_{{L^{2}_{A}}(\omega_{i})}\leq C\|u-\xi_{i}-\chi_{i}\|^{2}_{{\mathcal{E}}(\omega_{i})}$ (70) $$

for a constant C independent of $u-\xi_{i}-\chi_{i}$ and $\omega_{i}$.
We also have a similar inequality for boundary patches associated with Dirichlet data as $\xi_{i}=0$ and $u-\chi_{i}=0$ for points on $\partial\omega_{i}\cap\partial\Omega_{D}$. This observation together with ([69](#S10.E69)), ([28](#S2.E28)), ([29](#S2.E29)) and ([30](#S2.E30)) shows that ([15](#S2.E15)) implies ([16](#S2.E16)).

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 3
presentation: 3
questions: N/A
soundness: 2
strengths: - The experimental results seems to be good.
- The presentation is well.
summary: This paper present an oversampling method O2 by maximizing the probability synthetic point belonging to the minority class. Theoretical guarantees and numerical experimental results are attached.
weaknesses: The core idea of this paper feels somewhat basic, and there seems to be a gap between its contribution and the typical ICLR standard. The theorems, particularly Theorem 2, appear straightforward and lack deeper insights; for instance, the discrepancy guarantee only presents an inequality without explicit bounds. Additionally, the citation in the main body of the paper is limited, making it challenging to assess the originality of the work; this may be due to the paper's simplicity and a lack of relevant references. Lastly, the experimental results seem limited in significance, though I acknowledge that, as I am less specialized in experimental parts of traditional machine learning, my view here may be less comprehensive—other reviewers may offer a more positive perspective on this aspect.

---

confidence: 4
contribution: 1
initial_rating: 1
presentation: 2
questions: 1. Will there be any side effects of the proposed method?
2. Is it possible that the generated synthetic minority class introduce generative error to the model? How to ensure that the generated synthetic minority class are beneficial?
soundness: 2
strengths: 1. The problem studied in this paper is important in the literature.
2. The paper is easy to follow.
3. The proposed method is easy to implement.
summary: This paper introduces a new oversampling framework, optimized oversampling, which generates synthetic minority class points by maximizing the probability of belonging to the minority class to improve the classifier's performance for imbalanced datasets.
weaknesses: 1. The proposed method is a data-level approach, and it can be integrated with some classification algorithms such as CART and LR. However, this is very limited in practice as many other popular classifiers are not considered in this paper.
2. The methods compared in the numerical experiments are outdated, and generative model-based methods [1-2] should be considered.
3. The experimental section is slightly weak. The performance of the proposed methods is not significantly better than that of the compared methods, especially for the F1 score, which is the most important metric in imbalanced classification.

[1] Odena, Augustus, et al. "Conditional Image Synthesis with Auxiliary Classifier GANs." 2017

[2] Mariani, Giovanni, et al. "BAGAN: Data Augmentation with Balancing GAN." 2018

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: The optional filtering step using a binary classifier is an interesting addition. How does the choice of classifier for this step impact the overall performance of O2? Would different classifiers yield significantly different results?
soundness: 2
strengths: 1. The paper provides theoretical backing for the effectiveness of O2.
2. The authors benchmark O2 against a variety of state-of-the-art oversampling methods across multiple datasets.
3. The case study on predicting mortality risk in patients with blunt spleen trauma illustrates the practical implications of the proposed method.
summary: The paper introduces a novel oversampling framework called Optimized Oversampling (O2) aimed at addressing the issue of class imbalance in machine learning classification tasks. The authors propose a method that generates synthetic minority class points by maximizing the probability of belonging to the minority class, as estimated by a trained classification model. The paper presents theoretical guarantees that the points generated by O2 are more likely to belong to the minority class compared to those generated by existing methods. The authors benchmark O2 against state-of-the-art oversampling techniques on 16 publicly available imbalanced datasets using Classification Trees (CART) and Logistic Regression (LR) for downstream classification tasks. The results indicate that O2 outperforms existing methods, particularly in highly imbalanced datasets. Additionally, a case study is presented, demonstrating the application of O2 in predicting mortality risk in patients with blunt spleen trauma, highlighting improvements in both performance and interpretability.
weaknesses: 1. While the paper evaluates O2 on 16 datasets, the diversity of these datasets in terms of domain and characteristics could be expanded. For example, the experiments can have results on some large and standard long-tailed datasets such as ImageNet-LT and CIFAR-LT.
2. The paper could benefit from a more detailed exploration of how these penalty coefficients affect performance and whether they require dataset-specific tuning.
3. The theoretical results are based on certain assumptions that may not hold in all practical scenarios.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: See Weaknesses.
soundness: 2
strengths: 1. This paper proposes a framework for generating high-probability minority class samples and provides a theoretical analysis of its effectiveness.

2. The proposed method is simple and effective in terms of both process and outcome, which is commendable.

3. The authors applied the proposed method to medical research, effectively predicting the mortality of patients with spleen trauma.
summary: This paper introduces Optimized Oversampling ($O^2$), a new framework to address class imbalance in classification tasks. $O^2$ generates synthetic minority class points by maximizing their probability of belonging to the minority class, as estimated by a trained model. Theoretical analysis shows $O^2$ has an advantage over other methods, and it outperforms state-of-the-art techniques on 16 imbalanced datasets, especially with the CART classifier.$O^2$ offers a deterministic, optimization-based approach that enhances synthetic points from other methods, proving particularly effective on highly imbalanced, medium-sized datasets and in real-world applications.

This method approaches sample generation from the perspective of objective function optimization, but it may not be effective for the more current issues, such as image datasets with a long-tailed distribution. Furthermore, the method for generating new samples for binary imbalance problems is not novel.
weaknesses: 1. The research topic is outdated, as these issues were proposed 20 years ago. The algorithms being compared, except for LoRAS, are all fairly old.

2. The parameters likely influence the experimental results. The impact of parameters $\lambda_1$, $\lambda_2$, and $\lambda_3$ on the results needs further discussion.

3. The main issue with this paper is that it does not address the current problem of long-tailed image distributions in imbalanced datasets. For high-dimensional, multi-class image datasets such as CIFAR10-LT, CIFAR100-LT, and iNaturalist2018, the effectiveness of the proposed oversampling method remains to be validated. 

If the authors can demonstrate that their method is effective on these long-tailed datasets, I would change my opinion of this paper from negative to positive.

Meta-review: The paper concerns the problem of imbalanced classification. The Authors introduce a new oversampling method that maximizes the probability of belonging to minority class of the generated examples. 

The Reviewers identified several flaws in the submission, including a very limited theoretical contribution, outdated references, and insufficient empirical evaluation.

All ratings are below the bar and the Authors did not prepare any rebuttal.

Decision: Reject