Analysis and application of an overlapped FEM-BEM for wave propagation in unbounded and heterogeneous media
AAnalysis and application of an overlapped FEM-BEM forwave propagation in unbounded and heterogeneous media
V. Domínguez ∗ M. Ganesh † February 9, 2021
Abstract
An overlapped continuous model framework, for the Helmholtz wave propagation prob-lem in unbounded regions comprising bounded heterogeneous media, was recently intro-duced and analyzed by the authors (
J. Comput. Phys., , 109052, 2020 ). Thecontinuous Helmholtz system incorporates a radiation condition (RC) and our equivalenthybrid framework facilitates application of widely used finite element methods (FEM) andboundary element methods (BEM), and the resulting discrete systems retain the RC ex-actly. The FEM and BEM discretizations, respectively, applied to the designed interiorheterogeneous and exterior homogeneous media Helmholtz systems include the FEM andBEM solutions matching in artificial interface domains, and allow for computations of theexact ansatz based far-fields. In this article we present rigorous numerical analysis of a dis-crete two-dimensional FEM-BEM overlapped coupling implementation of the algorithm.We also demonstrate the efficiency of our discrete FEM-BEM framework and analysis us-ing numerical experiments, including applications to non-convex heterogeneous multipleparticle Janus configurations. Simulations of the far-field induced differential scatter-ing cross sections (DSCS) of heterogeneous configurations and orientation-averaged (OA)counterparts are important for several applications, including inverse wave problems. Ourrobust FEM-BEM framework facilities computations of such quantities of interest, withoutboundedness or homogeneity or shape restrictions on the wave propagation model.
AMS subject classifications:
Keywords:
Helmholtz, Heterogeneous, Unbounded, Wave Propagation, Finite Element Meth-ods, Integral Equations, Nyström Boundary Element Methods, Janus Configurations ∗ Department of Estadística, Informática y Matemáticas, Universidad Pública de Navarra, Tudela, Spain/Institute for Advanced Materials (INAMAT), Pamplona, Spain. Email: [email protected] † Department Applied Mathematics and Statistics Department, Colorado School of Mines, Golden, CO, USA.Email: [email protected] a r X i v : . [ m a t h . NA ] F e b Introduction
Simulation of scattered acoustic and electromagnetic fields, and hence understanding the im-pact of refractive indices of wave propagation media, are crucial for a large class of applica-tions [6, 24, 32, 34]. The term
Janus particles was mentioned in a Nobel Prize lecture [9] aboutthree decades ago, and since then additional interests include understanding the effect of a classof piecewise-continuous heterogeneous refractive indices induced Janus configurations. A sim-ple Janus particle is designed by combining two distinct homogeneous refractive indices, andJanus configurations in general may comprise multiple particles or complex structures withheterogeneous material properties. The example configuration Ω , of the type illustrated inFigure 1, has been investigated for synthesis and applications, see for example [29, 31, 38, 39].Figure 1: An incident wave u inc impinging on a given Janus-type heterogeneous multiple-particle configuration ( Ω ) in R . Our framework introduces two artificial boundaries Σ and Γ ,and hence two overlapped computational regions: (i) a bounded FEM-domain Ω (with boundary Σ ) and (ii) an unbounded BEM-region R \ Ω (exterior to the smooth interface Γ ). The FEMand BEM solutions are constrained to match in the overlapped region Ω :=( R \ Ω ) ∩ Ω . Developing numerical methods to understand scattering effects by Janus configurations(even with simple structures) has been an active area of research, see for example the Janusspherical acoustic configuration effect research [21, 34] published in 2020 and references therein.Motivated by details in the book [34], authors of [21] numerically investigated scattering effectsof homogeneous spheres with two standard (sound-soft and transmission) boundary proper-ties, leading to solving the Helmholtz equation with two distinct boundary conditions and theSommerfeld radiation condition (SRC), for the unknown scattered- and far-fields.As described in [21, 34], to understand the impact of Janus-type wave configurations, accu-rate simulations of certain quantities of interests (QoI) are important. The two key QoI are thefar-field intensity based differential scattering cross sections (DSCS) and also the configurationorientation effect QoI, the orientation-averaged (OA) DSCS. Simulation of the OA-DSCS isequivalent to simulating far-fields generated by the configuration, such as Ω in Figure 1, withthe incident wave u inc impinging Ω from hundreds of directions surrounding the configuration.The key numerical tool in [21] is the stable far-fields based T-matrix framework that was ana-lyzed (with a priori bounds ) and implemented in [13, 14]. For nonlinear inverse problems based2pplications that use the far-field (with phase or phase-less DSCS) data, efficient simulationof far-fields (such as at each Newton-type iterations) are important, see for example the 2019book [6] and extensive references therein.Far-field computations based numerical schemes for such QoI and tools, in general, donot allow for configurations with heterogeneous structures and unbounded regions. This ismainly because numerical partial differential equation (PDE) algorithms mainly apply eitherthe ubiquitous finite element method (FEM) or the boundary element method (BEM). TheFEM requires a bounded domain for finite number of tessellations, and the BEM requirementof the fundamental solution is practical mainly for a homogeneous medium PDE. The FEM isapplied to a variational formulation equation of the PDE that involves domain integrals, andthe BEM discretizes an equivalent boundary integral equation (BIE) for the boundary unknownin a chosen ansatz for the scattered field, in conjunction with the fundamental solution [6, 32].Accordingly for wave propagation models with a bounded heterogeneous medium ( Ω ), toapply the FEM, an artificial truncation of the unbounded region using, typically, a polygonalboundary Σ (as in Figure 1) is introduced and a wave absorbing boundary condition (ABC)on Σ is imposed, ignoring the SRC. The standard variational formulation on the boundarypolygonal domain (such as Ω in Figure 1) for the heterogeneous Helmholtz PDE with ABC on Σ is sign-indefinite [16,24] and this non-coercive restriction was removed recently by developing,analyzing, and implementing a new sign-definite variational formulation [17]. Another option toexactly incorporate the SRC is, for example, by ignoring the heterogeneity and consider insteadthe homogeneous model exterior to the polygonal boundary Σ or exterior to a smooth boundary Γ (as in Figure 1). The artificial smooth boundary choice is especially suitable for developinghigh-order spectrally accurate approximation of the scattered field in the homogeneous mediumusing high-order BEMs.For Janus-type configurations, it is important to avoid restrictions in either of the abovetwo options by using both the FEM and BEM. This can be achieved by appropriately couplingthe FEM and BEM solutions, depending on the choices of the artificial boundaries to computethe interior (heterogeneous) and exterior (homogeneous) media solutions. Such FEM-BEMcoupling mathematical frameworks have been developed and analyzed only by a few authors,but several researchers implemented the associated FEM-BEM computational frameworks.The widely used FEM-BEM coupling is obtained by choosing one artificial (FEM appropri-ate) polygonal boundary Σ , see for example the review article [36] and references therein. Inaddition to analysis difficulties [36], this approach introduces restricted regularity for the solu-tion exterior to Σ . The high-order regularity based (BEM appropriate) single smooth boundarychoice Γ was subsequently developed in [26, 27]. For recent implementations of the single arti-ficial boundary based FEM-BEM framework with high-order accuracy we refer to [15, 16, 19],and associated analysis issues are highlighted in [15, Section 6]. Using two artificial interfacesa mathematical framework was developed over four decades ago in [25] and was subsequentlyused in [8, 22]. A domain-overlapped framework [illustrated in Figure 1 with overlapped re-gion ( R \ Ω ) ∩ Ω ] was developed recently by the authors in [10] that takes advantage of apolygonal boundary ( Σ ) required for a wide-class FEM, and a smooth boundary ( Γ ) for high-order spectral BEM. Mathematical analysis of the continuous framework in [10] establishes theequivalence and regularity of the decomposed model.In this work we present rigorous numerical analysis of the FEM-BEM adaptive couplingframework introduced in [10] for solving the Helmholtz acoustic/electromagnetic wave propa-gation problem on the plane with a bounded heterogeneous region, and demonstrate efficiency3f the numerical algorithm for complex Janus-type heterogeneous configurations with compu-tational experiments. The algorithm works on a suitable partition of the plane defined from apolygonal domain containing a smooth curve that surround all heterogeneity. The overlappedpartition is made up of the bounded polygonal domain containing the heterogeneity, and theunbounded homogeneous region exterior to the smooth curve.On the bounded domain of the partition we approximate the solution by a FEM withclassical continuous piecewise polynomials on triangular grids, whereas a high-order NyströmBEM is used to compute the scattered wave in the unbounded region. Both solutions arecoupled by demanding to coincide in the two artificial boundaries that ensures the FEM andBEM solutions matching in the intersecting common domain of the partition. We prove that theconvergence of the scheme in natural norms is of the same order as the best approximations ofthe projection in the FEM and BEM finite dimensional spaces. In addition, since the Nyströmmethod is super-algebraically convergent, only a relatively few degrees of freedom are requiredin the BEM solver to keep the error of the same order as the FEM solution, and hence thealgorithm facilitates accurate far field, DSCS, and OA-DSCS computations.The rest of this article is organized as follows: In Section 2, we recall the Helmholtz modeland an equivalent decomposition decomposition framework, both at continuous level. In Sec-tion 3 we setup a discrete counterpart of the decomposition framework and the overlappedFEM-BEM algorithm. Rigorous numerical analysis of the FEM-BEM algorithm establishingoptimal order convergence, of the hybridized numerical solutions, in Section 4 forms the maintheoretical contribution of the article that include references to several results proven in the twoAppendix sections of this article. In Section 5, we computationally demonstrate the describedalgorithm and theoretically analysis using three distinct sets of experiments, in conjunction withthe theory and practical applicability, including implementation of the algorithm for multiple-particle Janus-type configurations with non-smooth solutions, and also for complex structuredheterogeneous regions. Throughout the article, let n be a piecewise-continuous refractive index with heterogeneityrestricted to a bounded domain Ω ⊂ R , and in the exterior we take n | Ω c0 ≡ so that theexterior Ω c0 (:= R \ Ω ) is a free-space unbounded medium.The scattered field u s is induced by an incident wave u inc (with wavenumber k > ) fromthe exterior region impinging on the heterogeneous medium Ω . It is convenient to assume thatthat incident field satisfies the homogeneous Helmholtz equation ∆ v + k v = 0 in all of the plane R , although it is sufficient to take it to satisfy outside of a compact set in R containing, say,a point-source. Physically appropriate incident waves such as the plane-wave and point-sourcehave these properties.We seek the total wave field u (:= u s + u inc ) ∈ H ( R ) , representing acoustic or electro-magnetic fields, satisfying the uniquely solvable problem governed by the variable coefficientHelmholtz equation and the SRC (with (cid:98) x := x / | x | ∈ S ): (cid:12)(cid:12)(cid:12)(cid:12) ∆ u + k n u = 0 , in R ,∂ (cid:98) x u s − i ku s = o ( | x | − / ) , as | x | → ∞ . (2.1)We note that in the heterogeneous medium Ω , u s is the interior unknown field and is the solu-tion of inhomogeneous Helmholtz equation with inhomogeneous term f = − (∆ u inc + k n u inc ) .4or the unknown scattered field u s exterior to Ω , in (2.1), the SRC lim | x |→∞ | x | / (cid:18) ∂u s ( x ) ∂ (cid:98) x − i ku s ( x ) (cid:19) = 0 . (2.2)holds uniformly in all directions (cid:98) x = x / | x | ∈ S . The scattered field u s is a radiating field and,as a consequence of the SRC, its behavior at infinity is captured by the far-field u ∞ ∈ L ( S ) ,where u ∞ ( (cid:98) x ) = lim | x |→∞ | x | / e − i k | x | u s ( x ) . (2.3)For the total field u induced by a plane wave with incident direction (cid:98) d ∈ S , it is appropriateto denote the associated far-field as u ∞ ( (cid:98) x ; (cid:98) d ) = u ∞ ( θ ; φ ) , with (cid:98) x = p ( θ ) = (cos θ, sin θ ) and (cid:98) d = p ( φ ) = (cos φ, sin φ ) , for some θ, φ ∈ [0 , π ) . The single-incident plane wave QoI DSCSand multiple-incident plane waves QoI OA-DSCS are then given by [21, 34] u DSCS ( θ ; φ ) = | u ∞ ( θ ; φ ) | ; u OADSCS ( θ ) = 12 π (cid:90) π | u ∞ ( θ ; φ ) | d φ. (2.4)Clearly, computation of the u OADSCS with high-accuracy requires discretization of the above inte-gral with over thousand discrete incident direction angles φ , leading to solving the same numberof Helmholtz model (2.1)-(2.2) for many distinct inputs. This motivates the difficulty of eval-uating u OADSCS
Janus-type configurations with piecewise-continuous refractive indices defined onnon-trivial geometries.For the given wave propagation problem (2.1), next we recall an equivalent decompositionframework introduced and analyzed in [10]. The framework introduce two artificial curves Γ and Σ with interior Ω and Ω , respectively, satisfying Ω ⊂ Ω ⊂ Ω ⊂ Ω , with the assumptionthat Γ is smooth and Σ is a polygonal boundary. A sketch of the different domains is displayedin Figure 1. We denote henceforth Ω c1 := R \ Ω . At continuous level, it is convenient toconsider the decomposition framework using operators defined on classical Sobolev spaces. Tothis end, for a general domain D ∈ R with boundary ∂D and for real m , let H m ( D ) denotethe classical Sobolev space.We also consider H s ( ∂D ) which is well defined for any s if D is smooth. Recall that in thiscase the trace operator γ ∂D : H s +1 / ( D ) → H s ( ∂D ) is continuous for any s > as consequenceof the Sobolev Trace Theorems, see for example [7] or [30, Ch 4].For Lipschitz domains Ω with boundary Σ = ∂ Ω , such as the chosen polygonal domain Ω , H s (Σ) is defined only for s ∈ [ − , [1, 20, 30]. We then commit a (slight) abuse of notationand set for s > H s (Σ) = (cid:8) γ Σ u : u ∈ H s +1 / ( D ) (cid:9) , (2.5)with D a open domain such that Σ ⊂ D , endowed with the image norm. It is known that thespace is independent of D , and that it is still the classical Sobolev for s ∈ (0 , , but does nothold for s = 1 . In other words, using the definitions, the trace operator γ Σ : H s +1 / ( D ) → H s (Σ) is still continuous for s ∈ (0 , ∪ (1 , ∞ ) .The decomposition framework starts with an interior and an exterior Helmholtz problem:• ( Interior Dirichlet Helmholtz problem in Ω with polygonal boundary Σ ):Given f Σ ∈ H / (Σ) , find ω int ∈ H (Ω ) so that (cid:12)(cid:12)(cid:12)(cid:12) ∆ ω int + k n ω int = 0 , in Ω ,γ Σ ω int = f Σ . (2.6)5 ( Exterior Dirichlet Helmholtz problem in Ω c1 with smooth boundary Γ ):Given f Γ ∈ H / (Γ) , find ω ext ∈ H (Ω c1 ) so that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∆ ω ext + k ω ext = 0 , in Ω c1 ,γ Γ ω ext = f Γ ,∂ (cid:98) x ω ext − i kω ext = o ( | x | − / ) . (2.7)The SRC in (2.7) ensures that the exterior Dirichlet problem is uniquely solvable [6, 32] for allwavenumbers k . For the interior Dirichlet Helmholtz problem (2.6) the well-posedness does nothold for all wavenumbers k [24], and throughout this article we assume that the wavenumber k is such that ω int is the unique solution of (2.6). The well-posedness assumption for the interiorDirichlet Helmholtz problems in Ω (and in the overlapped region Ω := Ω c ∩ Ω ) can be easilyavoided, for example, by modifying the artificial boundaries Σ and Γ [10, Section 2.1].It is convenient to use operators to describe the continuous decomposition framework.To this end, corresponding to the interior Dirichlet problem (2.6) we consider two operators K Ω Σ , K ΓΣ ; and associated with the exterior problem (2.7) we consider two operators K Ω c1 Γ , K ΣΓ .These pairs of operators are defined, using the unique solution ω int of (2.6) and ω ext of (2.7)and their traces respectively on Γ and Σ , as follows: K Ω Σ f Σ = ω int , K ΓΣ f Σ = γ Γ ω int ; K Ω c1 Γ f Γ = ω ext , K ΣΓ f Γ = γ Σ ω ext . (2.8)In [10], the authors proved that the unique solution u of (2.1) can be constructed in two stepsas follows:1. Find f Σ : Σ → C and f Γ : Γ → C so that (cid:12)(cid:12)(cid:12)(cid:12) f Σ − K ΣΓ f Γ = γ Σ u inc , − K ΓΣ f Σ + f Γ = − γ Γ u inc . (2.9a)2. Construct u = (cid:40) K Ω Σ f Σ , in Ω , K Ω c1 Γ f Γ + u inc , in Ω c1 . (2.9b)The boundary unknowns system (2.9a) can be written in matrix-valued operator form as ( I − K ) (cid:20) f Σ f Γ (cid:21) = (cid:20) γ Σ u inc − γ Γ u inc (cid:21) , K := (cid:20) K ΣΓ K ΓΣ (cid:21) . (2.10)( I is obviously the identity matrix operator.) It is not difficult to see that the off-diagonalblock K : H (Σ) × H (Γ) → H s (Σ) × H s (Γ) is continuous for any s ≥ [10].The following result summarizes the equivalence of the decomposition framework to theoriginal problem (2.1) at the continuous level. Theorem 2.1.
Assume that the only solution to problem (cid:12)(cid:12)(cid:12)(cid:12) ∆ v + k v = 0 , in Ω = Ω ∩ Ω c ,γ Γ v = 0 , γ Σ v = 0 (2.11) is the trivial one. Then I − K : H s (Σ) × H s (Γ) → H s (Σ) × H s (Γ) is invertible for any s ≥ .Therefore (2.9a) is uniquely solvable. Furthermore, u defined by (2.9b) is the unique solutionof the full Helmholtz problem (2.1) . roof. This result is proven in [10]. We will give here a sketch of the proof for the sake ofcompleteness. By Fredholm alternative it suffices to show that
I − K is one-to-one. We notethat any if ( f ∗ Σ , f ∗ Γ ) ∈ N ( I − K ) then ω := (K Ω Σ f ∗ Σ − K Ω c1 Γ f ∗ Γ ) | Ω is a solution for (2.11).By hypothesis, ω = 0 . Hence, u in (2.9b) is well defined for u inc = 0 . The well-posednessof the scattering problem (2.1), implies that u = 0 . In particular, K Ω c1 Γ f ∗ Γ | Ω c = 0 , and fromthe analytic continuation principle K Ω c1 Γ f (cid:63) Γ = 0 . Therefore, also K Ω Σ f ∗ Σ | Ω = 0 and hence ( f ∗ Σ , f ∗ Γ ) = . Remark . For the OA-DSCS calculations with a large number of incident plane waves u inc , amarked advantage of our equivalent formulation (2.9a)-(2.9b) compared to the original modelproblem (2.1) [with a large number of inhomogeneous terms f = − (∆ u inc + k n u inc ) in Ω ] isthat the interior and exterior homogeneous problems K Ω Σ f Σ , K Ω c1 Γ f Γ (and hence their traces K ΓΣ f Σ , K ΣΓ f Γ ) can be setup independent of the incident waves. Since the unknowns f Σ and f Γ are defined, respectively, only on boundary curves Σ and Γ , these unknowns can represented bya few boundary unknowns (or boundary basis functions) and hence the associated Helmholtzproblems setup K Ω Σ f Σ , K Ω c1 Γ f Γ is a naturally parallel process for computational purposes. The numerical discretization of the continuous decomposition framework (2.9) is essentiallyobtained by replacing the four continuous operators in (2.8) using appropriate FEM and BEMbased discrete counterparts. For K Ω Σ we will use the standard FEM with continuous piecewisepolynomial elements on a triangular conformal mesh of Ω [24].For the BEM certainly an extensive range of methods is at our disposal in the literature [6,32]. We will restrict ourselves to the spectrally accurate Nyström method [6, 28] .This schemeprovides a discretization of the key boundary integral operators of the associated Calderoncalculus; and in this article we will use only make use of the discrete Single- and Double-Layeroperators that converge super-algebraically, and it is not difficult to implement. A disadvantageof the Nyström method is that it requires an accurate differentiable parameterization of theboundary. This is because the method is based on splitting the integral operators into regularand singular parts for which appropriate decompositions and factorizations of the kernels ofthe operators are needed. This is not a severe restriction in our case since Γ is an auxiliaryuser-chosen artificial curve and therefore can be taken to be simple and smooth.In the next two subsections we recall the standard FEM and Nyström procedure and con-clude this section with the discrete FEM-BEM decomposition framework required for the mainfocus of this work on the numerical analysis of the FEM-BEM algorithm [10] . The FEM procedure
Let {T h } h be a sequence of regular triangular meshes with h denoting the discrete parameter,the diameter of the largest element of the grid. We then write h → to mean that the maximumof the diameters of the elements tends to 0. A technical mesh assumption, not restrictive inpractice, described in detail in Assumption 1 below, will be used in our proofs to ensure (i) afaster convergence of the FEM in stronger norms around Γ (see Theorem B.1 in Appendix);and (ii) the stability of the full method. 7n T h , we construct the finite dimensional FEM space P h,d := { v h ∈ C (Ω ) : u h | T h ∈ P d , ∀ T h ∈ T h } , where P d is the space of bivariate polynomial of degree d . Then given f h Σ ∈ γ Σ P h,d , we define K h Ω Σ : γ Σ P h,d → P h,d as K h Ω Σ f h Σ := u h that is the unique of the standard FEM discrete system: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u h ∈ P h,d b k,n ( u h , v h ) = 0 , ∀ v h ∈ P h,d ∩ H (Ω ) γ Σ u h = f h Σ , b k,n ( u, v ) := (cid:90) Ω ∇ u · ∇ v − k (cid:90) Ω n uv. (3.1)It is well known that K h Ω Σ [24] is well defined for any sufficiently fine mesh. Notice also that K h Ω Σ is defined on the discrete space γ Σ P h,d , the trace finite element space on the boundary Σ .Hence, with the help of Q h Σ : C (Σ) → γ Σ P h,d , (3.2)the nodal interpolation operator on γ Σ P h,d defined by the finite element space, we have K h Ω Σ Q h Σ ≈ K Ω Σ , providing an optimal approximation in P h,d for sufficiently smooth Dirichlet data f h Σ (see The-orem 4.1 in the next section). The BEM procedure
Let x = ( x ( t ) , x ( t )) : R → Γ (3.3)be a smooth π − periodic parameterization of Γ . We denote by SL k , DL k , the (parameterized)layer potentials, defined for z ∈ R \ Γ as, (SL k ϕ )( z ) := (cid:90) π Φ k ( z − x ( t )) ϕ ( t ) d t, (DL k g )( z ) := (cid:90) π (cid:0) ∇ y Φ k ( z − y ) (cid:1)(cid:12)(cid:12)(cid:12) y = x ( t ) · ν ( t ) g ( t ) d t, where Φ k = i4 H (1)0 ( k | · | ) is the fundamental solution for the constant coefficient Helmholtzoperator in R with wavenumber k ; ν ( t ) := ( x (cid:48) ( t ) , − x (cid:48) ( t )) is a non-normalized normal vector.Observe that | x (cid:48) ( t ) | > is then incorporated to the density in SL k and to the kernel in DL k . Wefollow the same convention for the single and double boundary operator defined for π -periodicdensities ϕ and g as (V k ϕ )( s ) := ( γ Γ SL k ϕ ) ( x ( s )) = (cid:90) π Φ k ( x ( s ) − x ( t )) ϕ ( t ) d t (3.4) (K k g )( s ) := ± g ( s ) + ( γ ∓ Γ DL k g )( s ) = (cid:90) π (cid:0) ∇ y Φ k ( x ( s ) − y ) (cid:1)(cid:12)(cid:12)(cid:12) y = x ( t ) · ν ( t ) g ( t ) d t. (3.5)8he Brakhage-Werner formulation, first introduced in [3] (see also [6, 32]) provides a robustrepresentation for the solution of the exterior Dirichlet solution for the Helmholtz problem(2.7): K Ω c1 Γ f Γ = (DL k − i k SL k )( I + K k − i k V k ) − ( f Γ ◦ x ) , (3.6)with I obviously being the identity operator. The above representation of the exterior scatteredfield, satisfying the SRC, provides an exact ansatz for the associated far-field as ( F ϕ )( (cid:98) z ) , (cid:98) z ∈ S defined, using the boundary density ϕ = ( I + K k − i k V k ) − ( f Γ ◦ x ) [6]: ( F ϕ ) ( (cid:98) z ) := (cid:114) k π exp (cid:0) − π i (cid:1) (cid:90) π exp( − i k ( (cid:98) z · x ( t ))) (cid:2) (cid:98) z · ( x (cid:48) ( t ) , − x (cid:48) ( t )) + 1 (cid:3) ϕ ( t ) d t. (3.7)Thus accurate computational approximations of ϕ provide spectrally accurate approximationsfor both the scattered and the far-field. For computing ϕ , the Nyström BEM solver makesuse of a decomposition of the single- and double- layer operator into logarithmic and regularparts: (V k ϕ )( s ) = (cid:90) π A ( s, t ) log sin s − t ϕ ( t ) d t + (cid:90) π B ( s, t ) ϕ ( t ) d t, (K k g )( s ) = (cid:90) π C ( s, t ) log sin s − t g ( t ) d t + (cid:90) π D ( s, t ) g ( t ) d t. Functions
A, B, C, D are smooth and π − biperiodic.Next with N being a positive integer BEM discretization parameter, we consider the grid { t j } j ∈ Z ⊂ R , t j := jπN , and the N -dimensional space of trigonometric polynomials defined by T N := span (cid:104) e (cid:96) : (cid:96) ∈ Z N (cid:105) , e (cid:96) ( t ) := exp(i (cid:96)t ) (3.8)where Z N = {− N + 1 , − N + 2 , . . . , N } . The interpolation operator for π − periodic functions ϕ T N (cid:51) Q N ϕ s.t. (Q N ϕ )( t j ) = ϕ ( t j ) , is known to be well defined.The Nyström method is based on the following approximations for the single- and double-layer operators: (V Nk ϕ )( s ) := (cid:90) π Q N ( A ( s, · ) ϕ (cid:1) ( t ) log sin s − t d t + (cid:90) π Q N ( B ( s, · ) ϕ (cid:1) ( t ) d t (K Nk g )( s ) := (cid:90) π Q N ( C ( s, · ) g (cid:1) ( t ) log sin s − t d t + (cid:90) π Q N ( D ( s, · ) g (cid:1) ( t ) d t. We stress that the integrals above can be computed exactly. Indeed, − π (cid:90) π log sin t e n ( t ) d t = − π (cid:90) π log sin t cos( nt ) d t = (cid:40) log 4 , n = 0 | n | n (cid:54) = 0 (cid:90) π (Q N g )( t ) d t = πN N − (cid:88) j =0 (Q N g )( t j ) = πN N − (cid:88) j =0 g ( t j ) , i.e., for the regular part the suggested approximation yields the trapezoidal/rectangular rulefor π − periodic functions.The evaluation of the potentials, as integral operators with smooth kernel, is carried out ina similar way: (cid:0) SL Nk ϕ (cid:1) ( z ) := (cid:90) π Q N (Φ k ( z − x ( · )) ϕ )( t ) d t (cid:0) DL Nk g (cid:1) ( z ) := (cid:90) π Q N (cid:0)(cid:0) ∇ y Φ k ( z − y ) (cid:1)(cid:12)(cid:12) y = x ( · ) · ν ( · ) g )( t ) d t. (3.9)We are ready to present the discrete version for K h Ω c1 Γ . Hence, in view of K Ω c1 Γ f Γ = (DL k − i k SL k ) L k ( f Γ ◦ x ) , with L k := ( I + K k − i k V k ) − , (3.10)we define for f = f Γ ◦ x K N Ω c1 Γ f := (DL Nk − i k SL Nk ) L Nk f , with L Nk := ( I + K Nk − i k V Nk ) − (3.11)Hence ϕ N := L Nk f is the spectrally accurate approximation to the density ϕ = L k f Γ . Inaddition , we can also introduce the associated far-field approximation, see (3.7), as (cid:0) F N ϕ N (cid:1) ( (cid:98) z ) := (cid:114) k π exp (cid:0) − π i (cid:1) πN N (cid:88) j = − N +1 exp( − i k ( (cid:98) z · x ( t j ))) (cid:2) (cid:98) z · ( x (cid:48) ( t j ) , − x (cid:48) ( t j )) + 1 (cid:3) ϕ N ( t j ) . (3.12) Remark . It is not difficult to show that ϕ N = L Nk f ⇒ Q N ϕ N = Q N L Nk Q N f . Since potentials and QoI uses in the computations only pointwise values of the density ϕ , wecan conclude that the true unknowns of the Nyström method are the values of the density atthe grid points { ϕ ( t j ) } j ∈ Z and that such values are computed by the algorithm using only thevalues of the right-hand-side f at the same grid. The BEM-FEM coupling method
Let K N ΣΓ := γ Σ K N Ω c1 Γ , K h ΓΣ := ( γ Γ K h Ω Σ ) ◦ x , be the discrete counterparts of K ΣΓ and K ΓΣ . That is, the first operator computes the BEMsolution and evaluate on the outer polygonal domain Σ , whereas the second one solves the FEMsolution and evaluate it on the interior and smooth Γ .Then our coupled FEM-BEM algorithm is:10 Solve ( f h Σ , f N ) ∈ γ Σ P h,d × T N s.t. (cid:18) I − (cid:20) Q h Σ K N ΣΓ Q N K h ΓΣ (cid:21)(cid:19) (cid:20) f h Σ f N (cid:21) = (cid:20) Q h Σ γ Σ u inc − Q N γ Γ ( u inc ◦ x ) (cid:21) (3.13a)• Construct w h = K h Ω Σ f h Σ , ω N = K N Ω c1 Γ f N , u h,N = (cid:40) w h , in Ω ,ω N + u inc , in Ω c1 . (3.13b) Remark . In view of Remark 3.1, for implementation of the above algorithm, the pointwisevalues of the numerical solution ( f h Σ , f N ) at the boundary nodes (of the FEM and BEM gridson Σ and Γ ) are the true unknowns. Hence (3.13a) leads to a relatively small algebraic system.We also note that, in view of the first-stage coupled and constrained solutions in (3.13a),the overlapped algorithm is designed in such a way that for sufficiently fine grids, the final-stageFEM and BEM parts of the algorithm in (3.13b) lead to two numerically coinciding solutionsin Ω = Ω c1 ∩ Ω . (We demonstrate this property using numerical experiments.) We analyze the above FEM-BEM scheme through rigorous derivation of the stability andconvergence estimates in appropriate Sobolev norms. For the polygonal region we keep usingthe standard space H s (Ω ) and H s (Σ) (with the convention (2.5) for s > ). For the smoothboundary Γ , since we switch to the parameterized spaced (via x ), we will work with π − periodicSobolev spaces. To this end, we define H s := { ϕ ∈ D (cid:48) ( R ) : ϕ = ϕ ( · + 2 π ) , (cid:107) ϕ (cid:107) H s < ∞} , with (cid:107) ϕ (cid:107) H s = | ϕ (0) | + (cid:88) n (cid:54) =0 | n | s | (cid:98) ϕ ( n ) | (4.1)and (cid:98) ϕ ( n ) being the n th Fourier coefficient: (cid:98) ϕ ( n ) := (cid:90) π ϕ ( t ) e − n ( t ) d t, It is a well established result that H s (Γ) and H s are isomorphic via the composition with x (see for example [28, Th. 8.13]). For the implementation as well as for the theoretical analysis, we need some projections onto thediscrete spaces defined on the boundaries Σ and Γ . In the first case we have already introducedthe Lagrange interpolation operator in P h,d (see (3.2)). For such operators it can be shown, asconsequence of well-known results, that (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch t − s Σ (cid:107) f Σ (cid:107) H t (Σ) (4.2a)11here C depends only on Σ , s ∈ [0 , and s ≤ t < d + 1 and t > / . Here and in what follows h Σ is the maximum of the diameters of elements in Σ induced by the mesh T h . The convergenceof order h d +1 − s Σ can be attained if we assume an extra regularity for f Σ : for any s ∈ [0 , and t > d + 1 , there exists C so that (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch d +1 − s Σ (cid:107) f Σ (cid:107) H t (Σ) . (4.2b)We refer the reader to Appendix A for proofs of such results. In this section we also introducea more flexible projection on γ Σ P dh , which will be required for analysis purposes. Roughlyspeaking this is a consequence of carrying out the analysis in H (Ω) and H / (Σ) norms which,although are the natural ones for the FEM method, contains discontinuous functions and forwhich the action of Q h Σ cannot be therefore considered. Hence, we claim that there exists P h Σ : H / (Σ) → γ Σ P h,d a projection on γ Σ P h,d satisfying (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch t − s Σ (cid:107) f Σ (cid:107) H t (Σ) , ≤ s < , s ≤ t < d + 1 , t ≥ / . (4.3a) (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch d +1 − s Σ (cid:107) f Σ (cid:107) H t (Σ) , ≤ s < , t > d + 1 . (4.3b)Proofs of the estimates (4.2a) and (4.3) are given in Proposition A.3 in Appendix A.The convergence of the trigonometric interpolation operator in the Sobolev frame is a wellknown result, see for example, [35, Th. 8.3.1]: (cid:107) Q N ϕ − ϕ (cid:107) H s ≤ CN s − t (cid:107) ϕ (cid:107) H t , (4.4)for any t ≥ s ≥ and t > / , where C depends only on s, t . For similar reasons, and again onlyfor our analysis, we use the L -projection, which turns out to be the H s -orthogonal projectionfor any s , on T N : T N (cid:51) P N ϕ , s.t. (cid:92) P N ϕ ( n ) = (cid:99) ϕ ( n ) , n ∈ Z N . It is straightforward to show that for any t ≥ s , (cid:107) P N ϕ − ϕ (cid:107) H s ≤ N s − t (cid:107) ϕ (cid:107) H t . (4.5)With the help of these two projections we define K h,N := (cid:20) Q h Σ K N ΣΓ P N Q N K h ΓΣ P h Σ (cid:21) (4.6)so that the (3.13a) can be set up as ( I − K h,N ) (cid:20) f Σ f (cid:21) = (cid:20) g Σ g (cid:21) for appropriate ( g Σ , g ) (cid:62) ∈ H / (Σ) × H / . Observe that in the case that the right-hand-side ( g Σ , g ) belongs to the discrete space γ Σ P h,d × T N (as in (3.13a)), so is the solution ( f Σ , f ) .In this case, P h Σ and P N in (4.6) are already acting on elements on the discrete space and cansafely removed. Thus the role of these projection is to facilite the analysis by setting up theequation in the continuous framework. 12 .2 Convergence for the FEM scheme We recall some classical convergence results as well as superconvergence phenomenon for theFEM solution in the following theorem.
Theorem 4.1.
For fixed g Σ ∈ H / (Σ) , g h Σ ∈ γ Σ P h,d , let u := K Ω Σ g Σ , u h := K h Ω Σ g h Σ , be the solution of the interior Helmholtz problem (2.6) and the approximation given by the FEM (3.1) . Then there exists C > independent of g Σ , g h Σ and h so that (cid:107) u − u h (cid:107) H (Ω) ≤ C (cid:104) inf v h ∈ P h,d (cid:107) u − v h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) H / (Σ) (cid:105) . (4.7) Furthermore, let D (cid:48) , D be domains with D (cid:48) ⊂ D (cid:48) ⊂ D ⊂ D ⊂ Ω \ Ω and {T h } h be asequence of regular grids with h → which are locally quasi-uniform in D . Then there exists δ ∈ (1 / , such that for any ε ∈ [0 , / and for any fine enough grid T h , (cid:107) u − u h (cid:107) H ε ( D (cid:48) ) ≤ C (cid:2) ( h δ h − εD + h − εD ) (cid:107) u − u h (cid:107) H (Ω) + h − εD (cid:107) g Σ − g h Σ (cid:107) L (Σ) + h d − εD (cid:107) u (cid:107) H d +1 ( D ) (cid:3) (4.8) with C > depending on ε , D and D (cid:48) , and h D being the maximum of the diameters of theelements of the grid contained in D .Proof. For (4.7) we refer to [37] (see also [11]). Estimate (4.8) can be derived using the su-perconvergence of the FEM solution in the interior of the computational domain where thesolution is smooth (actually analytic) cf. [33]. We give a proof of these results in Corollary B.3in Appendix B for the sake of completeness.Let us point out that the constant δ ∈ (1 / , in (4.8) depends on the regularity of the dualproblem with homogeneous Dirichlet condition and a right-hand-side in L (Ω) . Therefore, forconvex polygonal domains in R , we can take δ = 1 [20]. Assumption 1
Assume that for some open domain D ⊂ Ω \ Ω with Γ ⊂ D there exists ε > such that the sequence of grids {T h } h is quasi-uniform in D and satisfies h / h − ε D → where h D is the maximum of the diameters of the elements of the grid T h having non-emptyintersection with D . (cid:3) We note that this assumption allows locally refined grids but introduces a very weak re-striction on the ratio between the larger element in Ω and the smaller element in D . However,since the exact solution is smooth on D , it is reasonabe to expect that small elements are notgoing to be used in this subdomain. Lemma 4.2.
Let D and the sequence of grids {T h } be as in Assumption 1. Then, for any ε ∈ [0 , / there exists C > such that (cid:107) K ΓΣ f Σ − K h ΓΣ P h Σ f Σ (cid:107) H / ε ≤ Ch / h − εD (cid:2) inf v h ∈ P h,d (cid:107) K Ω Σ f Σ − v h (cid:107) H ( D (cid:48) ) + (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) (cid:3) + Ch d − εD (cid:107) f Σ (cid:107) L (Σ) . roof. Notice that (cid:107) f Σ − P h Σ f Σ (cid:107) L (Σ) = (cid:107) f Σ − P h Σ f Σ − P h Σ ( f Σ − P h Σ f Σ ) (cid:107) L (Σ) ≤ Ch / (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) . Take D (cid:48) so that Γ ⊂ D (cid:48) ⊂ D (cid:48) ⊂ D . We claim that there exist C independent of f Σ , so that (cid:107) K D Σ f Σ (cid:107) H d +1 ( D ) ≤ C (cid:107) f Σ (cid:107) L (Σ) . (4.9)This can be seen as consequence of that the differential equations in D becomes the homoge-neous Helmholtz equation and Γ is sufficiently far away from D .The continuity of the trace operator γ Γ : H ε ( D (cid:48) ) → H / ε (Γ) ∼ H / ε and a directapplication of (4.7-4.8) yield (cid:107) K ΓΣ f Σ − K h ΓΣ P h Σ f Σ (cid:107) H / ε ≤ C ( h δ h − εD + h − εD + h / − ε ) (cid:2) inf v h ∈ P h,d (cid:107) K Ω Σ f Σ − v h (cid:107) H (Ω ) + (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) (cid:3) + Ch d − εD (cid:107) f Σ (cid:107) L (Σ) . Using that δ ≥ / the result is proven.We are ready to prove the convergence, in operator norm, of the first off diagonal block in(4.6) to the corresponding one in (2.10) Proposition 4.3.
For any < ε ≤ ε and for any t ≥ there exists C > independent of f Σ , h and N such that (cid:107) K ΓΣ f Σ − Q N K h ΓΣ P h Σ f Σ (cid:107) H / ≤ C (cid:0) N − ε h − εD (cid:1) h / (cid:2) inf v h ∈ P h,d (cid:107) K Ω Σ f Σ − v h (cid:107) H (Ω ) + (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) (cid:3) + C ( N − ε h d − εD + h dD + N − t ) (cid:107) f Σ (cid:107) L (Σ) . (4.10) In particular, we have the convergence in operator norm: (cid:107) K ΓΣ − Q N K h ΓΣ P h Σ (cid:107) H / (Σ) → H / → as ( N, h ) → ( ∞ , . (4.11) Proof.
The identity K ΓΣ − Q N K h ΓΣ P h Σ = (cid:0) I − Q N (cid:1) K ΓΣ + (I − Q N ) (cid:0) K h ΓΣ P h Σ − K ΓΣ (cid:1) + (cid:0) K ΓΣ − K h ΓΣ P h Σ (cid:1) , (4.12)and the estimates (cid:107) (I − Q N )K ΓΣ f Σ (cid:107) H / ≤ CN − t (cid:107) K ΓΣ f Σ (cid:107) H t +1 / ≤ C (cid:48) N − t (cid:107) f Σ (cid:107) H (cid:107) (I − Q N )(K h ΓΣ P h Σ − K ΓΣ (cid:1) f Σ (cid:107) H / ≤ C ε N − ε (cid:107) (K h ΓΣ P h Σ − K ΓΣ (cid:1) f Σ (cid:107) H / ε , which are consequences of (4.4) and (4.9), yield (cid:107) K ΓΣ f Σ − Q N K h ΓΣ P h Σ f Σ (cid:107) H / ≤ C ε N − ε (cid:107) (K h ΓΣ P h Σ − K ΓΣ (cid:1) f Σ (cid:107) H / ε + C (cid:48) N − t (cid:107) f Σ (cid:107) H + (cid:107) (K h ΓΣ P h Σ − K ΓΣ (cid:1) f Σ (cid:107) H / Applying Lemma 4.2 twice (with ε = ε for the first term and ε = 0 for the third one) yield(4.10). Consequently, (4.11) follows. 14 .3 Convergence for the BEM scheme The inverse inequality (cid:107) ϕ N (cid:107) H t ≤ N t − s (cid:107) ϕ N (cid:107) H s , t ≥ s, ∀ ϕ N ∈ T N , (4.13)that is straightforward to derive from (3.8) and (4.1), will be used repeatedly in this section.The first result in this subsection summarizes the convergence of the BEM solver in a formatthat will be used later. This is based on the convergence in norm of the approximation operator L Nk to the continuous counterpart. We recall that L k : H s → H s is continuous for any s ∈ R . Theorem 4.4.
Fix t ≥ s > / . Then, there exists C > such that for any N large enoughand for any f ∈ H t , (cid:107)L k f − L Nk f (cid:107) H s ≤ CN s − t − α (cid:107) f (cid:107) H t , with α = min { s, } . (4.14) Therefore, L Nk : H s → H s is uniformly continous in N for any s > / .Moreover, for s, t ≥ , with t ≥ max { s − , } there exits C > such that, for N largeenough, (cid:107) Q N L k P N f − Q N L Nk P N f (cid:107) H s ≤ CN s − t − α (cid:107) P N f (cid:107) H t ≤ CN s − t − α (cid:107) f (cid:107) H t . (4.15) Proof.
The estimates, recall (3.10) and (3.11), (cid:107) K Nk − K k (cid:107) H t → H s + (cid:107) V Nk − V k (cid:107) H t → H s ≤ CN s − t − α with s, t, α as in the statement of the Theorem (see Chapter 12 and 13 in [28] or, for a moredetailed proof, [12, Th. 3.1]) proves that L Nk : H s → H s for s > / is well defined, for N large enough, and it is uniformly bounded. Estimate (4.14) follows now from the identity L k − L Nk = L Nk (cid:2)(cid:0) L Nk ) − − L − k (cid:3) L k = L Nk (cid:2) (K Nk − K k ) − i k (V Nk − V k ) (cid:3) L k . To prove the second estimate (4.15), we proceed in two steps: For s ∈ [1 , ∞ ) , t ≥ s − wehave (cid:107) Q N L Nk P N f − Q N L k P N f (cid:107) H s ≤ C (cid:107) ( L Nk − L k )P N f (cid:107) H s ≤ C (cid:48) N s − t − (cid:107) P N f (cid:107) H t +1 ≤ C (cid:48) N s − t − (cid:107) P N f (cid:107) H t . (4.16)(Notice that in the last step we have used the inverse inequality (4.13).) On the other hand,for s ∈ [0 , and t ≥ , we make use of the bound (cid:107) Q N g (cid:107) H ≤ C (cid:107) g (cid:107) H to derive (cid:107) Q N L Nk P N f − Q N L k P N f (cid:107) H s ≤ C (cid:107) ( L Nk − L k )P N f (cid:107) H ≤ C (cid:48) N − t − (cid:107) P N f (cid:107) H t +1 ≤ C (cid:48) N − t (cid:107) P N f (cid:107) H t . (4.17)Estimates (4.16) and (4.17) yield the desired result (4.15). Lemma 4.5.
For any domain D with D ∩ Γ = ∅ and any r, s there exists C > such that forany N and ϕ N ∈ T N , (cid:107) DL Nk ϕ N − DL k ϕ N (cid:107) H s ( D ) + (cid:107) SL Nk ϕ N − SL k ϕ N (cid:107) H s ( D ) ≤ CN − r (cid:107) ϕ N (cid:107) H . (4.18)15 roof. Since the kernels of the integral operators are smooth, the estimate follows from thealiasing effect of the trapezoidal rule for periodic functions. Indeed, for | n | ≤ N we have forany g smooth enough (cid:12)(cid:12)(cid:12)(cid:12) πN N − (cid:88) j =0 g ( t j ) e − n ( t j ) − (cid:90) π g ( t ) e − n ( t ) d t (cid:12)(cid:12)(cid:12)(cid:12) ≤ π (cid:88) (cid:96) (cid:54) =0 | (cid:98) g ( n + 2 (cid:96)N ) | , e n ( t ) := exp(i nt ) (see for instance [28, 35]). Hence, for any ϕ N ∈ T N (cid:12)(cid:12)(cid:12)(cid:12) πN N − (cid:88) j =0 g ( t j ) ϕ N ( t j ) − (cid:90) π g ( t ) ϕ N ( t ) d t (cid:12)(cid:12)(cid:12)(cid:12) ≤ π (cid:88) n ∈ Z N | (cid:98) ϕ ( n ) | (cid:88) (cid:96) (cid:54) =0 | (cid:98) g ( − n + 2 (cid:96)N ) | . The last term can be easily bounded using the Cauchy-Schwarz inequality: (cid:88) n ∈ Z N | ϕ N ( n ) | (cid:88) (cid:96) (cid:54) =0 | (cid:98) g ( − n + 2 (cid:96)N ) | ≤ (2 N ) − r (cid:107) ϕ N (cid:107) H (cid:34) (cid:88) n ∈ Z N (cid:32) (cid:88) (cid:96) (cid:54) =0 | (cid:96) − n/ (2 N ) | r (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) =: C r ( n/ (2 N )) × (cid:18) (cid:88) (cid:96) (cid:54) =0 | − n + 2 (cid:96)N | r | (cid:98) g ( − n + 2 (cid:96)N ) | (cid:19)(cid:35) / ≤ C r N − r (cid:107) ϕ N (cid:107) H (cid:107) g (cid:107) H r which is valid for any r > / . (We have used above that the series function C r ( z ) is boundedfor | z | ≤ / ) Corollary 4.6.
For any s ≥ , there exists C > so that for any N , (cid:107) Q N L Nk P N f (cid:107) H s ≤ C (cid:107) P N f (cid:107) H s ≤ C (cid:107) f (cid:107) H s . Proof.
The identity L k = 2I − M k , M k := 2 L k ( L − k − I) = 2 L k (K k − i k V k ) (4.19)and the mapping properties V k − i k K k : H s → H s +1 (see for instance, [35, Section 6.2]) yield (cid:107) Q N L k P N f − L k P N f (cid:107) H s = (cid:107) Q N M k P N f − M k P N f (cid:107) H s ≤ C (cid:48) N − (cid:107) M k P N f (cid:107) H s +1 ≤ CN − (cid:107) f (cid:107) H s (4.20)where we have used also (4.4) and (4.5). The result follows readily from the decomposition (cid:107) Q N L Nk P N f (cid:107) H s ≤ (cid:107)L k P N f (cid:107) H s + (cid:107) Q N L Nk P N f − Q N L k P N f (cid:107) H s + (cid:107) Q N L k P N f − L k P N f (cid:107) H s , and the estimates (4.20) and (4.15) in Theorem 4.4.We are ready to prove the convergence of the corresponding block in (4.6).16 roposition 4.7. For any t ≥ , there exists C > such that for any h , N and f ∈ H t , (cid:107) K ΣΓ f − Q h Σ K N ΣΓ P N f (cid:107) H / (Σ) ≤ C ( N − t (cid:107) f (cid:107) H t + h d +1 / (cid:107) f (cid:107) H ) . (4.21) In particular, we have the convergence in operator norm: (cid:107) K ΣΓ f − Q h Σ K N ΣΓ P N f (cid:107) H / → H / (Σ) → as ( N, h ) → ( ∞ , . (4.22) Proof.
For the purpose of this proof, we define R ΣΓ := γ Σ (DL k − i k SL k ) , R N ΣΓ := γ Σ (DL Nk − i k SL Nk ) . (4.23)Clearly R ΣΓ : H → H t (Σ) is continuous, for any t , and, from Lemma 4.5, (see also (3.10)) (cid:107) R ΣΓ Q N ϕ − R N ΣΓ Q N ϕ (cid:107) H m (Σ) ≤ C t N − t (cid:107) Q N ϕ (cid:107) H (4.24)for any t and m , with C independent of ϕ and N . We also have for any compact domain D far away from Γ and for any t ≥ , (cid:107) R ΣΓ ϕ (cid:107) H m ( D ) ≤ C (cid:107) ϕ (cid:107) H − t (4.25)which in particular implies cf. (4.2b) (cid:107) (I − Q h Σ )R ΣΓ ϕ (cid:107) H / (Σ) ≤ Ch d +1 / (cid:107) ϕ (cid:107) H − t , (cid:107) Q h Σ R ΣΓ ϕ (cid:107) H / (Σ) ≤ C (cid:107) ϕ (cid:107) H − t . (4.26)Write now K ΣΓ − Q h Σ K N ΣΓ P N = R ΣΓ L k − Q h Σ R N ΣΓ Q N L Nk P N = (I − Q h Σ )R ΣΓ L k + Q h Σ R ΣΓ L k (I − P N )+ Q h Σ R ΣΓ (I − Q N ) L k P N + Q h Σ R ΣΓ (Q N L k − Q N L Nk )P N + Q h Σ (R ΣΓ − R N ΣΓ )Q N L Nk P N =: L h + L h,N + L h,N + L h,N + L h,N . (4.27)Let us bound these terms now. From (4.26) and the continuity L k : H s → H s , (cid:107) L h f (cid:107) H / (Σ) ≤ Ch d +1 / (cid:107)L k f (cid:107) H ≤ C (cid:48) h d +1 / (cid:107) f (cid:107) H . Proceeding similarly we derive from (4.5), (cid:107) L h,N f (cid:107) H / (Σ) ≤ C (cid:107) (I − P N ) f (cid:107) H − t ≤ CN − t (cid:107) f (cid:107) H , and, by (4.4), (cid:107) L h,N f (cid:107) H / (Σ) ≤ C (cid:107) (I − Q N ) L k P N f (cid:107) H ≤ C (cid:48)(cid:48) N − t (cid:107) f (cid:107) H t . The fourth term is bounded using (4.15) in Theorem 4.4, (cid:107) L h,N f (cid:107) H / (Σ) ≤ C (cid:107) (Q N L k P N − Q N L Nk P N ) f (cid:107) H ≤ C (cid:48)(cid:48) N − t (cid:107) f (cid:107) H t . Finally, (4.26) again, (4.24) with Corollary 4.6 yield, (cid:107) L h,N f (cid:107) H / (Σ) ≤ CN − t (cid:107) Q N L Nk P N f (cid:107) H ≤ CN − t (cid:107) f (cid:107) H . Gathering these bounds yields the estimate (4.21), and consequently (4.22) holds.17 .4 Convergence of the full scheme
We are ready to prove the main result of this paper: stability and convergence for the FEM-BEM numerical algorithm.
Theorem 4.8.
For any N large enough and T h sufficiently fine satisfying Assumption 1, themapping I − K h,N : H / (Σ) × H / → H / (Σ) × H / is uniformly bounded, invertible and with inverse uniformly bounded.Moreover, if ( f Σ , f Γ ) is the solution of (2.9a) and ( f h Σ , f N ) that of (3.13a) , then for any r ≥ and < ε ≤ ε , with ε as in Assumption 1, we have the following estimate, with f = f Γ ◦ x , (cid:107) f Σ − f h Σ (cid:107) H / (Σ) + (cid:107) f − f N (cid:107) H / ≤ C (cid:104) (cid:107) γ Σ u inc − Q h Σ γ Σ u inc (cid:107) H / (Σ) + (cid:107) u inc ◦ x − Q N u inc ◦ x (cid:107) H / + (cid:0) h − εD N − ε + 1 (cid:1) h / (cid:2) inf v h ∈ P h,d (cid:107) K Ω Σ f Σ − v h (cid:107) H (Ω ) + (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) (cid:3) + ( N − ε h d − εD + h dD + N − r ) (cid:107) f Σ (cid:107) L (Σ) + ( N − t + h d +1 / ) (cid:107) f (cid:107) H t (cid:105) (4.28) with C independent of f Σ , f , h and N .Proof. From (4.11) and (4.22) in Propositions 4.3 and 4.7 we conclude (cid:107)K − K h,N (cid:107) H / (Σ) × H / → H / (Σ) × H / → , as ( h, N ) → (0 , ∞ ) which proves the first part of the theorem.On the other hand, for small enough c > , independent of h and N , it holds c (cid:0) (cid:107) f Σ − f h Σ (cid:107) H / (Σ) + (cid:107) f − f N (cid:107) H / (cid:1) ≤ (cid:13)(cid:13)(cid:13)(cid:13) ( I − K h,N ) (cid:20) f Σ − f h Σ f − f N (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) H / (Σ) × H / ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) γ Σ u inc − Q h Σ γ Σ u inc u inc ◦ x − Q N u inc ◦ x (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) H / (Σ) × H / + (cid:13)(cid:13)(cid:13)(cid:13) ( K − K h,N ) (cid:20) f Σ f (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) H / (Σ) × H / (4.29)and the estimate (4.28) follows from Propositions 4.3 and 4.7. Theorem 4.9.
Let ( u, ω ) = (K Ω Σ f Σ , K Ω c1 Γ f Γ ) be the exact solution of (2.9) , and let ( u h , ω N ) :=(K h Ω Σ f h Σ , K N Ω c1 Γ f N ) be that of the FEM-BEM system (3.13) . Then, for any compact set D ⊂ R d \ Ω , r ≥ , t > d + 1 / there exist C such that (cid:107) u − u h (cid:107) H (Ω ) + (cid:107) ω − ω N (cid:107) H r ( D ) ≤ C (cid:0) h d − εD N − ε + h d +1 / + N − t + h dD (cid:1) (cid:107) u inc (cid:107) H t +1 (Ω ) + C inf v h ∈ P h,d (cid:107) u − v h (cid:107) H (Ω ) . (4.30)18 roof. Notice that as consequence of Theorem 2.1, and with f = f ◦ x as before, (cid:107) f Σ (cid:107) H t +1 / (Σ) + (cid:107) f (cid:107) H t +1 / ≤ C t (cid:107) u inc (cid:107) H t +1 (Ω ) , and, cf. (4.2a) and (4.4), (cid:107) γ Σ u inc − Q h Σ γ Σ u inc (cid:107) H / (Σ) ≤ Ch d +1 / (cid:107) u inc (cid:107) H t +1 (Ω ) , (cid:107) u inc ◦ x − Q N u inc ◦ x (cid:107) H / ≤ CN − t (cid:107) u inc (cid:107) H t +1 (Ω ) . Since (cid:107) f Σ − P h Σ f Σ (cid:107) H / (Σ) ≤ C (cid:48) inf p h ∈ γ Σ P h,d (cid:107) f Σ − p h (cid:107) H / (Σ) ≤ C (cid:48) (cid:107) f Σ − f h Σ (cid:107) H / (Σ) , the estimate (4.28) yields (cid:107) f Σ − f h Σ (cid:107) H / (Σ) + (cid:107) f − f N (cid:107) H / ≤ c ( h, N ) (cid:107) f Σ − f h Σ (cid:107) H / (Σ) + C (cid:0) h d − εD N − ε + h dD + N − t + h d +1 / (cid:1) (cid:107) u inc (cid:107) H t +1 (Ω ) + C (cid:48) inf v h ∈ P h,d (cid:107) u − v h (cid:107) H (Ω ) . (4.31)with C ( h, N ) → as ( h, N ) → (0 , ∞ ) . (Actually C (cid:48) above can be shown to tend to zero as ( h, N ) → (0 , ∞ ) ). Taking h sufficiently small and N large enough, say such that c ( h, N ) < / ,the estimate for the error u − u h follows from (4.7) in Theorem 4.1.Regarding the other term, we notice that by Lemma 4.5 (cid:107) ω − ω N (cid:107) H r ( D ) = (cid:107) (DL k − i k SL k ) L k f − (DL Nk − i k SL Nk )Q N L Nk f N (cid:107) H r ( D ) ≤ C (cid:16) (cid:107)L k f − Q N L k f (cid:107) H + (cid:107) Q N L k f − Q N L Nk f N (cid:107) H (cid:17) . (4.32)Clearly, recall that M k , recall (4.19), is a pseudodifferencial operator of order − , (cid:107)L k f − Q N L k f (cid:107) H = (cid:107) (I − Q N )M k f (cid:107) H ≤ CN − t − / (cid:107) f (cid:107) H t +1 / . We claim that for the second term in (4.32) the following bound holds (cid:107) Q N L k f − Q N L Nk f N (cid:107) H ≤ C (cid:107) P N f − f N (cid:107) H / + N − t − / (cid:107) f (cid:107) H t +1 / (4.33)which, with estimate (4.31), should prove the result. Indeed, writing Q N L k f − Q N L Nk f N = Q N L k ( f − P N f ) + (Q N L k P N f − Q N L Nk P N f ) + Q N L Nk (P N f − f N ) , and using (cid:107) Q N L k ( f − P N f ) (cid:107) H ≤ (cid:107)L k ( f − P N f ) (cid:107) H + CN − (cid:107)L k ( f − P N f ) (cid:107) H ≤ C (cid:48) N − t − / (cid:107) f (cid:107) H t +1 / for the first term, (4.15) in Theorem 4.4 (with s = α = 0 ) for the second one and Corollary 4.6for the third term, (4.33) follows. 19 Numerical experiments
This section comprises three sets of numerical experiments to demonstrate our FEM-BEMalgorithm and analysis. As proved in Theorem 4.9, convergence of the FEM-BEM numericalsolution is dictated by the best approximation last term in (4.30). The best approximationaccuracy depends on the smoothness of the exact solution u of (2.1), induced by smoothnessof the refractive index n in the heterogeneous region.The first set of experiments is for the smooth solution case to observe the fast convergenceof our method under optimal conditions. Motivation for the second and third set of experi-ments are from the Janus particle configurations with non-smooth (only piecewise-continuous)refractive indices. For the second and third experiments, the total wave has limited regularityand belongs to H (Ω ) \ H / (Ω) . For the latter cases, we also demonstrate that our FEM-BEMalgorithm converges, with lower convergence rates, and good accuracy can be obtained usinghigh-order finite elements (such as the FEM space spanned by quadratic or cubic splines). Inparticular for the multi-particle Janus-type configuration experiments, we demonstrate thateven quadratic FEM is not sufficient, highlighting the difficulty associated with Janus con-figurations for wave propagation models and the need for efficient high-order FEM-BEM tocompute approximations for the practical QoI such as the DSCS and OA-DSCS. We recall thatthese QoI are defined in (2.4) and in the algorithm, the far-field is computed using the densitybased ansatz in (3.12).Basic setup of the three sets of experiments is similar. Starting from an initial, coarse mesh T H , we consider a set of meshes T H ⊂ T H/ ⊂ · · · T h ⊂ T h/ ⊂ · · · T h e obtained from successiveuniform refinement where each triangle is divided into four elements. For the BEM solver weproceed analogously, consider an initial, relatively low number N = N , and then we doublethe points: { N , N , . . . , N e } .One of aims of the numerical experiments is to demonstrate the error estimates proved inTheorem 4.9. This is carried out as follows: Let T h e be the finest of the chosen finite elementgrids and let N e be the largest of the chosen BEM discrete parameters. We then compute thesolution ( u h e / , ω N e ) using next refinement steps. We will use this pair as our reference exactsolution, i.e. with the notation of Theorem 4.9, ( u h e / , ω N e ) ≈ ( u, ω ) so that we compute (cid:107) u h − u h e / (cid:107) H (Ω ) ; (cid:107) Q h e / ( ω N − ω N e ) (cid:107) H ( D ) , (5.1)as the H − error for the finite element solution (total wave) and boundary element solution(scattered wave). Here, Q h e / is the Lagrange interpolation operator on P h e / ,d which means Q h e / u h e / = u h e / and Q h e / u h = u h for any u h in our list of experiments. Since both ω N and ω are smooth functions as long as D ⊂ ext(Γ) (the region exterior to Γ ), such a strategy givesa sufficiently good estimate for the true error of the numerical solution.In our computations we have taken D in (5.1) to be T h e / ∩ Ext( r Γ) ∩ Ω , with r = 1 . toensure no instabilities for evaluation of potentials in (3.11). Mesh generation in our experimentswere obtained using the open-source package GMSH [18], that is efficient for generation oftriangular meshes on complex geometries Ω . Our choice of a complex Ω for the third setof experiments is inspired by the complex structure of Baby-Yoda (apprentice Grogu fromDisney’s Mandalorian TV series), and we use this geometry as an illustration of the flexibilityand efficiency of our FEM-BEM approach, and the ability to approximate non-smooth fieldsfrom complex heterogeneous media. The total field induced by the Baby-Yoda based piecewisecontinuous refractive index illustrates the complexity of the wave propagation model (2.1).20e simulated numerical experiments using the quadratic ( P ) and cubic ( P ) spline elementsfor FEM, and for several different values for the BEM parameter N . Thus the number of BEMdegrees of freedom (DoF) for the Γ boundary unknown function is N and, throughout thissection, we denote L as the FEM DoF in Ω and M as the number of Dirichlet constrainednodes on Σ . Clearly M << L (since M ≈ √ L ), and because Γ is smooth boundary, thespectral accuracy of the Nyström BEM implies that N << M . Our FEM-BEM frameworkinvolving only relatively small algebraic linear systems (3.13a), for the ( N + M ) boundaryunknowns, were solved iteratively using GMRES with a very low tolerance of − , to ensurethat the reported errors corresponds to the method itself and not from an approximation ofthe algebraic system solutions.For a desired level of accuracy, choices of the discretization parameters crucially depend onthe wavelength of the problem, which in turn depends on the variable coefficient k n of themodel (2.1) and the size of computational regions for the main unknowns (that determine theDoF in algebraic systems). For the heterogeneous and unbounded region model problem (2.1),our efficient FEM-BEM framework reduces the computational regions to just bounded regions Ω and Σ . We recall that if Ω , diam denotes the diameter of the domain Ω and n max is themaximum value of the refractive index n , for our FEM-BEM framework, the FEM interiorproblem wavelength is ( k × n max × Ω , diam ) / (2 π ) ; and the BEM exterior problem wavelengthis ( k × Γ diam ) / (2 π ) , where Γ diam is the length of the artificial boundary curve Γ . Figure 2:
Experiment Γ and Σ . The octagonal shaped polygonaldomain Ω , with boundary Σ , is centered at origin and has circumradius 3. Γ is a smoothrounded-squared curve. The framework illustrate flexibility of artificial curve choices Γ and Σ . In the first set of experiments we choose the example refractive index induced by a Gaussiansmoothing window function in the radial variable r = | x | , x = ( x, y ) ∈ R : n ( x, y ) ≡ n ( r ) = (cid:40) . − r ) , r ∈ [0 , , , r > . r = 1 is less than (cid:15) (where (cid:15) is the machine-epsilon),the refractive index can be considered as numerically smooth. This example function and ournumerical experiments FEM-BEM decomposition framework with artificial curves are demon-strated in Figure 2. For the first set of experiments to demonstrate accuracy of the FEM-BEMapproximations, the incident wave is a plane wave with wavenumber k = 5 and direction (cid:98) d = (1 , . The rounded-square smooth curve Γ in Figure 2 is obtained using the π periodparametrization x ( t ) = 45 √ t ) cos t + (1 + sin t ) sin t, − (1 + cos t ) cos t + (1 + sin t ) sin t ) , t ∈ [0 , π ] . (5.2)The initial (coarsest) grid comprise 4,512 triangles for both P and P elements and we usedup to five uniform refinements. For N , the discrete parameter for the BEM part, the initialvalue N = 20 was refined five times (by doubling). Accuracy of the the first set of simulatedFEM and BEM solutions in the H -norm are displayed in Table 1. According to Theorem 4.9,for the Experiment H -error is dominated by the optimal rate h d in the P d element space for d = 2 , .For fixed and sufficiently high N (= 160) case, this estimated convergence can be observed fromthe last row in Table 1, as the FEM mesh size h is reduced by two with mesh refinement (andcorresponding L increase) the FEM solution errors decrease approximately by (1 / d -times, for d = 2 , . The spectral accuracy of the BEM solution can be observed from the last columnin Table 1 with high accuracy achieved using relative small BEM DoF. For the first set ofexperiments, the GMRES iterations convergence was attained using a small number ( m ) ofiterations with ≤ m ≤ . N / L FEM BEM FEM BEM FEM BEM FEM BEM FEM BEM
020 2.4e-01 6.4e-02 6.1e-02 5.5e-02 1.5e-02 5.4e-02 3.7e-03 5.4e-02 2.2e-05 5.4e-02040 2.4e-01 3.7e-02 6.1e-02 5.7e-03 1.5e-02 7.0e-04 3.7e-03 3.6e-04 2.1e-05 3.5e-04080 2.4e-01 2.8e-02 6.1e-02 4.3e-03 1.5e-02 4.6e-04 3.7e-03 6.7e-05 2.1e-05 4.8e-06160 2.4e-01 2.7e-02 6.1e-02 2.7e-03 1.5e-02 2.4e-04 3.7e-03 3.2e-05 2.1e-05 2.5e-06 N / L FEM BEM FEM BEM FEM BEM FEM BEM FEM BEM
020 1.2e-02 5.6e-02 1.5e-03 5.6e-02 2.6e-04 5.6e-02 2.1e-04 5.6e-02 2.0e-04 5.06-02040 1.2e-02 1.3e-03 1.4e-03 3.8e-04 1.7e-04 3.4e-04 2.1e-05 3.4e-04 1.1e-06 3.4e-04080 1.2e-02 1.1e-03 1.4e-03 1.3e-04 1.7e-04 6.9e-06 2.1e-05 9.2e-08 1.1e-06 2.7e-08160 1.3e-02 1.2e-03 1.4e-03 9.5e-05 1.7e-04 3.5e-06 2.1e-05 7.0e-08 1.1e-06 1.9e-08
Table 1:
Experiment k = 5 ). Estimated H − error for total wave in Ω (FEM-part)and for scattered wave away from Γ (BEM-part) for P (top), P (bottom) elements. For the case k = 5 , the Experiment . .We also performed higher frequency simulations, with k = 10 , , we observed similar FEM-BEM convergence rates similar to that in Table 1 for the smooth solution Experiment k = 5 , , cases, sufficientlyfine FEM mesh and also high-order finite elements are required to obtained good accuracy.22 .2 Experiment Figure 3:
Experiment Janus particles(each having a local piecewise-constant refractive index) of distinct sizes and non-smooth shapes.The FEM octagonal shaped polygonal domain Ω , with artificial boundary Σ , is centered at originand has circumradius 3.5. The BEM smooth artificial boundary Γ is an elliptical curve and theellipse circumscribes the multi-particle heterogeneous Janus configuration. For this set of experiments, a disconnected heterogeneous medium illustrated in Figure 3,is induced by a collection of non-smooth (polygonal) Janus particles of different sizes, and eachparticle is designed using two distinct materials/liquids [31]. As described in Section 1, thisis a model problem for the so-called Janus configuration which has received much attention inrecent years. For the Figure 3 Janus configuration based second set of numerical experiments,we used the following non-smooth (piecewise-continuous) refractive index function, defined inthe full unbounded wave propagation medium, for ( x, y ) ∈ R , as n ( x, y ) = .
333 ( x, y ) ∈ first halves of Janus particles , .
496 ( x, y ) ∈ second halves of Janus particles , Outside the collection of Janus particles . (5.3)In (5.3), we recall that the numbers . and . correspond to, respectively, the refractiveindex of water and toluene at 20C; such two chemicals mix have been used in the literatureto build Janus configurations [31]. Janus configuration based numerical investigations [21,29, 31, 34, 38, 39] have been mainly restricted to simple shapes, and results in this sectionfurther demonstrate that such (shape, size, and number of particles) restrictions can be removedusing efficient algorithms that can accommodate complex heterogeneous structures with efficientdecomposition framework. The simulated Experiment k = 5 , , (corre-sponding to, respectively, approximate interior problem wavelengths . , . , . ) with thecoarse mesh for FEM comprising 1401 triangles. For a fixed incident plane wave case, we tookthe direction to be (cid:98) d = (1 , . For simulating the OA-DSCS of the Janus configuration, we23hose one thousand equally-spaced incident directions (cid:98) d ( φ ) surrounding the configuration thatcorresponds to one thousand of the configuration orientations, with equally spaced orientationangles φ ∈ [0 , π ) , starting from the initial Janus configuration in Figure 3 (with φ = 0 ).The choice of equally spaced direction angles facilitates high-order approximations to theintegral in (2.4) to compute the OA-DSCS using the rectangle quadrature. We recall Remark 2.2to highlight that our discrete FEM-BEM algorithm is efficient for computing solutions with alarge number of incident waves. The first set of experiments for the Janus configuration isperformed with φ = 0 to ensure that high-order approximations of the DSCS integrand in (2.4)can be computed using our FEM-BEM algorithm. We recall that computation of the far-fieldin (2.4) is a smooth post-processing of the field density on Γ , and hence we expect that theaccuracy of far-field approximations should be better, after ensuring convergence and goodaccuracy of the Janus configuration based FEM-BEM solution.Our simulation results (for k = 5 , , ), demonstrating the FEM-BEM solution accuracyfor the ( φ = 0 ) Janus configuration in Figure 3, obtained using the P elements are in Table 2,and the counterpart P elements results are in Table 3. Because of the restricted regularity ofthe exact scattered field induced by only piecewise-continuity of the refractive index in (5.3),it can be observed from Tables 2-3, the limited (second-order) rate of convergence of the FEMsolution does not improve by using P elements compared to that for P elements. However,because of the restricted regularity, it is important to use high-order P elements especiallyfor the higher frequency ( k > ) cases. This is because, for example with k = 20 (despiteusing over one million FEM DoF L ), the P elements based FEM solutions with over errors observed in the last row of Table 2 are not acceptable, and for the same situation the P elements provide much better accurate solutions.In Figure 4 we visualize the total field in Ω and the scattered field in the overlapped region ( R \ Ω ) ∩ Ω , and in Figure 5, we demonstrate the constraint of (numerically) matching theFEM and BEM solutions in the overlapped region by plotting the error in the region in log-scale. These sets of experiments illustrate the computational difficulties associated with Janusconfigurations based wave models, especially taking into account that linear P elements arestandard for pure FEM algorithms (that do not satisfy the SRC) [24].Further advantages of our FEM-BEM algorithm to compute high-order approximationsto smooth far-fields, and hence the DSCS (with fixed incident direction angle φ = 0 ), aredemonstrated in Table 4. Recall that the far-field is a smooth function defined on S , and theuniform norm errors in the far-field were approximated by evaluating the DSCS at equallydistributed points in S . Because of the smooth post-processing of the FEM-BEM solutions tocompute far-fields, similar to the k = 5 results in Table 4, we observed higher accuracy for the k = 10 , cases compared to the FEM-BEM solutions accuracy. We conclude the Experiment k = 5 , , . 24 / L FEM BEM FEM BEM FEM BEM FEM BEM FEM BEM
020 1.6e+00 2.0e+00 3.2e-01 2.7e-01 7.8e-02 2.3e-01 1.9e-02 2.3e-01 5.2e-03 2.3e-01040 1.6e+00 1.9e+00 3.2e-01 1.3e-01 7.8e-02 1.3e-02 1.9e-02 1.2e-03 4.7e-03 6.1e-04080 1.6e+00 1.8e+00 3.1e-01 1.2e-01 7.7e-02 9.6e-03 1.9e-02 8.7e-04 4.7e-03 8.4e-05160 1.6e+00 1.8e+00 3.1e-01 1.1e-01 7.7e-02 7.6e-03 1.9e-02 7.9e-04 4.7e-03 5.4e-05 N / L FEM BEM FEM BEM FEM BEM FEM BEM FEM BEM
020 7.9e+01 7.8e+01 4.9e+01 4.4e+01 4.3e+01 4.0e+01 4.0e+01 4.0e+01 1.6e+01 4.0e+01040 7.2e+01 5.5e+01 1.7e+01 1.7e+01 5.3e+00 5.7e+00 5.6e-01 5.6e-01 5.4e-01 3.8e-02080 7.3e+01 5.5e+01 1.6e+01 1.7e+01 5.1e+00 5.4e+00 4.8e-01 4.8e-01 4.8e-01 3.2e-02160 7.3e+01 5.6e+01 1.6e+01 1.7e+01 4.8e+00 5.1e+00 4.5e-01 4.5e-01 4.5e-01 2.8e-02 N / L FEM BEM FEM BEM FEM BEM FEM BEM FEM BEM
020 6.0e+02 2.3e+02 1.5e+03 5.75e+02 1.1e+03 4.6e+02 1.6e+03 6.3e+02 1.4e+03 5.4e+02040 4.5e+02 1.5e+02 2.8e+02 1.35e+02 3.6e+02 1.7e+02 1.3e+03 7.2e+02 1.6e+03 8.5e+02080 3.3e+02 1.2e+02 2.3e+02 1.14e+02 4.1e+01 2.2e+01 1.2e+01 4.7e+00 1.2e+00 4.6e-01160 3.2e+02 1.2e+02 2.4e+02 1.17e+02 4.0e+01 2.2e+01 1.2e+01 4.5e+00 1.1e+00 4.2e-01
Table 2:
Experiment P elements with k = 5 (top), k = 10 (mid), and k = 20 (top). Estimated H − error for the total wave in Ω (FEM-part) and for the scattered waveaway from Γ (BEM-part). N / L FE BEM FE BEM FE BEM FE BEM FE BEM
020 1.4e-01 3.3e-01 1.7e-02 3.1e-01 3.1e-03 3.1e-01 2.1e-03 3.1e-01 2.9e-03 3.1e-01040 1.3e-01 9.4e-02 1.7e-02 5.2e-03 2.3e-03 8.2e-04 3.9e-04 8.0e-04 7.9e-05 8.0e-04080 1.3e-01 7.8e-02 1.7e-02 4.2e-03 2.3e-03 2.1e-04 3.9e-04 1.1e-05 7.9e-05 1.5e-06160 1.3e-01 7.8e-02 1.7e-02 3.7e-03 2.3e-03 2.0e-04 3.9e-04 1.2e-05 7.9e-05 1.5e-06 N / L FE BEM FE BEM FE BEM FE BEM FE BEM
020 4.9e+01 5.7e+01 4.3e+01 5.2e+01 4.3e+01 5.2e+01 4.3e+01 5.2e+01 4.3e+01 5.2e+01040 8.3e+00 1.2e+01 9.8e-01 1.3e+00 5.2e-02 5.8e-02 7.5e-03 9.6e-03 5.2e-03 8.9e-03080 7.9e+00 1.1e+01 1.0e+00 1.4e+00 4.1e-02 3.7e-02 4.3e-03 1.3e-03 5.8e-04 7.8e-05160 8.0e+00 1.1e+01 9.6e-01 1.3e+00 4.0e-02 3.5e-02 4.2e-03 6.5e-04 5.8e-04 5.0e-05 N / L FE BEM FE BEM FE BEM FE BEM FE BEM
020 1.5e+03 6.7e+02 2.5e+03 1.2e+03 1.4e+03 6.7e+02 1.4e+03 6.5e+02 1.4e+03 6.5e+02040 9.6e+02 5.3e+02 7.5e+02 4.1e+02 1.5e+03 1.0e+03 1.6e+03 1.1e+03 1.6e+03 1.1e+03080 6.3e+02 3.8e+02 6.3e+01 4.1e+01 1.4e+00 6.9e-01 7.9e-02 2.5e-02 1.4e-02 6.5e-03160 6.6e+02 3.9e+02 6.0e+02 3.9e+01 1.4e+00 7.3e-01 7.7e-02 2.4e-02 8.2e-03 9.9e-04
Table 3:
Experiment P elements with k = 5 (top), k = 10 (mid), and k = 20 (top). Estimated H − error for the total wave in Ω (FEM-part) and for the scattered waveaway from Γ (BEM-part). Experiment k = 20 ). Computed total wave (left) and scattered wave (right)using P -finite elements with , elements and , , nodes, and using a spectral BEMwith N = 80 . Figure 5:
Experiment ( R \ Ω ) ∩ Ω . / L N / L Table 4:
Experiment k = 5 , with P (top) and P (bottom) elements and φ = 0 (i.e., with d = (1 , in the incident plane wave). Figure 6: Orientation-Averaged DSCS for Experiment .3 Experiment
In this set of experiments, we demonstrate flexibility of our FEM-BEM algorithm with complexstructured heterogeneous regions. To this end, we consider a Baby-Yoda like domain, depictedin Figure 7 on which the refraction index function is defined. The function is taken to bepiecewise-constant at distinct parts of the domain with values indicated in Figure 7, leadingto the total field solution u with limited regularity. In Figure 7 we also show the curve Γ taken for our BEM computations, which is similar to the smooth curve in Experiment Σ is such that the FEM computational domain Ω is arectangle. Similar to the first two sets of experiments, we simulate the model with wavenumbers k = 5 , , , and for the configuration in Figure 7 these value respectively correspond to . , . , . interior wavelengths region Ω .Figure 7: Experiment
Similar to the first two sets of experiments, sample results in Table 5 demonstrate thepower of our overlapped FEM-BEM algorithm even for the case of simulation of non-smoothwave fields induced by only piecewise-continuous refractive index in complicated heterogeneousregions. Another marked advantage of our FEM-BEM algorithm is the low number of iterationsrequired for convergence of the GMRES to iteratively solve the resulting interface algebraiclinear systems (3.13a). In particular, as can be observed in Table 6, for each fixed wavenumber k the number of required GMRES iterations is independent of the numbers of FEM-BEM DoFs(level of discretizations). In addition, even as the frequency is doubled the required iterationsgrow only mildly (to achieve a fixed error tolerance), and because of the well-conditioning theinterface system (3.13a), even for the high-frequency case, less than GMRES iterations arerequired without using any preconditioner. We conclude the numerical experiments section withvisualizations in Figure 8 for (i) the total field in Ω ; (ii) the scattered field in the overlappedregion Ω = ( R \ Ω ) ∩ Ω ; and (iii) accurate matching of the FEM and BEM solutions in theoverlapped region in Figure 9. 28 / L FE BEM FE BEM FE BEM
020 3.1e-02 3.0e-01 2.0e-02 3.0e-01 1.6e-02 3.0e-01040 2.7e-02 2.2e-03 1.4e-02 2.0e-03 6.2e-03 2.0e-03080 2.7e-02 9.6e-04 1.4e-02 2.4e-04 6.2e-03 5.2e-05160 2.7e-02 9.5e-04 1.4e-02 2.4e-04 6.2e-03 5.2e-05 N / L FE BEM FE BEM FE BEM
020 2.1e+01 1.9e+01 2.1e+01 1.9e+01 2.1e+01 1.9e+01040 8.5e-02 2.8e-02 2.9e-02 5.7e-02 1.3e-02 5.1e-03080 8.5e-02 2.7e-02 2.9e-02 2.1e-03 1.3e-02 2.1e-04160 8.5e-02 2.6e-02 2.9e-02 2.0e-03 1.3e-02 2.1e-04 N / L FE BEM FE BEM FE BEM
020 1.7e+02 1.1e+02 1.7e+02 1.1e+02 1.7e+02 1.1e+02040 5.1e+00 4.3e+00 5.0e+01 4.4e+00 5.0e+01 4.4e+00080 1.0e+00 2.7e-01 1.3e-01 6.5e-03 3.0e-02 8.8e-04160 1.0e+00 2.7e-01 1.3e-01 6.2e-03 3.0e-02 8.6e-04
Table 5:
Experiment P elements with k = 5 (top), k = 10 (mid), and k = 20 (top). Estimated H − error for the total wave in Ω (FEM-part) and for the scattered waveaway from Γ (BEM-part). N / L N / L Table 6:
Experiment − ) using P (top) and P (bottom) elements. For each fixed k , the total number of requiredGMRES iterations is independent of N and L and grows mildly with respect to increase in thewavenumber k . Experiment k = 20 ). Computed total wave (left) and scattered wave (right) using P -finite elements with
478 464 elements and , , nodes, and using a spectral BEM with N = 80 . Figure 9:
Experiment Ω = ( R \ Ω ) ∩ Ω . eferences [1] R.A. Adams and J.J.F. Fournier. Sobolev spaces , volume 140 of
Pure and Applied Mathe-matics (Amsterdam) . Elsevier/Academic Press, Amsterdam, second edition, 2003.[2] S. Bertoluzza. The discrete commutator property of approximation spaces.
C. R. Acad.Sci. Paris Sér. I Math. , 329(12):1097–1102, 1999.[3] H. Brakhage and P. Werner. Über das Dirichletsche Aussenraumproblem für dieHelmholtzsche Schwingungsgleichung.
Arch. Math. , 16:325–329, 1965.[4] S.C. Brenner and L.R. Scott.
The mathematical theory of finite element methods , volume 15of
Texts in Applied Mathematics . Springer, New York, third edition, 2008.[5] P.G. Ciarlet.
The Finite Element Method for Elliptic Problems . Classics in Applied Math-ematics. Society for Industrial and Applied Mathematics, 2002.[6] D. Colton and R. Kress.
Inverse Acoustic and Electromagnetic Scattering Theory . Springer,4th edition, 2019.[7] M. Costabel. Boundary integral operators on Lipschitz domains: elementary results.
SIAMJ. Math. Anal. , 19(3):613–626, 1988.[8] J. Coyle and P. Monk. Scattering of time-harmonic electromagnetic waves by anisotropicin homogeneous scatterers or impenetrable obstacles.
SIAM J. Math. Anal. , 37:1590–1617,2004.[9] P. G. de Gennes. Soft Matter (Nobel Lecture).
Angew. Chem. Int. Ed. Engl. , 31:842–845,1992.[10] V. Domínguez, M. Ganesh, and F.J. Sayas. An overlapping decomposition framework forwave propagation in heterogeneous and unbounded media: Formulation, analysis, algo-rithm, and simulation.
J. Comput. Phys. , 403:109052, 2020.[11] V. Domínguez and F.-J. Sayas. Stability of discrete liftings.
C. R. Math. Acad. Sci. Paris ,337(12):805–808, 2003.[12] V. Domínguez and C. Turc. High order Nyström methods for transmission problems forHelmholtz equations. In
Trends in differential equations and applications , volume 8 of
SEMA SIMAI Springer Ser. , pages 261–285. Springer, [Cham], 2016.[13] M. Ganesh and S. C. Hawkins. Algorithm 975: TMATROM—a T-matrix reduced ordermodel software.
ACM Trans. Math. Softw. , 44:9:1–9:18, 2017.[14] M. Ganesh, S. C. Hawkins, and R. Hiptmair. Convergence analysis with parameter esti-mates for a reduced basis acoustic scattering T-matrix method.
IMA J. Numer. Anal. ,32:1348–1374, 2012.[15] M. Ganesh and C. Morgenstern. High-order FEM-BEM computer models for wave prop-agation in unbounded and heterogeneous media: application to time-harmonic acoustichorn problem.
J. Comput. Appl. Math. , 307:183–203, 2016.3116] M. Ganesh and C. Morgenstern. High-order FEM domain decomposition models for high-frequency wave propagation in heterogeneous media.
Comp. Math. Appl. (CAMWA) ,75:1961–1972, 2018.[17] M. Ganesh and C. Morgenstern. A coercive heterogeneous media Helmholtz model: for-mulation, wavenumber-explicit analysis, and preconditioned high-order FEM.
NumericalAlgorithms , 83:1441–1487, 2020.[18] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities.
International Journal for Numerical Methods inEngineering , 79:1309 – 1331, 2009.[19] A. Gillman, A. H. Barnett, and P.-G. Martinsson. A spectrally accurate direct solutiontechnique for frequency-domain scattering problems with variable media.
BIT NumericalMathematics , 55(1):141–170, 2015.[20] P. Grisvard.
Elliptic problems in nonsmooth domains , volume 69 of
Classics in AppliedMathematics . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2011.[21] S. C. Hawkins, T. Rother, and J. Wauer. Numerical study of acoustic scattering by Janusspheres.
J. Acoust. Soc. Am. , 147:4097, 2020.[22] C. Hazard and M. Lenoir. On the solutions of time-harmonic scattering problems formaxwell’s equations.
SIAM J. Math. Anal. , 27:1597–1630, 1996.[23] G. C. Hsiao and W. L. Wendland.
Boundary integral equations , volume 164 of
AppliedMathematical Sciences . Springer-Verlag, Berlin, 2008.[24] F. Ihlenburg.
Finite element analysis of acoustic scattering , volume 132 of
Applied Math-ematical Sciences . Springer-Verlag, New York, 1998.[25] A. Jami and M. Lenoir. A variational formulation for exterior problems in linear hydro-dynamics.
Comput. Methods Appl. Mech. Engrg. , 16:341–359, 1978.[26] A. Kirsch and P. Monk. Convergence analysis of a coupled finite element and spectralmethod in acoustic scattering.
IMA J. Numer. Anal. , 10(3):425–447, 1990.[27] A. Kirsch and P. Monk. An analysis of the coupling of finite-element and Nyström methodsin acoustic scattering.
IMA J. Numer. Anal. , 14(4):523–544, 1994.[28] R. Kress.
Linear integral equations , volume 82 of
Applied Mathematical Sciences . Springer,New York, third edition, 2014.[29] M. Lattuada and A. Hatton. Synthesis, properties and applications of janus nanoparticles.
Nano Today , 6:286–308, 2011.[30] W. McLean.
Strongly elliptic systems and boundary integral equations . Cambridge Uni-versity Press, Cambridge, 2000. 3231] T. M.Ruhland, A.H.Gröschel, N. Ballard, T. S. Skelhon, A.Walther, A. H. E. Müller, andS. A. F. Bon. Influence of Janus particle shape on their interfacial behavior at liquid-liquidinterfaces.
Langmuir , 29:1388–1394, 2013.[32] J.-C. Nédélec.
Acoustic and electromagnetic equations , volume 144 of
Applied MathematicalSciences . Springer-Verlag, New York, 2001. Integral representations for harmonic problems.[33] J. A. Nitsche and A. H. Schatz. Interior estimates for Ritz-Galerkin methods.
Math.Comp. , 28:937–958, 1974.[34] T. Rother.
Sound Scattering on Spherical Objects . Springer, New York, 2020.[35] J. Saranen and G. Vainikko.
Periodic integral and pseudodifferential equations with nu-merical approximation . Springer Monographs in Mathematics. Springer-Verlag, Berlin,2002.[36] F. J. Sayas. The validity of Johnson-Nédélec’s BEM-FEM coupling on polygonal interfaces.
SIAM Review , 55:131–146, 2013.[37] L.R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfyingboundary conditions.
Math. Comp. , 54(190):483–493, 1990.[38] M. Vafaeezadeh and W. R. Thiel. Janus interphase catalysts for interfacial organic reac-tions.
J. Mol. Liq. , 315:113735, 2020.[39] J. Zhang, A. Gryzbowski, and Granick. Janus particle synthesis, assembly and application.
Langmuir , 33:6964–6977, 2017.
A Sobolev convergence estimates for some pro jections
In this section we collect some useful results for projections on finite element spaces. The firstset of results is concerned with finite element spaces on polygonal compact closed boundaries in R . There finite element spaces are inherited by taking Dirichlet trace of finite element spaceson triangular meshes. We finish this section proving a commutation property for Scott-Zhangtype projections, which is required in Appendix B for deriving superconvergence results of thefinite element solution in stronger norms that are also valid in three dimensions. We think thatmost of the described results here belong to the folklore in finite element analysis. We includeresults below for the sake of completeness. A.1 Projections on polygonal boundary finite element spaces
Let Σ be a polygonal simply connected closed curve with interior Ω ⊂ R . We consider {T h } h ≥ asequence of regular triangular meshes of Ω and denote by { τ h } h ≥ that inherited on Σ . Withoutloss of generality we can assume that h Σ ≈ h , where h Σ and h are the maximum, respectively,of the diameters of the elements of the grids T h and τ h . In addition, let P h,d := { u h ∈ C (Ω) : u h | T h ∈ P d } , γ Σ P h,d Ω and its boundary Σ . We denote by Q h : C (Ω) → P h,d , Q h Σ : C (Σ) → γ Σ P h,d the corresponding nodal (Lagrange) interpolation operators. Note that γ Σ Q h u h = Q h Σ γ Σ u h , ∀ u h ∈ P h,d . Our objectives in this section are twofold: (a) derive convergence estimates in the Sobolevnorms (cid:107) · (cid:107) H s (Σ) for Q h Σ and s ≥ (see (2.5) for the definition we have taken for these spacesfor s > ); (b) define an alternative stable and convergent projection in weaker norms. Inparticular, we are interested in working with functions H / (Σ) , the trace space, for which theinterpolant cannot be defined since the space contains discontinuous functions. Next, we startwith the first objective: Proposition A.1.
For any s ∈ [0 , and s ≤ t < d + 1 with t > / there exists C > suchthat (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch t − s Σ (cid:107) f (cid:107) H t (Σ) . (A.1) Furthermore, for any t > d + 1 (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch d +1 − s Σ (cid:107) f (cid:107) H t (Σ) (A.2) with C > depending only on s ∈ [0 , and t .Proof. We start with a classical result: for s ∈ { , } and t > max { / , s } [5, Th. 3.1.6] (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ C (cid:48) h t − s Σ (cid:107) f Σ (cid:107) H t (Σ) where { Σ (cid:96) } (cid:96) are the edges of Σ and H t (Σ (cid:96) ) the corresponding classical Sobolev space, and (cid:107) f Σ (cid:107) H t (Σ) := (cid:88) (cid:96) (cid:107) f Σ (cid:107) H t (Σ (cid:96) ) . For t > , with t (cid:54) = 1 , , . . . , we then use that the trace operator is also continuous in this norm (cid:107) γ Σ u (cid:107) H t (Σ) ≤ C (cid:107) u (cid:107) H t +1 / (Ω) , ∀ u ∈ H t +1 / (Ω) (A.3)cf. [20, Th. 1.5.2.8] (the bound in (A.3) breaks precisely for t = 0 , , . . . ; see also [23, Th.4.2.7]) which proves (A.1) for the considered non-integer values of t . By interpolation of Sobolevspaces, we can extend the result for any t < d + 1 and s ∈ [0 , , and hence (A.1) holds.The bound (A.2) can be proven by starting from (cid:107) Q h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch d +1 − s Σ (cid:107) f Σ (cid:107) H d +1+ ε (Σ) , with ε ∈ (0 , and again using (A.3).We will focus next on constructing a projection P h Σ : H / (Σ) → γ Σ P h,d with same conver-gence rates in an extended Sobolev scale. A key fact in this construction is the existence of aScott-Zhang-type projection cf. [37] (see also [11]). Specifically, there exists a continuous linearmapping Π h,d : H s (Ω) → P h,d , with s > / , satisfying34. Π h is a projection: Π h v h = v h , ∀ v h ∈ P h,d . (A.4a)2. The image of any element with null trace has null trace as well: γ Σ Π h v = 0 , if γ Σ v = 0 . (A.4b)3. It is quasi-local: for any triangle K ∈ T h , Π h v | T = 0 , if v | S T = 0 , with S T := (cid:91) T (cid:48)∈T h T (cid:48) ∩ T (cid:54) = ∅ T (cid:48) . (A.4c)Furthermore, for any ≤ s ≤ t with s ∈ [0 , / and / < t ≤ d + 1 there exists C = C ( s, t ) so that (cid:107) Π h v − v (cid:107) H s ( T ) ≤ Ch t − sT (cid:107) v (cid:107) H t ( S T ) . (A.4d)The constant C ( s, t ) depends only on the chunkiness parameter of the grid. As conse-quence of (A.4d) we have (cid:34) (cid:88) T ∈T h h s − t ) T (cid:107) Π h v − v (cid:107) H s ( T ) (cid:35) / ≤ C (cid:107) v (cid:107) H t (Ω) . (A.4e)Since for s ∈ [0 , / ∪ [1 , / , we have (cid:107) v (cid:107) H s (Ω) ≤ C (cid:34) (cid:88) T ∈T h (cid:107) v (cid:107) H s ( T ) (cid:35) / , (A.5)the simpler estimate (cid:107) Π h v − v (cid:107) H s (Ω) ≤ Ch t − s (cid:107) v (cid:107) H t (Ω) . (A.6)can be derived from (A.4d), first for s ∈ [0 , / ∪ [1 , / , and then can be extended for s ∈ [1 / , by interpolation of Sobolev spaces. Remark
A.2 . The proof of (A.5) is based on working with the Slobodeckij form of the Sobolevnorm: for non-integer s > and for any Lipschitz domain D ⊂ R m (cid:107) u (cid:107) H s ( D ) = (cid:107) u (cid:107) H s ( D ) + (cid:88) | α | = s (cid:90) D (cid:90) D | ∂ α u ( x ) − ∂ α u ( y ) || x − y | m +2( s − s ) d x d y , where s is the largest integer less than s . Then, for t ∈ (0 , / , it is easy to derive the estimate (cid:107) v (cid:107) H t (Ω) ≤ C t, Ω (cid:34) (cid:88) T ∈T h (cid:107) v (cid:107) H t ( T ) + (cid:90) T | v ( x ) | ρ ( x , ∂T ) t d x (cid:35) where ρ ( x , ∂T ) = inf y ∈ ∂T | x − y | .
35 Hardy-type inequality (see [30, Lemma 3.32]) allows to bound the above integral term by (cid:90) T | v ( x ) | ρ ( x , ∂T ) t d x ≤ C (cid:107) v (cid:107) H t ( T ) . The constant C appearing above depends on t and on the chunkiness parameter of T . Thisinequality does not hold for t ∈ [1 / , since the last integral in the left hand side is expectedto be non-convergent. The result for t > is a simple extension of this argument.The last ingredient is a right inverse of the trace operator R ΣΩ : H s (Σ) → H s +1 / (Ω) for for s ∈ (0 , . For instance, one can take R ΣΩ g Σ := u, with u satisfying ∆ u = 0 , γ Σ u = g Σ , the Dirichlet solution for the Laplace operator. Such operator is continuous, not only for s ∈ (0 , , but it attains the end point s = 1 as well (cf. [30, Chapter 6]; see Theorem 6.12 andthe discussion following it).We are ready to define the desired projection on the finite element space on the boundary γ Σ P h,d . To this end, we first set P h Σ := γ Σ Π h R ΣΩ . Proposition A.3.
Let P h Σ : H s (Σ) → γ Σ P h,d as above. Then (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch t − s Σ (cid:107) f Σ (cid:107) H t (Σ) , ≤ s < , s ≤ t < d + 1 , t > , (A.7) where C is independent of f Σ and h . Furthermore, for any s ∈ [0 , and t > d + 1 , there exists C > such that (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ Ch d +1 − s Σ (cid:107) f (cid:107) H t (Σ) . (A.8) Proof.
By construction, Π h is a projection. Indeed, if Q h : C (Ω) → P h,d is the nodal d − Lagrangeinterpolation operator, and since γ Σ Q h = Q h Σ γ Σ , we have that γ Σ Π h R ΣΩ f h Σ − f h Σ = γ Σ (Π h − Q h )R ΣΩ f h Σ = γ Σ Π h (I − Q h )R ΣΩ f h Σ (cid:124) (cid:123)(cid:122) (cid:125) ∈ H (Ω) = 0 by (A.4b).On the other hand, for t ∈ (1 / , with t ≥ s , and s ∈ (0 , , and by the continuity of thetrace operator, (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) ≤ C (cid:107) Π h R ΣΩ f Σ − R ΣΩ f Σ (cid:107) H s +1 / (Ω) ≤ C (cid:48) h t − s Σ (cid:107) R ΩΣ f Σ (cid:107) H t +1 / (Ω) ≤ C (cid:48)(cid:48) h t − s Σ (cid:107) f Σ (cid:107) H t (Σ) . (A.9)To prove the estimate in H (Σ) = L (Σ) we recall the trace inequality cf [4, Th. 1.6.6] (cid:107) γ Σ u (cid:107) L (Σ) ≤ C (cid:0) (cid:107) u (cid:107) L (Ω) + (cid:107) u (cid:107) L (Ω) (cid:107) u (cid:107) H (Ω) (cid:1) ≤ C (1 + h − ) (cid:107) u (cid:107) L (Ω) + Ch Σ (cid:107) u (cid:107) H (Ω) which yields (cid:107) P h Σ f Σ − f Σ (cid:107) L (Σ) ≤ Ch / (cid:107) Π h R ΩΣ f Σ − R ΩΣ f Σ (cid:107) H (Ω) + Ch − / (cid:107) Π h R ΩΣ f Σ − R ΩΣ f Σ (cid:107) H (Ω) ≤ C (cid:48) h t Σ (cid:107) R ΩΣ f Σ (cid:107) H t +1 / (Ω) ≤ C (cid:48) h t Σ (cid:107) f Σ (cid:107) H t (Σ) .
36e have then proved (A.7) for ≤ s < and max { s, / } < t ≤ and therefore it only remainsto extend this result for t > . But, (cid:107) P h Σ f Σ − f Σ (cid:107) H s (Σ) = (cid:107) P h Σ ( f Σ − Q h Σ f ) − ( f Σ − Q h Σ f ) (cid:107) H s (Σ) ≤ Ch − s Σ (cid:107) f Σ − Q h Σ f Σ (cid:107) H (Σ) and the result follows from Proposition A.1.Note that, using the above arguments, (A.7) cannot be extended to s = 1 . A.2 Commutator properties for Scott-Zhang projections
We end this section showing a commutation property for Π h . We stress out that the result isalso valid in 3D for tetrahedral meshes with minor but direct modifications. Lemma A.4.
For any (cid:36) ∈ C ∞ (Ω) there exists C > so that for d ≥ and s ∈ [0 , / (cid:107) Π h ( (cid:36)w h ) − (cid:36)w h (cid:107) H s (Ω) ≤ Ch (cid:107) w h (cid:107) H s (Ω) , ∀ w h ∈ P h,d . For d = 1 (linear finite elements), we have instead (cid:107) Π h ( (cid:36)w h ) − (cid:36)w h (cid:107) H s (Ω) ≤ Ch min { , − s } (cid:107) w h (cid:107) H s (Ω) , ∀ w h ∈ P h, . Proof.
Consider Q h : C → P h,d the classical nodal interpolant on P h,d which satisfies (cid:107) u − Q h u (cid:107) H s ( T ) ≤ C s,t h t − sT (cid:107) u (cid:107) H t ( T ) , ∀ T ∈ T h (A.10)for any ≤ s < / , and s ≤ t ≤ d + 1 with t > for triangular meshes in bidimensionalpolygonal domains ( t > / for 3D polygonal domains). In (A.10), h T is the diameter of theelement T , and C s,t is a constant independent of T and u h .Recall also the inverse inequalities which hold locally on each triangle: (cid:107) w h (cid:107) H s ( T ) ≤ C s,t h t − sT (cid:107) w h (cid:107) H t ( T ) , ∀ w h ∈ P h,d (A.11)with C depending again only on d , s ≥ t and the chunkiness parameter of the grid. Then,locally on each element it holds (cid:107) (cid:36)w h (cid:107) H d +1 ( T ) ≤ C (cid:36) (cid:107) w h (cid:107) H d +1 ( T ) = C (cid:36) (cid:107) w h (cid:107) H d ( T ) ≤ C (cid:48) (cid:36),s h s − d − T (cid:107) w h (cid:107) H s ( T ) for any s ∈ [0 , d ] .We start with the s ∈ [0 , / case. Using (A.5), the convergence properties of Π h (A.4d),the interpolant operator Q h and the inverse inequality (A.11), (cid:107) Π h ( (cid:36)w h ) − (cid:36)w h (cid:107) H s (Ω) ≤ C s (cid:88) T ∈T h (cid:107) Π h ( (cid:36)w h − Q h ( (cid:36)w h )) − ( (cid:36)w h − Q h ( (cid:36)w h )) (cid:107) H s ( T ) ≤ C (cid:48) s (cid:88) T ∈T h h − sT (cid:107) (cid:36)w h − Q h ( (cid:36)w h ) (cid:107) H ( T ) ≤ C (cid:48) s (cid:88) T ∈T h h d +2 − sT (cid:107) (cid:36)w h (cid:107) H d +1 ( T ) ≤ C (cid:48)(cid:48) s,(cid:36) (cid:88) T ∈T h h T (cid:107) w h (cid:107) H s ( T ) ≤ C (cid:48) s,(cid:36) h (cid:107) w h (cid:107) H s (Ω) . s ∈ [1 , / , we proceed in a similar manner: Using the stability of the Scott-Zhangtype projection (cid:107) Π h ( (cid:36)w h ) − (cid:36)w h (cid:107) H s (Ω) ≤ C s (cid:107) Q h ( (cid:36)w h ) − (cid:36)w h (cid:107) H s (Ω) ≤ C (cid:48) s (cid:88) T ∈T h h d − sT (cid:107) (cid:36)w h (cid:107) H d +1 ( T ) ≤ C s,(cid:36) (cid:88) T ∈T h h d − sT (cid:107) w h (cid:107) H d ( T ) . The result follows now using either the direct bound (cid:107) w h (cid:107) H d ( T ) ≤ (cid:107) w h (cid:107) H s ( T ) , for d = 1 , or the inverse inequality (A.11) for d ≥ .We have then proven the result for s ∈ [0 , / ∪ [1 , / . For the intermediate values of s ∈ (1 / , we just invoke the Interpolation Theory of Sobolev spaces. Remark
A.5 . We note that in [2] the weaker estimate is proven: (cid:107) Π h ( (cid:36)v h ) − (cid:36)v h (cid:107) H s (Ω) ≤ Ch s (cid:107) v h (cid:107) H s (Ω) , ∀ v h ∈ P h,d for s ∈ (1 / , and for any projection on finite element spaces satisfying rather general assump-tions which include, in particular, our case. B Convergence/superconvergence of FEM in stronger norms
In this appendix Ω is a polygonal domain with boundary Σ in both R m with m = 2 , and b n, Ω ( u, v ) := ( ∇ u, ∇ v ) Ω − ( u, v ) n, Ω := (cid:90) Ω ∇ u · ∇ v − (cid:90) Ω n uv is the bilinear form associated to the Dirichlet problem (cid:12)(cid:12)(cid:12)(cid:12) ∆ u + nu = fγ Σ u = g Σ . Here − n is a L − function with compact support Ω in Ω . As before we consider a sequence ofregular grids made up of conformal triangular/tetrahedral elements T h . The parameter h refersagain to the maximum of the diameters of the elements in such a way that we write h → tomean that the diameters of the elements tend to zero. On T h we construct a continuous finiteelement space of the elements which are polynomials of degree d on each T ∈ T h and denote itby P h,d .We will work with the numerical solution given by the usual FEM scheme: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u h ∈ P h,d b n, Ω ( u h , v h ) = − ( f, v h ) Ω , ∀ v h ∈ P h,d ∩ H (Ω) γ Σ u = g h Σ . Here, γ Σ P h,d (cid:51) g h Σ ≈ g Σ . Some preliminary results can be listed at this point. Provided that g Σ ∈ H / (Σ) , f ∈ L (Ω) then u ∈ H (Ω) with (cid:107) u (cid:107) H (Ω) ≤ C (cid:0) (cid:107) f (cid:107) L (Ω) + (cid:107) g Σ (cid:107) H / (Σ) (cid:1) .
38e recall the standard error result [11, 37] (cid:107) u − u h (cid:107) H (Ω) ≤ C (cid:104) inf v h ∈ P h,d (cid:107) u − v h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) H / (Ω) (cid:105) which gives us a convergence estimate in the natural norm. We however want to exploreconvergence estimates in stronger norms in the subdomains where the solution u is more regular.The following theorem presents the main result of this appendix. The proof of this result usessimilar ideas as those presented in [33] for proving superconvergence in H for a more generalclass of partial differential equations. We include it in this work for the sake of completeness. Theorem B.1.
Let D (cid:48) , D be two domains with D (cid:48) ⊂ D (cid:48) ⊂ D ⊂ D ⊂ Ω \ Ω and let {T h } bea sequence of regular grids with h → which are locally quasi-uniform in D (cid:48) . Then for any ε ∈ [0 , / , u ∈ H r +1 ( D (cid:48) ) with r ∈ [ ε, d ] , there exists C > depending on ε , r , D , D (cid:48) so thatfor any grid fine enough T h it holds (cid:107) u − u h (cid:107) H ε ( D (cid:48) ) ≤ C (cid:2) h − εD (cid:107) u − u h (cid:107) H ( D ) + h r − εD (cid:107) u (cid:107) H r +1 ( D ) + h − εD (cid:107) u − u h (cid:107) H ( D ) (cid:3) (B.1) where h D is the maximum of the diameters of the triangles contained in D .Proof. Let Π h be a Scott-Zhang type projection satisfying (A.4), and hence the commutatorproperty stated in Lemma A.4. Take a smooth cut-off function (cid:36) satisfying (cid:36) | D (cid:48) ≡ , supp (cid:36) ⊂ D, Π h ( (cid:36)v h ) | D (cid:48) = v h | D (cid:48) , ∀ v h ∈ P h,d , Π h ( (cid:36)v ) | D ∈ H ( D ) , ∀ v ∈ H ( D ) . (B.2)We consider the H ( D ) inner product β D ( u, v ) := ( ∇ u, ∇ v ) D + ( u, v ) D , and set e h := u − u h , (cid:101) e h := Π h ( (cid:36)u ) − Π h ( (cid:36)u h ) ∈ P h,d , (cid:101) e h := (cid:36)u − Π h ( (cid:36)u ) . Using (B.2), (cid:36)e h | D (cid:48) = ( (cid:101) e h + (cid:101) e h ) | D (cid:48) . (B.3)Then (cid:107) e h (cid:107) H ε ( D (cid:48) ) ≤ (cid:107) (cid:36)e h (cid:107) H ε ( D ) ≤ (cid:107) (cid:101) e h (cid:107) H ε ( D ) + (cid:107) (cid:101) e h (cid:107) H ε ( D ) , (B.4)For bounding the first term we start applying the inverse inequality (cid:107) (cid:101) e h (cid:107) H ε ( D ) ≤ C ε h − εD (cid:107) (cid:101) e h (cid:107) H ( D ) = C ε h − εD β D ( (cid:101) e h , v h ) , with v h = 1 (cid:107) (cid:101) e h (cid:107) H ( D ) (cid:101) e h ∈ P h,d ∩ H ( D ) , (B.5)and continue applying the decomposition (cid:101) e h = ( (cid:36)u − (cid:36)u h ) + (Π h − I)( (cid:36)u − (cid:36) Π h u ) + (Π h − I)( (cid:36) Π h u − (cid:36)u h ) . (cid:107) (cid:101) e h (cid:107) H ε ( D ) ≤ C ε h − εD (cid:16) β D ( (cid:36)u − (cid:36)u h , v h ) + β D ((Π h − I)( (cid:36)u − (cid:36) Π h u ) , v h )+ β D ((Π h − I)( (cid:36) Π h u − (cid:36)u h ) , v h ) (cid:17) ≤ C ε h − εD (cid:16) β D ( (cid:36)e h , v h ) + (cid:107) (Π h − I)( (cid:36)u − (cid:36) Π h u ) (cid:107) H ( D ) + (cid:107) (Π h − I)( (cid:36) Π h u − (cid:36)u h ) (cid:107) H ( D ) (cid:17) ≤ C ε h − εD (cid:0) β D ( (cid:36)e h , v h ) + (cid:107) (cid:36)u − (cid:36) Π h u (cid:107) H ( D ) + h D (cid:107) (Π h u − u ) + ( u − u h ) (cid:107) H ( D ) (cid:1) ≤ C ε h − εD (cid:0) β D ( (cid:36)e h , v h ) + (cid:107) Π h u − u (cid:107) H ( D ) + h D (cid:107) e h (cid:107) H ( D ) (cid:1) , (B.6)where to bound the last term, we used the the commutator property, cf. Lemma A.4 with w h = Π h u − u h . Next, we consider the first term in (B.6) and obtain a bounding estimate.Observe that β D ( (cid:36)e h , v h ) = (cid:90) D e h ∇ (cid:36) · ∇ v h − (cid:90) D v h ∇ e h · ∇ (cid:36) + β D ( e h , (cid:36)v h )= 2 (cid:90) D e h ∇ (cid:36) · ∇ v h + (cid:90) D v h e h ∆ (cid:36) + β D ( e h , (cid:36)v h − Π h ( (cid:36)v h )) + β D ( e h , Π h ( (cid:36)v h )) . Thus, and after applying again Lemma A.4 to the third term above, β D ( (cid:36)e h , v h ) ≤ C (cid:107) e h (cid:107) L ( D ) + Ch D (cid:107) e h (cid:107) H ( D ) + β D ( e h , Π h ( (cid:36)v h )) . (B.7)Finally, using the orthogonality relation for the Galerkin solution, β D ( e h , Π h ( (cid:36)v h )) = (cid:90) Ω ∇ e h · ∇ Π h ( (cid:36)v h ) − (cid:90) Ω n e h Π h ( (cid:36)v h ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + (cid:90) D (1 + n ) e h Π h ( (cid:36)v h ) ≤ C (cid:107) e h (cid:107) L ( D ) (cid:107) (cid:36)v h (cid:107) L (Ω) ≤ C (cid:107) e h (cid:107) L ( D ) . (B.8)Plugging (B.8) in (B.7) and next (B.7) in (B.6) yield (cid:107) (cid:101) e h (cid:107) H ε ( D (cid:48) ) ≤ Ch − εD (cid:0) (cid:107) e h (cid:107) L ( D ) + (cid:107) Π h u − u (cid:107) H ( D ) + h D (cid:107) e h (cid:107) L ( D ) (cid:1) (B.9)Therefore, (B.9) with (B.4) prove (cid:107) e h (cid:107) H ε ( D (cid:48) ) ≤ C ε (cid:0) h − εD (cid:107) e h (cid:107) L ( D ) + h − εD (cid:107) e h (cid:107) H ( D ) + h − εD (cid:107) Π h u − u (cid:107) H ( D ) + (cid:107) Π h ( ωu ) − ωu (cid:107) H ε ( D ) (cid:1) , (B.10)for any ε ∈ [0 , / . The result follows readily from the convergence estimates for Π h .A convergence estimate in L (Ω) is needed to obtain a convergence result of the finite elementsolution in D . This is done next using a variant of the classical Aubin-Nische argument40 emma B.2 (Aubin-Nitsche trick) . Then for δ ∈ (1 / , it holds (cid:107) u − u h (cid:107) L (Ω) ≤ C (cid:104) h δ (cid:107) u − u h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) L (Σ) (cid:105) where h denotes the maximum of the diameters of the elements of the mesh T h of Ω .Proof. Let us denote as before e h := u − u h and take w the (unique) solution of (cid:12)(cid:12)(cid:12)(cid:12) w ∈ H (Ω)∆ w + n w = − e h or in variational form , (cid:12)(cid:12)(cid:12)(cid:12) w ∈ H (Ω) b Ω ,n ( w, v ) = ( e h , v ) Ω , ∀ v ∈ H (Ω) . It is known that there exists δ ∈ (1 / , such that (cid:107) w (cid:107) H δ (Ω) ≤ C δ (cid:107) e h (cid:107) L (Ω) with C δ independent of e h .Then (cid:107) e h (cid:107) L (Ω) = ( e h , e h ) Ω = − (∆ w + n w, e h ) Ω = b Ω ,n ( w, e h ) + (cid:90) Σ ∂ n w e h = b Ω ,n ( w − Π h w, e h ) + (cid:90) Σ ∂ n we h where we have made use of the Galerkin orthogonality of the discrete solution and the fact that Π h w ∈ P h,d ∩ H (Ω) for w ∈ H (Ω) . Notice that from the trace theorem (cid:107) ∂ n w (cid:107) L (Σ) ≤ (cid:107) γ Σ ∇ w (cid:107) L (Σ) ≤ C (cid:107) w (cid:107) H δ (Ω) . Therefore, (cid:107) e h (cid:107) L (Ω) ≤ C (cid:2) (cid:107) w − Π h w (cid:107) H (Ω) (cid:107) e h (cid:107) H (Ω) + (cid:107) w (cid:107) H δ (Ω) (cid:107) e h (cid:107) L (Σ) (cid:3) ≤ C (cid:48) (cid:2) h δ (cid:107) e h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) L (Σ) (cid:3) (cid:107) w (cid:107) H δ (Ω) ≤ C (cid:48)(cid:48) (cid:2) h δ (cid:107) e h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) L (Σ) (cid:3) (cid:107) e h (cid:107) L (Ω) . The proof is now completed.The parameter δ in Lemma B.2 is nothing but the extra regularity grade, from a Sobolevpoint of view, of the homogeneous Dirichlet problem ∆ w + n w = g, γ Σ w = 0 . It is a very well established result, see for instance [20], that δ ∈ (1 / , and it attains to be for convex polygons in R . Corollary B.3.
Under the same assumptions as those taken in Theorem B.1 and Lemma B.2,we have that for any ε ∈ [0 , / (cid:107) u − u h (cid:107) H ε ( D (cid:48) ) ≤ C (cid:0) ( h δ h − εD + h − εD ) (cid:107) u − u h (cid:107) H (Ω) + h − εD (cid:107) g Σ − g h Σ (cid:107) L (Σ) + h r − εD (cid:107) u (cid:107) H r +1 ( D ) (cid:1) . In particular, if h δ h − εD → and g h Σ is taken satisfying h − εD (cid:107) g Σ − g h Σ (cid:107) L (Σ) ≤ C ( h ) (cid:107) g Σ − g h Σ (cid:107) H / (Σ) with C ( h ) → as h → , there exists η ( h ) with η ( h ) → as h → such that (cid:107) u − u h (cid:107) H ε ( D (cid:48) ) ≤ η ( h ) (cid:0) inf v h ∈ P h,d (cid:107) u − v h (cid:107) H (Ω) + (cid:107) g Σ − g h Σ (cid:107) H / (Σ) (cid:1) ..