Planewave density interpolation methods for 3D Helmholtz boundary integral equations

TL;DR

This paper presents planewave density interpolation methods that regularize singular and nearly singular integral kernels in 3D Helmholtz boundary integral equations, enabling accurate numerical evaluation regardless of target position.

Contribution

The paper introduces novel planewave density interpolation techniques for regularizing integral kernels in 3D Helmholtz equations, compatible with Nyström and boundary element methods.

Findings

Methods effectively handle singular and nearly singular integrals.

Numerical examples demonstrate accuracy in acoustic scattering problems.

Compatible with standard quadrature rules for improved computation.

Abstract

This paper introduces planewave density interpolation methods for the regularization of weakly singular, strongly singular, hypersingular and nearly singular integral kernels present in 3D Helmholtz surface layer potentials and associated integral operators. Relying on Green's third identity and pointwise interpolation of density functions in the form of planewaves, these methods allow layer potentials and integral operators to be expressed in terms of integrand functions that remain smooth (at least bounded) regardless the location of the target point relative to the surface sources. Common challenging integrals that arise in both Nystr\"om and boundary element discretization of boundary integral equation, can then be numerically evaluated by standard quadrature rules that are irrespective of the kernel singularity. Closed-form and purely numerical planewave density interpolation procedures are presented in this paper, which are used in conjunction with Chebyshev-based Nystr\"om and Galerkin boundary element methods. A variety of numerical examples---including problems of acoustic scattering involving multiple touching and even intersecting obstacles, demonstrate the capabilities of the proposed technique.

Tables

Table 1. Table 1: Far- and near-field errors in the solution of the Neumann problem ( 2 ) produced by the Nyström discretization of the BM integral equation ( 7 ) using the numerical PWDI procedure of Section  4.2 for M = 2 , 3 𝑀 2 3 M=2,3 and various grid sizes  N 𝑁 N .

Nyström method — BM integral equation
	$M=2$		$M=3$
$N$	Direct	Regularized	Direct	Regularized
	Far-field
4	$8.55\cdot {10}^{-2}$	$2.01\cdot {10}^{-1}$	$8.89\cdot {10}^{-2}$	$6.66\cdot {10}^{-2}$
8	$1.73\cdot {10}^{-3}$	$8.44\cdot {10}^{-3}$	$3.55\cdot {10}^{-4}$	$2.03\cdot {10}^{-3}$
16	$3.46\cdot {10}^{-4}$	$9.46\cdot {10}^{-4}$	$5.31\cdot {10}^{-5}$	$2.60\cdot {10}^{-4}$
32	$4.13\cdot {10}^{-5}$	$1.21\cdot {10}^{-4}$	$7.22\cdot {10}^{-6}$	$3.22\cdot {10}^{-5}$
	Near-field
4	$1.61\cdot {10}^{-0}$	$9.37\cdot {10}^{-1}$	$1.19\cdot {10}^{-0}$	$2.48\cdot {10}^{-0}$
8	$7.32\cdot {10}^{-3}$	$1.63\cdot {10}^{-2}$	$3.31\cdot {10}^{-3}$	$6.72\cdot {10}^{-3}$
16	$1.86\cdot {10}^{-3}$	$1.53\cdot {10}^{-3}$	$1.79\cdot {10}^{-4}$	$8.35\cdot {10}^{-4}$
32	$2.27\cdot {10}^{-5}$	$2.02\cdot {10}^{-4}$	$2.08\cdot {10}^{-5}$	$8.29\cdot {10}^{-5}$

Table 2. Table 2: Far- and near-field relative errors in the solution of the Dirichlet problem ( 1 ) produced by the boundary element discretization of the BW integral equation ( 6 ) using the analytical PWDI technique for different interpolation orders M 𝑀 M and mesh sizes  h ℎ h . The surface Γ Γ \Gamma considered in this example is a sphere of unit radius and centered at the origin.

Boundary element method
		Far-field		Near-field
$h$	DoF	$M=0$	$M=1$	$M=0$	$M=1$
		BW integral equation
$5.34\cdot {10}^{-1}$	114	$4.48\cdot {10}^{-2}$	$4.50\cdot {10}^{-2}$	$7.11\cdot {10}^{-2}$	$5.98\cdot {10}^{-2}$
$2.72\cdot {10}^{-1}$	400	$1.29\cdot {10}^{-2}$	$1.32\cdot {10}^{-2}$	$2.70\cdot {10}^{-2}$	$1.65\cdot {10}^{-2}$
$1.36\cdot {10}^{-1}$	1507	$3.54\cdot {10}^{-3}$	$3.54\cdot {10}^{-3}$	$2.76\cdot {10}^{-2}$	$4.14\cdot {10}^{-3}$
$6.87\cdot {10}^{-2}$	6009	$8.98\cdot {10}^{-4}$	$8.96\cdot {10}^{-4}$	$4.81\cdot {10}^{-3}$	$1.05\cdot {10}^{-3}$
		BM integral equation
$5.34\cdot {10}^{-1}$	114	$5.55\cdot {10}^{-2}$	$5.38\cdot {10}^{-2}$	$8.80\cdot {10}^{-2}$	$8.29\cdot {10}^{-2}$
$2.72\cdot {10}^{-1}$	400	$1.65\cdot {10}^{-2}$	$1.61\cdot {10}^{-2}$	$3.02\cdot {10}^{-2}$	$1.94\cdot {10}^{-2}$
$1.36\cdot {10}^{-1}$	1507	$3.97\cdot {10}^{-3}$	$3.95\cdot {10}^{-3}$	$3.04\cdot {10}^{-2}$	$5.49\cdot {10}^{-3}$
$6.87\cdot {10}^{-2}$	6009	$1.02\cdot {10}^{-3}$	$1.01\cdot {10}^{-3}$	$4.71\cdot {10}^{-3}$	$1.38\cdot {10}^{-3}$

Equations

$$\{\begin{array}{rl}\Delta u_D^s+k^2{u}_D^s=& 0\mathrm{in}{\mathbb{R}}^3\setminus \Omega ,\\ u_D^s+u^{\mathrm{inc}}=& 0\mathrm{on}\Gamma ,\\ \underset{|\mathbf{r}|\to \infty }{\lim }|\mathbf{r}|\left(\frac{\partial u_D^s}{\partial |\mathbf{r}|}-ik{u}_D^s\right)=& 0,\end{array}$$

$$\{\begin{array}{rl}\Delta u_N^s+k^2{u}_N^s=& 0\mathrm{in}{\mathbb{R}}^3\setminus \Omega ,\\ \frac{\partial u_N^s}{\partial \mathbf{n}}+\frac{\partial u^{\mathrm{inc}}}{\partial \mathbf{n}}=& 0\mathrm{on}\Gamma ,\\ \underset{|\mathbf{r}|\to \infty }{\lim }|\mathbf{r}|\left(\frac{\partial u_N^s}{\partial |\mathbf{r}|}-ik{u}_N^s\right)=& 0,\end{array}$$

$$\left(\mathcal{S}\phi \right)(\mathbf{r}):={\int }_{\Gamma }G(\mathbf{r},q)\phi (q)\mathrm{d}s(q)\text{mbox}and\left(\mathcal{D}\phi \right)(\mathbf{r}):={\int }_{\Gamma }\frac{\partial G(\mathbf{r},q)}{\partial n(q)}\phi (q)\mathrm{d}s(q),$$

$$\begin{array}{rl}\left(S\phi \right)(p):={\int }_{\Gamma }G(p,q)\phi (q)\mathrm{d}s(q),& \left(K^{\prime }\phi \right)(p):={\int }_{\Gamma }\frac{\partial G(p,q)}{\partial n(p)}\phi (q)\mathrm{d}s(q),\\ \left(K\phi \right)(p):={\int }_{\Gamma }\frac{\partial G(p,q)}{\partial n(q)}\phi (q)\mathrm{d}s(q),& \left(N\phi \right)(p):=f.p.{\int }_{\Gamma }\frac{{\partial }^{2}G(p,q)}{\partial n(p)\partial n(q)}\phi (q)\mathrm{d}s(q),\end{array}$$

$$u^s(\mathbf{r})=(\mathcal{D}\phi )(\mathbf{r})-\mathrm{i}\eta (\mathcal{S}\phi )(\mathbf{r})\mathbf{r}\in {\mathbb{R}}^3\setminus \Gamma .$$

$$(\mathrm{BW})\frac{1}{2}{\phi }_D(p)+(K{\phi }_D)(p)-\mathrm{i}\eta (S{\phi }_D)(p)=-u^{\mathrm{inc}}(p),p\in \Gamma ,$$

$$(\mathrm{BM})\frac{\mathrm{i}\eta }{2}{\phi }_N(p)-\mathrm{i}\eta (K^{\prime }{\phi }_N)(p)+(N{\phi }_N)(p)=-\frac{\partial u^{\mathrm{inc}}(p)}{\partial n(p)},p\in \Gamma .$$

$${\mathbf{e}}_i(p):=\frac{\partial \mathbf{x}}{\partial x_i}(p),i=1,2.$$

$$g_{ij}:=⟨{\mathbf{e}}_i,{\mathbf{e}}_j{⟩}_p={\mathbf{e}}_i\cdot {\mathbf{e}}_j,1\le i,j\le 2,$$

$$\mathrm{grad}\phi =(g^{11}{\partial }_1\phi +g^{12}{\partial }_2\phi ){\mathbf{e}}_1+(g^{21}{\partial }_1\phi +g^{22}{\partial }_2\phi ){\mathbf{e}}_2={\partial }_1\phi {\mathbf{e}}^1+{\partial }_2\phi {\mathbf{e}}^2.$$

$$\mathrm{Hess}(\phi ):T_p\Gamma \to T_p\Gamma ,\mathrm{Hess}(\phi )(Y)={\nabla }_Y\mathrm{grad}\phi ,Y\in T_p\Gamma ,$$

$${\nabla }_{{\mathbf{e}}_i}{\mathbf{e}}_j=\sum _{\ell =1}^2{\Gamma }_{ik}^{\ell }{\mathbf{e}}_{\ell },$$

$$\mathrm{Hess}(\phi )=\sum _{i,j=1}^2\left({\partial }_i{\partial }_j\phi -\sum _{\ell =1}^2{\Gamma }_{ij}^{\ell }{\partial }_{\ell }\phi \right){\mathbf{e}}^i\otimes {\mathbf{e}}^j,$$

$${\partial }^{\alpha }\phi :={\partial }_1^{{\alpha }_1}{\partial }_2^{{\alpha }_2}\phi =\frac{{\partial }^{|\alpha |}\phi }{\partial {x_1}^{{\alpha }_1}\partial {x_2}^{{\alpha }_2}}$$

$$\phi ({\exp }_p(v))=\phi (p)+v^{\top }\mathrm{grad}\phi (p)+\frac{1}{2}{v}^{\top }\mathrm{Hess}(\phi )(p)v+\mathcal{O}(|v{|}^3),p\in \Gamma ,v\in T_p(\Gamma )$$

$$\phi (q)=\phi (p)+v^{\top }\mathrm{grad}\phi (p)+\frac{1}{2}{v}^{\top }\mathrm{Hess}(\phi )(p)v+\mathcal{O}(|p-q{|}^3),p,q\in \Gamma ,$$

$$\phi (q)=\phi (p)+\sum _{j=1}^M{\nabla }_{\Gamma }^j\phi (p)[v\otimes \cdots \otimes v]+\mathcal{O}(|p-q{|}^{M+1}),q,p\in \Gamma ,$$

$${\nabla }_{\Gamma }^j\phi (Y)={\nabla }_Y({\nabla }_{\Gamma }^{j-1}\phi ),j\ge 2,{\nabla }_{\Gamma }\phi :=\mathrm{grad}\phi ,Y\in T_p\Gamma ,$$

$${\Delta }_{\mathbf{r}}\Phi (\mathbf{r},p)+k^2\Phi (\mathbf{r},p)=0,\mathbf{r}\in {\mathbb{R}}^3\mathrm{for}\mathrm{all}p\in \Gamma .$$

$$\Phi (q,p):=\underset{\epsilon \to 0}{\lim }\Phi (q+\epsilon n(q),p)\text{mbox}and{\Phi }_n(q,p):=\underset{\epsilon \to 0}{\lim }\nabla \Phi (q+\epsilon n(q),p)\cdot n(q),$$

$$1_{\Omega }(\mathbf{r})\Phi (\mathbf{r},p)=-{\int }_{\Gamma }\frac{\partial G(\mathbf{r},q)}{\partial n(q)}\Phi (q,p)\mathrm{d}s(q)+{\int }_{\Gamma }G(\mathbf{r},q){\Phi }_n(q,p)\mathrm{d}s(q)$$

$$\begin{array}{r}{u}^s(\mathbf{r})=-1_{\Omega }(\mathbf{r})\Phi (\mathbf{r},p)+{\int }_{\Gamma }\frac{\partial G(\mathbf{r},q)}{\partial n(q)}\left\{\phi (q)-\Phi (q,p)\right\}\mathrm{d}s(q)\\ -{\int }_{\Gamma }G(\mathbf{r},q)\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}\mathrm{d}s(q),\end{array}$$

$$\frac{1}{2}\left\{\phi (p)-\Phi (p,p)\right\}+{\int }_{\Gamma }\frac{\partial G(p,q)}{\partial n(q)}\left\{\phi (q)-\Phi (q,p)\right\}\mathrm{d}s(q)$$

$$-{\int }_{\Gamma }G(p,q)\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}\mathrm{d}s(q)$$

$$\underset{q\to p}{\lim }{\partial }^{\alpha }\left\{\phi (q)-\Phi (q,p)\right\}=0\text{mbox}forall|\alpha |\le M,\text{mbox}and$$

$$\underset{q\to p}{\lim }{\partial }^{\alpha }\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}=0\text{mbox}forall|\alpha |\le M,$$

$$\left|\phi (q)-\Phi (q,p)\right|\lesssim |q-p{|}^{M+1}\text{mbox}and\left|\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right|\lesssim |q-p{|}^{M+1},$$

$$\left|\frac{\partial G(p,q)}{\partial n(q)}\left\{\phi (q)-\Phi (q,p)\right\}\right|\lesssim |q-p{|}^M\text{mbox}and\left|G(p,q)\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}\right|\lesssim |q-p{|}^M.$$

$$\begin{array}{rl}\frac{1}{2}\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}+{\int }_{\Gamma }\frac{{\partial }^{2}G(p,q)}{\partial n(p)\partial n(q)}\left\{\phi (q)-\Phi (q,p)\right\}\mathrm{d}s(q)& \\ -{\int }_{\Gamma }\frac{\partial G(p,q)}{\partial n(p)}\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}\mathrm{d}s(q)=-\frac{\partial u^{\mathrm{inc}}(p)}{\partial n(p)}& \text{mbox}forallp\in \Gamma ,\end{array}$$

$$\begin{array}{rl}\left|\frac{{\partial }^{2}G(p,q)}{\partial n(p)\partial n(q)}\left\{\phi (q)-\Phi (q,p)\right\}\right|\lesssim & |q-p{|}^{M-2}\text{mbox}and\\ \left|\frac{\partial G(p,q)}{\partial n(p)}\left\{\mathrm{i}\eta \phi (q)-{\Phi }_n(q,p)\right\}\right|\lesssim & |q-p{|}^M,\end{array}$$

$$\left|\frac{\partial G(\mathbf{r},q)}{\partial n(q)}\{\phi (q)-\Phi (q,p^{\ast })\}\right|\lesssim \frac{|q-p^{\ast }{|}^{M+1}}{|q-\mathbf{r}{|}^2}\le |q-p^{\ast }{|}^{M-1}$$

Full text

**Keywords: Helmholtz equation, integral equations, Nyström methods, boundary element methods

**

AMS subject classifications: 65N38, 35J05, 65T40, 65F08.

1 Introduction

Challenging weakly singular, strongly singular, hypersingular, and nearly singular surface integrals are ubiquitous to boundary integral equation (BIE) formulations of linear partial differential equations (PDEs). A plethora of numerical and semi-analytical procedures including, for instance, singularity subtraction111This technique has also been referred to as singularity extraction by some authors. [6, 18, 24, 23, 36, 38], Duffy-like transformations [14, 32, 31, 29, 21], polar singularity cancelation [4, 22, 34], singularity extraction [33, 35], among others techniques, have been proposed in the literature for the evaluation of these difficult integrals in the context of both Nyström and boundary element methods. Despite all these significant efforts and the compelling advantages that BIE methods offer over standard volume discretization techniques such as finite element and finite difference methods—especially in handling unbounded domains and seamlessly incorporating radiation conditions at infinity for time-harmonic wave scattering—they still face criticism of being difficult to implement. From the authors’ viewpoint, the main source of practical difficulties arises from the significant effort researchers and practitioners have to invest into understanding and implementing a specific set of techniques tailored to handle the various integration scenarios concerning the target point location relative to the surface sources. We hereby address this issue for both Nyström and boundary element discretizations of three-dimensional Helmholtz BIEs by introducing a universal semi-analytical procedure capable of regularizing all the aforementioned challenging surface integrals at the continuous level, i.e., prior to numerical integration. For the sake conciseness, we specifically consider combined-field BIE formulations of sound-soft (Dirichlet) and sound-hard (Neumann) scattering problems leading to the well-known Brakhage-Werner [2] and Burton-Miller [5] integral equations, respectively, which feature all four boundary integral operators of Calderón calculus.

As mentioned above, there is extensive literature on the subject. We refer the reader to [28] for a thorough review concerning Nyström methods and various approaches to deal with nearly singular integrals. Regarding BEMs specifically, two main groups of techniques can be distinguished. On one hand we have semi-analytical techniques [6, 18, 24, 23, 36, 38] whereby singular terms are extracted from the kernel to be integrated in closed form, while the remaining smoother part is integrated numerically by means of standard quadrature rules. And, on the other hand, we have techniques based on regularizing coordinate transformations [14, 32, 31, 29, 21, 22, 34] whereby specialized changes of variables are utilized to turn singular integrands into regular (analytic) integrands to which standard quadrature rules can be directly applied to achieve any desired accuracy. Although effective at dealing with the specific classes of integrands (on polygonal surface meshes) and basis functions for which they have been designed, none of the aforementioned techniques handle nearly singular integrals arising when target points lying off the surface are close to the surface sources.

In detail, this paper presents planewave density interpolation (PWDI) methods for the regularization of weakly singular, strongly singular, hypersingular and nearly singular integral kernels present in Helmholtz layer potentials and associated boundary integral operators. Relying on Green’s third identity and a certain Taylor-like interpolation of the surface density in terms of homogeneous solutions of the underlaying PDE (planewaves in this case), density interpolation methods [27, 28] allow layer potentials and operators to be expressed in terms of integrand functions that are smooth (at least bounded) regardless the target point location. The resulting surface integrals can then be numerically evaluated by means of standard off-the-shelf quadrature rules that are irrespective of the singularity of the associated integral kernels. As such, kernel-regularized layer potentials and operators can be directly evaluated at target points that are arbitrarily close to their surface sources enabling, in particular, the straightforward Nyström or boundary element method (BEM) discretization of BIEs involving multiple obstacles that are close, touching, or even intersecting each other. Indeed, we demonstrate through numerical experiments that BIEs posed on the surface of composite obstacles, i.e., obstacles that can be expressed as unions of geometrically simpler intersecting obstacles, can be recast as BIEs posed on the union of the boundaries of the simpler domains, which, upon application of the proposed PWDI kernel-regularization technique, can be directly solved using BEM retaining the expected order of convergence. This aspect of the proposed technique may significantly simplify the numerical solution of many real-world problems involving intricate obstacles, as it effectively allows bypassing the often involved task of meshing complex surfaces.

The structure of this paper is as follows: The theoretical basis of Taylor interpolation on regular surfaces, and of density interpolation methods in general, are established in Section 3. Two PWDI procedures are next introduced in Section 4. One amounts to the non-trivial extension to three dimensions of the low-order closed-form analytic procedure put forth in [27] (Section 4.1), while the other is a purely numerical procedure for the construction of arbitrarily high-order planewave density interpolants (Section 4.2). Section 5 then provides the details on the discretization of kernel-regularized layer potentials and integral operators by means of a Chebyshev-based Nyström method (Section 5.1) and a Galerkin BEM (Section 5.2). (The compatibility of the proposed approach with fast methods is addressed in Appendix C.) Section 6, finally, presents a variety of numerical examples that validate and demonstrate the various capabilities of the PWDI technique in the context of both Nyström and boundary element methods.

2 Preliminaries

For the sake of definiteness we focus in this paper on scattering problems related to acoustic sound-soft and sound-hard scatterers, e.g., either Dirichlet or Neumann boundary conditions. We thus seek scattered fields that are solutions of the following exterior Dirichlet and Neumann boundary value problems

[TABLE]

and

[TABLE]

respectively, where $\Omega\subset\mathbb{R}^{3}$ is a bounded obstacle whose boundary $\Gamma$ is a piecewise smooth, oriented and closed surface. (The incident fields $u^{\mathrm{inc}}$ in equations (1) and (2) are assumed to be solutions of the Helmholtz equation in all of $\mathbb{R}^{3}$.)

The Dirichlet (1) and Neumann (2) scattering problems can be formulated via well-posed boundary integral equations by means of the combined field approach introduced by Brakhage-Werner [2] and Burton-Miller [5], respectively. The combined field approach relies on the use of Helmholtz single- and double-layer potentials, hereby denoted as

[TABLE]

for $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$, respectively, where $G(\boldsymbol{r},\boldsymbol{r}^{\prime}):=(4\pi)^{-1}{\operatorname{e}^{\operatorname{i}\!k|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}}/{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}$ is the outgoing free-space Green function for the Helmholtz equation in $\mathbb{R}^{3}$ with wavenumber $k>0$. (In what follows we utilize the symbol $\boldsymbol{r}$ to denote points that do not lie of the surface $\Gamma$ while the symbols $p$ and $q$ are used exclusively to refer to points on the surface $\Gamma$.)

Interior and exterior Dirichlet/Neumann traces of the single- and double-layer potentials give rise to the four boundary integral operators of the Calderón calculus associated with the Helmholtz equation. Specifically, the Helmholtz single-layer ($S$), double-layer ($K$), adjoint double-layer ($K^{\prime}$) and hypersingular ($N$) operators are defined as

[TABLE]

for $p\in\Gamma$, where $\mathbb{n}(q)$ denotes the outward pointing unit normal to $\Gamma$ at $q\in\Gamma$. As usual, the initials f.p. in the definition of the hypersingular operator $N$ stand for Hadamard finite-part integral.

The combined field approach consists of looking for a scattered field $u^{s}=u^{s}_{D}$ (resp. $u^{s}=u^{s}_{N}$) in the form

[TABLE]

where $\varphi=\varphi_{D}:\Gamma\to\mathbb{C}$ (resp. $\varphi=\varphi_{N}:\Gamma\to\mathbb{C}$) is an unknown density function and $\eta\in\mathbb{R}$ is the coupling parameter. The enforcement of Dirichlet and Neumann boundary conditions on $\Gamma$ leads to the following combined field boundary integral equations

[TABLE]

and respectively

[TABLE]

Both boundary integral equations (6) and (7) are well-posed in appropriate functional spaces provided that $\eta\in\mathbb{R},\eta\neq 0$ [9].

As is well known, one of the main challenges in the numerical discretization of boundary integral equations (6) and (7) is posed by the singular character of the kernels of the boundary integral operators defined in equations (4) as the integration point $q$ approaches the target point $p$. Indeed, for a sufficiently regular surface $\Gamma\subset\mathbb{R}^{3}$ the operators $S$, $K$, and $K^{\prime}$ feature kernels with weak (integrable) singularities of type $\mathcal{O}(|p-q|^{-1})$, while the operators $N$ feature hyper-singular kernels of type $\mathcal{O}(|p-q|^{-3})$ as $\Gamma\ni q\to p\in\Gamma$. The numerical evaluation of the layer potentials (3), on the other hand, faces the significant challenge of dealing with the nearly singular character of the integral kernels at points $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ lying near the boundary at which, although smooth, the integrands exhibit large derivatives that ultimately hinder the accuracy of standard integration procedures.

In what follows we present a density interpolation method aimed at expressing the boundary integral operators (4) and layer potentials (3) in terms of surface integrands of prescribed regularity. For presentation simplicity and without loss of generality, instead of treating each one of the integral operators (4) and layer potentials (3) separately, we focus on the combined field integral operators of the BW (6) and BM (7) integral equations and the associated combined field potential (5).

3 Kernel regularization via density interpolation

Before we briefly embark on the presentation of the proposed density interpolation method, we first state some useful results of the differential geometry of surfaces that will provided the theoretical basis and the notation for the derivations presented below in this section. The main result of the next section is summarized in Remark 3.1.

3.1 Taylor series on smooth surfaces

We assume throughout this section that $\Gamma$ is a regular surface. First, given a system of coordinates around $p\in\Gamma$ with $\mathbf{x}(x_{1},x_{2})=p,\ \mathbf{x}:V\subset\mathbb{R}^{2}\to\Gamma$, we define the covariant basis of the tangent space $T_{p}\Gamma$ of $\Gamma$ at a point $p$ as

[TABLE]

We will follow the usual convention of not using the argument $p$ whenever there is no possibility of confusion. Using the Riemannian metric tensor

[TABLE]

we define the contravariant basis as $\mathbf{e}^{i}:=\sum_{j=1}^{2}g^{ij}\mathbf{e}_{j}$, $i=1,2,$ in terms of the inverse of the metric tensor $(g^{ij})=(g_{ij})^{-1}$. We also denote by $g$ the determinant of the metric tensor $(g_{ij})$, that is $g=g_{11}g_{22}-g_{12}^{2}$. With these notations in place, we have that the unit normal at $p\in\Gamma$ is given by $\mathbb{n}={\mathbf{e}_{1}\wedge\mathbf{e}_{2}}/{\sqrt{g}}.$

Given a function $\varphi:\Gamma\to\mathbb{C}$ we define its tangential gradient (or contravariant gradient) by the formula $\langle d\varphi,X\rangle_{p}={\rm grad}\ \varphi\cdot X$, for all $X\in T_{p}\Gamma,$ where $d\varphi=\partial_{1}\varphi\ \mathbf{e}_{1}+\partial_{2}\varphi\ \mathbf{e}_{2}$ is a $1-$form. An explicit formula for ${\rm grad}\ \varphi$ is given by

[TABLE]

We also define the Hessian of $\varphi$, ${\rm Hess}(\varphi)$ at $p\in\Gamma$ as the linear operator

[TABLE]

where $\nabla$ is the Riemannian connection on $\Gamma$. The latter can be expressed as

[TABLE]

in terms of the Christoffel symbols defined by $\displaystyle\Gamma^{\ell}_{ij}:=\frac{\partial\mathbf{e}_{i}}{\partial x_{j}}\cdot\mathbf{e}^{\ell}.$

It can be shown that ${\rm Hess}(\varphi)$ can be also viewed as a symmetric bilinear form on $T_{p}\Gamma$ given by ${\rm Hess}(\varphi)(X,Y)=\langle{\rm Hess}(\varphi)X,Y\rangle_{p}$, $X,Y\in T_{p}\Gamma$. The expression of ${\rm Hess}(\varphi)$ can be computed explicitly in the form

[TABLE]

where $\mathbf{e}^{i}\otimes\mathbf{e}^{j}=\mathbf{e}^{i}\ (\mathbf{e}^{j})^{\top}$. For a scalar function $\varphi:\Gamma\to\mathbb{C}$ and a multi-index $\alpha=(\alpha_{1},\alpha_{2}),$ $\alpha_{j}\in\mathbb{Z}$, $\alpha_{j}\geq 0$, $j=1,2$, we denote by

[TABLE]

where $|\alpha|=\alpha_{1}+\alpha_{2}$.

Finally, we need to make use of the exponential map on $\Gamma$. This map is defined on an open neighborhood $\mathcal{U}$ of the origin in $T_{p}\Gamma$, that is $\exp_{p}:\mathcal{U}\subset T_{p}\Gamma\to\Gamma$ such that, for $v\in\mathcal{U}$ with $|v|$ small enough, $\exp_{p}(v)$ is defined as the point on $\Gamma$ which is distance $|v|$ away on the geodesic originating at $p$ and having velocity $v/|v|$ at $p$. With these notations in place, we are in the position to state Taylor’s formula in the form

[TABLE]

as $|v|\to 0$, or equivalently as

[TABLE]

as $|p-q|\to 0$, where $v=\exp_{p}^{-1}(q)\in T_{p}\Gamma$ in the case when $\varphi$ is a smooth function defined on $\Gamma$. Taylor’s formula can be carried to higher-order terms in the form

[TABLE]

as $|p-q|\to 0$, where $v=\exp_{p}^{-1}(q)\in T_{p}\Gamma$. and where the $j$-th tensor $(\nabla_{\Gamma}^{j})\varphi:\underbrace{T_{p}\Gamma\times\cdots\times T_{p}\Gamma}_{j\ {\rm times}}\to\mathbb{R}$ is defined recursively as

[TABLE]

in terms of the Riemannian connection $\nabla$ on $\Gamma$. Clearly, as the surface gradient (9) and the Hessian (10), the higher order terms $\nabla_{\Gamma}^{j}\varphi$, $j\geq 2$, can be expressed as a linear combination of tensor products of the form $\mathbf{e}^{i_{1}}\otimes\cdots\otimes\mathbf{e}^{i_{j}},\ i_{\ell}\in\{1,2\}$, whose coefficients, in turn, can be expressed as a linear combination of $\partial^{\alpha}\varphi$ for all $|\alpha|\leq j$.

Remark 3.1**.**

The main take away message of this section is that Taylor’s formula (14) implies that if two smooth density functions, say $\varphi$ and $\psi$, are such that $\partial^{\alpha}\varphi(p)=\partial^{\alpha}\psi(p)$ for some $p\in\Gamma$ and for all $|\alpha|\leq M$, where the derivatives are taken with respect to any local parametrization of the surface around the point $p$, then $\varphi(p)=\psi(q)+\mathcal{O}(|p-q|^{M+1})$ as $\Gamma\ni q\to p\in\Gamma$. In the next section we will use that result produce a suitable Taylor-like interpolation of the density that will be used to regularize the boundary integrals.

3.2 Kernel-regularized boundary integral operators and layer potentials

Our density interpolation method relies on use of certain families of smooth functions $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$ that are solutions of the Helmholtz equation

[TABLE]

Letting

[TABLE]

for any given $p\in\Gamma$, denote the Dirichlet and Neumann traces of such functions, respectively, we have that an application of the Green’s third identity [9, 26] leads to

[TABLE]

for all $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ and $p\in\Gamma$, where $\boldsymbol{1}_{\Omega}$ denotes the characteristic function of the domain $\Omega$, i.e., $\boldsymbol{1}_{\Omega}=1$ in $\Omega$ and $\boldsymbol{1}_{\Omega}=0$ in $\mathbb{R}^{3}\setminus\overline{\Omega}$. Therefore, combining the layer potential (5) with formula (16) we obtain the following equivalent expression for the combined field potential (5):

[TABLE]

which is valid for all $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ and $p\in\Gamma$.

Letting then $\boldsymbol{r}=p+\varepsilon\mathbb{n}(p)$, $\varepsilon>0$, and taking the limit of both sides of equation (17) as $\varepsilon\to 0^{+}$ we obtain the following reformulation of the BW boundary integral equation (6)

[TABLE]

where we have utilized the standard jump conditions of the single- and double-layer operators [9, 26].

The scope of the proposed density interpolation technique is to explicitly and efficiently construct a family of functions $\Phi(\boldsymbol{r},p)$ such that the integrands that enter in equation (18) are regular (at least bounded) as $\Gamma\ni q\to p\in\Gamma$. To this end, for a given $\eta\in\mathbb{R},\eta\neq 0$, and a scalar function $\varphi:\Gamma\to\mathbb{C}$ which is assumed to be $(M+1)$-times continuously differentiable at $p\in\Gamma$, we say that a family of functions $\Phi(\boldsymbol{r},p)$ defined above satisfies Taylor-like interpolation conditions of order $M\geq 0$ at $p\in\Gamma$ if its Dirichlet and Neumann traces defined in equations (15), satisfy

[TABLE]

respectively, where all the derivatives are taken with respect to $q$ on the surface. In the light of the Taylor’s formula (14), it is clear that

[TABLE]

hold for all $p\in\Gamma$ at which the Taylor-like interpolation conditions (19b) are satisfied, regardless of the surface parametrization underlaying (19b) (see Remark 3.1). These estimates imply, in turn, that

[TABLE]

Therefore, from the estimates above we conclude that the proposed procedure effectively regularizes the singularities of the kernels of the boundary integral operators in equation (18) provided $\Phi$ satisfies the Taylor-like interpolation conditions (19b) for $M\geq 0$.

Similarly, taking the exterior normal derivative of the expression (17) we obtain that the BM integral equation (7) can be equivalently expressed as

[TABLE]

where the integrands satisfy

[TABLE]

and are at least bounded provided $\Phi$ satisfies the Taylor-like interpolation conditions (19b) at $p\in\Gamma$ for $M\geq 2$.

Finally, we apply the proposed density interpolation technique to the combined field potential (5) at observation points $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ near the boundary $\Gamma$. Letting $p=p^{*}={\rm arg}\min_{q\in\Gamma}|\boldsymbol{r}-q|\in\Gamma$ in the formula (17) for the combined field potential, we obtain that the corresponding integrands in (17) satisfy

[TABLE]

and

[TABLE]

provided $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$ interpolates $\varphi$—in the sense of the conditions in (19b)—at a nearly singular point $p=p^{*}\in\Gamma$. Clearly, for sufficiently large interpolation orders $M$, not only the integrands vanish at $p^{*}$, also their derivatives do. Finally, we mention that kernel-regularized expressions for the gradient of the combined field potential can be obtained by direct differentiation of (17).

As we will see in the next section (and in numerical results presented Section 6), expressions for the potential and its normal derivative stemming from (17) with $p=p^{*}={\rm arg}\min_{q\in\Gamma}|\boldsymbol{r}-q|$, can be exploited to produce kernel-regularized operators for problems involving multiple obstacles that are close or even intersecting each other.

Remark 3.2**.**

Derivations similar to the ones presented above can be carried out to produce kernel-regularized expressions for all four integral operators of Calderón calculus (4). In fact, such expressions for the double-layer and hypersingular operators are given by the left-hand-side of (18) and (21), respectively, that result from setting $\eta=0$ where $\Phi$ must satisfy (19b) with the corresponding $\eta=0$ value. Similarly, kernel-regularized expressions for the single-layer and adjoint double-layer operators are given by the left-hand-side of (18) and (21), respectively, that result from dividing them by $-\operatorname{i}\!\eta$ and taking the limit $\eta\to\infty$ where $\Phi$ must satisfy (19b) with the corresponding $(-\operatorname{i}\!\eta)^{-1}=0$ value.

Remark 3.3**.**

Maue’s formula [9, Theorem 2.23] provides an alternative expression for the hypersingular operator. In fact, for a sufficiently regular surface $\Gamma$ and density function $\varphi$, the hypersingular operator can be equivalently expressed as

[TABLE]

where tangential rotational operator at a point $p\in\Gamma$ is defined as $\vv{\operatorname{curl}}^{p}_{\Gamma}=-\mathbb{n}(p)\wedge{\rm grad}^{p}_{\Gamma}$ in terms of the surface gradient (9) at $p\in\Gamma$ which is here denoted as ${\rm grad}_{\Gamma}^{p}$. As usual the initials p.v. in front of the integral sign stand for principal value integral. Using Maue’s formula (22) we hence obtain that the BM integral equation (7) can be alternatively expressed as

[TABLE]

where the most singular integrand satisfies

[TABLE]

and is at least bounded provided $M\geq 2$. Both (21) and (23) forms of the BM integral equation are considered in the numerical examples presented in Section 6 below.

Remark 3.4**.**

Yet another expression for the hypersingular operator can be easily derived from (22) by “moving” the operator $\vv{\operatorname{curl}}_{\Gamma}^{p}$ outside the surface integral, i.e.,

[TABLE]

where the scalar rotational operator on the right-hand-side is defined as $\operatorname{curl}_{\Gamma}^{p}=\mathbb{n}(p)\cdot{\rm grad}_{\Gamma}^{p}$. It thus follows from (24) that the hypersingular operator can be evaluated by applying the proposed technique to the single-layer operator alone, although separate Taylor-like interpolants for each one of the three components of $\vv{\operatorname{curl}}_{\Gamma}^{q}\varphi(q)$ and $\mathbb{n}(q)\varphi(q)$ are needed. This approach is utilized in Section 6 to produce accurate BEM discretizations of the BM integral equation using $M=0$ and $1$.

3.3 Multiple-scattering approach to scattering by composite surfaces

Let $\Omega_{j}$, $j=1,2$, be open and simply connected domains with smooth boundaries $\Gamma_{j}=\partial\Omega_{j}$. Suppose $\Omega\subset\mathbb{R}^{3}$ is given by the union $\Omega=\Omega_{1}\cup\Omega_{2}$ where $\Omega_{1}\cap\Omega_{2}\neq\emptyset$. For the sake of conciseness and simplicity we focus here on the exterior Dirichlet problem (1) which we proceed to formulate as a multiple scattering problem encompassing the two obstacles $\Omega_{1}$ and $\Omega_{2}$. This formulation is advantageous in many practical applications where suitable discrete representations of the “combined surface” $\Gamma=\partial(\Omega_{1}\cup\Omega_{2})$ (in terms of surfaces meshes in the case BEM or manifold representations in terms of coordinate patches in the case Nyström methods) are difficult to produce, but separate discretizations of its component parts, $\Gamma_{1}=\partial\Omega_{1}$ and $\Gamma_{2}=\partial\Omega_{2}$, are easy to generate. An important example in this regard are Van der Waals molecular surfaces which are given by union of a typically large number of spherical atoms (e.g. [7]).

Instead of considering the BW integral equation (6) on $\Gamma=\partial(\Omega_{1}\cup\Omega_{2})$, we pose it on $\tilde{\Gamma}=\Gamma_{1}\cup\Gamma_{2}$. Letting $\tilde{\varphi}:\tilde{\Gamma}\to\mathbb{C}$ be a density function we look for the scattered field in the form of the combined field potential $u^{s}(\boldsymbol{r})=(\tilde{\mathcal{D}}\tilde{\varphi})(\boldsymbol{r})-\operatorname{i}\!\eta(\tilde{\mathcal{S}}\tilde{\varphi})(\boldsymbol{r})$ where $\tilde{\mathcal{D}}$ and $\tilde{\mathcal{S}}$ are the double- and single-layer potentials in (3), but defined in terms of surface integrals over $\tilde{\Gamma}$. The enforcement of the Dirichlet boundary condition on $\tilde{\Gamma}$ yields the integral equation

[TABLE]

on $\tilde{\Gamma}$ for the unknown density function $\tilde{\varphi}$, where $\varphi_{j}=\tilde{\varphi}|_{\Gamma_{j}}$ and $f_{j}|_{\Gamma_{j}}=-u^{\mathrm{inc}}|_{\Gamma_{j}}$, $j=1,2$. The operators $S_{ij}$ and $K_{ij}$, $i,j=1,2,$ in (25) are the single- and double-layer operators in (4) but defined in terms of boundary integrals over $\Gamma_{i}$ and target points $p\in\Gamma_{j}$. Note that the Dirichlet data $(f_{1},f_{2})$ in (25) requires the incident field $u^{\mathrm{inc}}$ to be defined on $\tilde{\Gamma}\setminus\Gamma$. Typically, $u^{\mathrm{inc}}$ is given by an explicit expression that can be directly evaluated almost everywhere in $\mathbb{R}^{3}$ including $\tilde{\Gamma}$. If that is not the case, $(f_{1},f_{2})$ can be defined by simply extending $u^{\mathrm{inc}}$ to $\tilde{\Gamma}\setminus\Gamma$ by zero. Uniqueness of solutions for the integral equation (25) is stablished in Appendix A.

In order to evaluate the integral operators in (25) we apply the proposed technique to each one of the operators involving integration over the closed surfaces $\Gamma_{1}$ and $\Gamma_{2}$. This is achieved by regularizing $I/2+K_{ij}-\operatorname{i}\!\eta S_{ij}$ for $i=j$ as an integral operator acting on $\Gamma_{i}$ using (18), and regularizing $K_{i,j}-\operatorname{i}\eta S_{ij}$ for $i\neq j$ as a layer potential that involves integration over $\Gamma_{i}$ and evaluation at target points $\boldsymbol{r}\in\Gamma_{j}$ using (17).

The effectiveness of this approach is demonstrated by numerical examples based on the BEM presented in Section 6. A more extensive study of the multiple-scattering approach to scattering by composite surface—that, in particular, will include the BM integral equation (7)—will be presented in a future contribution.

In the next section we present an explicit construction of families of functions $\Phi$ such that the conditions (19b) are satisfied for any given smooth function $\varphi$.

4 Interpolating functions

This section is devoted to the construction of the functions $\Phi$ introduced in the previous section. Specifically, we look for expressions of the kind

[TABLE]

where the derivatives $\partial^{\alpha}\varphi$ are defined in (11) and where the expansion functions $\Phi^{(1)}_{\alpha}$ and $\Phi^{(2)}_{\alpha}$ are taken to be linear combinations of planewaves:

[TABLE]

where $W_{\ell}(\boldsymbol{r}-p)=\exp\left\{{\rm i}k\boldsymbol{d}_{\ell}\cdot\left(\boldsymbol{r}-p\right)\right\}$, $\ell=1,\ldots,L$, have distinct directions $\boldsymbol{d}_{\ell}$ ($|\boldsymbol{d}_{\ell}|=1$). The directions $\boldsymbol{d}_{\ell}$ may depend on $p\in\Gamma$, in which case we make the dependence explicit in the notation $\boldsymbol{d}_{\ell}=\mathbb{d}_{\ell}(p)$. Clearly, any linear combination of the form (27) amounts to a homogeneous solution of the Helmholtz equation in the variable $\boldsymbol{r}$.

The following lemma, whose proof follows directly from Taylor’s theorem in two-dimensions, establishes simple point conditions on the traces of the expansion functions (27) that guarantee that the interpolation requirements (19b) are satisfied.

Lemma 4.1**.**

Let $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$ be given by (26), where $\Phi^{(j)}_{\alpha}:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$, $|\alpha|\leq M$, $j=1,2$, are the linear combinations of planewaves defined in (27). Then, sufficient conditions for $\Phi$ to satisfy (19b) at $p\in\Gamma$, are that the Dirichlet trace $\Phi^{(j)}_{\alpha}(q,p)$ and the Neumann trace $\Phi^{(j)}_{n,\alpha}(q,p)=\lim_{\epsilon\to 0}\nabla_{\boldsymbol{r}}\Phi^{(j)}_{\alpha}(q+\epsilon\mathbb{n}(q),p)\cdot\mathbb{n}(q)$, where $q\in\Gamma$ and $j=1,2$, satisfy

[TABLE]

for all sub-indices $\beta\in\mathbb{Z}^{2}_{+}$ such that $|\beta|\leq M$.

The following two sections address the problem of finding explicit expressions for $\Phi^{(1)}_{\alpha}$ and $\Phi^{(2)}_{\alpha}$, $|\alpha|\leq M$, by utilizing the point conditions (28b).

4.1 Closed-form planewave expansion functions in the case $M=1$

In this section we find closed-form expressions for the families of functions $\{\Phi^{(1)}_{\alpha}\}_{|\alpha|\leq M}$ and $\{\Phi^{(2)}_{\alpha}\}_{|\alpha|\leq M}$, defined in (27), whose traces satisfy the requirements in Lemma 4.1 for the interpolation order $M=1$.

We thus search for functions $\{\Phi^{(1)}_{\alpha}\}_{|\alpha|\leq 1}$ and $\{\Phi^{(2)}_{\alpha}\}_{|\alpha|\leq 1}$ that are linear combinations of planewaves whose directions $\boldsymbol{d}_{\ell}$ depend on $p\in\Gamma$, that is $W_{\ell}(\boldsymbol{r}-p)=\exp\left\{{\rm i}k\mathbb{d}_{\ell}(p)\cdot\left(\boldsymbol{r}-p\right)\right\}$, $\ell=1,\ldots,L$. The planewave directions $\mathbb{d}_{\ell}(p)$ are expressed in terms of the basis $\{{\boldsymbol{\tau}}_{1}(p),{\boldsymbol{\tau}}_{2}(p),\mathbb{n}(p)\}$, that is

[TABLE]

where the unitary contravariant vectors ${\boldsymbol{\tau}}_{j}(p),j=1,2,$ are defined by

[TABLE]

in terms of the Riemann metric tensor $(g_{ij})$ in (8) and its determinant $g$. Note that we dropped the dependence on $p$ of all the quantities in equations (29), as there is no risk of confusion.

We begin by defining, for $p\in\Gamma$ and a direction $\mathbb{d}(p)$ (which may not be unitary), the following functions

[TABLE]

for $\boldsymbol{r}\in\mathbb{R}^{3}$. Clearly, for unitary directions $\mathbb{d}_{1}(p)$ and $\mathbb{d}_{2}(p)$, the products

[TABLE]

[TABLE]

are linear combinations of planewaves of the form (27) with (unitary) directions $\frac{1}{\sqrt{2}}(\mathbb{d}_{1}(p)+\mathbb{d}_{2}(p))$ and $\frac{1}{\sqrt{2}}(\mathbb{d}_{1}(p)-\mathbb{d}_{2}(p))$ provided $\mathbb{d}_{1}(p)\cdot\mathbb{d}_{2}(p)=0$.

The following lemma introduces the sought linear combinations of planewaves:

Lemma 4.2**.**

Let $S$ and $C$ be the functions defined in (30) and denotes by $L=-\mathbf{e}_{1}\cdot\partial_{1}\mathbb{n}=\partial^{2}_{1}\mathbb{x}\cdot\mathbb{n}$, $M=-\mathbf{e}_{1}\cdot\partial_{2}\mathbb{n}=-\mathbf{e}_{2}\cdot\partial_{1}\mathbb{n}=\partial_{1}\partial_{2}\mathbb{x}\cdot\mathbb{n},$ and $N=-\mathbf{e}_{2}\cdot\partial_{2}\mathbb{n}=\partial^{2}_{2}\mathbb{x}\cdot\mathbb{n}$ the second fundamental form coefficients at $p\in\Gamma$. Then, the Dirichlet and Neumann traces of

[TABLE]

and

[TABLE]

satisfy the requirements (28b) and (28a) of Lemma 4.1, respectively, for $M=1$.

Proof.

The proof follows directly from the computation of the tangential derivatives of $\Phi^{(1)}_{\alpha}$ and $\Phi^{(2)}_{\alpha}$ at $\boldsymbol{r}=p\in\Gamma$ and the use of the identities [13]

[TABLE]

where $(g^{ij})$ denotes the inverse of the metric tensor $(g_{ij})$. ∎

The extension of the trigonometric ansatz technology utilized above to the construction of closed-form families of functions $\{\Phi^{(1)}_{\alpha}\}_{|\alpha|\leq M}$ and $\{\Phi^{(2)}_{\alpha}\}_{|\alpha|\leq M}$ that satisfy conditions in (28b) for $M\geq 2$ is challenging. In particular, additional distinct planewave directions ought to be incorporated in the ansatz. Because of the aforementioned difficulties, we advocate for the algebraic approach presented in following section to construct high-order planewave expansion functions in the case $M\geq 2$.

4.2 Higher-order planewave expansion functions

In this section we develop a purely algebraic algorithm to construct expansion functions $\{\Phi^{(1)}_{\alpha}\}_{|\alpha|\leq M}$ and $\{\Phi^{(2)}_{\alpha}\}_{|\alpha|\leq M}$, at a given (regular) point $p\in\Gamma$. Unlike the analytical approach presented in the previous section, we now select a collection of planewave directions $\{\boldsymbol{d}_{\ell}\}_{\ell=1}^{\ell=L}$, that are independent of $p\in\Gamma$. While the desired interpolation order $M$ and the number $L$ of planewave directions are parameters in our algorithm, the planewave directions themselves can be selected either randomly or uniformly from the unit sphere in three dimensions.

It is clear from Lemma 4.1 that in order to find the desired expansion coefficients $\{a_{\ell,\alpha}\}_{\ell=1,|\alpha|\leq M}^{\ell=L}$ (resp. $\{b_{\ell,\alpha}\}_{\ell=1,|\alpha|\leq M}^{\ell=L}$) that determine the expansion functions $\{\Phi^{(1)}_{\alpha}\}_{|\alpha|\leq M}$ (resp. $\{\Phi^{(2)}_{\alpha}\}_{|\alpha|\leq M}$) at $p\in\Gamma,$ one has to impose the $N=(M+1)(M+2)$ independent conditions (28a) (resp. (28b)) which have to be satisfied exactly. Consequently, a solvable linear system for the coefficients could be produced provided the number of planewave directions satisfies $L\geq N$. In order to form such a linear system, we proceed to sort the $N/2$ indices $\alpha=(\alpha_{1},\alpha_{2})$ satisfying $\alpha_{1}+\alpha_{2}\leq M$ by introducing a bijective mapping $J:\{|\alpha|\leq M\}\to\{1,\ldots,N/2\}$. Therefore, letting $\boldsymbol{\delta}_{j}$, $j=1,\ldots,N,$ denote the canonical vectors of $\mathbb{R}^{N}$ we have that conditions (28a) lead to the linear system

[TABLE]

for the coefficient vector $\mathbb{a}_{j}=[a_{1,J^{-1}(j)},\ldots,a_{L,J^{-1}(j)}]^{T}\in\mathbb{C}^{L}$, while conditions (28b) yield the system

[TABLE]

for the coefficient vector $\mathbb{b}_{j}=[b_{1,J^{-1}(j)},\ldots,b_{L,J^{-1}(j)}]^{T}\in\mathbb{C}^{L}$, where $\boldsymbol{C}(p)$ is a $N\times L$ complex-valued matrix. Note that we have assumed in these derivations that the first $N/2$ rows of $\boldsymbol{C}(p)$ correspond to the conditions on the Dirichlet traces, while the remaining $N/2$ rows correspond to the conditions on the Neumann traces.

In the case when $L>N$, the solution of the linear systems (34) and (35) must be understood in the least-squares sense. Letting $\boldsymbol{C}^{\dagger}(p)=\left[{\mathbb{c}}^{\dagger}_{1}(p),\dots,{\mathbb{c}}^{\dagger}_{{N}}(p)\right]\in\mathbb{C}^{L\times N}$ denote the Moore-Penrose pseudoinverse of $\boldsymbol{C}(p)$, that is $\boldsymbol{C}(p)\boldsymbol{C}^{\dagger}(p)=\boldsymbol{I}$ with $\boldsymbol{I}\in\mathbb{R}^{N\times N}$ being the identity matrix and $\boldsymbol{C}(p)$ being full-rank, the unknown vectors $\mathbb{a}_{j}$ and $\mathbb{b}_{j}$ are respectively given by

[TABLE]

Note that under the assumption that $\boldsymbol{C}(p)$ is full-rank, $\boldsymbol{C}^{\dagger}(p)$ can be computed explicitly via the formula $\boldsymbol{C}^{\dagger}(p)=\boldsymbol{C}^{*}(p)\left(\boldsymbol{C}(p)\boldsymbol{C}^{*}(p)\right)^{-1}$, where $\boldsymbol{C}^{*}$ denotes the Hermitian transpose of $\boldsymbol{C}$.

We present in Appendix B a recursive approach to compute the entries of the matrices $\boldsymbol{C}(p)$ in the case $M=3$.

Remark 4.3**.**

We have found difficult to prove whether for a given set of distinct planewave directions $\{\boldsymbol{d}_{\ell}\}_{\ell=1}^{\ell=L}$ the matrix $\boldsymbol{C}(p)\in\mathbb{C}^{N\times L}$ is full-rank for any $p\in\Gamma$ and will be left for future work. In practice, however, planewave directions selected from a “uniform” spherical grid give rise to numerically invertible matrices $\boldsymbol{C}(p)\boldsymbol{C}^{*}(p)$. In detail, the planewave directions for the construction of numerical PWDI interpolants used throughout this paper are given by $(\cos\theta_{m}\sin\phi_{n},\sin\theta_{m}\sin\phi_{n},\cos\phi_{n})$ where $\theta_{m}=2\pi(m-1/2)/L_{\theta}$ for $m=1,\ldots,L_{\theta}$ and $\phi_{n}=\pi(n-1/2)/L_{\phi}$ for $n=1,\ldots,L_{\phi}$, with $L=L_{\theta}\times L_{\phi}=2\times 2,4\times 3,5\times 4,6\times 5$ for interpolation orders $M=0,1,2,3$, respectively.

5 Numerical evaluation of integral operators and layer potentials

This section presents Nyström and BEM discretization schemes based on standard quadrature rules, for the numerical evaluation of the kernel-regularized integral operators and layer potentials associated to the combined field integral equations (6) and (7).

5.1 Chebyshev-based Nyström method

We here briefly describe the 3D boundary integral equation method introduced in our previous contribution [28]. The surface $\Gamma$ is represented as the union $\Gamma=\bigcup_{k=1}^{N_{p}}\overline{\mathcal{P}^{k}}$ of non-overlapping patches $\mathcal{P}^{k}$, $k=1,\dots,N_{p}$, where $\mathcal{P}^{k}\cap\mathcal{P}^{l}=\emptyset$ if $k\neq l$. It is assumed (throughout this section) that each surface patch $\mathcal{P}^{k}$ has associated a bijective $\mathcal{C}^{\infty}$ coordinate map $\mathbb{x}^{k}:\mathcal{H}\mathcal{\to}\overline{\mathcal{P}^{k}}$,

[TABLE]

where $\mathcal{H}=[-1,1]\times[-1,1]\subset\mathbb{R}^{2}$. Furthermore, the coordinate maps (37) are selected in such a way that the unit normal $\mathbb{n}^{k}(\boldsymbol{\xi})={\partial_{1}\mathbb{x}^{k}(\boldsymbol{\xi})\wedge\partial_{2}\mathbb{x}^{k}(\boldsymbol{\xi})}/{|\partial_{1}\mathbb{x}^{j}(\boldsymbol{\xi})\wedge\partial_{2}\mathbb{x}^{j}(\boldsymbol{\xi})|}$ at the point $\mathbb{x}^{k}(\boldsymbol{\xi})\in\mathcal{P}^{k}$ points outward to the surface $\Gamma$. The surface integral of a sufficiently regular function $F:\Gamma\to\mathbb{R}$—such as the integrands in (18), (21) and (17)—can then be expressed as

[TABLE]

In order to numerically evaluated the integral above with high-precision, we utilize open Chebyshev grids in the parameter space $\mathcal{H}$. Accordingly, $\mathcal{H}$ is discretized by means of the so-called Fejér’s first quadrature rules [12] which yields the approximation

[TABLE]

where the quadrature points $t_{j}$ are the Chebyshev zero points

[TABLE]

and the Fejér quadrature weights are given by

[TABLE]

A key feature of this discretization scheme is that the quadrature rule (38) yields spectral (super-algebraic) accuracy for integration of smooth $C^{\infty}(\mathcal{H})$ functions. As expected, however, slower convergence rates are achieved for less regular integrands (such as surface density functions associated to problems of scattering by piecewise smooth obstacles). Yet another important feature of this discretization scheme is that all the partial derivatives of the coordinate maps $\mathbb{x}^{k}$, unit normals $\mathbb{n}^{k},$ and functions $\varphi(\mathbb{x}^{k}(\boldsymbol{\xi}))$ that are needed for the construction of the planewave density interpolant $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$, can be efficiently and accurately computed at the grid points $\left(t_{i},t_{j}\right)$, $1\leq i,j\leq N$, by means FFT differentiation. More details can be found in [28].

Finally, the proposed procedure for the numerical evaluation of the BW combined field operator $K-{\rm i}\eta S$, using the method of Section 4.2, is summarized in Algorithm 1. A completely analogous procedure can be followed for evaluation of the BM combined field operator $N-{\rm i}\eta K^{\prime}$.

5.2 Galerkin boundary element method

This section concerns the use of the proposed planewave density interpolation method within the standard second-order BEM discretization using triangular surface meshes.

To fix ideas, we consider once again the BW boundary integral equation (6) which upon constructing an appropriate density interpolation function $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$, can be equivalently expressed (in strong form) as (18). Throughout this section we assume that $\Omega\subset\mathbb{R}^{3}$ is a bounded Lipschitz polyhedral domain. Therefore, associated to the surface $\Gamma=\partial\Omega$ there is a triangulation $\mathcal{T}_{h}$ such that $\Gamma=\overline{\bigcup_{T\in\mathcal{T}_{h}}T}$. Note that both the single- and double-layer operators are bounded on $H^{1/2}(\Gamma)$ and that the Green’s identities used in the derivation of (18) still hold true for $\Gamma$ being Lipschitz [25]. Therefore, assuming that $u^{\mathrm{inc}}|_{\Gamma}\in H^{1/2}(\Gamma)$ we readily have that the variational formulation of (18)—or equivalently (6)—reads as: Find $\varphi\in H^{1/2}(\Gamma)$ such that:

[TABLE]

for all test functions $\psi\in H^{-1/2}(\Gamma)$222We assume here the continuous extension of the standard real pairing $(u,v)=\int_{\Gamma}uv\,\mathrm{d}s$ for $u,v\in L^{2}(\Gamma)$ to the dual pairing $\langle\cdot,\cdot\rangle_{H^{1/2}(\Gamma)\times H^{-1/2}(\Gamma)}$..

In order to find an approximate weak solution $\varphi_{h}\in H^{1/2}(\Gamma)$ of (41), we resort to a Galerkin BEM for which we consider the (finite-dimensional) subspace $W_{h}=\{v\in C(\Gamma):v|_{T}\mbox{ is a linear function }\forall T\in\mathcal{T}_{h}\}\subset H^{1/2}(\Gamma)\subset H^{-1/2}(\Gamma)$. Clearly, the set $\{v_{j}\}_{j=1}^{N}$ of linear polynomials supported on $\bigcup_{p_{j}\in T}T$ satisfying the condition $v_{j}(p_{i})=\delta_{i,j}$ for all $i,j=1\ldots,N$, where $p_{j}$, $j=1,\ldots,N$ are the mesh nodes, forms a basis of the subspace $W_{h}$. In detail, letting $p_{\tau_{1}},p_{\tau_{2}}$ and $p_{\tau_{3}}$ with $\tau_{1},\tau_{2},\tau_{3}\in\{1,\ldots,N\},$ denote the vertices of a triangle $T\in\mathcal{T}_{h}$, and defining the unit vectors (see Figure 1)

[TABLE]

[TABLE]

we have that the sought approximate density function $\varphi_{h}\in W_{h}$ is given by

[TABLE]

where the basis functions are

[TABLE]

with $h_{j}=(p_{\tau_{i}}-p_{\tau_{j}})\cdot\boldsymbol{\nu}^{T}_{j}>0,$ $i\neq j$. The discrete variational formulation can thus be expressed as: Find $\varphi_{h}\in W_{h}$ such that:

[TABLE]

for all basis functions $v_{j}\in W_{h}$, $j=1,\ldots,N$.

Consider now the term in (44) associated with the single-layer operator; that is

[TABLE]

Clearly, both the outer and inner integrals in (45) have to be evaluated by means of quadrature rules for which typically the vertices of the triangules $T\in\mathcal{T}_{h}$ are used as quadrature points. Suppose the vertices $p_{\tau_{\ell}}$, $\ell=1,2,3$, are used as quadrature points for evaluation of the outer integral (over $T$). In order to regularize the integrand of the inner integral (over $K$) then, the normal derivative $\Phi_{n}$ of the planewave interpolant $\Phi$ has to approximate to sufficiently high-order the density function $\operatorname{i}\!\eta\varphi_{h}$ at the vertices $p_{\tau_{\ell}}$, $\ell=1,2,3$. Unfortunately, the problem with this numerical integration scheme is that the construction of $\Phi$ requieres the surface unit normal and surface tangent vectors to be properly defined at the vertices $p_{\tau_{\ell}}$, $\ell=1,2,3,$ which does not typically happen for general polyhedral surfaces.

In order to circumvent this issue we propose to use of a second-order Gauss quadrature rule for triangles that makes use of quadrature points $\tilde{p}_{\ell}$, $\ell=1,2,3,$ that lie in the interior $T^{\circ}$ of the triangle $T$ [11, 15]. This quadrature rule yields the approximation

[TABLE]

where the quadrature points are given by

[TABLE]

Since the quadrature points $\tilde{p}_{\ell}$, $\ell=1,2,3$, lie in the interior of the triangle $T$ (see Figure 1), the surface unit normal and the surface tangent vectors—which are required in the construction of the planewave interpolant—are uniquely defined at those points.

The inner integral, on the other hand, can be approximated by means of any sufficiently high-order quadrature rule. Using the standard node-based quadrature rule, for instance, we obtain

[TABLE]

where $p_{\kappa_{m}}$, $m=1,2,3,$ are the vertices of the triangle $K$ (see Figure 1). A completely analogous approach can be followed to evaluate the term in (44) involving the double-layer operator.

Remark 5.1**.**

Note that the terms inside the double sum on the right-hand side of (48) are always well defined even when the triangles $T$ and $K$ coincide. The effect of the PWDI technique lies then in the regularization of the nearly-singular integral kernels that arise when the triangles $T$ and $K$ coincide or are close to each other. In the latter case, however, the effectiveness of the proposed technique is affected by the limited (piecewise planar) global regularity assumed on the surface parametrization. It is thus not worth to pursue interpolation orders $M\geq 2$ in the context of the proposed BEM for piecewise planar surface representations.

In what follows we describe in some detail the construction of the closed-form planewave interpolant for $M=1$ in the BEM context. As discussed in Section 4.1 above, the construction of the planewave interpolant $\Phi:\mathbb{R}^{3}\times\Gamma\to\mathbb{C}$ requires the knowledge of a local smooth parametrization of the surface $\Gamma$ at and around the interpolation point $p\in\Gamma$. By construction, the interpolation point $p$ lies always in the interior of some triangle $T\in\mathcal{T}_{h}$. Therefore, the local surface parametrization has constant tangent vectors $\mathbb{e}^{T}_{1}$ and $\mathbb{e}^{T}_{2}$ that can be computed directly from the node data. In fact, selecting $\mathbb{e}^{T}_{1}=\boldsymbol{\nu}^{T}_{1}$ and $\mathbb{e}^{T}_{2}=\boldsymbol{\tau}^{T}_{1}$, for instance, we have that the expressions for the planewave interpolants (31)-(32) simplify significantly due to the fact that metric tensor becomes the identity, i.e.,  $g_{i,j}=\delta_{i,j}$, and the second fundamental form coefficients vanish, i.e., $L=M=N=0$, at $p\in T^{\circ}$.

Having defined the tangent vectors $\mathbb{e}^{T}_{1}$ and $\mathbb{e}^{T}_{2}$, the surface derivatives of $\varphi_{h}\in W_{h}$ at quadrature points $\tilde{p}_{\ell}\in T^{\circ}$, $\ell=1,2,3,$ can be computed by direct differentiation of (42) which in turn involves the derivatives

[TABLE]

of the basis functions (43). An algorithmic description of the numerical evaluation of $(v_{j},(K-\operatorname{i}\eta S)\varphi_{h})$ for $j=1,\ldots,N$, is given in Algorithm 2.

In order to tackle the BM integral equation (7), in turn, we apply the closed-form density interpolation technique to both single- and double-layer operators separately. In detail, we first resort to identities (22) and (24) to express the hypersingular operator in terms of single-layer operators. Then, upon integration by parts, the discrete variational formulation for the BM integral equation reads as: Find $\varphi_{h}\in W_{h}$ such that

[TABLE]

for all basis functions $v_{j}\in W_{h}$, $j=1,\ldots,N$, where $(\cdot,\cdot)$ denotes the standard real pairing $(\varphi,\psi)=\int_{\Gamma}\varphi(q)\cdot\psi(q)\,\mathrm{d}s(q)$. Noting that $(v_{j},K^{\prime}\varphi_{h})=(Kv_{j},\varphi_{h})$ we hence conclude that it suffices to apply the proposed technique to both $S$ and $K$ separately (see Remark 3.2).

Finally, in order to produce accurate evaluations of the combined field potential (5) at target points $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ near the surface $\Gamma$, we resort once again to the interior quadrature points (47). Indeed, in the context of the BEM the combined field potential at a point $\boldsymbol{r}\in\mathbb{R}^{3}\setminus\Gamma$ can be expressed as

[TABLE]

where the integrals over $T$ are approximated as follows

[TABLE]

6 Numerical examples

This section presents a variety of numerical experiments that illustrate different aspects of the proposed methodology.

6.1 Validation of the density interpolation procedures

Our first numerical example is devoted to the validation of the two density interpolation procedures introduced above in Section 4. We start of by taking $\Gamma$ as the (smooth) boundary of the bean-shaped obstacle displayed in Figure 2(a), on which we define the density function

[TABLE]

where $\varphi$ is a given smooth density and $\Phi$ is the planewave interpolant at $p^{*}\in\Gamma$. Note that, by construction, $\rho$ and its first $M$ tangential derivatives vanish at $p^{*}\in\Gamma$. As was discussed in Section 5.1, the surface $\Gamma$ is here represented by means of six non-overlapping rectangular patches each of which is discretized using Chebyshev grids consisting of $50\times 50$ points. Figure 2(a) (resp. 2(d)) displays the real part of $\rho$ produced by the closed-form (resp. algebraic) procedure. Figures 2(b) and 2(c) (resp. 2(e) and 2(f)), in turn, display slices of the real part of $\partial^{\alpha}\rho$, $|\alpha|=M$, at $p^{*}$ in the parameter space, obtained using the closed-form (resp. algebraic) procedure with $M=1$ (resp. $M=3$). The density function $\varphi$ utilized here is selected as the Dirichlet trace of the field produced by a point source at the point $\boldsymbol{r}_{0}=(0.1,-0.1,0.25)$ placed inside $\Gamma$. The wavenumber and the coupling parameter considered in this example are $k=10$ and $\eta=k$, respectively. These results demonstrate that the prescribed Taylor interpolation order $M$ is achieved by the proposed procedures. Similar results are obtained for the imaginary part of $\rho$ as well as for $\rho_{n}(q)={\rm i}\eta\varphi(q)-\Phi_{n}(q,p^{*})$, which, for the sake brevity, are not displayed here.

6.2 Nyström and Boundary Element methods

This section illustrates the capabilities of the density interpolation method for the regularization of the combined field potential and associated BW and BM integral operators.

6.3 Simple surfaces

In our first example we let $\Gamma$ be a unit sphere at the origin. (This simple surface has the advantage that can be easily represented using both quadrilateral patches and triangular meshes, allowing us to tackle the same problems using Nyström and BEM methods.) In order to assess the numerical errors, we consider an exact solution for both Dirichlet (1) and Neumann (2) problems, which is given by $u^{s}_{\rm exact}(\boldsymbol{r})=\operatorname{e}^{\operatorname{i}\!k|\boldsymbol{r}-\boldsymbol{r}_{0}|}/|\boldsymbol{r}-\boldsymbol{r}_{0}|-\operatorname{e}^{\operatorname{i}\!k|\boldsymbol{r}-\boldsymbol{r}_{1}|}/|\boldsymbol{r}-\boldsymbol{r}_{1}|$ where the source points $\boldsymbol{r}_{0}=(0.2,0.1,0.1)$ and $\boldsymbol{r}_{1}=(-0.1,0.3,-0.1)$ lie inside the unit sphere $\Gamma$. The real part of the Dirchlet trace of $u^{s}_{\rm exact}$ on $\Gamma$ (for $k=1$) is plotted (in colors) in Figures 3(a) and 3(b). The far-field errors in the numerical solution $u^{s}$ are measured by

[TABLE]

where $\Xi$ is the spherical grid (of radius $r=10$) displayed in Figure 3(a). The near-field errors, in turn, are measured using (53) but with $\Xi$ being the grid of the unit cube displayed in Figure 3(b). Note that $\Gamma$ touches the cube surface grid at the center of its six faces.

The accuracy of the Chebyshev-based Nyström method is assessed first. Figures 4(a) and 4(a) display the far- and near-field errors, respectively, in the approximate Dirichlet solution obtained from the kernel-regularized BW integral equation (18) for $k=\eta=1$. These figures display the errors obtained using the (closed-form) analytical (A-PWDI) procedure as well as the (algebraic) numerical (N-PWDI) density interpolation procedure introduced in Sections 4.1 and 4.2, respectively, for various discretization sizes $N$ and density interpolation orders $M$. The surface $\Gamma$ is here represented using six quadrilateral patches and each one of them is discretized using a Chebyshev grid of $N\times N$ points. The near fields, in particular, were computed using the kernel-regularized combined field potential (17). As can be observed in these results, the proposed technique yields third-order convergence of the far fields (as the grid size $N$ increases) for interpolation orders $M=0$ and $1$ and fifth-order convergence for $M=2$ and $3$. In the near field, on the other hand, third-order convergence is observed for all orders with the only exception of $M=3$ for which fifth order is achieved as $N$ increases. It is worth mentioning that the associated linear systems were solved iteratively by means of GMRES [30] with an error tolerance of $10^{-8}$. A nearly constant number of iterations ($\sim\!10$) was needed in all the examples considered in Figure 4. As in the case of the Laplace equation in 2D using the trapezoidal rule [28, Remark 6.1], the fact that the interpolation orders $M=0$ and $1$, and also $M=2$ and $3$, render the same order of convergence can be explained by the fact that Fejér’s quadrature rule in this case integrates certain odd singular terms exactly.

The geometric setup of Figure 3 is next used to assess the accuracy of the Chebyshev-based Nyström method when dealing with the more challenging BM integral equation (7) for the solution of the Neumann problem (2). The relevant numerical results are summarized in Table 1 where it can be clearly seen that, despite the overall smaller errors obtained for $M=3$, both interpolation orders $M=2$ and $M=3$ yield far- and near-field errors that exhibit the same nearly third-order convergence rate as $N$ increases. Direct (21) and regularized (23) versions of the hyper-singular operator were considered in these examples. A number of GMRES iterations as large as 200 was needed to achieved the desired accuracy in some of the examples considered in Table 1 due to the known unfavorable spectral properties of the hypersingular operator present in the BM integral equation. As is well known, the number of GMRES iterations can be drastically reduced by considering appropriate preconditioners for the hypersingular operator (cf. [1, 3]).

We next consider once again the Dirichet and Neumann problems posed in the exterior of the unit sphere, but now utilizing BEM discretizations of the associated BW (6) and BM (7) integral equations. The numerical results are summarized in Table 2. The close-form analytical density interpolation procedure of Section 5.2 is used in all the examples included in this table. The discrete variational formulations corresponding to the BW and BM integral equations are given in (44) and (50), respectively. We recall that the latter is here discretized as indicated in Section 5.2—by expressing it in terms of kernel-regularized single- and double-layer operators. As expected, these results demonstrate that far-field errors exhibit second-order convergence rates for both boundary integral equations and interpolation order $M=0$ and 1, as the mesh size $h=\max_{T\in\mathcal{T}_{h},i,j=1,2,3}|p_{\tau_{i}}-p_{\tau_{j}}|$ decreases. In fact, the errors obtained using the interpolation orders $M=0$ and $1$ are almost identical. The closeness of the errors observed might be explained by the possible dominance of the Galerkin-BEM $\mathcal{O}(h^{2})$ errors over the errors introduced by the numerical integration procedure. The near-field errors, on the other hand, exhibit nearly second-order convergence rates for both orders $M=0$ and $1$, with significantly smaller errors obtained for $M=1$.

In order to demonstrate the accuracy of the Nyström method when dealing with more complex geometries, we consider the scattering of a planewave $u^{\mathrm{inc}}(\boldsymbol{r})=\operatorname{e}^{\operatorname{i}\!k\boldsymbol{r}\cdot\boldsymbol{d}}$, in the direction $\boldsymbol{d}=(\cos\frac{\pi}{3},-\sin\frac{\pi}{3},0)$, that impinges on the three sound-soft obstacles shown in Figure 5. The resulting scattered field is solution of the exterior Dirichlet problem (1) that is here solved by means of the Nyström method applied to the kernel-regularized BW integral equation (18). The far-field errors reported in Figure 5 were produced by (53) with the set $\Xi$ being the spherical grid displayed in Figure 3(a). The reference solution $u^{s}_{\rm exact}$ in (53) was generated using a fine discretization of the BW integral equation consisting of Chebyshev grids comprising $36\times 36$ points per surface patch. Both (algebraic) numerical (with $M=2$) and (closed-form) analytical (with $M=1$) density interpolation procedures are utilized in this example. Third- and fifth-order convergence rates of the far-field errors are observed for $M=1$ and $2$, respectively, for the smooth surface cases (bean and ellipsoid). Only second-order convergence is achieved in the non-smooth surface case (cube) for both interpolation orders $M=1$ and $2$. The poor convergence rate observed in the latter case is explained by the singular behavior of the integral equation solution $\varphi$ along the edges of the cube (cf. [10]). In fact, for the interior point source problem described above in this section—in which case $\varphi$ is smooth up to the edges of the cube—third- and fifth-order convergence rates are attained. The total field solution of a higher-frequency scattering problem for the bean-shaped obstacle—whose diameter is $10\lambda$ ($k=2\pi/\lambda=10\pi$)—is shown in Figure 6. The near fields displayed in Figure 6 are accurate to at least four decimal places, everywhere including near and on the surface of the bean obstacle.

To finalize this section, we present examples aiming at demonstrating the capability of the BEM solver of handling complex geometries of engineering relevance. To this end, we consider a triangular mesh representation of a Falcon airplane produced by Gmsh [17], which is used in the solution of two Dirichlet problems (1) with different incident fields. In the first example we validate our BEM solver for this challenging geometry by considering an incident field given by two point sources placed inside the airplane’s fuselage. The numerical solution is then compared with the exact solution. For the wavenumber $k=0.5\pi$ and the mesh size $h=1.62$ we obtain a relative error (53) of $8.1\cdot 10^{-2}$ at a sphere containing the airplane. The near-field error at two different planes intersecting the airplane is displayed on the first row of Figure 7. Finally, our second example considers a planewave incident field in the direction $\boldsymbol{d}=(\cos\frac{\pi}{4},-\sin\frac{\pi}{4},0)$ for the same wavenumber ($k=0.5\pi$). Two views of the real part of total field are displayed in the second row of Figure 7 where it can be clearly seen the acoustic shadow under the airplane.

6.4 Composite surfaces

In this final section we apply the multiple-scattering approach put forth in Section 3.3, to the solution of the Dirichlet problem (1), posed in the exterior of the composite domain $\Omega=\Omega_{1}\cup\Omega_{2}$ with boundary $\Gamma=\partial(\Omega_{1}\cup\Omega_{2})$, where $\Omega_{1}$ is a sphere or radius $0.5$ and $\Omega_{2}$ is a hemisphere of radius $1.5$. The incident field used in this example is a planewave $u^{\mathrm{inc}}(\boldsymbol{r})=\operatorname{e}^{\operatorname{i}k\boldsymbol{d}\cdot\boldsymbol{r}}$ in the direction $\boldsymbol{d}=(\cos\frac{\pi}{4},0,-\sin\frac{\pi}{4})$ and $k=\eta=1$. The multiple-scattering BW integral equation (25), posed on $\tilde{\Gamma}=\partial\Omega_{1}\cup\partial\Omega_{2}$, is discretized using the BEM detailed in Section 5.2 with $M=1$. Figure 8 presents the far-field errors for various mesh sizes $h$. The error is defined here as in (53) with $\Xi$ being the spherical grid shown in Figure 3(a) and the reference solution $u^{s}_{\rm exact}$ being produced using a fine mesh discretization, with $h=0.11$, of the surfaces $\Gamma_{1}$ and $\Gamma_{2}$. Three different $\Omega$ configurations, shown in inset plots in Figure 8, are considered, including one (on the left-hand-side) in which the two obstacles are touching at a single point. Clearly, second-order convergence is achieved, as $h$ decreases, in all three configurations. The real part of the total field solution of the problem of scattering (for $k=\eta=5$) together with the absolute value of the error obtained using the multiple-scattering approach and the standard approach, are shown in Figure 9. The reference solution for the error estimation is produced using a fine-grid discretization of $\Gamma$ using $h=0.1$.

Appendix A Uniqueness of multiple scattering integral equation formulation for composite surfaces

This appendix is devoted to establishing uniqueness of solution of the boundary integral equation system (25). We present the following result whose proof asumes that $\Omega_{1}\cap\Omega_{2}\neq\emptyset$ is open:

Lemma A.1**.**

The integral equation system (25) has at most one solution.

Proof.

Assume that (25) admits a non-trivial solution $\tilde{\varphi}$ when $(f_{1},f_{2})=(0,0)$. The potential $u^{s}(\boldsymbol{r})=(\tilde{\mathcal{D}}\tilde{\varphi})(\boldsymbol{r})-\operatorname{i}\eta(\tilde{\mathcal{S}}\tilde{\varphi})(\boldsymbol{r})$ is then a homogeneous solution of the Helmholtz equation in $\mathbb{R}^{3}\setminus(\Gamma_{1}\cup\Gamma_{2})$ whose limits from outside vanish on both $\Gamma_{1}$ and $\Gamma_{2}$. Since the uniqueness of the exterior Dirichlet problem implies that $u^{s}=0$ in $\mathbb{R}^{3}\setminus(\Omega_{1}\cup\Omega_{2})$ (cf. [9]), we have from the integral equation and the jump conditions of the single- and double-layer potentials on $\Gamma_{j}$, that $[u^{s}]=\varphi_{j}$ and $[\partial_{n}u^{s}]=\operatorname{i}\!\eta\varphi_{j}$ on $\Gamma_{j}$, $j=1,2$. Therefore, $u^{s}$ is a homogeneous solution of the Helmholtz equation in $\Omega_{1}\cap\Omega_{2}$ that satisfies the Robin boundary condition $\partial_{n}u^{s}-\operatorname{i}\!\eta u^{s}=0$ on $\partial(\Omega_{1}\cap\Omega_{2})\subset\Gamma_{1}\cup\Gamma_{2}$. By uniqueness of this interior boundary value problem (cf. [9]) we have that $u^{s}=0$ in ${\Omega_{1}\cap\Omega_{2}}$. Similarly, we get that $u^{s}=0$ in ${\Omega_{1}\setminus\Omega_{2}}$ and in ${\Omega_{2}\setminus\Omega_{1}}$, and thus $u^{s}=0$ in all of ${\Omega_{1}\cup\Omega_{2}}$. Therefore, since $u^{s}=0$ in $\mathbb{R}^{3}\setminus(\Gamma_{1}\cup\Gamma_{2})$ we conclude that $[u^{s}]=\varphi=0$ on $\Gamma_{1}\cup\Gamma_{2}$, which is a contradiction. The proof is now complete. ∎

Appendix B Explicit expression of the matrix $\boldsymbol{C}(p)$ for $M=3$

This appendix is devoted to the derivation of the matrix $\boldsymbol{C}$ used in the algebraic procedure introduced in Section 4.2 for the construction of planewave interpolant in the case $M=3$.

In order to ease derivations we define

[TABLE]

where $\boldsymbol{d}_{\ell}$ ($|\boldsymbol{d}_{\ell}|=1$) is a given planewave direction that does not depend on the expansions point $p$. Evaluating all the partial derivatives of $w_{\ell}$ and $w_{n,\ell}$ up to third order, we get

[TABLE]

where the indices $\alpha,\beta,\gamma\in\mathbb{Z}_{+}^{2}$ satisfy $|\alpha|=|\beta|=|\gamma|=1$.

Selecting the ordering

[TABLE]

for the multi-indices, and enforcing the point conditions by evaluating the expressions (54b) at $q=p$ we obtain that the entries of $\boldsymbol{C}(p)$, that are denoted by $c_{n,\ell}$ for $n=1,\ldots,N$ and $\ell=1,\ldots,L$, are given by the following (recursive) relations

[TABLE]

[TABLE]

Remark B.1**.**

Since the identities $\partial^{\alpha}\mathbb{x}=0$ for $|\alpha|>1$ and $\partial^{\alpha}\mathbb{n}=0$ for $|\alpha|>0$ hold true for any planar coordinate patch, we have that letting $\tau_{1}={\rm i}k\left(\boldsymbol{d}_{\ell}\cdot\partial^{(1,0)}\mathbb{x}\right)$, $\tau_{2}={\rm i}k\left(\boldsymbol{d}_{\ell}\cdot\partial^{(0,1)}\mathbb{x}\right)$ and $\tau_{3}={\rm i}k\left(\boldsymbol{d}_{\ell}\cdot\boldsymbol{n}\right)$, the $\ell$-th column of $\boldsymbol{C}$ corresponding to a given planewave direction $\boldsymbol{d}_{\ell}$ is simply given the vector $[\mathbb{c},\tau_{3}\mathbb{c}]^{T}\in\mathbb{C}^{20}$ where $\mathbb{c}=[1,\tau_{1},\tau_{2},\tau_{1}^{2},\tau_{1}\tau_{2},\tau_{2}^{2},\tau_{1}^{3},\tau_{1}^{2}\tau_{2},\tau_{1}\tau_{2}^{2},\tau_{2}^{3}]$.

Appendix C Compatibility of the PWDI method with fast methods

This appendix discusses the compatibility of the PWDI technique with standard fast multipole methods (FMMs) for the acceleration of the solution of the combined field integral equation (6) by means of iterative linear algebra solvers such as GMRES [30]. In a nutshell, FMMs rely on a certain low-rank approximation of the relevant integral kernel $K:\Gamma\times\Gamma\to\mathbb{C}$ in the form

[TABLE]

which is typically valid for observation points $p$ that are sufficiently far away from the source point $q$. In the case of the BW combined field integral operator, for example, the kernel is given by $K(p,q)=\partial G(p,q)/\partial n(q)-{\rm i}\eta G(p,q),$ where $G$ is the free-space Green function defined in (2). When a kernel expansion of the form (58) is available, a FMM can be devised to effectively reduce the computational cost of the numerical evaluation of the integral operator

[TABLE]

from $\mathcal{O}(N_{t}^{2})$ to nearly $\mathcal{O}(N_{t})$ operations, where $N_{t}$ denotes the total number of quadrature points (resp. mesh nodes) used in the discretization of the boundary integral over $\Gamma$ by the Nyström (resp. boundary element) method. At low-frequencies, when the kernel $K$ is non-oscillatory (i.e., for obstacles $\Omega$ of just a few wavelengths $\lambda=2\pi/k$ in diameter), such approximation can be achieved, for example, by means of multipole/planewave expansions [20, 8, 19], polynomial interpolation [16], or projection methods based on equivalent sources [37].

Unfortunately, however, unlike the non-regularized combined field integral equation (6), the integral equation (18) resulting from application of the density interpolation technique requires evaluation of boundary integrals of the form

[TABLE]

for $p\in\Gamma$, where not only $G$ depends on $p$ and $q$ but also the Dirichlet and Neumann traces $u$ and $u_{n}$ of the interpolating function $U$ in (26). Therefore, in principle, any FMM for the numerical evaluation of (60) would require low-rank approximations of two additional kernels, namely

[TABLE]

Fortunately, direct use of existing FMMs for the evaluation of (60) is still possible when the algorithm introduced in Section 4.2 is used for the construction of the interpolating function $U$. Indeed, since the directions $\boldsymbol{d}_{\ell}$ of the planewaves

[TABLE]

are independent of the both observation and source points, the resulting interpolating functions is trivially separable, i.e,

[TABLE]

where the expansion coefficients are given by

[TABLE]

in terms of $a_{\ell,\alpha}$ and $b_{\ell,\alpha}$ being given by solution of the linear systems (34) and (35), respectively. Therefore, upon integration over $\Gamma$ we obtain

[TABLE]

for all $p\in\Gamma$. We thus conclude from here that the kernel regularized integral operator (59) can in fact be evaluated by direct application of FMM-accelerated double- and single-layer operators, to each one of the planewaves $w_{\ell}$, $\ell=1,\ldots,L$. Furthermore, since the coefficients $a_{\ell,\alpha}$ and $b_{\ell,\alpha}$ for $|\alpha|\leq M$ and $\ell=1,\ldots,L$, the multipole moments (or equivalent sources), and the local expansions coefficients associated to each one of the planewave $w_{\ell}$, are independent of the density function $\varphi$, they can be precomputed, stored, reused at each iteration of the linear algebra solver (GMRES) thus further reducing the computational cost associated to a forward map evaluation. It can be easily checked, finally, that the overall computational cost of one forward map evaluation using the proposed FMM-accelerated algorithm is $\mathcal{O}(N_{t}\log N_{t})$ for the high-order Nyström method of Section 5.1, and $\mathcal{O}(N_{t})$ for the low-order BEM of Section 5.2.