An overlapping decomposition framework for wave propagation in heterogeneous and unbounded media: Formulation, analysis, algorithm, and simulation

TL;DR

This paper introduces an overlapping decomposition framework that accurately models wave propagation in heterogeneous, unbounded media by combining finite element and boundary element methods, avoiding approximation of the radiation condition.

Contribution

The work presents a novel overlapping FEM-BEM algorithm that exactly incorporates the radiation condition for wave propagation in complex media, enhancing simulation accuracy and efficiency.

Findings

Efficient FEM-BEM algorithm for 2D wave propagation.

Exact incorporation of radiation condition in simulations.

Effective handling of complex heterogeneous media.

Abstract

A natural medium for wave propagation comprises a coupled bounded heterogeneous region and an unbounded homogeneous free-space. Frequency-domain wave propagation models in the medium, such as the variable coefficient Helmholtz equation, include a faraway decay radiation condition (RC). It is desirable to develop algorithms that incorporate the full physics of the heterogeneous and unbounded medium wave propagation model, and avoid an approximation of the RC. In this work we first present and analyze an overlapping decomposition framework that is equivalent to the full-space heterogeneous-homogenous continuous model, governed by the Helmholtz equation with a spatially dependent refractive index and the RC. Our novel overlapping framework allows the user to choose two free boundaries, and gain the advantage of applying established high-order finite and boundary element methods (FEM and BEM) to simulate an equivalent coupled model. The coupled model comprises auxiliary interior bounded heterogeneous and exterior unbounded homogeneous Helmholtz problems. A smooth boundary can be chosen for simulating the exterior problem using a spectrally accurate BEM, and a simple boundary can be used to setup a high-order FEM for the interior problem. Thanks to the spectral accuracy of the exterior computational model, the resulting coupled system in the overlapping region is relatively very small. Using the decomposed equivalent framework, we develop a novel overlapping FEM-BEM algorithm for simulating the acoustic or electromagnetic wave propagation in two dimensions. Our FEM-BEM algorithm for the full-space model incorporates the RC exactly. Numerical experiments demonstrate the efficiency of the FEM-BEM approach for simulating smooth and non-smooth wave fields, with the latter induced by a complex heterogeneous medium and a discontinuous refractive index.

Tables

Table 1. Table 1: Experiment 1: ℙ 3 subscript ℙ 3 \mathbb{P}_{3} Finite element space and k = π / 4 , π , 4 ​ π 𝑘 𝜋 4 𝜋 4 𝜋 k=\pi/4,\pi,4\pi (top, middle, bottom tables). In the first row and the first column, L 𝐿 {L} and 2 ​ N 2 𝑁 2N are the number of degrees of freedom used to compute the FEM and BEM solutions, respectively. The number of GMRES iterations required for solving the system, with a residual tolerance of 10 − 8 superscript 10 8 10^{-8} , is given within the parenthesis. Estimated (relative) uniform errors in the far-field are given in columns two to five.

$N/L$	7,999		31,657		125,953		502,465		2,007,169
010	3.1e-03	(012)	6.6e-05	(012)	2.2e-06	(012)	1.2e-06	(012)	1.2e-06	(012)
020	3.1e-03	(012)	6.5e-05	(012)	2.0e-06	(012)	2.5e-10	(012)	4.7e-11	(012)
040	3.1e-03	(012)	6.5e-05	(012)	2.0e-06	(012)	1.8e-10	(012)	1.4e-11	(012)
080	3.1e-03	(012)	6.4e-05	(012)	2.0e-06	(012)	1.5e-10	(012)	9.0e-12	(012)

Table 2. Table 2: Experiment 1: ℙ 4 subscript ℙ 4 \mathbb{P}_{4} Finite element space and k = π / 4 , π , 4 ​ π 𝑘 𝜋 4 𝜋 4 𝜋 k=\pi/4,\pi,4\pi (top, middle, bottom tables). In the first row and the first column, L 𝐿 {L} and 2 ​ N 2 𝑁 2N are the number of degrees of freedom used to compute the FEM and BEM solutions, respectively. The number of GMRES iterations required for solving the system, with a residual tolerance of 10 − 8 superscript 10 8 10^{-8} , is given within the parenthesis. Estimated (relative) uniform errors in the far-field are given in columns two to five.

$N/L$	14,145		56,129		223,617		892,673		3,567,105
010	3.9e-04	(012)	9.4e-06	(012)	1.4e-06	(012)	1.2e-06	(012)	1.2e-06	(012)
020	3.9e-04	(012)	8.9e-06	(012)	2.5e-07	(012)	6.9e-10	(012)	8.4e-11	(012)
040	3.9e-04	(012)	8.9e-06	(012)	2.5e-07	(012)	7.0e-10	(012)	1.0e-10	(012)
080	3.9e-04	(012)	8.9e-06	(012)	2.5e-07	(012)	7.0e-10	(012)	9.9e-11	(012)

Table 3. Table 3: Experiment 2: ℙ 3 subscript ℙ 3 \mathbb{P}_{3} Finite element space and k = π / 4 , π , 4 ​ π 𝑘 𝜋 4 𝜋 4 𝜋 k=\pi/4,\pi,4\pi (top, middle, bottom tables). In the first row and the first column, L 𝐿 {L} and 2 ​ N 2 𝑁 2N are the number of degrees of freedom used to compute the FEM and BEM solutions, respectively. The number of GMRES iterations required for solving the system, with residual tolerance of 10 − 8 superscript 10 8 10^{-8} , is given within the parenthesis. Estimated (relative) uniform errors in the far-field are given in columns two to five.

$N/L$	39,085		69,381		622,573		2,488,441
010	2.8e-03	(015)	2.8e-03	(015)	2.8e-03	(015)	2.8e-03	(015)
020	5.8e-05	(015)	8.4e-07	(015)	8.4e-07	(015)	8.4e-07	(015)
040	5.3e-05	(015)	1.0e-07	(015)	6.1e-09	(015)	6.9e-10	(015)
080	5.8e-05	(015)	7.4e-08	(015)	6.5e-09	(015)	4.4e-10	(015)

Table 4. Table 4: Experiment 2: ℙ 4 subscript ℙ 4 \mathbb{P}_{4} Finite element space and k = π / 4 , π , 4 ​ π 𝑘 𝜋 4 𝜋 4 𝜋 k=\pi/4,\pi,4\pi (top, middle, bottom tables). In the first row and the first column, L 𝐿 {L} and 2 ​ N 2 𝑁 2N are the number of degrees of freedom used to compute the FEM and BEM solutions, respectively. The number of GMRES iterations required for solving the system, with residual tolerance of 10 − 8 superscript 10 8 10^{-8} , is given within the parenthesis. Estimated (relative) uniform errors in the far-field are given in columns two to five.

$N/L$	69,381		276,905		1,106,385		4,423,073
010	2.8e-03	(015)	2.8e-03	(015)	2.8e-03	(015)	2.8e-03	(015)
020	1.3e-06	(015)	8.4e-07	(015)	8.4e-07	(015)	8.4e-07	(015)
040	1.3e-06	(015)	1.6e-07	(015)	6.8e-10	(015)	6.8e-10	(015)
080	1.3e-06	(015)	1.6e-07	(015)	6.9e-10	(015)	6.8e-10	(015)

Equations

$$\left|\begin{array}[]{rcl}\Delta u+k^{2}n^{2}\>u&=&0,\quad\text{in }\mathbb{R}^{m},\\ \partial_{r}u^{s}-\mathrm{i}k{u}^{s}&=&o(|r|^{\frac{m+1}{2}}),\quad\text{as }|r|\to\infty.\end{array}\right.$$

$$\left|\begin{array}[]{rcl}\Delta w+k^{2}n^{2}\>w&=&0,\quad\text{in }\Omega_{2},\\ \gamma_{\Sigma}w&=&f_{\Sigma}^{\rm inp}.\end{array}\right.$$

$${\mathrm{K}}_{D\Sigma }{f}_{\Sigma }^{\mathrm{inp}}:=w{|}_D.$$

$$\left|\begin{array}[]{rcl}\Delta\widetilde{\omega}+k^{2}\widetilde{\omega}&=&0,\quad\text{in }\Omega_{1}^{\rm c},\\ \gamma_{\Gamma}\widetilde{\omega}&=&f_{\Gamma}^{\rm inp},\\ \partial_{r}\widetilde{\omega}-{\rm i}k\widetilde{\omega}&=&o(|r|^{(m-1)/2}).\end{array}\right.$$

$${\mathrm{K}}_{D\Gamma }{f}_{\Gamma }^{\mathrm{inp}}:=\tilde{\omega }{|}_D,$$

$$\left|\begin{array}[]{ccccrcl}\lx@intercol(f_{\Sigma},f_{\Gamma})\in H^{1/2}(\Sigma)\times{H^{1/2}(\Gamma)}\hfil\lx@intercol\\[4.73611pt] f_{\Sigma}&-&\mathrm{K}_{\Sigma\Gamma}f_{\Gamma}&=&\gamma_{\Sigma}u^{\rm inc}\\ -\mathrm{K}_{\Gamma\Sigma}f_{\Sigma}&+&f_{\Gamma}&=&-\gamma_{\Gamma}u^{\rm inc}\end{array}\right.$$

$$u:=\{\begin{array}{ll}{\mathrm{K}}_{{\Omega }_2\Sigma }{f}_{\Sigma },& \text{in Ω2},\\ {\mathrm{K}}_{{\Omega }_1^{\mathrm{c}}\Gamma }{f}_{\Gamma }+u^{\mathrm{inc}},& \text{in Ω1c}.\end{array}$$

$${\mathrm{K}}_{{\Omega }_{12}\Sigma }{f}_{\Sigma }=u^{\mathrm{inc}}{|}_{{\Omega }_{12}}+{\mathrm{K}}_{{\Omega }_{12}\Gamma }{f}_{\Gamma },$$

$$\left|\begin{array}[]{rcl}\Delta v+k^{2}v&=&0,\quad\text{in }\Omega_{12},\\ \gamma_{\Gamma}v&=&0,\quad\gamma_{\Sigma}v\,=\,0\end{array}\right.$$

$$\Vert f{\Vert }_{H^s(D)}^2:=\sum _{|\mathbf{\alpha }|\le s}{\int }_D|{\partial }_{\mathbf{\alpha }}f{|}^2$$

$$\Vert f{\Vert }_{H^s(D)}^2:=\Vert f{\Vert }_{H^{s_0}(D)}^2+\sum _{|\mathbf{\alpha }|\le s_0}{\int }_D{\int }_D\frac{|{\partial }_{\mathbf{\alpha }}f(\mathbf{x})-{\partial }_{\mathbf{\alpha }}f(\mathbf{y}){|}^2}{|\mathbf{x}-\mathbf{y}{|}^{m+2\beta }}\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y}.$$

$$H^s(D):=\{f\in L^2(D):\Vert f{\Vert }_{H^s(D)}<\infty \},$$

$$\Vert \phi {\Vert }_{H^s(\partial D)}^2:=\sum _{j=1}^J\Vert ({\mu }^j\phi )\circ {\mathbf{x}}^j{\Vert }_{H^s({\mathbb{R}}^{m-1})}^2.$$

$$H^s(\partial D):=\{\phi \in L^2(\partial D):\Vert \phi {\Vert }_{H^s(\partial D)}<\infty \}.$$

$$H^s(\partial D):=\{{\gamma }_{\partial D}u:u\in H^{s+1/2}(D)\}$$

$$\Vert \phi {\Vert }_{H^s(\partial D)}:=\underset{\genfrac{}{}{0ex}{}{0\ne u\in H^{s+1/2}(D)}{{\gamma }_{\Sigma }u=\phi }}{\inf }\Vert u{\Vert }_{H^{s+1/2}(D)}.$$

$${\Phi }_k(\mathbf{x},\mathbf{y}):=\{\begin{array}{ll}\frac{\mathrm{i}}{4}{H}_0^{(1)}(kr),& \mathbf{x},\mathbf{y}\in {\mathbb{R}}^2,\\ \frac{1}{4\pi r}\exp (\mathrm{i}kr),& \mathbf{x},\mathbf{y}\in {\mathbb{R}}^3,\end{array}$$

$$({\mathrm{SL}}_k\phi )(\mathbf{x}):={\int }_{\Gamma }{\Phi }_k(\mathbf{x}-\mathbf{y})\phi (\mathbf{y})\mathrm{d}{\sigma }_{\mathbf{y}},({\mathrm{DL}}_{k}g)(\mathbf{x}):={\int }_{\Gamma }{\partial }_{\mathbf{\nu }(\mathbf{y})}{\Phi }_k(\mathbf{x}-\mathbf{y})g(\mathbf{y})\mathrm{d}{\sigma }_{\mathbf{y}},\mathbf{x}\in {\mathbb{R}}^m\setminus \Gamma ,$$

$${\mathrm{V}}_k\phi$$

$${\mathrm{K}}_{k}g$$

$$\frac{1}{2}\mathrm{I}+{\mathrm{K}}_k-\mathrm{i}k{\mathrm{V}}_{k,\sigma }:H^s(\Gamma )\to H^s(\Gamma ).$$

$${\cal L}_{k,\sigma}:=\Big{(}\tfrac{1}{2}{\rm I}+\mathrm{K}_{k}-\mathrm{i}k\,\mathrm{V}_{k,\sigma}\Big{)}^{-1}{:H^{s}(\Gamma)\to H^{s}(\Gamma)}$$

$${\mathrm{K}}_{{\Omega }_1^{\mathrm{c}}\Gamma }=({\mathrm{DL}}_k-\mathrm{i}k{\mathrm{SL}}_{k,\sigma }){\mathcal{L}}_{k,\sigma }.$$

$${\mathrm{K}}_{{\Omega }_1^{\mathrm{c}}\Gamma }:H^s(\Gamma )\to H_{\mathrm{loc}}^{s+1/2}({\Omega }_1^{\mathrm{c}})$$

$${\mathrm{K}}_{D\Gamma }:H^s(\Gamma )\to H^r(D).$$

$${\mathrm{K}}_{\Sigma \Gamma }:H^s(\Gamma )\to H^r(\Sigma )$$

$${\mathrm{DL}}_k:H^s(\Gamma )\to H_{\mathrm{loc}}^{s+1/2}({\Omega }_1^{\mathrm{c}}),{\mathrm{SL}}_k:H^{s-1}(\Gamma )\to H_{\mathrm{loc}}^{s+1/2}({\Omega }_1^{\mathrm{c}}),$$

$$\Vert {\mathrm{K}}_{{\Omega }_2\Sigma }{f}_{\Sigma }^{\mathrm{inp}}{\Vert }_{H^1({\Omega }_2)}\le C\Vert f_{\Sigma }^{\mathrm{inp}}{\Vert }_{H^{1/2}(\Sigma )},$$

$$\Vert {\mathrm{K}}_{{\Omega }_2\Sigma }{f}_{\Sigma }^{\mathrm{inp}}{\Vert }_{H^{s+1/2}({\Omega }_2)}\le C\Vert f_{\Sigma }^{\mathrm{inp}}{\Vert }_{H^s(\Sigma )}.$$

$${\mathrm{K}}_{D\Sigma }:H^0(\Sigma )\to H^r(D).$$

$${\mathrm{K}}_{\Gamma \Sigma }:H^0(\Sigma )\to H^r(\Gamma )$$

Full text

The coupled model comprises auxiliary interior bounded heterogeneous and exterior unbounded homogeneous Helmholtz problems. A smooth boundary can be chosen for simulating the exterior problem using a spectrally accurate BEM, and a simple boundary can be used to setup a high-order FEM for the interior problem. Thanks to the spectral accuracy of the exterior computational model, the resulting coupled system in the overlapping region is relatively very small. Using the decomposed equivalent framework, we develop a novel overlapping FEM-BEM algorithm for simulating the acoustic or electromagnetic wave propagation in two dimensions. Our FEM-BEM algorithm for the full-space model incorporates the RC exactly. Numerical experiments demonstrate the efficiency of the FEM-BEM approach for simulating smooth and non-smooth wave fields, with the latter induced by a complex heterogeneous medium and a discontinuous refractive index.

AMS subject classifications: 65N30, 65N38, 65F10, 35J05

Keywords: Heterogeneous, Unbounded, Wave Propagation, Finite/Boundary Element Methods.

1 Introduction

Wave propagation simulations, governed by the Helmholtz equation, in bounded heterogeneous and unbounded homogenous media are fundamental for numerous applications [13, 33, 39].

Finite element methods (FEM) are efficient for simulating the Helmholtz equation in a bounded heterogeneous medium, say, $\Omega_{0}\subset\mathbb{R}^{m}$ ($m=2,3$). The standard (non-coercive) variational formulation of the variable coefficient Helmholtz equation in $H^{1}(\Omega_{0})$ [33] has been widely used for developing and analyzing the sign-indefinite FEM, see for example [7, 3, 12, 27, 29, 40]. The open problem of developing a coercive variational formulation for the heterogeneous Helmholtz model was solved recently in [28], and an associated preconditioned sign-definite high-order FEM was also established using direct and domain decomposition methods in [28].

For a large class of applications the wave propagation occurs in the bounded heterogeneous medium and also in its complement, $\mathbb{R}^{m}\setminus\Omega_{0}$, the exterior unbounded homogeneous medium. Using the fundamental solution, the constant coefficient Helmholtz equation exterior to $\Omega_{0}$ can be reformulated as an integral equation (IE) on the boundary of $\Omega_{0}$. Algorithms for simulating the boundary IE (BIE) are known as boundary element methods (BEM). Several coercive and non-coercive BIE reformulations [13, 39] of the exterior Helmholtz model have been used to develop algorithms for the exterior homogeneous Helmholtz models, see for example the acoustic BEM survey articles [34, 11], respectively, by mathematical and engineering researchers, each with over 400 references.

The exterior wave propagation BEM models lead to dense complex algebraic systems, and the standard variational formulation based interior wave FEM models lead to sparse complex systems with their eigenvalues in the left half of the complex plane [26, 38]. Developing efficient preconditioned iterative solvers for such systems has also dominated research activities over the last two decades [19], in conjunction with efficient implementations using multigrid and domain decomposition techniques, see [25, 27] and references therein.

For applications that require solving both the interior heterogenous and exterior homogeneous problems, various couplings of the FEM and BEM algorithms with appropriate conditions on polygonal interfaces have also been investigated in the literature[5, 32, 6]. The review article [43] describes some theoretical validations of the coupling approaches considered in the earlier literature and delicate choices of the coupling interface. The coupling methods in [5, 32, 6, 31, 43] lead to very large algebraic systems with both dense and sparse structures. For wave propagation models, given the complexity involved in even separately solving the FEM and BEM algebraic systems, it is efficient to avoid large combined dense and sparse structured systems arising from the coupling methods in [5, 32, 6, 31, 43].

Such complicated-structured coupled large-scale systems can be avoided, for the Helmholtz PDE interior and exterior problems, using the approach proposed in [35] and recently further explored in [24] using high-order elements for a class of applications with complex heterogeneous structures. The FEM-BEM algorithms in [35, 24] are based on the idea of using a non-overlapping smooth interface to couple the interior and exterior solutions. As described in [24, Section 6], there are several open mathematical analysis problems remain to be solved in the coupling and FEM-BEM framework of [35, 24].

The choice of smooth interface in the FEM-BEM algorithms of [35, 24] is crucial because the methods require solving several interior and exterior wave problems to setup the interface condition. In particular, the number of FEM and BEM problems to be solved is twice the number of degrees of freedom required to approximate the unknown interface function. The interface function can be approximated by a few degrees of freedom only on smooth interfaces. Efficient spectrally accurate BEM algorithms have been developed for simulating scattered waves exterior to smooth boundaries in two and three dimensional domains [8, 9, 13, 20]. However for standard interior FEM algorithms, it is desirable to have simple polygonal/polyhedral boundaries, and in particular those with right angles, which facilitate the development and implementation of high-order FEM algorithms.

To this end, we develop an equivalent framework for the heterogeneous and unbounded region wave propagation model with two artificial interfaces. In particular, our novel FEM-BEM framework is based on an interior smooth interface $\Gamma$ for simulating scattered exterior waves using a spectrally accurate Nyström BEM, and an exterior simple polygonal/polyhedral interface $\Sigma$ for the efficient high-order FEM simulation of the absorbed interior waves. In Figure 1, we sketch the resulting overlapped decomposition of a heterogeneous and unbounded medium in which the absorbed and scattered waves are induced by an input incident wave $u^{\rm inc}$.

The decomposition facilitates the application of efficient high-order FEM algorithms in the interior polygonal/polyhedral domain $\Omega_{2}$, that contains the heterogeneous region $\Omega_{0}\subset\overline{\Omega}_{1}$. The unbounded exterior region $\mathbb{R}^{m}\setminus\overline{\Omega}_{1}$ does not include the heterogeneity and has a smooth boundary $\Gamma$. It therefore supports spectrally accurate BEM algorithms to simulate exterior scattered waves, and also exactly preserves the radiation condition (RC), even in the computational model.

In addition, the decomposition framework provides an analytical integral representation of the far-field using the scattered field, and hence our high-order FEM-BEM model provides relatively accurate approximations of the far-field arising from the heterogeneous model. For inverse wave models, accurate modeling of the far-field plays a crucial role in the identification of unknown wave propagation configuration properties from far-field measurements [2, 13, 23].

Our approach in this article is related to some ideas presented in [10, 16, 15]. The choice of two artificial boundaries leads to two bounded domains $\overline{\Omega}_{0}\subset\Omega_{1}\subset\Omega_{2}$ and an overlapping region between $\Omega_{1}^{\rm c}=\mathbb{R}^{m}\setminus\overline{\Omega}_{1}$ and $\Omega_{2}$. We prove that, under appropriate restrictions of the scattered and absorbed fields in the overlapping region $\Omega_{12}(:=\Omega_{1}^{\rm c}\cap\Omega_{2})$, our decomposed model is equivalent to the original Helmholtz model in the full space $\mathbb{R}^{m}$. The unknowns in our decomposed framework, which exactly incorporates the RC, are: (a) the trace of the scattered wave on $\Gamma$ that will yield the solution in the unbounded domain $\Omega_{1}^{\rm c}$, through a boundary layer potential ansatz of the scattered field; (b) the trace of the total wave in the boundary $\Sigma$ of $\Omega_{2}$, that will provide the Dirichlet data to determine the total absorbed wave in the bounded domain $\Omega_{2}$. These properties will play a crucial role in designing and implementing our high-order FEM-BEM algorithm.

The FEM-BEM numerical algorithm can be discerned at this point: It comprises approximating the absorbed wave field in a finite dimensional space using an FEM spline ansatz in the bounded domain $\Omega_{2}$, and by a BEM ansatz for the scattered field in the unbounded region, exterior to $\Gamma$, and these fields are constrained to (numerically) coincide on the overlapping domain $\Omega_{12}$, and hence on the interface boundaries. Since these artificial boundaries can be freely chosen, we can ensure a bounded simple polygonal/polyhedral domain, more suitable for high-order FEM, and an unbounded region with a smooth boundary for spectrally accurate BEM. In particular, the framework brings the best of the two numerical (FEM and BEM) worlds to compute the fields accurately for the full heterogeneous model problem, without the need to truncate the unbounded wave propagation region and approximate the RC.

The algorithmic construction and solving of the interface linear system, which determines key unknowns of the model on the interface boundaries (that is, the ansatz coefficients of the trace of the FEM and BEM solutions), is challenging. However, important properties of the continuous problem, such as a compact perturbation of the identity, are inherited by the numerical scheme. Consequently, the system of linear equations for the interface unknowns is very well conditioned. Such properties, in conjunction with a cheaper matrix-vector multiplication for the underlying matrix, support the use of iterative solvers such as GMRES [42, 41] to compute the ansatz coefficients. Major computational aspects of our high-order FEM and BEM discretizations in the framework are independent and hence the underlying linear systems can be solved, a priori, by iterative Krylov methods. We show that the number of GMRES iterations, to solve the interface system, is independent of various levels of discretization for a chosen frequency of the model. For increasing frequencies, we also demonstrate that the growth of the number of GMRES iterations is lower than the frequency growth.

Instead of using an iterative scheme for the interface system arising in our algorithm, one may also consider the construction and storage of the matrix and a direct solver for the system. The advantage of the latter is that the interface problem matrix can be reused for numerous incident input waves that occur in many practical applications, for example, to compute the monostatic cross sections, and also for developing appropriate reduced order model (ROM) [21] versions of our algorithm. The matrix arising in our interface system is relatively small because of the spectral accuracy of the BEM algorithm, and because the system involves only unknowns on the artificial interface boundaries. Hence post-processing of the computed fields, such as for the evaluation of the far-field, can be done quickly and efficiently. The far-field output also plays a crucial role in developing stable ROMs for wave propagation models [22, 21].

The paper is organized as follows. In Section 2 we present the decomposition framework and prove that, under very weak assumptions, the decomposition is well-posed and is equivalent to the full heterogeneous and unbounded medium wave propagation model. In Section 3 we present a numerical discretization for the two dimensional case, combining high-order finite elements with spectrally accurate convergent boundary elements [36] and describe the algebraic and implementation details. In Section 4 we demonstrate the efficiency of the FEM-BEM algorithm for simulating wave propagation in two distinct classes of (smooth and non-smooth) heterogenous media.

2 Decomposition framework and well-posedness analysis

Let $\Omega_{0}\subset\mathbb{R}^{m}$, $m=2,3$, be a bounded domain. The ratio of the speed of wave propagation inside the heterogeneous (and not necessarily connected) region $\Omega_{0}$ and on its free-space exterior $\Omega_{0}^{\rm c}:=\mathbb{R}^{m}\setminus\overline{\Omega}_{0}$ is described through a refractive index function $n$ that we assume in this article to be piecewise smooth with $1-n$ having compact support in $\overline{\Omega}_{0}$ (i.e, $n|_{\Omega_{0}^{\rm c}}\equiv 1$).

The main focus of this article is to study the wave propagation in $\mathbb{R}^{m}$, induced by the impinging of an incident wave $u^{\rm inc}$, say, a plane wave with wavenumber $k>0$. More precisely, the continuous wave propagation model is to find the total field $u(:=u^{s}+u^{\rm inc})\in H^{1}_{\rm loc}(\mathbb{R}^{m})$ that satisfies the Helmholtz equation and the Sommerfeld RC:

[TABLE]

It is well known that (2.1) is uniquely solvable [36]. (Later in this section, we introduce the classical Sobolev spaces $H^{s}$, for $s\geq 0$, with appropriate norms.)

2.1 A decomposition framework

The heterogeneous-homogeneous model problem (2.1) is decomposed by introducing two artificial curves/surfaces $\Gamma$ and $\Sigma$ with interior $\Omega_{1}$ and $\Omega_{2}$ respectively satisfying $\overline{\Omega}_{0}\subset\Omega_{1}\subset\overline{\Omega}_{1}\subset\Omega_{2}$. We assume from now on that $\Gamma$ is smooth and $\Sigma$ is a polygonal/polyhedral boundary. A sketch of the different domains is displayed in Figure 1. Henceforth, $\Omega_{i}^{\rm c}:=\mathbb{R}^{m}\setminus\overline{\Omega}_{i},~{}i=0,1,2$.

We introduce the following decomposed heterogeneous and homogeneous media auxiliary models:

• •

For a given function $f_{\Sigma}^{\rm inp}\in H^{1/2}(\Sigma)$, we seek a propagating wave field $w$ so that $w$ and its trace $\gamma_{\Sigma}w$ on the boundary $\Sigma$ satisfy

[TABLE]

Throughout the article, we assume that this interior problem is uniquely solvable. We introduce the following operator notation for the heterogeneous auxiliary model: For any Lipschitz $m$- or $(m-1)$-dimensional (domain or manifold) $D\subset\Omega_{2}$, we define the solution operator $\mathrm{K}_{D\Sigma}$ associated with the auxiliary model (2.2) as

[TABLE]

Two cases will be of particular interest for us: $\mathrm{K}_{\Omega_{2}\Sigma}{f_{\Sigma}^{\rm inp}}$, which is nothing but $w$ satisfying (2.2), and $\mathrm{K}_{\Gamma\Sigma}f_{\Sigma}^{\rm inp}=\gamma_{\Gamma}w$, the trace of the solution ${w}$ of (2.2) on $\Gamma{\subset\Omega_{2}}$.

• •

In the exterior unbounded homogeneous medium $\Omega_{1}^{\rm c}:=\mathbb{R}^{m}\setminus\Omega_{1}$, for a given function $f_{\Gamma}^{\rm inp}{\in H^{1/2}(\Gamma)}$ we seek a scattered field $\widetilde{\omega}$ satisfying

[TABLE]

Unlike problem (2.2), (2.4) is always uniquely solvable [36]. We define the associated solution operator $\mathrm{K}_{D\Gamma}$ as

[TABLE]

with special attention to $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}f_{\Gamma}^{\rm inp}$ and $\mathrm{K}_{\Sigma\Gamma}f_{\Gamma}^{\rm inp}$, namely the scattered field $\widetilde{\omega}$ satisfying (2.4) and its trace $\gamma_{\Sigma}\widetilde{\omega}$.

The decomposition framework that we propose for the continuous problem is the following:

[TABLE]

Solve the interface boundary integral system to find $(f_{\Sigma},f_{\Gamma})$, using data $(\gamma_{\Sigma}u^{\rm inc},\gamma_{\Gamma}u^{\rm inc})$ :

[TABLE] 2. 2.

Construct the total field for the model problem (2.1) using the solution $(f_{\Sigma},f_{\Gamma})$ of (2.6a), by solving the auxiliary models (2.2) and (2.4):

[TABLE]

We claim that, provided (2.6a) is solvable, the decomposed framework-based field $u$ defined in (2.6b) is the solution of (2.1). Notice that we are implicitly assuming in (2.6b) that

[TABLE]

where we recall the notation $\Omega_{12}=\Omega_{1}^{\rm c}\cap\Omega_{2}$. Indeed, in view of (2.6a), both functions in (2.7) agree on $\Sigma\cup\Gamma$ (the boundary of $\Omega_{12}$). Assuming, as we will do from now on, that the only solution to the homogeneous system

[TABLE]

is the trivial one and noticing that $n|_{\Omega_{12}}\equiv 1$ which implies that $\mathrm{K}_{\Omega_{12}\Sigma}f_{\Sigma}$ and $\mathrm{K}_{\Omega_{12}\Gamma}f_{\Gamma}$ are solutions of the Helmholtz equation in $\Omega_{12}$, we can conclude that (2.7) holds. Since $u$ defined in (2.6b) belongs to $H^{1}_{\rm loc}(\mathbb{R}^{m})$, it is simple to check that this function is the solution of (2.1).

We remark that the hypothesis we have taken on the artificial boundaries/domains, i.e. the well-posedness of problems (2.2) and (2.8), are not very restrictive in practice: $\Sigma$ or $\Gamma$ can be modified if needed. Alternatively, one can consider different boundary conditions on $\Gamma$ and $\Sigma$ (such as Robin conditions), redefining $\mathrm{K}_{D\Sigma}$ and $\mathrm{K}_{D\Gamma}$ accordingly, which will lead to a variant of the framework that we analyze in this article. In a future work we shall explore other boundary conditions on the interfaces and analysis of the resulting variant models.

2.2 Well-posedness of the decomposed continuous problem

The aim of this subsection is to prove that the system of equations (2.6a), under the above stated hypothesis, has a unique solution. Consequently, we can conclude that the decomposition for the exact solution presented in (2.6b) exists and is unique. To this end, we first derive some regularity results related to the operators $\mathrm{K}_{D\Sigma}$ and $\mathrm{K}_{D\Gamma}$ in Sobolev spaces. For the topic of Sobolev spaces, we refer the reader to [1, 37].

2.2.1 Functional spaces

Let $D\subset\mathbb{R}^{m}$ be a Lipschitz domain. For any non-negative integer $s$, we denote

[TABLE]

the Sobolev norm, where the summation uses the standard multi-index notation in $\mathbb{R}^{m}$. For $s=s_{0}+\beta$ with $s_{0}$ a non-negative integer and $\beta\in(0,1)$, we set

[TABLE]

The Sobolev space $H^{s}(\Omega)$ ($s\geq 0$) can be defined as,

[TABLE]

endowed with the above natural norm.

If $\partial D$ denotes the boundary of $D$, we can introduce $H^{s}(\partial D)$ with a similar construction using local charts: Let $\{\partial D^{j},\mu^{j},{\bf x}^{j}\}_{j=1}^{J}$ be an atlas of $\partial D$, that is, $\{\partial D\}_{j}$ is an open covering of $\partial D$, $\{\mu^{j}\}$ a subordinated Lipschitz partition of unity on $\partial D$, and ${\bf x}^{j}:\mathbb{R}^{m-1}\to\partial D$ being Lipschitz and injective with $\partial D^{j}\subset\mathop{\rm Im}{\bf x}^{j}$, then we define

[TABLE]

We note that $(\mu^{j}\varphi)\circ{\bf x}^{j}$ can be extended by zero outside of the image of ${\bf x}^{j}$. We then set

[TABLE]

The space $H^{s}(\partial D)$ is well defined for $s\in[0,1]$: Any choice of $\{\partial D^{j},\mu^{j},{\bf x}^{j}\}$ gives rise to an equivalent norm (and inner product). If $\partial D$ is a ${\cal C}^{m}$-boundary, such as $\Gamma$ in Figure 1, this construction can be set up for $s\in[0,m]$ by taking $\{{\bf x}^{j},{\omega}^{j}\}$ to be in ${\cal C}^{m}$ as well. In particular, if $\partial D$ is smooth we can define $H^{s}(\partial D)$ for any $s\geq 0$. Further, the space $H^{-s}(\partial D)$ can be defined as the realization of the dual space of $H^{s}(\partial D)$ when the integral product is taken as a representation of the duality pairing.

It is a classical result that the trace operator $\gamma_{\partial D}u:=u|_{\partial D}$ defines a continuous onto mapping from $H^{s+1/2}(D)$ into $H^{s}(\partial D)$ for any $s\in(0,1)$. Actually, if $\partial D$ is smooth then $s\in(0,\infty)$. In these cases, we can alternatively define

[TABLE]

endowed with the image norm:

[TABLE]

We will use this definition to extend $H^{s}(\partial D)$ for $s>1$ in the Lipschitz case. Notice that with this definition, the trace operator from $H^{s+1/2}(D)$ into $H^{s}(\partial D)$ is continuous for any $s>0$.

2.2.2 Boundary potentials and integral operators

Let $\Phi_{k}$ be the fundamental solution for the two- or three-dimensional constant coefficient Helmholtz operator ($\Delta+k^{2}I$) equation, defined for $\mathbf{x},\mathbf{y}\in\mathbb{R}^{m}$ with $r:=|\mathbf{x}-\mathbf{y}|$ as

[TABLE]

where $H^{(1)}_{n}$ denotes the first kind Hankel function of order $n$. For a smooth curve/surface $\Gamma$, with outward unit normal ${\bm{\nu}}$ and normal derivative at ${\bf y}\in\Gamma$ denoted by $\partial_{{\bm{\nu}}({\bf y})}$, let

[TABLE]

denote the single- and double-layer potentials, with density functions $\varphi$ and $g$, respectively.

The single- and double-layer boundary integral operators are then given, via the well-known jump relations [13] for the boundary layer potentials, by

[TABLE]

where $\gamma^{-}_{\Gamma}$ and $\gamma^{+}_{\Gamma}$ are trace operators on $\Gamma$, respectively, from the interior $\Omega_{1}$ and exterior $\Omega_{1}^{\rm c}$. Given a real non-vanishing smooth function $\sigma:\Gamma\to\mathbb{R}$, and $\mathrm{V}_{k,\sigma}\phi:=\mathrm{V}_{k}(\sigma\phi)$ for any $\phi\in H^{s}(\Gamma)$, we consider the combined field acoustic layer operator

[TABLE]

Throughout this article, ${\rm I}$ denotes the identity operator. The standard combined field operator used in the literature [13] is based on the choice $\sigma\equiv 1$. In this article, we do not restrict ourselves to the usual choice for reasons which will be fully explained later. Since $\Gamma$ is smooth and that $\mathrm{K}_{k},\mathrm{V}_{k,\sigma}:H^{s}(\Gamma)\to H^{s+1}(\Gamma)$ are continuous, the operator in (2.13) is invertible as a consequence of the Fredholm alternative and the injectivity of (2.13), which follows from a very simple modification of the classical argument in [13, Th 3.33]).

Thus the inverse of the combined field integral operator

[TABLE]

is well defined. Further, using (2.11)-(2.12) and with $\mathrm{SL}_{k,\sigma}\phi:=\mathrm{SL}_{k}(\sigma\phi)$ for any $\phi\in H^{s}(\Gamma)$, we can write the solution operator occurring in the construction (2.6b) as

[TABLE]

The above solution operator, a variant of the Brakhage-Werner formulation (BWF) [4, 13], will be used in this article for both theoretical and computational purposes. The choice $\sigma\equiv 1$ reduces to the standard BWF [4, 13].

2.2.3 Well-posedness analysis of the interface model

In this subsection, we first develop two key results before proving well-posedness of the boundary integral system (2.6a).

Lemma 2.1**.**

The operator

[TABLE]

is continuous for any $s\in[0,\infty)$. Further, for any bounded Lipschitz domain/manifold $D\subset\Omega_{1}^{\rm c}$ with $\overline{D}\cap\Gamma=\emptyset$, the solution operator $\mathrm{K}_{D\Gamma}$ in  (2.5), for the homogeneous media problem (2.4), satisfies the following mapping property for any $s,r\geq 0$

[TABLE]

In particular,

[TABLE]

is continuous and compact, for $s,r\in\mathbb{R}$.

Proof.

The first desired property follows from the identities (2.15), (2.14) and the well known mapping properties

[TABLE]

see for instance [37, Th. 6.12]. If $\overline{D}\cap\Gamma=\emptyset$, the kernels in the boundary potentials in $\mathrm{DL}_{k}$ and $\mathrm{SL}_{k}$ are smooth functions in $\overline{D}\times\Gamma$ and hence the properties (2.17) and (2.18) hold. ∎

Next we consider the heterogeneous media model solution operator $\mathrm{K}_{\Omega_{2}\Sigma}$, as defined in (2.2)- (2.3). We recall the well known classical estimate [33]

[TABLE]

with $C>0$ being a constant independent of ${f_{\Sigma}^{\rm inp}}$. Below, we generalize this to obtain a higher regularity, using boundary layer potentials and boundary integral operators, defined in this case on barely Lipschitz curves/surfaces to improve the estimate for domains $D$ with $\overline{D}\subset\Omega_{2}\setminus\overline{\Omega}_{1}$.

Lemma 2.2**.**

There exists a constant $C=C(k,n,\Omega_{2})$ so that for any $s\in[0,1]$ and $f_{\Sigma}^{\rm inp}\in{H^{s}(\Sigma)}$,

[TABLE]

Furthermore, if $D\subset\overline{D}\subset\Omega_{2}\setminus\overline{\Omega}_{1}$ the following solution operator mapping property holds for any $r\in\mathbb{R}$

[TABLE]

Consequently,

[TABLE]

is continuous and compact, for any $r\in\mathbb{R}$.

Proof.

Throughout this proof we let $s\in[0,1]$ and, for notational convenience, we denote $v:=\mathrm{K}_{\Omega_{2}\Sigma}f_{\Sigma}^{\rm inp}$. Since, by definition,

[TABLE]

By the third Green identity (see for instance [37, Th. 6.10]) we have the representation

[TABLE]

with $\mathop{\rm supp}{g^{v}_{n}}\subset\Omega_{0}$, where we have used the notation

[TABLE]

In the expression above $\mathrm{SL}_{k,\Sigma}$ and $\mathrm{DL}_{k,\Sigma}$ denote respectively the single- and double-layer potential from the corresponding densities, defined on $\Sigma$, associated with the constant coefficient Helmholtz operator $\Delta+k^{2}{\rm I}$. Next we prove that

[TABLE]

To this end, we start from the decomposition $v=v_{1}+v_{2}$, where the harmonic $v_{1}$ and the interior wave-field $v_{2}$ are solutions of

[TABLE]

Classical potential theory results, see [37, Th 6.12] and the discussion which follows it (see also references therein), show that there exists $C>0$ so that

[TABLE]

for any $f_{\Sigma}^{\rm inp}\in H^{s}(\Sigma)$. On the other hand, following  [30, Ch. 4] or [14] there exists $\varepsilon>0$ and $C_{\varepsilon}>0$ such that

[TABLE]

By the trace theorem (applied to $\nabla v_{2}$),

[TABLE]

Combining these estimates with (2.23) we conclude that

[TABLE]

Notice also that if $D\subset\overline{D}\subset\Omega_{2}\setminus\overline{\Omega}_{1}$, because the kernels of the potentials operators and the Newton potential are smooth in the corresponding variables, we gain from the extra smoothing properties of the underlying operators in (2.23) to derive

[TABLE]

where the constants $C$ and $C^{\prime}$ are independent of $f_{\Sigma}^{\rm inp}$.

∎

For deriving the main desired result of this section, it is convenient to define the following off-diagonal operator matrix

[TABLE]

Then (2.6a) can be written in operator form

[TABLE]

where ${\cal I}$ denotes the $2\times 2$ block identity operator. A simple consequence of Lemmas 2.1 and 2.2 is that

[TABLE]

is continuous for any $s\geq 0$. Next we prove that this operator is indeed an isomorphism:

Theorem 2.3**.**

For any $s\geq 0$,

[TABLE]

is an invertible compact perturbation of the identity operator.

Proof.

The continuity of ${\cal K}:H^{0}(\Sigma)\times H^{0}(\Gamma)\to H^{s}(\Sigma)\times H^{s}(\Gamma)$ for any $s\geq 0$ has already been established in the two preceding lemmas. In particular, ${\cal K}$ is compact. Moreover, the null space ${\cal I}-{\cal K}$ consists of smooth functions. For any $(g_{\Sigma},g_{\Gamma})\in N({\cal I}-{\cal K})$, we construct

[TABLE]

Note that $w:=(v-\vartheta)$ defined, in principle, in $\Omega_{12}={\Omega_{2}\cap\Omega_{1}^{\rm c}}$ satisfies

[TABLE]

By the well-posedness of problem (2.8), we have $w=0$ in $\Omega_{12}$. We define $u$ on $\mathbb{R}^{m}$ as

[TABLE]

Note that $u$ is well defined in $\Omega_{12}$, and it is a solution of (2.1) with incident wave $u^{\rm inc}=0$. Therefore, $u=0$ which implies that $\vartheta=0$ in $\Omega_{2}^{\rm c}$. The principle of analytic continuation yields that $\vartheta=0$ also in $\Omega_{1}^{\rm c}$ and therefore $g_{\Gamma}=\gamma_{\Gamma}\vartheta=0$. Finally,

[TABLE]

and hence the desired result follows. ∎

3 A FEM-BEM algorithm for the decomposed model

In this section we consider the numerical discretizations on the proven equivalent decomposed system (2.6). In this article, we restrict to the two-dimensional (2-D) case. [The 3-D algorithms and analysis for (2.6) will be different to the 2-D case, and in a future work we shall investigate a 3-D FEM-BEM computational model.] Briefly, the approach consists of replacing the continuous operators $\mathrm{K}_{\Omega_{2}\Sigma}$ and $\mathrm{K}_{D\Gamma}$ with suitable high-order FEM and BEM procedures-based discrete operators. The stability of such a discretization depends on the numerical methods chosen in each case.

For discretization of the differential operator $\mathrm{K}_{\Omega_{2}\Sigma}$ based on the heterogeneous domain model, we could consider a standard FEM with triangular, quadrilateral or even more complex elements. We will choose the first case, for the sake of simplicity, and we expect the analysis developed in this case could cover these other types of elements, with appropriate minor modifications.

The BEM procedure, for discretizing the exterior homogeneous medium associated $\mathrm{K}_{D\Gamma}$ through boundary integral operators, is more open since an extensive range of methods is available in the literature. We will restrict ourselves to the spectral Nyström method [36] (see also [17]). This scheme provides a discretization of the four integral operators of the associated Calderon calculus, and has exponential rate of convergence. In this article, we will make use of high-order discretizations of the single- and double-layer operators that are easy to implement.

A key restriction of the standard Nyström method to achieve spectrally accurate convergence is the requirement of a smooth diffeomorphic parameterization of the boundary. This is because the method starts from appropriate decompositions and factorizations of the kernels of the operators to split these into regular and singular parts. This is not a severe restriction in our case since $\Gamma$ is an auxiliary user-chosen smooth curve and can therefore be easily constructed as detailed as required.

Next we briefly consider these two known numerical procedures and hence describe our combined FEM-BEM algorithm and implementation details.

3.1 The FEM procedure

Let $\{\mathcal{T}_{h}\}_{h}$ be a sequence of regular triangular meshes where $h$ is the discrete mesh parameter, the diameter of the largest element of the grid. Hence we write $h\to 0$ to mean that the maximum of the diameters of the elements tends to 0. Using ${\cal T}_{h}$, we construct the finite dimensional spline approximation space

[TABLE]

where $\mathbb{P}_{d}$ is the space of bivariate polynomials of degree $d$. We define the FEM approximation $\mathrm{K}_{\Omega_{2}\Sigma}^{h}$ to $\mathrm{K}_{\Omega_{2}\Sigma}$ as follows: The FEM operator

[TABLE]

for $f_{\Sigma,h}^{\rm inp}\in\gamma_{\Sigma}\mathbb{P}_{h,d}$, is constructed as $u_{h}:=\mathrm{K}_{\Omega_{2}\Sigma}^{h}f_{\Sigma,h}^{\rm inp}$, where $u_{h}\in\mathbb{P}_{h,d}$ is the solution of the discrete FEM equations:

[TABLE]

The discrete FEM operator $\mathrm{K}_{\Omega_{2}\Sigma}^{h}$ is well defined for sufficiently small $h$.

3.2 The BEM procedure

Let

[TABLE]

be a smooth $2\pi-$periodic regular parameterization of $\Gamma$. We denote by the same symbol ${\rm SL}_{k}$, ${\rm DL}_{k}$, ${\rm V}_{k}$ and ${\rm K}_{k}$ the parameterized layer potentials and boundary layer operators:

[TABLE]

where $\bm{\mu}(t):=(x_{2}^{\prime}(t),-x_{1}^{\prime}(t))=|{\bf x}^{\prime}(t)|\>{\bm{\nu}}\circ{\bf x}(t)$. Observe that $|{\bf x}^{\prime}(t)|$ is incorporated into the density in ${\rm SL}_{k}$ and to the kernel in ${\rm DL}_{k}$. We follow the same convention for the single- and double-layer weakly singular boundary integral operators. For high-order approximations, it is important to efficiently take care of the singularities. In particular, for the spectrally accurate Nyström BEM solver, we use the following representations of the layer operators with smooth $2\pi$ bi-periodic kernels $A,\ B,\ C,\ D$ [13]:

[TABLE]

The Nyström method, based on a discrete positive integer parameter $N$, starts with setting up a uniform grid

[TABLE]

and the space of trigonometric polynomials of degree at most $N$

[TABLE]

with $\mathbb{Z}_{N}:=\{-N+1,-N+2,\ldots,N\}.$ We next introduce the interpolation operator $\mathrm{Q}_{N}$

[TABLE]

to define discretizations of the single and double layer operators:

[TABLE]

We stress that the above integrals can be computed exactly using the identities:

[TABLE]

and for $g_{N}\in\mathbb{T}_{N}$,

[TABLE]

which are based on properties of the trapezoidal/rectangular rule for $2\pi-$periodic functions.

The high-order approximation evaluation of the potentials is achieved in a similar way:

[TABLE]

leading to the rectangular rule approximation as in (3.6)

Now we are ready to describe the discrete operator $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}^{N}$ that is a high-order approximation to the exterior homogeneous model continuous operator $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}$. First, we introduce the parameterized counterpart of the continuous operator in (2.13),

[TABLE]

(which corresponds to $\sigma\circ{\bf x}=\frac{1}{|{\bf x}|}$). Then we define

[TABLE]

We remark that the definition of $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}^{N}$ requires only evaluation of input functions at the grid points. In particular, it is well defined on continuous functions. Indeed, we have

[TABLE]

and since the discrete boundary layer operators only use pointwise values of the density at the grid points (i.e., ${\rm Q}_{N}\varphi$), evaluation of $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}^{N}g$ requires only values of $g$ at the grid points. So we can replace, when necessary,

[TABLE]

The discrete operator $\mathrm{K}_{\Sigma\Gamma}^{N}g$ is defined accordingly by taking the trace of $\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}^{N}g$ on $\Sigma$. Thus our algorithm is based on the idea of taking the trace of FEM and BEM solutions on $\Gamma$ and $\Sigma$ respectively.

3.3 The FEM-BEM computational model

In addition to the discrete operators defined above, we need one last discrete operator to describe the FEM-BEM algorithm. Let

[TABLE]

denote the usual Lagrange interpolation operator on $\gamma_{\Sigma}\mathbb{P}_{h,d}$, the inherited finite element space on $\Sigma$. Our full FEM-BEM algorithm is:

[TABLE]

• •

Step 1: Solve the finite dimensional system

[TABLE]

• •

Step 2: Construct the FEM-BEM solution

[TABLE]

Remark 3.1*.*

We have committed a slight abuse of notation in the right-hand-side of (3.12a) by writing

[TABLE]

instead of the correct, but more complex, {\rm Q}_{N}\big{(}(\gamma_{\Gamma}u^{\rm inc})\circ{\bf x}\big{)}. Similarly,

[TABLE]

should be read in the lower extra-diagonal block of the matrix in (3.12a). Indeed, this is equivalent to replacing a space on $\Gamma$ with that obtained via the parameterization (3.2). Since both spaces are isomorphic, being strict in the notation for description of these operators is not absolutely necessary. In particular, we avoid complicated notation and use a compact way to describe the algorithm and associated theoretical results.

Remark 3.2*.*

Complete numerical analysis of the FEM-BEM algorithm is beyond the scope of this article. In a future work, we shall carry out a detailed numerical analysis of the FEM-BEM algorithm. Below we give the main results. In summary, the analysis is based on the following assumption on the mesh-grid:

Assumption 1

There exists $\varepsilon_{0}>0$ such that the sequence of grids $\{{\cal T}_{h}\}_{h}$ satisfies

[TABLE]

where $D\subset\Omega_{2}\setminus\overline{\Omega}_{0}$ is an open neighborhood of $\Gamma$, and $h_{D}$, the maximum of the diameters of the elements of the grid ${\cal T}_{h}$ with non-empty intersection with $D$.

We note that this assumption allows locally refined grids, but introduces a very weak restriction on the ratio between the largest element in $\Omega_{2}$ and the smallest element in $D$. However, since the exact solution is smooth on $D$, the partial differential equation in this domain is just the homogeneous Helmholtz equation, and it is reasonable to expect that small elements are not going to be used in this subdomain.

Using Assumption 1, in a future work we shall prove the well-posedness of the discrete system (• ‣ 3.3) and optimal order of convergence of the FEM-BEM solution. In particular, after deriving convergence of the individual FEM and BEM approximations, we shall prove the following convergence result: For any region $\Omega_{R}\subset\Omega_{1}^{\rm c}=\mathbb{R}^{2}\setminus\overline{\Omega}_{1}$, $0<\epsilon\leq\epsilon_{0}$, $r\geq 0$, $t\geq d+3/2$,

[TABLE]

where $h_{D}$ is as in (3.13) and $h_{\Sigma}$ is the maximum distance between any two consecutive Dirichlet/constrained nodes in $\mathcal{T}_{h}$; $(u,\omega)=(\mathrm{K}_{\Omega_{2}\Sigma}f_{\Sigma},\mathrm{K}_{\Omega_{1}^{\rm c}\Gamma}f_{\Gamma})$ is the exact solution of (2.6); and $(u_{h},\omega_{N})$ is the unique solution of the numerical method (• ‣ 3.3).

Next we describe algebraic details required for implementation of the algorithm, followed by numerical experiments in Section 4 to demonstrate the efficiency of the FEM-BEM algorithm to simulate wave propagation in the heterogeneous and unbounded medium.

3.4 FEM-BEM algebraic systems and evaluation of wave fields

Simulation of approximate interior and exterior wave fields $u_{h,N}$ using the solution of (3.12a) and the representation in (3.12b) requires: (i) computing the interior solution $u_{h}$ by once solving the finite element system (3.1) using the Dirichlet data $f_{\Sigma}^{h}$; and (ii) the exterior solution ${\omega_{N}}$ in $\Omega_{1}^{\rm c}$ by evaluating the layer potential value $(\mathrm{DL}_{k}^{N}-\mathrm{i}\mathrm{SL}_{k}^{N}){\cal L}_{k}^{N}f^{N}_{\Gamma}$, using the representation in (3.7).

Since ${\cal L}_{k}^{N}f^{N}_{\Gamma}\in\mathbb{T}_{N}$ and that the dimension of $\mathbb{T}_{N}$ is $2N$, using (3.4)–(3.7), the degrees of freedom (DoF) required to compute the exterior solution ${\omega_{N}}$ is equal to the number of interpolatory uniform grid points $t_{j},~{}j=-N+1,\dots,N$ in (3.3) that determine the interpolatory operator $\mathrm{Q}_{N}$ in (3.5). The linear algebraic system corresponding to the Dirichlet problem (3.1) for $u_{h}\in\mathbb{P}_{h,d}$ is obtained by using an ansatz that is a linear combination of the basis functions spanning $\mathbb{P}_{h,d}$. Coefficients in the $u_{h}$ ansatz are values of $u_{h}$ at the nodes that determine $\{\mathcal{T}_{h}\}_{h}$. The nodes include constrained/boundary Dirichlet nodes on $\Sigma$ and free/interior non-Dirichlet nodes in $\Omega_{2}$.

Henceforth, for a chosen mesh for the bounded domain $\Omega_{2}$, we use the notation $M$ and $L$ to denote the number of Dirichlet- and free-nodes nodes in the mesh, respectively. The FEM system (3.1) to compute the solution $u_{h}$ leads to an $L$-dimensional linear system for the unknown vector ${\bf u}_{L}$ (that are values of ${u}_{h}$ at the interior nodes). The system is governed by a real symmetric sparse matrix, say, ${\bf A}_{L}$. The matrix ${\bf A}_{L}$ is obtained by eliminating the row and column vectors associated at the boundary nodes. Let ${\bf D}_{L,M}$ be the $L\times M$ matrix that is used to move the Dirichlet condition to the right-hand-side of the system. Thus for a given Dirichlet data vector $\widehat{{\bf f}}_{M}$, we may theoretically write ${\bf u}_{L}={\bf A}^{-1}_{L}{\bf D}_{L,M}\widehat{{\bf f}}_{M}$. Let ${\bf T}_{2N,L}$ be the $2N\times L$ sparse matrix so that ${\bf T}_{2N,L}{\bf u}_{L}(={\bf T}_{2N,L}{\bf A}^{-1}_{L}{\bf D}_{L,M}\widehat{{\bf f}}_{M})$ is the trace of the finite element solution ${w}_{h}$ of (3.1) at the $2N$ interior points $\mathbf{x}(t_{j})\in\Gamma,j=-N+1,\dots,N$ that are the BEM grid points.

For describing the full FEM-BEM system, using the above representation, it is convenient to define the $2N\times M$ matrix

[TABLE]

The matrix ${\bf A}^{-1}_{L}$ in (3.15), in general, should not be computed in practice. We may consider instead a ${\bf L}_{L}{\bf D}_{L}{\bf L}_{L}^{\top}$ factorization [18] (for example, implemented in the Matlab command ldl), where ${\bf D}_{L}$ is a block diagonal matrix with $1\times 1$ or $2\times 2$ blocks and ${\bf L}_{L}$ is a block (compatible) unit lower triangular matrix. Hence, each multiplication by ${\bf A}^{-1}_{L}$ is reduced to solving two (block) triangular and one $2\times 2$ block diagonal system which can be efficiently done, leading to evaluation of ${\bf\widetilde{K}}_{2N,M}$ on $M$-dimensional vectors. Of course the ${\bf L}_{L}{\bf D}_{L}{\bf L}_{L}^{\top}$ factorization is a relatively expensive process, but worthwhile in our method to simulate the complex heterogeneous and unbounded region model. (We further quantify this process using numerical experiments in Section 4.)

The ansatz for the unknown density $f^{N}_{\Gamma}\in\mathbb{T}_{N}$ is a linear combination of $2N$ known basis functions ${\exp({\rm i}\ell t)},~{}\ell=-N+1,\dots,N$ in (3.4) that span $\mathbb{T}_{N}$. The $2N$-dimensional BEM system for the unknown vector $\widetilde{\bf f}_{2N}$ (that are values of the unknown density at the Nyström node points $t_{j},~{}j=-N+1,\dots,N$) is governed by a complex dense matrix and an input $2N$-dimensional vector $\widetilde{{\bf f}}_{2N}$ determined by the Dirichlet data on $\Gamma$ in the exterior homogeneous model (2.4) evaluated at $t_{j},~{}j=-N+1,\dots,N$. We may write

[TABLE]

where ${\bf B}_{2N}$ is the $2N\times 2N$ Nyström matrix corresponding to the discrete boundary integral operator in (3.9). Similar to ${\bf T}_{2N,L}$, let ${\bf P}_{M,2N}$ be the matrix representation of the (discrete) combined potential generated by a density at the $M$ Dirichlet nodes of ${\cal T}_{h}$. That is, ${\bf P}_{M,2N}\bm{\varphi}_{2N}$ is the vector form of ${\rm Q}_{h}^{\Sigma}\gamma_{\Sigma}(\mathrm{DL}_{k}^{N}-{\rm i}k\mathrm{SL}_{k}^{N})\varphi$, following the BEM representation (3.9) for evaluation of the exterior field at the $M$ Dirichlet nodes on $\Sigma$. Similar to the interior problem based matrix in (3.15), corresponding to the exterior field it is convenient to introduce the $M\times 2N$ matrix

[TABLE]

Obviously, $M<<L$ (since $M\sim L^{1/2}$ in the 2D case for quasi-uniform grids) and, thanks to the choice of the smooth boundary $\Gamma$, the standard Nyström BEM is spectrally accurate, which further implies that $2N<<M$. (We will quantify this substantially smaller “$<<$” claim using numerical experiments in Section 4.) Thus the cost of setting up an ${\bf LU}$ decomposition of the dense matrix ${\bf B}_{2N}$ is negligible and consequently the matrix $\widehat{\bf K}_{M,2N}$ product with any $2N$-dimensional vector can be efficiently evaluated.

The implementation procedure described above to compute the interior and exterior fields using (3.12b) requires the $M$-dimensional vector $\widehat{{\bf f}}_{M}$ with the values of the unknown at the Dirichlet nodes on $\Sigma$ and the $2N$-dimensional ${\widetilde{{\bf f}}_{2N}}$ at the $2N$ uniform grid points ${\bf x}(t_{j}),~{}j=-N+1,\dots,N$ on $\Gamma$. Since $\Sigma$ and $\Gamma$ are artificial boundaries for the decomposition of the original model, the vectors $\widehat{{\bf f}}_{M},\widetilde{{\bf f}}_{2N}$ are unknown. The interface system (3.12a), that uses the data $u^{\rm inc}$ in the original model, completes the process to compute $\widehat{{\bf f}}_{M},\widetilde{{\bf f}}_{2N}$. In particular, for the matrix-vector form description of (3.12a), we obtain input data vectors, say $\widehat{{\bf u}}_{M}^{\rm inc}$ and $\widetilde{{\bf u}}_{2N}^{\rm inc}$, using the vector form representations of ${\rm Q}^{h}_{\Sigma}\gamma_{\Sigma}u^{\rm inc}$ and ${\rm Q}_{N}\gamma_{\Gamma}u^{\rm inc}$, respectively.

More precisely, using (3.15)–(3.17), the matrix-vector algebraic system corresponding to  (3.12a) takes the form

[TABLE]

where ${\bf I}_{M},{\bf I}_{2N}$ are, respectively, the $M\times M$ and $2N\times 2N$ identity matrices.

In our implementation, instead of solving the full linear system in (3.18) we work with the Schur complement

[TABLE]

After solving for $\widetilde{\bf f}_{2N}$ in (3.19a), the main computational cost for finding $\widehat{\bf f}_{M}$ involves only the matrix-vector multiplication ${\bf\widehat{K}}_{M,2N}\widetilde{\bf f}_{2N}$. The latter requires solving a BEM system, which can be carried out using a direct solve because $2N$ is relatively small.

4 Numerical experiments

In this section we consider two sets of numerical experiments to demonstrate the overlapping decomposition framework based FEM-BEM algorithm. In the first set of experiments, the heterogeneous domain $\Omega_{0}$ has non-trivial curved boundaries and the refractive index function $n$ is smooth; and in the second set of experiments $\Omega_{0}$ is a complex non-smooth structure and $n$ is a discontinuous function. For these two sets of experiments, we consider the $\mathbb{P}_{d}$ Lagrange finite elements with $d=2,3,4$ for the interior FEM model with mesh values $h$, and several values of the Nyström method parameter $N$ to achieve spectral accuracy and to make the BEM errors less than those in the FEM discretizations. The reported CPU times in the section are based on serial a implementation of the algorithm in Matlab (2017b) on a desktop with a 10-core Xeon E5-2630 processor and $128$GB RAM.

In our numerical experiments to compute $\widetilde{\bf f}_{2N}$ in (3.19a), we solve the linear system using: (i) the iterative GMRES method with the (relative) residual set to $10^{-8}$ in all the cases; and (ii) the direct Gaussian elimination solve which requires the full matrix ${\bf A}_{\rm Sch}$ in (3.19a). Both approaches are compared for the numerical experiments in Section 4.3. As an error indicator of our full FEM-BEM algorithm, we analyze the widely used quantity of interest (QoI) in numerous wave propagation applications: the far-field arising from both the interior and exterior fields induced by the incident field impinging from a particular direction. For a large class of inverse wave models [13], the far-field measured at several directions is fundamental to infer various properties of the wave propagation medium.

To computationally verify the quality of our FEM-BEM algorithm in Section 4, we analyze the numerical far-field error at thousands of direction unit vectors $\bf z$. Using (3.16), we define a spectrally accurate approximation to the QoI as

[TABLE]

The exact representation of the QoI is [13]

[TABLE]

Using the angular representation of the direction vectors ${\bf z}$, we compute approximate far-fields at $1,000$ uniformly distributed angles. We report the QoI errors for various grid parameter sets $(h,N)$, and demonstrate high-order convergence of our FEM-BEM algorithm. The maximum of the estimated errors in the approximate QoI, using the values at the $1,000$ uniform directions, are used below to validate the efficiency and high-order accuracy of the FEM-BEM algorithm.

4.1 Star-shaped domain with five-star-pointed refractive index

In Experiment 1 set, we choose $\Omega_{0}$ to be the star-shaped region sketched in the interior of the disk $\Omega_{1}$ in Figure 2, and the refractive index function is defined using polar coordinates as

[TABLE]

with

[TABLE]

Notice that $\widetilde{\chi}(x)$ is a smooth cut-off function with $\mathop{\rm supp}\chi=(-\infty,1]$. Therefore, the function $\chi$ is smooth and also symmetric around $1/2$: $\chi(1-x)=1-\chi(x)$ for any $x$.

For this example, $\Omega_{2}$ is the rectangle $[-6,6]\times[-8,8]$ with boundary $\Sigma$, so that the diameter of the interior domain is $20$. Thus, for a chosen wavenumber $k$, the interior heterogeneous model is of wavelength $10k/\pi$. For our numerical experiments we choose three wavenumbers $k=\pi/4,\pi,4\pi$, to simulate the problems with acoustic characteristic size of $2.5,10,40$ wavelengths, respectively. The smooth boundary $\Gamma$ for this example is a circle centered at zero and radius $3.5$

For the interior FEM model, the initial coarse grid consists of $2,654$ triangles, which is refined up to four times, in the usual way. We show the simulated far-field error results in Tables 1 and 2 using $\mathbb{P}_{3}$ and $\mathbb{P}_{4}$ elements, respectively. In these tables estimates of the (relative) maximum errors in computing the QoI far-fields are presented as well as the number (given within parentheses) of GMRES iterations needed to achieve convergence with the residual tolerance $10^{-8}$. Next we discuss some key aspects of the computed results in Tables 1-2.

To compute the errors for a set of discretization parameters, as exact/truth solutions we used the FEM-BEM algorithm solutions obtained with $N=640$ and the next level of FEM mesh refinement to these in the tables. The fast spectrally accurate convergence of the Nyström BEM, after achieving a couple of digits of accuracy, can be observed by following the far-field maximum errors in the last columns in Tables 1-2. In particular the last columns results, for the FEM spline degree $d=3,4$ cases, demonstrate that relatively small DoF $2N$ is required for the Nyström BEM solutions accuracy to match that of the FEM solutions, especially compared to the FEM DoF $L$. The last rows in Tables 1-2 clearly demonstrate that higher values of $N$ are not useful because of the stagnation of the errors due to limited accuracy of the FEM discretizations. Further, a closer analysis of the results in Tables 1-2 shows that the the computed far-fields exhibit superconvergence, with ${\cal O}(h^{2d})$ errors. In addition, in Figure 7 we demonstrate the faster convergence of the (Experiment 1) smooth total field solutions in the $H^{1}$- norm, and compare with the rate of convergence for a non-smooth solution (Experiment 2) case.

In the Experiment 1 set, with a smooth heterogeneous region $\Omega_{0}$ and a smooth refractive index function $n$, it can be shown that the exact near-field solution for the model problem is smooth. However, this fact alone is not sufficient to explain in detail the superconvergence of the computed far-fields. We may conjecture that some faster convergence is occurring in the background for the near-field in some weak norms, and that the calculation of the far-fields is benefitting from this to achieve the superconvergence. In a future work, we shall explore the numerical analysis our FEM-BEM algorithm.

In Figure 3, we illustrate the convergence of the GMRES iterations and show that as the frequency is increased four-fold, the number of required iterations for the solutions to converge with the $10^{-8}$ residual tolerance increases at (a slightly) slower rate.

Next we consider how the size of the overlapped FEM-BEM region $\Omega_{12}$ affects the speed of convergence of the GMRES iterations. To this end, we have run a set of additional experiments for the star-shaped (Experiment 1) problem with $k=\pi$, using several choices of $\Gamma$, to obtain larger to smaller diameter overlapped regions $\Omega_{12}$. In particular, we chose several BEM smooth boundaries $\Gamma$ to be circles centered at the origin with radii spanning from 2.625 (closer to the heterogeneity) to 5.856 (closer to the FEM boundary $\Sigma$), yielding several $\Omega_{12}$, respectively, with larger to smaller sizes. For all these simulation cases, we fixed the BEM DoF to be $2N=160$, and the fixed $\mathbb{P}_{3}$ elements were obtained using $445,440$ triangles with the number of free-nodes (FEM DoF) to be $L=1,106,385$. We present the corresponding results in Figure 4.

In the left panel of Figure 4, we can see a sample of the curves $\Gamma$ used for the set of experiments with varying size $\Omega_{12}$, and correspondingly in the right panel of Figure 4, we present the number of GMRES iterations required to converge with, again, the residual tolerance $10^{-8}$. Results in Figure 4 clearly demonstrate that the number of GMRES iterations increases as the size of the overlapped region $\Omega_{12}$ decreases. This can be explained as follows: At the continuous level, the interacting operators $\mathrm{K}_{\Sigma\Gamma}$ and $\mathrm{K}_{\Gamma\Sigma}$ tend to lose the compactness property, as the overlapped region becomes thinner. (We shall explore this observation theoretically in a future work.). On the other hand, it is interesting to note from these experiments that the choice of $\Gamma$ being very close to the heterogeneity does not affect the convergence of the GMRES iterations. We could conjecture that this might happen for the considered set of experiments because the exact solution for Experiment 1 problem is smooth. However, we have noticed a similar behavior for the next Experiment 2 problem, with a complex non-smooth heterogeneous region, for which regularity of the total wave field is limited.

4.2 Pikachu-shaped domain with piecewise smooth refractive index

In Experiment 2 set of experiments, we consider a more complicated non-smooth heterogeneous region shown in the interior of the curved domain $\Omega_{1}$ in Figure 5. The region $\Omega_{0}$ is set to be a polygonal Pikachu-shaped domain with the discontinuous refractive index function

[TABLE]

where $r=\sqrt{(x+0.18)^{2}+(y+0.6)^{2}},\theta=\arctan\!2((y+0.6),(x+0.18))$. The grids used in our computation, are adapted to the region $\Omega_{0}$, in such a way that any triangle $\tau\in{\cal T}_{h}$ is either contained or has empty intersection with $\Omega_{0}$. As the boundary of $\Omega_{1}$ and for the smooth curve $\Gamma$ for the exterior model, we choose

[TABLE]

For the interior FEM model, we choose $\Omega_{2}$ to be a polygonal domain as in Figure 5 with boundary $\Sigma$. We then proceed as in the previous experiment, using an initial coarse grid with $8,634$ triangles which is refined up to four times. The solution $u$ of the model is not smooth in $\Omega_{0}$ and $\overline{\Omega}_{0}^{c}$, because of the non-smoothness of the region $\Omega_{0}$ and the jump in the refractive index function. One may consider the use of a graded mesh around the boundary of $\Omega_{0}$ to obtain faster convergence. Based on the size of $\Omega_{2}$, the choices $k=\pi/4,\pi,4\pi$ lead to approximately $2.5$, $10$, and $40$ wavelengths interior FEM model, respectively, for simulations in Experiment 2.

We observe from the integer numbers (within in parentheses) in Tables 3-4 that the number of GMRES iterations grow, slower than the quadruple growth of the three frequencies considered in Experiment 2. The estimated (relative) maximum far-field errors for the non-smooth Experiment 2 model are given in Tables 3-4, demonstrating high-order accuracy of our FEM-BEM model as the finite element space degree, grid size, and the BEM DoF are increased. In Figure 7, for $d=2,3,4$, we compare convergence of the total field in the $H^{1}$-norm for the smooth (Experiment 1) and non-smooth (Experiment 2) simulations.

In Figure 7 we depict the simulated wave field solution for $k=\pi$, with $\mathbb{P}_{4}$ finite elements on a grid with $138,144$ triangles and $L=1,106,385$ free-nodes for the FEM solution, and $2N=320$ for the BEM solution. Specifically, we plot the simulated absorbed and scattered field numerical solution $u_{h,N}$ inside $\Omega_{2}$ in Figure 7.

4.3 Direct solver implementation and comparison with iterative solver

In this subsection we discuss the direct solver implementation of our method and compare its performance with the iterative approach we have used for simulating results described earlier in the section. When computing the matrix in (3.19a), the main issue is concerned with the matrix ${\bf\widetilde{K}}_{2N,M}$, which comprises the calculation of finite element solution followed by its evaluation at the nodes of the BEM. Because of the spectral accuracy of the Nyström BEM approximation, the DoF $2N$ is expected to be smaller, in practice, even compared to the number $M$ of FEM boundary Dirichlet (constrained) nodes (that is, $M>2N$). Accordingly, in our implementation we use instead the representation

[TABLE]

where we recall that ${\bf A}_{L}={\bf L}_{L}{\bf D}_{L}{\bf L}_{L}^{\top}$ is symmetric. This representation requires solving $2N$ (independent) finite element problems, one for each column of $\widetilde{\bf K}_{2N,M}^{\top}$, and a (sparse) matrix-vector multiplication. The first process, consumes the bulk of computation time (but is a naturally parallel task w.r.t. $N$) and can be carried out with wall-clock time similar to solving one FEM problem [24, Section 5.1.5].

The common CPU time for the direct and iterative solver amounts to the assembly of the finite element matrices ${\bf A}_{L}$ and ${\bf D}_{L,M}$, the ${\bf L}{\bf D}{\bf L}^{\top}$ factorization of the former, the boundary element matrix ${\bf B}_{2N}$ and the auxiliary matrices ${\bf T}_{2N,L}$ and ${\bf P}_{M,2N}$. Consequently the major difference in computation between the two approaches is: (i) the construction and storage of the matrix in (3.19a), followed by exactly solving the linear system for the direct method; versus (ii) the setting up of the system (3.19a) for matrix-vector multiplication and approximately solving the linear system with the GMRES iterations. The former approach is faster especially if the number of GMRES iterations is not very low (in single-digits) because of modern fast multi-threaded implementation of the direct solver. However, the latter approach is memory efficient and needed especially for large scale 3-D models.

Using a desktop machine, with a $10$-core processor and $128$GB RAM, we were able to apply the direct solver to simulate the example 2-D models in Experiment 1 and 2, even with millions of FEM (sparse) DoF within our FEM-BEM framework . For one of the largest cases reported in Table 2, with $\mathbb{P}_{4}$ elements for the wavenumber $k=4\pi$ ($40$ wavelengths case), with

[TABLE]

the GMRES approach system setup CPU time was $172$ seconds; and the direct approach setup CPU time was $332$ seconds. Because of requiring $102$ GMRES iterations, the solve time to compute a converged iterative solution was 586 seconds. However, because of the very efficient multi-threaded direct solvers (in Matlab) the direct solve time to compute the exact solution was only 0.014 seconds.

The size of the interface linear system for the experiment is only $160\times 160$ and hence our algorithm can be very efficiently used for a large number of incident waves $u^{{\rm inc}}$, that occur only in the small interface system. Thus we conclude that our FEM-BEM framework provides options to apply direct or iterative approaches to efficiently simulate wave propagation in heterogeneous and unbounded media. For 2-D low and medium frequency models with sufficient RAM, it seems to be efficient even to use the direct solver, and for higher frequency cases iterative solvers are efficient because of the demonstrated well-conditioned property of the system.

5 Conclusions

In this article we developed a novel continuous and discrete computational framework for an equivalent reformulation and efficient simulation of an absorbed and scattered wave propagation model, respectively, in a bounded heterogeneous medium and an unbounded homogeneous free-space. The model is governed by the Helmholtz equation and a decay radiation condition at infinity. The decomposed framework incorporates the radiation condition exactly and is based on creating two overlapping regions, without truncating the full space unbounded propagation medium. The overlapping framework has the advantage of choosing a smooth artificial boundary for the unbounded region of the reformulation, and a simple polygonal/polyhedral boundary for the bounded part of the two regions. The advantage facilitates the application of a spectrally accurate BEM for approximating the scattered wave, and setting up a high-order FEM for simulating the absorbed wave. We prove the equivalence of the decomposed overlapping continuous framework and the given model. The efficiency of our two-dimensional FEM-BEM computational framework was demonstrated in this work using two sets of numerical experiments, one comprising a smooth and the other a non-smooth heterogeneous medium.

Acknowledgement

Víctor Domínguez thanks the support of the project MTM2017-83490-P. Francisco-Javier Sayas was partially supported by the NSF grant DMS-1818867.