Computational homogenization of time-harmonic Maxwell's equations
11 Computational homogenization of time-harmonic Maxwell’sequations ∗ Patrick Henning and Anna Persson March 4, 2020
Abstract
In this paper we consider a numerical homogenization technique for curl-curl-problemsthat is based on the framework of the Localized Orthogonal Decomposition and which wasproposed in [
D. Gallistl, P. Henning, B. Verf¨urth. SIAM J. Numer. Anal. 56-3:1570–1596,2018 ] for problems with essential boundary conditions. The findings of the aforementionedwork establish quantitative homogenization results for the time-harmonic Maxwell’s equationsthat hold beyond assumptions of periodicity, however, a practical realization of the approachwas left open. In this paper, we transfer the findings from essential boundary conditionsto natural boundary conditions and we demonstrate that the approach yields a computablenumerical method. We also investigate how boundary values of the source term can effect thecomputational complexity and accuracy. Our findings will be supported by various numericalexperiments, both in 2 D and 3 D . Key words : Maxwell’s equations, multiscale, localized orthogonal decomposition, finite elementmethod, a priori analysis
AMS subject classifications:
In recent years, the interest in so called metamaterials has been growing rapidly. These materialsare typically engineered on a nano- or microscale to exhibit unusual properties, such as band gapsand negative refraction [26, 40, 27]. A common example is photonic crystals, which are composedof periodic cell structures interacting with light.To model photonic crystals and other metamaterials properly, a numerical method that canhandle such heterogeneous (or multiscale) media efficiently is needed. It is well known that classi-cal finite element methods struggle to produce good approximations in this case, unless the meshwidth is sufficiently small to resolve the microscopic variations intrinsic to the problem. Indeed,the mesh typically needs to be fine enough to resolve all features in the medium/metamaterial.This leads to issues with computational cost and available memory, which can be often onlyovercome by so-called multiscale methods. Multiscale methods are based on solving local de-coupled problems on the microscale whose solutions are used to correct the original modeland turn it into a homogenized macroscale model that can be solved cheaply. In the contextof electromagnetic waves in multiscale media, corresponding methods were proposed in, e.g.,[8, 24, 41, 35, 42, 20, 10, 25, 9].In [17] a method based on the Localized Orthogonal Decomposition (LOD) framework wasproposed for multiscale curl-curl-problems with essential boundary conditions. The LOD tech-nique was first introduced in [32] for elliptic equations in the H -space and was later developed ∗ The authors acknowledge the support of the Brummer&Partners MathDataLab and P. Henning also acknowl-edges the support of the Swedish Research Council (grant 2016-03339) and the G¨oran Gustafsson foundation. Department of Mathematics, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden. a r X i v : . [ m a t h . NA ] M a r to include many other types of equations and spaces, see, e.g., [2, 7, 18, 21, 29, 30, 31, 37, 33]and the references therein. The technique relies on a splitting of the solution space into a coarseand a fine part, using a suitable (quasi-) interpolation operator. The coarse part is typically aclassical finite element space and the fine part contains all the details not captured by the coarsespace. The basis functions of the coarse space are then enhanced by computing so called cor-rectors in the detail space. Due to an exponential decay of these correctors, they are quasi-localand can be computed efficiently by solving small problems. The resulting space (consisting of thecorrected coarse basis functions) can be used in a Galerkin approach to find approximations thatexhibit superior approximation properties in the energy norm when compared with the originalcoarse space. These results hold without assumptions on the structure of the problem, such asperiodicity, scale separation or symmetry.Despite the wide range of successful applications, a construction of the LOD for H (curl)-conforming N´ed´elec finite element spaces turned out to be very challenging. The reason was theabsence of suitable projection operators that are at the heart of the method (since correctionsare always constructed in the kernel of such operators). The required projections need to becomputable, local, H (curl)-stable, and they need to commute with exterior derivatives. Withthe groundbreaking work by Falk and Winther [14], such a projection was finally found andpaved the way for extending the LOD framework to H (curl)-spaces [17]. However, so far, theconstruction by Falk and Winther is only explicitly available for de Rham complexes with naturalboundary conditions, whereas the theoretical results in [17] are obtained for curl-curl-problemswith essential boundary conditions. This discrepancy (which turns out to be a non-trivial issue)prevented a practical implementation of the H (curl)-LOD until now. Hence, the main objectiveof this paper is to explain how the method can indeed be used for practical computations and todemonstrate its performance in numerical experiments. Here we have several goals. First, we showhow the H (curl)-LOD can be extended from essential boundary conditions to natural boundaryconditions. This enables us use the explicitly known Falk-Winther projection constructed in [14],which we used in our implementation. For the arising method we prove that linear convergenceis obtained in the H (curl)-norm whenever the source term fulfills f · n = 0. For general righthand sides f ∈ H (div) we prove convergence of order √ H . These reduced convergence rates forgeneral f originate from the trace of the error on the domain boundary, which dominates theoverall accuracy of the numerical method. This contribution can only be eliminated if the sourceterm has a vanishing normal trace on ∂ Ω. In order to restore a full linear convergence for generalsources f we propose the usage of so-called source term correctors that can be used to make theerror component on ∂ Ω sufficiently small. As it turns out, such source term correctors can be alsoa helpful tool to solve the curl-curl-problem with essential boundary conditions since it allowsthe usage of the same Falk-Winther projection as for the case of natural boundary conditions.With this, we are able to practically solve the same types of problems as considered in [17]. Allthe proposed methods were implemented and we include numerical examples in both 2 D and 3 D to support our theoretical findings. To the authors’ knowledge this is also the first time thatthe edge-based Falk-Winther projection was implemented and the code for the projection (usingthe FEniCS software [28]) is available at [1]. Finally, we shall also demonstrate numerically thatif we use a projection that is locally not commuting with exterior derivatives, then the arising H (curl)-LOD suffers from reduced convergence rates. This shows that the required features ofthe projection are not only artifacts of the proof-technique but are central to obtain a convergingmethod.Another method based on the LOD framework was recently proposed in [38]. A differentprojection is constructed by utilizing the discrete Helmholtz decomposition of the fine detailspace. Numerical experiments show that the proposed projection has the localization property,but an analytical proof of this property remains open.The paper is organized as follows. In Section 2, the curl-curl-problem is introduced and inSection 3 the method for natural boundary conditions is described. Essential boundary conditionsare briefly discussed in Section 4. Finally, numerical experiments are presented in Section 5. Theconstruction of the Falk-Winther projection and a discussion of how to implement it is left to theappendix. Given a computational domain Ω ⊂ R , we consider the curl-curl-problem of the following form.Find u : Ω → C such that curl( µ curl u ) + κ u = f , in Ω ,µ curl u × n = 0 , on Γ , (1)where n is the outward unit normal to Γ := ∂ Ω. In the context of the time-harmonic Maxwell’sequations, u = E describes the (unknown) electric field and the right hand side f is a source termrelated to current densities. The coefficients µ and κ are typically complex-valued and scalarparameter functions that describe material properties. In particular, they depend on electricpermittivity and the conductivity of the medium through which the electric field propagates.We indirectly assume that the propagation medium is a highly variable “multiscale medium”,i.e., the coefficients µ and κ are rapidly varying on a fine scale. In this work we mainly focuson natural boundary conditions (or Neumann-type boundary condition) which are of the form µ curl u × n = 0 on Γ. In the context of Maxwell’s equations this corresponds to a boundarycondition for (a 90 degree rotation of) the tangential component of the magnetic field.In order to state a variational formulation of problem (1), we introduce a few spaces. Forthat purpose, let G ⊆ R denote an arbitrary bounded Lipschitz domain. We define the space of H (curl)-functions on G by H (curl , G ) := { v ∈ L ( G, C ) | curl v ∈ L ( G, C ) } . The space is equipped with the standard inner product( v , w ) H (curl ,G ) := (curl v , curl w ) L ( G ) + ( v , w ) L ( G ) where ( · , · ) L ( G ) is the complex L -inner product. Furthermore, we define the space of H (div)-functions by H (div , G ) := { v ∈ L ( G, C ) | div v ∈ L ( G, C ) } with corresponding inner product( · , · ) H (div ,G ) := (div v , div w ) L ( G ) + ( v , w ) L ( G ) . The restriction of H (div , G ) to functions with a zero normal trace is given by H (div , G ) := { v ∈ H (div , G ) | v · n | ∂G = 0 } . If we drop the dependency on G , we always mean the corresponding space with G = Ω. Withthis we consider the following variational formulation of problem (1): find u ∈ H (curl) such that B ( u , v ) := ( µ curl u , curl v ) L (Ω) + ( κ u , v ) L (Ω) = ( f , v ) L (Ω) , for all v ∈ H (curl) . (2)Observe that the natural boundary condition µ curl u × n = 0 on Γ is intrinsically incorporatedin the variational formulation.For problem (2) to be well-defined, we shall also make a set of assumptions on the data, whichread as follows.(A1) Ω ⊂ R is an open, bounded, contractible polyhedron with a Lipschitz boundary.(A2) The coefficients fulfill µ ∈ L ∞ (Ω , R × ) and κ ∈ L ∞ (Ω , C × ).(A3) The source term fulfills f ∈ H (div).(A4) The coefficients are such that the sesquilinear form B : H (curl) × H (curl) → C given by (2)is elliptic on H (curl), i.e., there exists α > |B ( v , v ) | ≥ α (cid:107) v (cid:107) H (curl) for all v ∈ H (curl) . We assume that this property is independent of the computational domain Ω (i.e., we canrestrict the sesquilinear form to subdomains of Ω and have still ellipticity).Under the above assumptions (A1)-(A4), we have that (2) is well-defined by the Lax-Milgram-Babuˇska theorem [5]. Note that none of the assumptions is very restrictive and that they arefulfilled in many physical setups (cf. [17, 16] for discussions and examples). In [41] it wasalso demonstrated that the coercivity can be relaxed to a frequency-dependent inf-sup-stabilitycondition for indefinite problems. However, in order to avoid the additional technical details thatwould come with the indefinite setting, we will only focus on the simpler H (curl)-elliptic case asgiven by (A4). In this section, we follow the arguments presented in [17] for the case of essential boundaryconditions (i.e., u × n = 0 on Γ) to transfer them to the case of natural boundary conditions asgiven in (1). Let us start with some basic notation that is used for a conventional discretization of (1). Thissetting will be also used to exemplify why such a standard discretization of multiscale curl-curl-problems is problematic.Throughout this paper, we make the following assumptions on the discretization.(B1) T H is a shape-regular and quasi-uniform partition of Ω, in particular we have ∪T H = Ω.(B2) The elements of T H are (closed) tetrahedra and are such that any two distinct elements T, T (cid:48) ∈ T H are either disjoint or share a common vertex, edge or face.(B3) The mesh size (i.e., the maximum diameter of an element of T H ) is denoted by H .On T H we consider lowest-order H -, H (div)- and H (curl)-conforming elements. Here, S ( T H ) ⊂ H (Ω) denotes the space of T H -piecewise affine and continuous functions (i.e., scalar-valued first-order Lagrange finite elements) and we set ˚ S ( T H ) := S ( T H ) ∩ H (Ω) as the subset with a zerotrace on ∂ Ω. The lowest order N´ed´elec finite element (cf. [34, Section 5.5]) is denoted by N ( T H ) := { v ∈ H (curl) |∀ T ∈ T H : v | T ( x ) = a T × x + b T with a T , b T ∈ C } . Later on, we will also need the lowest-order Raviart–Thomas finite element space, which is givenby RT ( T H ) := { v ∈ H (div) |∀ T ∈ T H : v | T ( x ) = a T x + b T with a T ∈ C , b T ∈ C } . With this, a standard H (curl)-conforming discretization of problem (2) is to find u H ∈ N ( T H )such that B ( u H , v H ) = ( f , v H ) L (Ω) , for all v ∈ N ( T H ) . (3)As a Galerkin approximation, u H is a quasi-best approximation in the H (curl)-norm and fulfillsthe classical estimate (cid:107) u − u H (cid:107) H (curl) ≤ C inf v H ∈N ( T H ) (cid:107) u − v H (cid:107) H (curl) ≤ CH (cid:0) (cid:107) u (cid:107) H (Ω) + (cid:107) curl u (cid:107) H (Ω) (cid:1) . (4)A serious issue is faced if µ and κ are multiscale coefficients that are rapidly oscillating with acharacteristic length scale of the oscillations that is of order δ (cid:28)
1. In this case, both (cid:107) u (cid:107) H (Ω) and (cid:107) curl u (cid:107) H (Ω) will blow up, even for differentiable coefficients µ and κ . Even worse, if µ and κ are discontinuous, then the required regularity u , curl u ∈ H (Ω , C ) might not be available atall (cf. [11, 12, 6]). In such a setting the estimate (4) becomes meaningless, as convergence ratescan become arbitrarily slow. On top of that, convergence can only be observed in the fine scaleregime H < δ . This imposes severe restrictions on the mesh resolution and the costs for solving(3) can become tremendous. It is worth to mention that the picture does not change if we look atthe error in the L -norm. This is because (cid:107) u − u H (cid:107) L (Ω) and (cid:107) u − u H (cid:107) H (curl) are typically of thesame order, due to the large kernel of the curl-operator. In particular, for large mesh sizes H thespace N ( T H ) does not contain any good L -approximations of u . This is because fast variationsin κ cause equally fast oscillations in u with an amplitude of order O (1). These oscillationscannot be captured on coarse meshes and due to their non-negligible amplitude they are crucialin the L -sense. For periodic coefficients, the appearance of such oscillations can be explainedwith analytical homogenization theory (cf. [20]). An example for a function u that illustrates thisaspect is given in Figure 2 below. Note that this is a crucial difference to H -elliptic problems,where good L -approximations on coarse meshes are at least theoretically available (though theyare typically not found by standard Galerkin methods). A very good review about this topic for H -elliptic problems is given in [36]. As discussed in the previous subsection, a conventional discretization of the curl-curl-problem(2) is problematic as it requires very fine meshes T H and hence very high-dimensional discretespaces. Motivated by this issue, the question posed in [17] is whether it is possible to construct a quasi-local sub-scale correction operator K that acts on H (curl) and which only depends on B ( · , · )(but not on the source term f ), such that the original multiscale problem can be replaced by anumerically homogenized equation. By “numerically homogenized equation” we mean an equationthat can be discretized with a Galerkin method in the standard (coarse) space N ( T H ), such thatthe corresponding (homogenized) solution u H ∈ N ( T H ) is a good coarse scale approximationof u (in the dual space norm (cid:107) · (cid:107) H (div) (cid:48) , cf. [17] for details) and such that u H + K ( u H ) (i.e.,homogenized solution plus corrector) yields a good approximation of u in the H (curl)-norm. The(numerically) homogenized equation takes the form B (( I + K ) u H , ( I + K ) v ) = ( f , ( I + K ) v ) L (Ω) for all v ∈ N ( T H )and provided that K is available, this can be solved cheaply in a coarse (and hence low-dimensional)space N ( T H ). In [17], such a homogenized problem was obtained for essential boundary conditionsusing the technique of Localized Orthogonal Decompositions (LOD).The central tool in the construction of a suitable corrector operator K is the existence of localand stable projection operators that commute with the exterior derivative. For that we applythe theory of finite element exterior calculus (FEEC), where we refer to the excellent surveypapers [3, 4] for an introduction to the topic. A cochain complex of vector spaces is essentiallya sequence of homomorphisms between these spaces with the property that the image of each ofthe homomorphisms is a subset of the kernel of the next homomorphism. If the vector spaces arethe spaces of differential p -forms on some smooth manifold and if the homomorphisms are givenby the exterior derivative, we call the cochain complex a de Rham complex . In the following, weconsider the L de Rham complex with natural boundary conditions which is of the form H (Ω) grad −→ H (curl , Ω) curl −→ H (div , Ω) div −→ L (Ω) . Obviously, we have for the exterior derivative curl ◦ grad = and div ◦ curl = 0, which shows thecentral property of a cochain complex (that is the inclusion of the image of each homomorphisminto the kernel of the next). We also consider the corresponding finite element subcomplex (oflowest order spaces) S ( T H ) grad −→ N ( T H ) curl −→ RT ( T H ) div −→ P ( T H ) , where P ( T H ) ⊂ L (Ω) is the space of T H -piecewise constant functions. In this setting (and inmore general form), Falk and Winther [14] proved the existence of stable and local projections π H between any Hilbert space of the de Rham complex and the corresponding discrete space ofthe finite element subcomplex, so that the following diagram commutes: H (Ω) grad −→ H (curl , Ω) curl −→ H (div , Ω) div −→ L (Ω) − → π VH − → π EH − → π FH − → π TH S ( T H ) grad −→ N ( T H ) curl −→ RT ( T H ) div −→ P ( T H ) . Here, V stands for vertex, E stands for edges, F for faces and T for tetrahedra. In particular(and for our setting very essential) is the commutation π EH ( ∇ v ) = ∇ π VH ( v ) for all v ∈ H (Ω) . (5)The operators are local (i.e., local information can only spread in small nodal environments) andthey admit local stability estimates, for T ∈ T H , of the form (cid:107) π H ( v ) (cid:107) L ( T ) (cid:46) (cid:107) v (cid:107) L ( U ( T )) + H | v | H Λ( U ( T )) and | π H ( v ) | H Λ( T ) (cid:46) | v | H Λ( U ( T )) (6)for all v ∈ H Λ(Ω), where H Λ(Ω) stands for either of the spaces H (Ω), H (curl , Ω) or H (div , Ω),the semi-norm is given by | v | H Λ = (cid:107) dv (cid:107) L where d stands for the exterior derivative on H Λ(Ω)and π H stands for the corresponding projection operator from the commuting diagram. Theconstant C > U ( K ) is a small environment of K thatconsists of elements from T H and which has a diameter of order H . We refer the reader to [14,Theorem 2.2] for further details on this aspect.The idea proposed in [17] is now to correct/enrich the conventional discrete spaces Λ H =image( π H ) by information from the kernel of the projections. Due to the direct and stabledecomposition H Λ(Ω) = Λ H (Ω) ⊕ W (Ω) , where W (Ω) := kern( π H ), we know that the exact solution u to the curl-curl-problem (2) canbe uniquely decomposed into a coarse contribution in Λ H (Ω) = N ( T H ) and a fine (or detail)contribution in W := W (Ω) = kern( π EH ), i.e., u = π EH ( u ) + ( u − π EH ( u )) . (7)In order to identify and characterize the components of this decomposition, a Green’s correctoroperator G : H (curl) (cid:48) → W can be introduced. For F ∈ H (curl) (cid:48) it is defined by the equation B ( G ( F ) , w ) = F ( w ) for all w ∈ W . (8)It is well-defined by the Lax-Milgram-Babuˇska theorem due to assumption (A4) and since W isa closed subspace as the kernel of a H (curl)-stable projection. The main feature of the Green’scorrector is that it allows us to characterize the H (curl)-stable decomposition of the exact solution u ∈ H (curl) to problem (2) as u = u H − ( G ◦ L )( u H ) + G ( f ) , (9)where u H := π EH ( u ) ∈ N ( T H ) is the coarse part. The operator L is the differential operatorassociated with B ( · , · ), i.e., L ( v ) := B ( v , · ).To see that (9) holds, we start from the unique decomposition (7) and insert it into (2), wherewe restrict the test functions to W ⊂ H (curl). This yields B ( u − π EH ( u ) , w ) = −B ( π EH ( u ) , w ) + ( f , w ) L (Ω) = −B (( G ◦ L )( u H ) , w ) + B ( G ( f ) , w )for all w ∈ W . Since both u − π EH ( u ) and G ( f ) − ( G ◦ L )( u H ) are elements of W , they must beidentical, proving (9). In the last subsection we saw that the exact solution can be decomposed into u = u H − ( G ◦L )( u H ) + G ( f ), where u H := π EH ( u ). This motivates to define the ideal corrector operator K : N ( T H ) → W as K := −G ◦ L and to consider a numerical homogenized equation of the form: find u H ∈ N ( T H ) such that B (( I + K ) u H , ( I + K ) v ) = ( f , ( I + K ) v ) L (Ω) for all v ∈ N ( T H ) . (10)Since this is a Galerkin approximation, C´ea’s lemma implies that ( I + K ) u H is an H (curl)-quasibest approximation of u in ( I + K ) N ( T H ) and hence (cid:107) u − ( I + K ) u H (cid:107) H (curl) (cid:46) (cid:107) u − ( I + K ) u H (cid:107) H (curl) (9) = (cid:107)G ( f ) (cid:107) H (curl) . In fact, if µ and κ are self-adjoint, it can be shown that u H = π EH ( u ) and hence the error has theexact characterization u − ( I + K ) u H = G ( f ) (this can be proved analogously to [17, Theorem10]).As we just saw, in order to quantify the discretization error e H = u − ( I + K ) u H , all we haveto do is to estimate G ( f ). Since Ω is a contractible domain, we know that G ( f ) ∈ H (curl) admitsa regular decomposition (cf. [22, 23]) of the form G ( f ) = z + ∇ θ, where z ∈ [ H (Ω)] and θ ∈ H (Ω) , (11)and with (cid:107) z (cid:107) H (Ω) + (cid:107) θ (cid:107) H (Ω) (cid:46) (cid:107)G ( f ) (cid:107) H (curl) . (12)Exploiting the commuting property of the Falk-Winther projections, we obtain that we can write G ( f ) ∈ W as G ( f ) = G ( f ) − π EH ( G ( f )) (5) = ( z − π EH ( z )) + ∇ ( θ − π VH ( θ )) . (13)By the Bramble-Hilbert lemma together with the local stability estimates (6) we easily see thatfor any T ∈ T H it holds (cid:107) z − π EH ( z ) (cid:107) L ( T ) ≤ CH (cid:107) z (cid:107) H ( U ( T )) (14)and analogously (cid:107) θ − π VH ( θ ) (cid:107) L ( T ) ≤ CH (cid:107) θ (cid:107) H ( U ( T )) . (15)Again, U ( T ) is a generic environment of T that consists of few (i.e., O (1)) elements from T H .Observe that this yields a local regular decomposition of interpolation errors in the spirit ofSch¨oberl (cf. [39, Theorem 1] and [17, Lemma 4]).Next, we use the identity (13) in the definition of the Green’s corrector G to obtain α (cid:107)G ( f ) (cid:107) H (curl) ≤ B ( G ( f ) , G ( f )) = ( f , z − π EH ( z )) + ( f , ∇ ( θ − π VH ( θ ))) L (Ω) (16)= ( f , z − π EH ( z )) − (div f , θ − π VH ( θ )) L (Ω) + ( f · n , θ − π VH ( θ )) L (Γ) . Using the interpolation error estimates (14) and (15) for the Falk-Winther projections, we imme-diately have ( f , z − π EH ( z )) (cid:46) H (cid:107) f (cid:107) L (Ω) (cid:107) z (cid:107) H (curl) (12) (cid:46) H (cid:107) f (cid:107) L (Ω) (cid:107)G ( f ) (cid:107) H (curl) . (17)and − (div f , θ − π VH ( θ )) L (Ω) (cid:46) H (cid:107) div f (cid:107) L (Ω) (cid:107) θ (cid:107) H (Ω) (12) (cid:46) H (cid:107) div f (cid:107) L (Ω) (cid:107)G ( f ) (cid:107) H (curl) . (18)It remains to estimate the last term. For that, let F H, Γ denote the set of faces (from our trian-gulation) that are located on the domain boundary Γ. We obtain with the local trace inequality | ( f · n , θ − π VH ( θ )) L (Γ) | ≤ (cid:88) S ∈F H, Γ (cid:107) f · n (cid:107) L ( S ) (cid:107) θ − π VH ( θ ) (cid:107) L ( S ) ≤ (cid:107) f · n (cid:107) L (Γ) (cid:88) S ∈F H, Γ (cid:107) θ − π VH ( θ ) (cid:107) L ( S ) / (cid:46) (cid:107) f · n (cid:107) L (Γ) (cid:88) T ∈T H T ∩ Γ ∈F H, Γ H − (cid:107) θ − π VH ( θ ) (cid:107) L ( T ) + H (cid:107)∇ θ − ∇ π VH ( θ ) (cid:107) L ( T ) / , (15) (cid:46) (cid:107) f · n (cid:107) L (Γ) (cid:88) T ∈T H T ∩ Γ ∈F H, Γ H (cid:107)∇ θ (cid:107) L ( U ( T )) + H (cid:107)∇ θ (cid:107) L ( U ( T )) / (cid:46) √ H (cid:107) f · n (cid:107) L (Γ) (cid:107)G ( f ) (cid:107) H (curl) . (19)Here we used the regularity of the mesh T H and that each T -neighborhood U ( T ) contains only O (1) elements. Combining the estimates (16), (17), (18) and (19), we obtain the followingconclusion. Conclusion 3.1.
Let u ∈ H (curl) denote the exact solution to curl-curl-problem (3) with naturalboundary conditions and let u H ∈ N ( T H ) denote the solution to ideal numerically homogenizedequation (10), then it holds the error estimate (cid:107) u − ( I + K ) u H (cid:107) H (curl) (cid:46) H (cid:107) f (cid:107) H (div) + √ H (cid:107) f · n (cid:107) L (Γ) . In particular, if f ∈ H (div) (i.e., it has a vanishing normal trace), then we have optimal orderconvergence with (cid:107) u − ( I + K ) u H (cid:107) H (curl) (cid:46) H (cid:107) f (cid:107) H (div) . From Conclusion 3.1 we see that in general we expect the error to be dominated by the term √ H (cid:107) f · n (cid:107) L (Γ) . Note that in the case of essential boundary conditions, the problematic term( f · n , θ − π VH ( θ )) L (Γ) will always vanish since the corresponding regular decomposition guaranteesthat θ = 0 on Γ (cf. [22, 39]).Before we discuss how to restore the linear convergence for source terms with arbitrary normaltrace, we need to discuss how to localize the corrector operator K . From a computational point of view, the ideal multiscale method (10) is not feasible, as it requiresthe computation of global correctors K ( ψ E ) for each nodal basis function ψ E of the N´ed´elec space N ( T H ). The equation that describes a corrector K ( ψ E ) ∈ W reads B ( K ( ψ E ) , w ) = −B ( ψ E , w ) for all w ∈ W . This problem has the same computational complexity as the original curl-curl-problem. Solvingit for every basis function would hence multiply the original computing costs. To avoid this issue,we first split the global corrector K into a set of element correctors K T . For a tetrahedron T ∈ T H ,the corresponding element corrector K T ( v H ) ∈ W is defined as B ( K T ( v H ) , w ) = −B T ( v H , w ) for all w ∈ W , (20)where B T ( v , w ) := ( µ curl v , curl w ) L ( T ) + ( κ v , w ) L ( T ) . Obviously, we have K ( v H ) = (cid:88) T ∈T H K T ( v H ) . (21)The advantage of this formulation is that we can now exploit the strong localization of theelement correctors K T . In fact, they show an extremely fast decay to zero away from the element T . The decay can be quantified and is exponential in units of H . This is a consequence of thelocal support of the source term (i.e., B T ( · , · )) and the fact that K T maps into the kernel ofthe Falk-Winther projection. This allows to restrict the computational domain in (20) to smallenvironments N m ( T ) of T that have a diameter of order m · H . The decay was first proved in[17, Theorem 14] for the case of essential boundary conditions (i.e., for K ( ψ E ) × n = 0 on Γ).For natural boundary conditions we can proceed similarly (with minor modifications) to obtainthe following result. Lemma 3.2.
Given T ∈ T H , we define the (open) m -layer neighborhood recursively by N ( T ) := int( T ) and N m ( T ) := int (cid:16)(cid:91) { K ∈ T H | K ∩ N m − ( T ) (cid:54) = ∅} (cid:17) . On these patches (for m > ) we define the truncated element correctors K T,m ( v H ) ∈ W ( N m ( T ) ) := { w ∈ W | w ≡ \ N m ( T ) } as solutions to B ( K T,m ( v H ) , w ) = −B T ( v H , w ) for all w ∈ W ( N m ( T ) ) (22) For v H ∈ N ( T H ) , an approximation of the global corrector is given (in the spirit of (21) ) by K m ( v H ) := (cid:88) T ∈T H K T,m ( v H ) . The error between the truncated approximation K m and the ideal corrector K can be bounded by (cid:107)K ( v H ) − K m ( v H ) (cid:107) H (curl) ≤ Cρ m (cid:107) v H (cid:107) H (curl) , where < ρ < and C > are generic constants that are independent of m , H or the speed ofthe variations of µ and κ . Remark 3.3 (Boundary conditions of the truncated element correctors) . Observe that the trun-cated element correctors K T,m ( v H ) given by (22) have mixed boundary conditions. On ∂ N m ( T ) \ Γthey fulfill essential boundary conditions, i.e., K T,m ( v H ) × n = 0, and on ∂ N m ( T ) ∩ Γ they ful-fill natural boundary conditions, i.e., µ curl K T,m ( v H ) × n = 0. This is in contrast to the idealelement correctors given by (20) which only involve natural boundary conditions.From Lemma 3.2 we see that that we can approximate the global corrector K efficiently bycomputing a set of element correctors K T,m with a small local support N m ( T ). With this wecan formulate the following main theorem on the overall accuracy of localized multiscale scheme(which belongs to the class of LOD methods, cf. [17, 19, 32]). Theorem 3.4.
Let u ∈ H (curl) denote the exact solution to the curl-curl-problem with naturalboundary conditions. Furthermore, let K m denote the corrector approximation from Lemma 3.2and let u H,m ∈ N ( T H ) denote the unique solution to the numerically homogenized equation B (( I + K m ) u H,m , ( I + K m ) v H ) = ( f , ( I + K m ) v H ) L (Ω) for all v H ∈ N ( T H ) . Then the following error estimates holds (cid:107) u − ( u H,m + K m ( u H,m )) (cid:107) H (curl) (cid:46) ( H + ρ m ) (cid:107) f (cid:107) H (div) + √ H (cid:107) f · n (cid:107) L (Γ) . (23) In particular, if m (cid:38) | log( H ) | / | log( ρ ) | and f ∈ H (div) , then we have optimal order convergencein H , i.e., it holds (cid:107) u − ( u H,m + K m ( u H,m )) (cid:107) H (curl) (cid:46) H (cid:107) f (cid:107) H (div) . Proof.
The proof is analogous to [17, Conclusion 15]. Let u H = π EH ( u ) ∈ N ( T H ) be the Falk-Winther projection of the exact solution. Then we have with Lemma 3.2 (cid:107) u − u H − K m ( u H ) (cid:107) H (curl) (9) ≤ (cid:107)G ( f ) (cid:107) H (curl) + (cid:107)K ( u H ) − K m ( u H ) (cid:107) H (curl) ≤ (cid:107)G ( f ) (cid:107) H (curl) + Cρ m (cid:107) π EH ( u ) (cid:107) H (curl) (6) ≤ (cid:107)G ( f ) (cid:107) H (curl) + Cρ m (cid:107) u (cid:107) H (curl) ≤ (cid:107)G ( f ) (cid:107) H (curl) + Cρ m (cid:107) f (cid:107) H (div) . Conclusion 3.1 together with the quasi-best approximation property of Galerkin methods finishesthe proof.1Theorem 3.4 shows that we can still achieve the same order of accuracy as for the idealmultiscale method (10) if the number of layers m is selected proportional to | log( H ) | . Practicallythis means that in order to compute a sufficiently accurate corrector K m , we need to solvelocal problems on small patches with a diameter of order H | log( H ) | . The number of theselocal problems equals the number of all tetrahedra times the number of edges of a tetrahedron(which is equal to the number of basis functions of N ( T H ) that have support on a tetrahedron T ∈ T H ). Hence, the total number of local problems is 6 · (cid:93) T H . All these problems are small,cheap to solve and they are fully independent from each other, which allows for parallelization ina straightforward way.There are two more issues to be discussed, first, we need a computational representation ofthe Falk-Winther projection operator π EH . As this summarizes the results of the previous works[14, 17] and in particular [15], we postpone it to the appendix. The second issue that remains isthe question if and how we can restore a full linear convergence for the error in (23), if f does nothave a vanishing normal trace. We will investigate this question in the next subsection. In Theorem 3.4 we have seen that we can in general not expect a full linear convergence in themesh size H , unless f ∈ H (div). The estimate (23) suggests a convergence of order H / , whichwe can also confirm numerically (cf. Section 5.2). Hence, as the full linear rate is typicallynot achieved for arbitrary source terms, this poses the question if we can modify the methodaccordingly.One popular way to handle the influence of dominating source terms or also general (inhomo-geneous) boundary conditions is to compute source and boundary corrections (cf. [19] for detailsin the context of H -elliptic problems). Here we recall decomposition (9), which says that we canwrite the exact solution as u = u H + K ( u H ) + G ( f ) . Practically, we can approximate G ( f ) in the spirit of Lemma 3.2 by a corrector G m ( f ) based onlocalization. This means, that for each T ∈ T H , we can solve for an element source corrector G T,m ( f ) ∈ W ( N m ( T ) ) with B ( G T,m ( f ) , w ) = ( f , w ) L ( T ) for all w ∈ W ( N m ( T ) ) . The global source corrector is then given by G m ( f ) := (cid:88) T ∈T H G T,m ( f ) . After this, we can solve the corrected homogenized problem which is of the following form: find u corr H,m ∈ N ( T H ) such that B (( I + K m ) u corr H,m , ( I + K m ) v H ) = ( f , ( I + K m ) v H ) L (Ω) − B ( G m ( f ) , ( I + K m ) v H )for all v H ∈ N ( T H ). The final approximation is given by u ms H,m := ( I + K m ) u corr H,m + G m ( f )and it is possible to prove that the error converges exponentially fast in m to zero: (cid:107) u − u ms H,m (cid:107) H (curl) (cid:46) ρ m (cid:107) f (cid:107) L (Ω) , for a generic 0 < ρ < . The error estimate can be proved using similar arguments as in the proof of Theorem 3.4. Herewe recall that ρ is independent of m and H . By selecting m (cid:38) k | log( H ) | / | log( ρ ) | , we can achievethe convergence rate H k for any k ≥
1. For very large m (i.e N m ( T ) = Ω) the method is exact.2 Remark 3.5.
Note that this strategy can be also valuable for handling general non-homogenousboundary conditions. For example, considering the boundary condition µ curl u × n = g on Γ(which can be straightforwardly incorporated into the variational formulation), it is possible tocompute boundary correctors to increase the accuracy. This works analogously to the case ofsource correctors, i.e., one needs to compute element boundary correctors G E ( g ) ∈ W for eachboundary edge E ⊂ Γ. The edge boundary correctors solve problems of the form B ( G E ( g ) , w ) = ( g , w ) L ( E ) for all w ∈ W . Though the above strategy can potentially yield any degree of accuracy, it involves a com-putational overhead, since it additionally requires to compute the whole set of element sourcecorrectors G T,m ( f ). One way to reduce this overhead is to only compute correctors G T,m ( f ) forwhich T is close to the boundary. For example, we can set T in H := { T ∈ T H | N m ( T ) ∩ Γ (cid:54) = ∅} and T out H := T H \ T in H . In this case we have for the error u − (cid:0) u H + K m ( u H ) + (cid:88) T ∈T out H G T,m ( f ) (cid:1) = ( K m − K ) u H + (cid:88) T ∈T in H ( G T,m − G T ) f + (cid:88) T ∈T out H ( G T,m − G T ) f − (cid:88) T ∈T in H G T,m ( f ) . Here, the first three terms show an exponential convergence with order ρ m , whereas the last term isnot polluted by the boundary influence and converges with order H , since ( (cid:80) T ∈T in H G T,m ( f )) | Γ = 0. In the following we want to consider the curl-curl-problem with essential boundary conditions asoriginally done in [17], i.e., we seek u : Ω → C such thatcurl( µ curl u ) + κ u = f , in Ω , u × n = 0 , on ∂ Ω . (24)For the variational formulation we denote the restriction of H (curl) to functions with a vanishingnormal trace on Γ by H (curl) := { v ∈ H (curl , Ω) | v × n | ∂ Γ = 0 } . Hence, the weak formulation of (24) becomes: find u ∈ H (curl) with B ( u , v ) = ( f , v ) L (Ω) for all v ∈ H (curl) . Assuming again (A1)-(A4), this problem is well-posed.To discretize the problem, let the subspace of N ( T H ) which contains functions with a ho-mogenous tangential trace to be given by˚ N ( T H ) := N ( T H ) ∩ H (curl) . With this a multiscale approximation analogously to the case of natural boundary conditions(cf. Section 3) can be obtained by considering Falk-Winther projections between the spaces withvanishing traces (cf. [17, 14]), i.e., local and stable projections˚ π EH : H (curl) → ˚ N ( T H ) and ˚ π VH : H (Ω) → ˚ S ( T H ) := S ( T H ) ∩ H (Ω)3with the commuting property ∇ ˚ π VH ( v ) = ˚ π EH ( ∇ v ) for all v ∈ H (Ω). In this setting, the followingresult can be proved (see [17]):Recall the notation from Section 3.4. For T ∈ T H and m ∈ N > , let the local detail space begiven by ˚ W ( N m ( T ) ) := { w ∈ kern(˚ π EH ) | w ≡ \ N m ( T ) } ⊂ H (curl)and the element correctors ˚ K T,m ( v H ) ∈ ˚ W ( N m ( T ) ) be the solutions to B (˚ K T,m ( v H ) , w ) = B T ( v H , w ) for all w ∈ ˚ W ( N m ( T ) ) . (25)We define the global corrector as usual by ˚ K m ( v H ) := (cid:80) T ∈T H ˚ K T,m ( v H ). With this, the finalmultiscale approximation is given by ˚ u H,m ∈ ˚ N ( T H ) with B (( I + ˚ K m )˚ u H,m , ( I + ˚ K m ) v ) = ( f , ( I + ˚ K m ) v ) for all v ∈ ˚ N ( T H ) (26)and it holds the error estimate (cid:107) u − ( I + ˚ K m )˚ u H,m (cid:107) H (curl) (cid:46) ( H + ρ m ) (cid:107) f (cid:107) H (div) , where 0 < ρ < f . This is due to the regular decomposition of functions v ∈ H (curl),which is of the form v = z + ∇ θ , where z ∈ [ H (Ω)] and θ ∈ H (Ω). Hence, the problematicterm ( f · n , θ ) L (Γ) in (16) (which cause the degeneration of convergence rates) is equal to zero.The main issue with essential boundary conditions is the availability of the required Falk-Winther projection ˚ π EH (in the de Rham complex with vanishing traces). So far, the projectionhas not yet been explicitly constructed, which is also the main reason why the approach of [17]was not yet numerically validated. It is worth to mention that a corresponding constructionwould be similar to the construction in the case of natural boundary conditions (see AppendixA), with the difference that all degrees of freedom on the boundary need to be set to zero in theconstruction (cf. the descriptions in [14] and in [13, Section 7.1]). To our understanding, thisleads to local equations (as part of the construction of the projection) which are formulated on acomplex with mixed boundary conditions, which we however had problems to solve numerically.Hence, an implementation of ˚ π EH remains open, which raises the question if it is possible to workwith the projection π EH on the full space instead.In [17, Section 8] it was suggested to start from π EH and then set all weights to zero thatbelong to edge-based basis functions ψ E for a boundary edge E ⊂ Γ. This means the multiscalemethod is based on the kernel of the operator π ,EH := I ,EH ◦ π EH : H (curl) → ˚ N ( T H ) , where I ,EH : H (curl) → ˚ N ( T H ) is the standard edge-based N´ed´elec interpolation operator. Thischoice works indeed well if f ∈ H (div) has a vanishing normal trace, but it fails for general f (i.e.,if f · n (cid:54) = 0 on Γ). The reason is that π ,EH does no longer fulfill the very important commutingdiagram property in a “one-coarse layer”-environment of Γ. This leads to a local loss of stability.In the arising multiscale method, it can be compensated if f · n = 0. However, in any other casethe stability issues on the boundary cause a pollution of the numerical solution. We added acorresponding numerical experiment in Section 5.4 which shows that the convergence rate is nolonger linear but only O ( H / ).As an alternative we suggest to simply keep on using the same π EH as for the case of natural4boundary conditions and proceed otherwise analogously to (26). Practically, this means that weselect the local detail space for the computation of K T,m (analogously to (25)) as W ( N m ( T ) ) := { w ∈ H (curl) | π EH ( w ) = 0 and w ≡ \ N m ( T ) } . Otherwise, nothing changes and the multiscale approximation is given by u H,m ∈ ˚ N ( T H ) with B (( I + K m ) u H,m , ( I + K m ) v ) = ( f , ( I + K m ) v ) for all v ∈ ˚ N ( T H ) . (27)Since typically π EH ( v ) (cid:54)∈ ˚ N ( T H ) for w ∈ H (curl), the modified formulation introduces an errorat the boundary which reduces the overall convergence rate. This is similar to the case of naturalboundary conditions with f · n (cid:54) = 0. Indeed, the experiments in Section 5.3 and Section 5.4 showthat we get full linear convergence in the case when f · n = 0 and convergence of order √ H when f · n (cid:54) = 0. Again, a full linear rate can be reconstructed by using source term corrections close tothe boundary as described in Section 3.5. This is investigated numerically in Section 5.4. Herewe stress that the use of source correctors to handle the case of essential boundary conditionswith f · n (cid:54) = 0 is currently the only available option for optimal convergence rates. A preferablealternative would be a multiscale construction based on a computable Falk-Winther projectionfor essential boundary conditions as this would not require additional source correctors. However,as explained above, such a projection is currently not available. In this section we present some numerical examples in both 2 D and 3 D to verify the theoreticalfindings. For simplicity, we choose the computational domain to be either the unit square or theunit cube, depending on the dimension of the problem. In 2 D , the considered curl-curl-problemis defined analogously to the 3 D problem given by (2), however, with the difference that the2-dimensional curl-operator maps vector fields to scalar values via curl( v ) := ∂ x v − ∂ x v . TheFalk-Winther projection in 2 D and 3 D is discussed in Appendix A. All other modifications inthe problem and method formulation are straightforward. Even though there are no analyticalresults for the 2 D case, our numerical experiments will show that the convergence behavior isanalogous to the theoretically supported 3 D setting.To compute the error of the method in the examples below, we compute a reference solutionon a finer (uniform) mesh, of size h = √ · − in 2 D and h = √ · − in 3 D , using classical FEMin the lowest order N´ed´elec space N ( T h ). The relative errors are reported in the energy norminduced by the sesquilinear form B ( · , · ). That is, if e denotes the difference between the full LODmultiscale approximation u H,m + K m ( u H,m ) and the reference solution u h , then the relative erroris computed as (cid:32) B ( e , e ) B ( u h , u h ) (cid:33) / = (cid:32) ( µ curl e , curl e ) + ( κ e , e )( µ curl u h , curl u h ) + ( κ u h , u h ) (cid:33) / . To achieve a multiscale character of the problem, we choose the coefficients, µ and κ , to becheckerboards. The values in the cubes (squares in 2 D ) are set to 1 (black) and 0 .
001 (white),respectively. See Figure 1 for a visualization. The size of the checkerboard (i.e., the size of eachconstant block) is always chosen to be twice the size of the reference mesh. This means that thereference solution always resolves the coefficients. A typical reference solution in the 2 D case isplotted in Figure 2, where we clearly see the multiscale structure from the checkerboard.All implementations are made using the FEniCS software [28]. The code for the Falk-Wintherprojection in both 2 D and 3 D can be obtained from [1].5 (a) Two dimensions (b) Three dimensions Figure 1: Checkerboard coefficient for two and three dimensions.Figure 2: Reference solution in 2 D with f ( x, y ) = [sin(2 πx ) , sin(2 πy )] and natural boundaryconditions. In this example we show that we obtain linear convergence of the method for problems that fulfill f · n = 0, as predicted in Theorem 3.4.In 2 D we choose f ( x, y ) = [sin(2 πx ) , sin(2 πy )] and we compute LOD approximations onuniform meshes of size H = √ · − j , j = 0 , ...,
5. Furthermore, we let the number of layers in thepatch N m ( T ) range through m = 1 , , , , , H decreases. For each value of H , the error forthe LOD method is compared to the classical finite element approximation in the N´ed´elec space N ( T H ).In 3 D we choose f ( x, y, z ) = [ − x ( x − z − , , z ( z − x − H = √ · − j , j = 0 , .., T H .6 -1 -2 -1 (a) Two dimensions -1 -2 -1 (b) Three dimensions Figure 3: Relative errors for the LOD (blue ∗ ) and the FEM (red ◦ ) approximations of Example1 plotted against the mesh size H . The dashed line is cH . Here we show that the convergence is of order H / if f · n (cid:54) = 0, in accordance with Theorem 3.4.In 2 D we choose f ( x, y ) = [1 , m , we use the ideal method. That is, we do not use any localization when computing theLOD approximations. In Figure 4 we clearly see that the H / -term dominates the convergence( H / marked by a dashed line). In the plot the error for iteration j = 0 is omitted for the LODapproximation, because the error in this particular setting is close to zero. We also note that,despite the reduced convergence rate, the method still outperforms the classical FEM in terms ofaccuracy relative to the number of degrees of freedom in the solution space.To improve the error we test the source term corrections proposed in Section 3.5. We presentresults in 2 D for two cases: 1. when the correction is computed for all triangles and 2. whenit is only computed for triangles that intersect the outer boundary Γ. In Table 1 we see howthis affects the error for some values of H and m . As expected, computing the corrections for alltriangles in the domain decreases the error drastically. Note that for the ideal method (withoutlocalization) the error would vanish completely. Furthermore, we can see that computing thecorrections only for triangles intersecting the boundary still reduces the error significantly.Essentially, source corrections allow for an arbitrary accuracy depending on the size of thelocalization patches (i.e., the value of m ) and the number of coarse elements on which the sourcecorrection is computed. From our experiments we see that one layer of coarse elements aroundthe boundary is often sufficient to obtain a reasonable accuracy.Table 1: Relative errors for the LOD approximation of Example 2 with and without correctionfor the source term. H = √ · − , m = 2 H = √ · − , m = 3LOD without correction 0.1731 0.1235Corrections for triangles at the boundary 0.101 0.0828Corrections for all triangles 0 . · − . · − -1 Figure 4: Relative errors for the LOD without localization (blue ∗ ) and the FEM (red ◦ ) approx-imations of Example 2 plotted against the mesh size H . The dashed line is cH / . In this example we test the convergence rate for the problem with essential boundary conditions,i.e., problem (24), as described in Section 4 when f · n = 0 is fulfilled. We choose the same settingas in Section 5.1, except that we impose homogeneous essential boundary conditions instead. InFigure 5, we clearly see that linear order convergence is achieved in 2 D . In 3 D linear convergenceis achieved except for the coarsest mesh size. The dashed line is plotted from the second point.The experiment demonstrates that if f · n = 0, then the same Falk-Winther projection can beused for solving both the curl-curl-problem with essential boundary conditions and the curl-curl-problem with natural boundary conditions. In both cases, an optimal linear convergence rate isobtained. -1 -2 -1 (a) Two dimensions -1 -2 -1 (b) Three dimensions Figure 5: Relative errors for the LOD (blue ∗ ) and the FEM (red ◦ ) approximations of Example3 plotted against the mesh size H . The dashed line is cH . In this example we consider again essential boundary conditions, i.e., problem (24). Otherwisewe put ourselves in the same setting as in Section 5.2.8First, we conclude that the convergence rate is of order H / if f · n (cid:54) = 0, see Figure 6. Herewe have used the LOD method without localization to show that the reduced convergence rateis not due to the localization parameter. We also note that the method still outperforms theclassical FEM (in terms of accuracy relative to the number of degrees of freedom), despite thereduced convergence rate.As in Section 5.2 we also try to improve the error by using the source term correctors. Theresults are displayed in Table 2 for some values of H and m . As we can see, the error is improveddrastically when computing the correction for all triangles. If the computations are restricted tothe triangles intersecting the boundary, the error is still improved. -2 -1 -2 -1 Figure 6: Relative errors for the LOD without localization (blue ∗ ) and the FEM (red ◦ ) approx-imations of Example 4 plotted against the mesh size H . The dashed line is cH / .Table 2: Relative errors for the LOD approximation of Example 4 with and without correctionfor the source term. H = √ · − , m = 2 H = √ · − , m = 3LOD without correction 0.186 0.134Corrections for triangles at the boundary 0.117 0.0964Corrections for all triangles 3 . · − . · − Finally, we comment on the operator π ,EH = I ,EH ◦ π EH , see Section 4. This is the operatorwhich puts all weights to zero for the edge-basis functions ψ E corresponding to a boundary edge.In Figure 7, we see that this approach still shows reduced convergence of order √ H , which issimilar to the current approach of just keeping π HE . We conclude that it is not enough to simplyput all weights to zero to achieve the Falk-Winther projection with essential boundary conditions˚ π EH that is required for the results of [17] to hold true (i.e., a linear order convergence in H independent of the normal trace of f ). Acknowledgements.
The authors thank Barbara Verf¨urth for the many fruitful discussions onthe general methodology and Joachim Sch¨oberl for the helpful and valuable comments on theFalk-Winther operator.9 -2 -1 -2 -1 Figure 7: Relative errors for the LOD approximation of Example 4 without localization (blue ∗ )based on π ,EH plotted against the mesh size H . The dashed line is cH / . References [1] Falk-Winther projection code source available at GitHub. https://github.com/annaper3/falk-winther-projection .[2] A. Abdulle and P. Henning. Localized orthogonal decomposition method for the wave equa-tion with a continuum of scales.
Math. Comp. , 86(304):549–587, 2017.[3] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homologicaltechniques, and applications.
Acta Numer. , 15:1–155, 2006.[4] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus: from Hodgetheory to numerical stability.
Bull. Amer. Math. Soc. (N.S.) , 47(2):281–354, 2010.[5] I. Babuˇska. Error-bounds for finite element method.
Numer. Math. , 16:322–333, 1971.[6] A. Bonito, J.-L. Guermond, and F. Luddens. Regularity of the Maxwell equations in het-erogeneous media and Lipschitz domains.
J. Math. Anal. Appl. , 408(2):498–512, 2013.[7] D. L. Brown, D. Gallistl, and D. Peterseim. Multiscale Petrov-Galerkin method for high-frequency heterogeneous Helmholtz equations. In
Meshfree methods for partial differentialequations VIII , volume 115 of
Lect. Notes Comput. Sci. Eng. , pages 85–115. Springer, Cham,2017.[8] L. Cao, Y. Zhang, W. Allegretto, and Y. Lin. Multiscale asymptotic method for Maxwell’sequations in composite materials.
SIAM J. Numer. Anal. , 47(6):4257–4289, 2010.[9] P. Ciarlet, S. Fliss, and C. Stohrer. On the approximation of electromagnetic fields by edgefinite elements. Part 2: A heterogeneous multiscale method for Maxwell’s equations.
Comput.Math. Appl. , 73(9):1900–1919, 2017.[10] P. Ciarlet Jr. and C. Stohrer. Finite element heterogeneous multiscale method for the classicalHelmholtz equation.
C. R. Acad. Sci. , 352:755–760, 2014.[11] M. Costabel. A remark on the regularity of solutions of Maxwell’s equations on Lipschitzdomains.
Math. Methods Appl. Sci. , 12(4):365–368, 1990.0[12] M. Costabel, M. Dauge, and S. Nicaise. Singularities of Maxwell interface problems.
M2ANMath. Model. Numer. Anal. , 33(3):627–649, 1999.[13] A. Demlow. Convergence and quasi-optimality of adaptive finite element methods for har-monic forms.
Numer. Math. , 136(4):941–971, 2017.[14] R. S. Falk and R. Winther. Local bounded cochain projections.
Math. Comp. , 83(290):2631–2656, 2014.[15] R. S. Falk and R. Winther. Double complexes and local cochain projections.
Numer. MethodsPartial Differential Equations , 31(2):541–551, 2015.[16] P. Fernandes and M. Raffetto. Existence, uniqueness and finite element approximation of thesolution of time-harmonic electromagnetic boundary value problems involving metamaterials.
COMPEL , 24(4):1450–1469, 2005.[17] D. Gallistl, P. Henning, and B. Verf¨urth. Numerical homogenization of H (curl)-problems. SIAM J. Numer. Anal. , 56(3):1570–1596, 2018.[18] F. Hellman, P. Henning, and A. M˚alqvist. Multiscale mixed finite elements.
Discrete Contin.Dyn. Syst. Ser. S , 9(5):1269–1298, 2016.[19] P. Henning and A. M˚alqvist. Localized orthogonal decomposition techniques for boundaryvalue problems.
SIAM J. Sci. Comput. , 36(4):A1609–A1634, 2014.[20] P. Henning, M. Ohlberger, and B. Verf¨urth. A new heterogeneous multiscale method fortime-harmonic Maxwell’s equations.
SIAM J. Numer. Anal. , 54(6):3493–3522, 2016.[21] P. Henning and A. Persson. A multiscale method for linear elasticity reducing Poissonlocking.
Comput. Methods Appl. Mech. Engrg. , 310:156–171, 2016.[22] R. Hiptmair. Finite elements in computational electromagnetism.
Acta Numer. , 11:237–339,2002.[23] R. Hiptmair and J. Xu. Nodal auxiliary space preconditioning in H ( curl ) and H (div) spaces. SIAM J. Numer. Anal. , 45(6):2483–2509, 2007.[24] M. Hochbruck, B. Maier, and C. Stohrer. Heterogeneous multiscale method for Maxwell’sequations.
Multiscale Model. Simul. , 17(4):1147–1171, 2019.[25] M. Hochbruck and C. Stohrer. Finite element heterogeneous multiscale method for time-dependent Maxwell’s equations. In
Spectral and high order methods for partial differentialequations—ICOSAHOM 2016 , volume 119 of
Lect. Notes Comput. Sci. Eng. , pages 269–281.Springer, Cham, 2017.[26] A. Lamacz and B. Schweizer. A negative index meta-material for Maxwell’s equations.
SIAMJ. Math. Anal. , 48(6):4155–4174, 2016.[27] R. Lipton and B. Schweizer. Effective Maxwell’s equations for perfectly conducting split ringresonators.
Arch. Ration. Mech. Anal. , 229(3):1197–1221, 2018.[28] A. Logg, K.-A. Mardal, G. N. Wells, et al.
Automated Solution of Differential Equations bythe Finite Element Method . Springer, 2012.[29] A. M˚alqvist and A. Persson. A generalized finite element method for linear thermoelasticity.
ESAIM Math. Model. Numer. Anal. , 51(4):1145–1171, 2017.1[30] A. M˚alqvist and A. Persson. Multiscale techniques for parabolic equations.
Numer. Math. ,138(1):191–217, 2018.[31] A. M˚alqvist, A. Persson, and T. Stillfjord. Multiscale differential Riccati equations for linearquadratic regulator problems.
SIAM J. Sci. Comput. , 40(4):A2406–A2426, 2018.[32] A. M˚alqvist and D. Peterseim. Localization of elliptic multiscale problems.
Math. Comp. ,83(290):2583–2603, 2014.[33] A. M˚alqvist and D. Peterseim. Generalized finite element methods for quadratic eigenvalueproblems.
ESAIM Math. Model. Numer. Anal. , 51(1):147–163, 2017.[34] P. Monk.
Finite element methods for Maxwell’s equations . Numerical Mathematics andScientific Computation. Oxford University Press, New York, 2003.[35] M. Ohlberger and B. Verf¨urth. A new heterogeneous multiscale method for the Helmholtzequation with high contrast.
Multiscale Model. Simul. , 16(1):385–411, 2018.[36] D. Peterseim. Variational multiscale stabilization and the exponential decay of fine-scale cor-rectors. In
Building bridges: connections and challenges in modern approaches to numericalpartial differential equations , volume 114 of
Lect. Notes Comput. Sci. Eng. , pages 341–367.Springer, [Cham], 2016.[37] D. Peterseim. Eliminating the pollution effect in Helmholtz problems by local subscalecorrection.
Math. Comp. , 86(305):1005–1036, 2017.[38] X. Ren, A. Hannukainen, and A. Belahcen. Homogenization of multiscale eddy currentproblem by localized orthogonal decomposition method.
IEEE Transactions on Magnetics ,55(9):1–4, Sep. 2019.[39] J. Sch¨oberl. A posteriori error estimates for Maxwell equations.
Math. Comp. , 77(262):633–649, 2008.[40] B. Schweizer. Resonance meets homogenization: construction of meta-materials with aston-ishing properties.
Jahresber. Dtsch. Math.-Ver. , 119(1):31–51, 2017.[41] B. Verf¨urth. Numerical homogenization for indefinite H(curl)-problems. In J. U. K. Mikula,D. Sevcovic, editor,
Proceedings of Equadiff 2017 conference , pages 137–146, 2017.[42] B. Verf¨urth. Heterogeneous multiscale method for the Maxwell equations with high contrast.
ESAIM Math. Model. Numer. Anal. , 53(1):35–61, 2019.
A Construction of the Falk-Winther operator
In the following we describe the implementation of the Falk-Winther projection operator π EH : H (curl) → N ( T H ) as constructed in [15, 14]. Here we note that a discussion of its implementationcan be already partially found in [17], which however contains some minor mistakes, which arefixed in the following descriptions. Our implementation is discussed for the 3 D case, but we alsocomment on what needs to be changed for the 2 D case.Adapting the notation from the aforementioned references, we let ∆ ( T H ) denote the set ofvertices of T H and ∆ ( T H ) is the set of edges. The space N ( T H ) is spanned by the edge-orientedbasis functions ( ψ E ) E ∈ ∆ ( T H ) that are uniquely defined for any E ∈ ∆ ( T H ) through the property ˆ E ψ E · t E ds = 1 and ˆ E (cid:48) ψ E · t E ds = 0 for all E (cid:48) ∈ ∆ ( T H ) \ { E } . t E denotes the unit tangent to the edge E with a globally fixed sign. In the following wedenote by y ( E ) and y ( E ) the endpoints of E , where we make the orientation convention that t E = ( y ( E ) − y ( E )) / | E | , where | E | := length( E ) . The (edge-based) Falk-Winter projection π EH : H (curl) → N ( T H ) is of the form π EH u := S u + (cid:88) E ∈ ∆ ( T H ) ˆ E (cid:0) ( I − S ) Q E u (cid:1) · t E ds ψ E . (28)Here, S is constructed in such a way that π EH commutes with the nodal Falk-Winther projection(i.e., π EH ( ∇ v ) = ∇ π VH ( v ) for all v ∈ H (Ω)). However, since S is not a projection, the second terminvolving the operator Q E needs to be introduced. This ensures that π EH is invariant on N ( T H ).We shall describe the construction of the various constituents in the following subsections, wherewill also comment on their practical realization. A.1 Construction of π VH For a better understanding of the construction, we start with describing the nodal Falk-Wintherprojection π VH : H (Ω) → S ( T H ). For that, we denote for each vertex y ∈ ∆ ( T H ) the associatednodal patch by ω y := int (cid:16) (cid:91) { T ∈ T H : y ∈ T } (cid:17) . With this, we also define the piecewise constant function z y by z y = (cid:40) (meas( ω y )) − in ω y \ ω y (29)Furthermore, the restriction of S ( T H ) to ω y is denoted by S ( T H ( ω y )) (i.e., scalar-valued first-orderLagrange finite elements on ( T H ) | ω y ).The construction of π VH is simply based on a local H -projection on each patch ω y . Moreprecisely, we let Q y : H ( ω y ) → S ( T H ( ω y )) be defined through the solution to a local Neumannproblem, where Q y ( v ) ∈ S ( T H ( ω y )) solves ˆ ω y ∇ Q y ( v ) · ∇ w = ˆ ω y ∇ v · ∇ w for all w ∈ S ( T H ( ω y ))under the constraint ´ ω y Q y ( v ) dx = 0. Since Q y does not preserve constants, we need to restorethem through a local averaging operator M : L (Ω) → S ( T H ), which is defined by M v := (cid:88) y ∈(cid:52) ( T H ) (cid:32) ˆ ω y z y v dx (cid:33) λ y . Here, λ y ∈ S ( T H ) is the standard nodal basis function (hat function) associated with the vertex y , i.e., it fulfills λ y (˜ y ) = δ y ˜ y for all ˜ y ∈ ∆ ( T H ).With this, the node-based Falk-Winther projection is given by π VH ( v ) := M ( v ) + (cid:88) y ∈(cid:52) ( T H ) ( Q y ( v ))( y ) λ y . The operator is designed in the same way in 2 D . Though the construction of π VH is not requiredfor the construction of π EH , it can be helpful for validating the implementation through checkingthe commuting diagram property.3 A.2 Construction of S In the next step, we discuss the construction of S in (28). We start with defining the extendededge patch, for an edge E ∈ ∆ ( T H ), as the union of the two nodal patches that are associatedwith the corners of the edge. More precisely, for an edge E = conv { y , y } ∈ ∆ ( T H ) with vertices(corners) y , y ∈ ∆ ( T H ), the extended edge patch is given by ω ext E := ω y ∪ ω y . Furthermore, recalling the definition of z y from (29), we define for any edge E = conv { y , y } ∈ ∆ ( T H ) with vertices y , y ∈ ∆ ( T H ) the piecewise constant function ( δz ) E ∈ L (Ω) by( δz ) E := z y − z y . Before we can proceed, we need some additional notation. For that, we let T H ( ω ext E ) be therestriction of the mesh T H to the extended edge-patch ω ext E . We define the space of lowest-orderN´ed´elec finite elements over T H ( ω ext E ) and with a zero tangential trace on ∂ω ext E by˚ N ( T H ( ω ext E )) := { v H ∈ N ( T H ( ω ext E )) | v H × n = 0 on ∂ω ext E } . Similarly, we let ˚ RT ( T H ( ω ext E )) denote the space of lowest-order Raviart–Thomas finite elementsover T H ( ω ext E ) and with a zero normal trace on ∂ω ext E , i.e.,˚ RT ( T H ( ω ext E )) := { v H ∈ RT ( T H ( ω ext E )) | v H · n = 0 on ∂ω ext E } . Given E ∈ ∆ ( T H ), we search for a representation of ( δz ) E as the divergence of a ˚ RT -function.More precisely, we seek for each ( δz ) E a function z E ∈ ˚ RT ( T H ( ω ext E ))such that div z E = − ( δz ) E and ( z E , curl τ ) = 0 for all τ ∈ ˚ N ( T H ( ω ext E )) . This problem is well-posed, since ( δz ) E is T H -piecewise constant with zero average over ω ext E .Existence and uniqueness then follow from the exactness of the corresponding discrete complexwith vanishing traces (cf. [14]). A motivation for the crucial role of z E for the commutingproperty of the projection π EH is given in Section A.3 below. Remark A.1.
In practice, z E can be found by solving the following saddle point problem. Find( z E , v ) ∈ ˚ RT ( T H ( ω ext E )) × ˚ N ( T H ( ω ext E )) such that(div z E , div v ) + ( w , curl v ) = ( − ( δz ) E , div w ) for all w ∈ ˚ RT ( T H ( ω ext E ))( z E , curl τ ) = 0 for all τ ∈ ˚ N ( T H ( ω ext E )) . Using the functions z E , we define the operator M : L (Ω; C ) → N ( T H ) that maps anyfunction u ∈ L (Ω; C ) to M u := (cid:88) E ∈ ∆ ( T H ) ˆ ω ext E u · z E dx ψ E . As a last constituent of S , we define the operator Q y, − : H (curl , ω y ) → S ( T H ( ω y ))4again through the solution of a local discrete Neumann problem. More precisely, for u ∈ H (curl , ω y ), the image Q y, − u ∈ S ( T H ( ω y )) is the solution to( u − ∇ Q y, − u , ∇ v ) = 0 for all v ∈ S ( T H ( ω y )) ˆ ω y Q y, − u dx = 0 . Note the close relation of Q y, − to the operator Q y from Section A.2.With all the parts available, we can define the operator S : H (curl , Ω) → ˚ N ( T H ) via S u := M u + (cid:88) y ∈ ∆ ( T H ) ( Q y, − u )( y ) ∇ λ y . (30)From a practical point of view, it is useful to rewrite the operator S in terms of the edge-basedbasis functions ψ E of N ( T H ). Given a vertex y ∈ ∆ ( T H ), we can expand λ y in terms of thebasis functions ( ψ E ) E ∈ ∆ ( T H ) . We obtain ∇ λ z = (cid:88) E ∈ ∆ ( T H ) ˆ E ∇ λ z · t E ds ψ E = (cid:88) E ∈ ∆ ( z ) sign( t E · ∇ λ z ) ψ E where ∆ ( z ) ⊆ ∆ ( T H ) is the set of all edges that contain z . Thus, we can rewrite the operator S from (30) as S u := M u + (cid:88) E ∈ ∆ ( T H ) (cid:0) ( Q y ( E ) , − u )( y ( E )) − ( Q y ( E ) , − u )( y ( E )) (cid:1) ψ E . (31)At this point we have the commuting property, but S is not a projection yet. Remark A.2.
For the practical realization of the operator S , we seek a matrix representation, S , of S : N ( T h ) → N ( T H ), where T h is a refinement of T H that resolves the multiscale variationsof the problem. For a given edge E ∈ ∆ ( T H ), define the vector M E M E := (cid:34) ˆ ω ext E z E · ψ e dx (cid:35) e ∈T h ( ω ext E ) . To compute Q y, − , the mean constraint is enforced by a classical Lagrange multiplier, which leadsto a saddle point system. For a node y ∈ ∆ ( T H ) define the matrices A y := (cid:34) ˆ ω y ∇ λ z · ∇ λ w dx (cid:35) z,w ∈ ∆ ( T H ( ω y )) B y := (cid:34) ˆ ω y ∇ λ z · ψ e dx (cid:35) z ∈ ∆ ( T H ( ω y )) e ∈ ∆ ( T h ( ω y )) and the vector C y = (cid:34) ˆ ω y λ z dx (cid:35) z ∈ ∆ ( T H ( ω y )) . Solve the following system for Q y (cid:20) A y C ∗ y C y (cid:21) (cid:20) Q y V y (cid:21) = (cid:20) B y (cid:21) and let Q j = Q y [ y j ( E ) , :], for j = 1 ,
2. Finally, the matrix representation of S can be expressedas S ( E, [ e , ..., e N ]) = M E + Q − Q , where e , ...., e N denotes the edges in the patch T h ( ω ext E ).5 Remark A.3.
In the 2 D case, the computation of z E needs to be slightly changed. The lowest-order N´ed´elec space with zero tangential trace on ∂ω ext E , ˚ N ( T H ( ω ext E )), should be replaced withthe first-order (scalar-valued) Lagrange finite element space with zero trace on ∂ω ext E , that is,˚ S ( T H ( ω ext E )) := { v ∈ S ( T H ) | v = 0 on ∂ω ext E } . Except for curl having a different meaning in 2 D , this is the only change in the interpolation π EH that needs to be made for the 2 D case. A.3 Commuting property ∇ ◦ π VH = S ◦ ∇ For a better understanding of the construction of S , we will demonstrate that it indeed impliesthe commuting property ∇ π VH ( v ) = S ( ∇ v ) for all v ∈ H (Ω) . We start with applying S to a gradient ∇ v for some v ∈ H (Ω). First, we easily observe Q y, − ( ∇ v ) = Q y ( v ) . Next, for an edge E = { y , y } we consider the tangential component of ∇ M ( v ) (i.e., the localdegree of freedom) which is given by ˆ E ∇ M ( v ) · t E ds = | E | − ˆ E ∇ M ( v ) · ( y − y )= | E | − ˆ E (cid:32)(cid:32) ˆ ω y z y v dx (cid:33) ∇ λ y + (cid:32) ˆ ω y z y v dx (cid:33) ∇ λ y (cid:33) · ( y − y )= | E | − ˆ E (cid:32) ˆ ω y z y v dx − ˆ ω y z y v dx (cid:33) ∇ λ y · ( y − y )= (cid:32) ˆ ω y z y v dx − ˆ ω y z y v dx (cid:33) = ˆ ω ext E ( δ z ) E v dx = − ˆ ω ext E div z E v dx = ˆ ω ext E z E ∇ v dx. Since ∇ M ( v ) ∈ N ( T H ) is uniquely determined by these degrees of freedom, we have ∇ M ( v ) = (cid:88) E ∈ ∆ ( T H ) (cid:32) ˆ ω ext E z E ∇ v dx (cid:33) ψ E = M ( ∇ v ) . Combining the equations Q y, − ( ∇ v ) = Q z ( v ) and ∇ M ( v ) = M ( ∇ v ) for all v ∈ H (Ω), weobtain ∇ π VH ( v ) = ∇ M ( v ) + (cid:88) z ∈ ˚ (cid:52) ( T H ) ( Q z v )( z ) ∇ λ z = M ( ∇ v ) + (cid:88) y ∈ ˚ (cid:52) ( T H ) Q y, − ( ∇ v )( y ) ∇ λ y (30) = S ( ∇ v )as desired. Observe that it can be helpful for the validation of the Falk-Winther implementationto check the commuting property ∇ π VH v = S ( ∇ v ) numerically.6 A.4 Construction of Q E As mentioned before, the operator S commutes with vertex-based Falk-Winther operator π VH ,but it is not a projection yet. In order to modify S accordingly, another operator Q E : H (curl , ω ext E ) → N ( T H ( ω ext E ))needs to be introduced for all edges E ∈ ∆ ( T H ). Given an edge E and some some u ∈ H (curl , ω ext E ) the image Q E ( u ) ∈ N ( T H ( ω ext E )) is defined by the system( u − Q E u , ∇ τ ) = 0 for all τ ∈ S ( T H ( ω ext E ))(curl( u − Q E u ) , curl v ) = 0 for all v ∈ N ( T H ( ω ext E )) . Again, existence and uniqueness follow from the exactness of the discrete complex (cf. [3, 4, 14]).With Q E available, it is now possible to compute the edge-based Falk-Winther projectionaccording to (28) by π EH u = S u + (cid:88) E ∈ ∆ ( T H ) ˆ E (cid:0) ( I − S ) Q E u (cid:1) · t E ds ψ E , where I denotes the identity operator. Remark A.4.