An Approach to the Implementation of Laplace and a Broadband Helmholtz Fast Multipole Method as an Application Independent Library

Sanjay Velamparambil

Ansys Inc., 1995 N. 57th Court, Suite 100, Boulder, CO-80301, USA
sanjay.velamparambil@ansys.com

Submitted On: August 9, 2023; Accepted On: July 8, 2024

ABSTRACT

In this paper, we propose an approach to develop an application independent library of Laplace and Helmholtz fast multipole method (FMM) that can be used in different applications. For this purpose, we consider a generalized problem and a corresponding canonical problem (defined below). In the first main contribution, we show that it is possible to capture the essential characteristics of the canonical summation from sampling the values of certain potentials or signature functions. In the second main contribution, we show that partial derivatives of arbitrary orders acting on the far field can be represented as product of sparse matrices within the library, transparent to the user. Combining the two ideas, we show that once the FMM is configured to compute the canonical summation, the same setup can be used to work with a much wider, general class of problems.

Index Terms: Fast multipole methods, integral equations.

I. INTRODUCTION

In this paper, we describe an approach to develop a stand-alone, application independent library of the fast multipole method (FMM) [1, 2, 3] to be used in different physical applications and developed by multiple teams with differing backgrounds. In particular, we are concerned with the evaluation of summations of the form:

	$y_m$	$={\displaystyle \sum _{n=1}^N}{A}_{mn}{x}_n,1\le m\le M,$		(1a)
with:
	$A_{mn}$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{g}_m(\overrightarrow{x}){\mathrm{ℳ}}_{\mathrm{\sqcap }}{\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}}{\mathrm{\{}}_{\mathrm{\backslash }}\mathrm{\Leftarrow }{\overrightarrow{\mathrm{\S }}}^{\prime }\mathrm{\Rightarrow }{\mathrm{ℒ}}_{\mathrm{⊑}}^{\prime }\mathrm{𝒢}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\iff }{\overrightarrow{\mathrm{\S }}}^{\prime }\mathrm{\Rightarrow }\mathrm{\lceil }{\mathrm{\Omega }}^{\prime }\mathrm{\lceil }\mathrm{\Omega }\mathrm{\iff }$		(1b)
		$\mathrm{\hspace{0.17em}}u,v=x,y,z;\overrightarrow{x},{\overrightarrow{x}}^{\prime }\in {ℝ}^3,$
where ${\mathrm{ℳ}}_{\mathrm{\sqcap }}$ and ${\mathrm{ℒ}}_{\mathrm{⊑}}$ are linear partial differential operators with constant coefficients and:
	$G(\overrightarrow{x},{\overrightarrow{x}}^{\prime })$	$=\{\begin{array}{cc}\frac{1}{\parallel \overrightarrow{x}-{\overrightarrow{x}}^{\prime }\parallel }\hfill & \text{for the Laplace Equation}\hfill \\ \frac{e^{ik\parallel \overrightarrow{x}-{\overrightarrow{x}}^{\prime }\parallel }}{\parallel \overrightarrow{x}-{\overrightarrow{x}}^{\prime }\parallel }\hfill & \text{for the Helmholtz Equation}.\hfill \end{array}$		(1c)

Note that, in general, ${\mathrm{\Omega }}^{\prime }\ne \mathrm{\Omega }$. Furthermore, following the terminology used by the method of moments (MoM) practitioners in electromagnetics, we shall refer to $f_j(\overrightarrow{x})$ and $g_i(\overrightarrow{x})$ as the basis and testing functions, respectively, with the understanding that they can be Dirac $\delta$-functions to represent “point” particles.

Such computations can arise from a wide variety of applications, such as electromagnetics, acoustics and elastic wave scattering. Since its development in the 1980s [4, 5] for the Laplace equation to the elegant extension to the Helmholtz and Maxwell equations in the 1990s [2, 3, 6], FMM and its multilevel variants have transformed computational physics, especially, computational electromagnetics. It is beyond the scope of this paper to provide a comprehensive overview of these remarkable developments. Instead, we refer the reader to the recent review by He et al. [7].

Although the theory of FMM has been well established and many high quality implementations exist, both in the public domain and in proprietary products, the development of an application independent library is still a challenging problem due to various practical constraints such as:

Different formulations and physics lead to different choices for ${\mathrm{ℳ}}_{\mathrm{\sqcap }}$ and ${\mathrm{ℒ}}_{\mathrm{⊑}}$.

Different domain and geometry; the integrals could be defined over surfaces or volumetric regions, discretized with surface elements such as flat or curvilinear triangles and/or quadrilaterals and volume elements such as tetrahedra. It is also possible to have mixed formulations involving both surface elements and volume elements. The quality of such a mesh can often be poor, with highly non-uniform elements.

Different basis and testing procedures in reducing the integral equations to a matrix equation; note that each basis and testing function can be supported by one or more mesh elements.

Familiarity/background knowledge of the library user. Even though one can expect the users to have a sound understanding of the integral equation and the underlying physics, it is not necessary that they have a sufficient understanding of FMM to develop the components needed to integrate their code with a library.

We illustrate this with an example. In an implementation of Laplace FMM that uses spherical harmonics, the library user may have to provide the multipole coefficients of basis/testing functions. Although standard implementations of spherical harmonics are available, it may be tricky to ensure consistency between third party implementation of spherical harmonics and the ones used in the FMM implementation.

One of the ways to address these issues is to replace the basis and testing functions with weighted (point) “particles” as shown in Fig. 1. Although this approach simplifies the implementation of the library, it has several drawbacks. First, as illustrated in Fig. 1, replacing mesh elements and basis/testing functions with weighted particles generally leads to incompatible partitioning of the geometry, requiring error-prone bookkeeping, complicated data interchange between the application and the library, and messy corrections in the near field. In addition, given that the number of quadrature points is generally larger than the number of basis/testing functions, the matrix as seen by the library tends to be larger as well, potentially affecting the performance.

Figure 1: The two triangles are near to each other, but some of the quadrature points can still be well separated (shown in ellipses) leading to tricky bookkeeping and near field corrections.

Therefore, the objective of this paper is to present an approach that can handle a wide class of problems through the use of good abstractions to capture the physics, the formulation and the discretization. The main contributions of this work are:

Abstractions to characterize the geometry of the elements in a way that allows consistent handling of non-uniform meshes.

A technique to compute required multipole coefficients of the shape/basis and testing functions from sampled values of their corresponding potentials or signature functions which are easy to define and simple to evaluate, without involving any special functions.

We show that the partial derivatives of the far field, represented by multipole or plane wave expansions, can be represented as products of sparse matrices independent of the basis and testing functions allowing the user to compute these derivatives without having to directly work with multipole expansions.

Finally, we demonstrate that once the library is integrated into an existing code to compute what we call the canonical problem, defined in equation (4), the same implementation can be used to evaluate summations in equation (1) without modifying the library and without adding new subroutines to the application code.

II. PRELIMINARIES

A. General class of problems

We assume that the reader is familiar with the basic mathematical formulation and the octree based computational structure of the Laplace and Helmholtz multi-level fast multipole methods (MLFMM). Except when specificity is needed, we shall use the abbreviation MLFMM to refer to either of the two versions.

In the context of MLFMM, we begin by noting that if $g_m$ and $f_n$ in equation (1) belong to well separated boxes, the Green’s function in equation (1) does not involve any singularities. Since:

${\displaystyle \frac{\partial G}{\partial u}}$

$=-{\displaystyle \frac{\partial G}{\partial u^{\prime }}},u=x,y,z,$

(2)

whenever $\overrightarrow{x}\ne {\overrightarrow{x}}^{\prime }$, the differential operator ${\mathrm{ℒ}}_{\mathrm{⊑}}^{\prime }$ appearing on the source point ${\overrightarrow{x}}^{\prime }$ can be transferred to the observation point $\overrightarrow{x}$ when $g_m$ and $f_n$ are geometrically well separated. Therefore, for far field computations, it is possible to transform equation (1b) to a linear combination of integrals of the form:

	$$A_{mn}={\int }_{\mathrm{\Omega }}{g}_m(\overrightarrow{x}){\mathrm{𝒟}}_{\mathrm{\surd }\mathrm{\coprod }\mathrm{\nabla }}{\int }_{{\mathrm{\Omega }}^{\prime }}{\mathrm{\{}}_{\mathrm{\backslash }}\mathrm{\Leftarrow }{\overrightarrow{\mathrm{\S }}}^{\prime }\mathrm{\Rightarrow }\mathrm{𝒢}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\iff }{\overrightarrow{\mathrm{\S }}}^{\prime }\mathrm{\Rightarrow }\mathrm{\lceil }{\mathrm{\Omega }}^{\prime }\mathrm{\lceil }\mathrm{\Omega }\mathrm{\iff }$$		(3a)
	$$\text{where}{\mathrm{𝒟}}_{\mathrm{\surd }\mathrm{\coprod }\mathrm{\nabla }}\mathrm{\Im }\frac{{\partial }^{\mathrm{\surd }\mathrm{\Downarrow }\mathrm{\coprod }\mathrm{\Downarrow }\mathrm{\nabla }}}{\partial {\mathrm{\S }}^{\mathrm{\surd }}\partial {\mathrm{†}}^{\mathrm{\coprod }}\partial {\mathrm{‡}}^{\mathrm{\nabla }}}\mathrm{\iff }\overrightarrow{\mathrm{\S }}\ne {\overrightarrow{\mathrm{\S }}}^{\prime }\mathrm{\iff }$$		(3b)
with $p,q$ and $r$ being non-negative integers.

One of the key results presented in this work is that the differential operator ${\mathrm{𝒟}}_{\mathrm{\surd }\mathrm{\coprod }\mathrm{\nabla }}$ can be exactly represented as a product of certain sparse matrices in the calculation of far interactions. Anticipating this result, we define a new summation:

	$y_m$	$={\displaystyle \sum _{n=1}^N}{P}_{mn}{x}_n,1\le m\le M,$		(4a)
	$\text{where}{P}_{mn}$	$={\displaystyle {\int }_{\mathrm{\Omega }}}{g}_m(\overrightarrow{x}){\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}}{f}_n({\overrightarrow{x}}^{\prime })G(\overrightarrow{x},{\overrightarrow{x}}^{\prime })\mathit{d}{\mathrm{\Omega }}^{\prime }\mathit{d}\mathrm{\Omega },$		(4b)
and refer to it as the canonical problem.

B. Some notations and definitions

If $r=\parallel \overrightarrow{x}\parallel$ and $\widehat{\mathbf{x}}=\frac{\overrightarrow{x}}{\parallel \overrightarrow{x}\parallel }$ for $\overrightarrow{x}\in ℝ$, we shall often write $\overrightarrow{x}=(r,\widehat{\mathbf{x}})=(r,\theta ,\varphi )$ in spherical coordinates. We shall denote the unit sphere centered at the origin by ${\mathrm{\Omega }}_0=\{\overrightarrow{r}\in {ℝ}^3:|\overrightarrow{r}|=1\}$ and will typically use the symbol $\widehat{\mathbf{s}}=(1,\theta ,\varphi )$ to denote a vector on ${\mathrm{\Omega }}_0$.

Following [8, 9], we define the associated Legendre function $P_n^m(x)$ of degree $n>0$ and order $-n\le m\le n$ and by the formula:

$P_n^m(x)$

$={(-1)}^m{\displaystyle \frac{{(1-x^2)}^{m/2}}{2^{n}n!}}{\displaystyle \frac{d^{n+m}}{d^{n+m}x}}{(x^2-1)}^n.$

(5a)

Following Messiah [9], we define the spherical harmonic $Y_n^m(\theta ,\varphi )$ of degree $n$ and order $m$ by the relation:

$Y_n^m(\theta ,\varphi )$

$={\left[{\displaystyle \frac{2n+1}{4\pi }}\right]}^{1/2}{\left[{\displaystyle \frac{(n-m)!}{(n+m)!}}\right]}^{1/2}{P}_n^m(\cos \theta )e^{im\varphi }.$

(5b)

We have:

	$Y_n^{-m}(\theta ,\varphi )={(-1)}^m{[Y_n^m(\theta ,\varphi )]}^{\ast },$			(5c)
	${\displaystyle {\int }_0^{\pi }}{\displaystyle {\int }_0^{2\pi }}{Y}_n^m(\theta ,\varphi )Y_{n^{\prime }}^{m^{\prime }}{}^{\ast }(\theta ,\varphi )\sin \theta d\theta d\varphi$	$={\delta }_{n{n}^{\prime }}{\delta }_{m{m}^{\prime }},$		(5d)

where ${}^{\ast }$ denotes the complex conjugate.

Let $\overrightarrow{x}=r\widehat{\mathbf{x}}=(r,\theta ,\varphi )\in ℝ$. We define inner and outer multipole functions for the Laplace equation by the formulae [10]:

	$$O_n^m(\overrightarrow{x})=\sqrt{\frac{4\pi }{2n+1}}\frac{{(-1)}^n{i}^{-m}}{A_n^m}\frac{Y_n^m(\theta ,\varphi )}{r^{n+1}},$$		(6a)
	$$I_n^m(\overrightarrow{x})=\sqrt{\frac{4\pi }{2n+1}}{i}^m{A}_n^m{r}^n{Y}_n^m(\theta ,\varphi ),$$		(6b)
	$$\text{where}{A}_n^m=\frac{{(-1)}^n}{\sqrt{(n-m)!(n+m)!}}.$$		(6c)

Similarly, for the Helmholtz equation, we define the inner and outer functions, denoted by ${\tilde{I}}_n^m(\overrightarrow{x})$ and ${\tilde{O}}_n^m(\overrightarrow{x})$, respectively, by the formulae [10]:

	${\tilde{I}}_n^m(\overrightarrow{x})$	$=j_n(kr)Y_n^m(\widehat{\mathbf{x}}),$		(7a)
	${\tilde{O}}_n^m(\overrightarrow{x})$	$=h_n(kr)Y_n^m(\widehat{\mathbf{x}}),$		(7b)
where $j_n(kr)$ and $h_n(kr)$ are, respectively, spherical Bessel and Hankel functions of the first kind and of order $n$.

III. ABSTRACTIONS FOR THE FAR FIELD COMPUTATIONS

We first develop the abstractions needed to evaluate the canonical summation in equation (4). Assuming the existence of an octree structure, let $B_s$ and $B_r$ be the boxes in which the basis function $f_j$ and the testing function $g_i$, respectively, reside¹. Referring to Fig. 2, let:

	$B_s\subset D_s$	$\stackrel{\text{def}}{=}\{\overrightarrow{x}\in {ℝ}^3||\overrightarrow{x}-{\overrightarrow{c}}_s|<R_s\},$		(8)
	$B_r\subset D_r$	$\stackrel{\text{def}}{=}\{\overrightarrow{x}\in {ℝ}^3||\overrightarrow{x}-{\overrightarrow{c}}_r|<R_r\}.$		(9)

We assume that $\mathrm{supp}{f}_j\subset D_s$ and $\mathrm{supp}{g}_i\subset D_r$ and that ${\overline{D}}_s\cap {\overline{D}}_r=\mathrm{\varnothing }$.

Figure 2: The basis function, $f_j$, and the testing function, $g_i$, (“generalized particles”) are enclosed in spheres $D_s$, with center at ${\overrightarrow{c}}_s$ and $D_r$, with center at ${\overrightarrow{c}}_r$, respectively, under the condition that $D_s\cap D_r=\mathrm{\varnothing }$.

To represent the matrix element in equation (4), we first let:

	${\psi }_j(\overrightarrow{x})$	$={\displaystyle {\int }_{\mathrm{supp}{f}_j}}{f}_j({\overrightarrow{x}}^{\prime })G(\overrightarrow{x},{\overrightarrow{x}}^{\prime })\mathit{d}{\overrightarrow{x}}^{\prime },$		(10a)
	$P_{ij}$	$={\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){\psi }_j(\overrightarrow{x})\mathit{d}\overrightarrow{x}.$		(10b)

For a given accuracy $ϵ>0$, let $N\ge 0$ be given. Let ${𝒪}_n^m(\overrightarrow{x})$ and ${ℐ}_n^m(\overrightarrow{x})$, respectively, denote the outer and inner multipoles of either Laplace or Helmholtz kernel. Then, we recall the following foundational results from the theory of MLFMM [2, 10]:

The potential ${\psi }_j(\overrightarrow{x})$ can be approximated by a finite outer multipole series valid outside $D_s$:

${\psi }_j(\overrightarrow{x})$

$\approx {\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=-n}^n}{\alpha }_n^m{𝒪}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_s),\overrightarrow{x}\notin {\overline{D}}_s.$

(11a)

The potential ${\psi }_j(\overrightarrow{x})$ can be approximated as an inner multipole series valid inside $D_r$:

${\psi }_j(\overrightarrow{x})$

$\approx {\displaystyle \sum _{n=0}^{N^{\prime }}}{\displaystyle \sum _{m=-n}^n}{\beta }_n^m{ℐ}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r),\overrightarrow{x}\in D_r.$

(11b)

There exists a linear operator, $\mathrm{𝒯}$, called a translation operator, relating the coefficients $\{{\alpha }_n^m\}$ and $\{{\beta }_n^m\}$ by:

${\beta }_n^m$

$={\displaystyle \sum _{p=0}^N}{\displaystyle \sum _{q=-p}^q}\mathrm{𝒯}\mathrm{\Leftarrow }\mathrm{\backslash }\mathrm{\iff }\mathrm{\Updownarrow }\mathrm{\varnothing }\mathrm{\surd }\mathrm{\iff }\mathrm{\coprod }\mathrm{\Rightarrow }{\alpha }_{\mathrm{\surd }}^{\mathrm{\coprod }}\mathrm{\swarrow }$

(11c)

To evaluate $P_{ij}$, we substitute equation (11b) into equation (10b), and obtain:

	$P_{ij}$	$\approx {\displaystyle \sum _{n=0}^{\min \{N,N^{\prime }\}}}{\displaystyle \sum _{m=-n}^n}{\beta }_n^m{\rho }_n^m,$		(11d)
	${\rho }_n^m$	$\stackrel{\text{def}}{=}{\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){ℐ}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r)\mathit{d}\overrightarrow{x}.$		(11e)

Borrowing the terminology from Laplace FMM parlance, we shall refer to ${\alpha }_n^m$ as the “Q2M” coefficients and ${\rho }_n^m$ as the “L2P” coefficients.

Expressing the Q2M and L2P coefficients in terms of multipole functions presents a challenge to implementation when multiple teams are involved. For example, a team developing a stand-alone library of FMM, to be used by multiple product teams, may not be deeply familiar with the applications themselves. On the other hand, a product team, although experts in, say, integral equations, may have little or no familiarity with FMM. Given that it is these coefficients that couple the application side (integral equation) with the FMM implementation, it is imperative that they are evaluated correctly and efficiently.

An obvious way to do this is to let the FMM library be aware of the basis and testing functions, and thus let the library evaluate these functions. This has the advantage that if the library decides to change, say, the definitions of inner functions, it will remain transparent to the application. However, this approach has the severe disadvantage that the library is deeply tied to a particular class of basis and testing functions. It is then left with the task of potentially adding a new class of basis and testing functions for every new application. This is inefficient and will eventually make the code unmaintainable.

A second approach is to let the application evaluate the coefficients $({\alpha }_n^m,{\rho }_n^m)$ using formulae similar to those in equation (11d) and then hand-over the results to the FMM kernel. This has the great advantage that the FMM library remains completely decoupled from the application. However, this approach has the disadvantage that the application team must have some familiarity with FMM and multipoles. This may not always be the case.

Therefore, our goal is to evaluate the coefficients ${\alpha }_n^m$ and ${\rho }_n^m$ within the library without exposing the functions $({𝒪}_n^m,{ℐ}_n^m)$ to the application.

We now present an approach that eliminates these difficulties and it is inspired by ideas from [11, 12, 13, 14].

A. Evaluation of Q2M and L2P coefficients for the Laplace equation

We begin by noting the basis function $f_j$ is defined by the application and that the application “knows” how to evaluate ${\psi }_j(\overrightarrow{x})$ in equation (10a) at any point $\overrightarrow{x}\in {ℝ}^3$. Moreover, if $\overrightarrow{x}\notin \mathrm{supp}{f}_j$, the integral can be accurately evaluated using simple quadrature.

Consider a sphere $S\stackrel{\text{def}}{=}\{\overrightarrow{x}\in {ℝ}^3||\overrightarrow{x}-{\overrightarrow{c}}_s|=R>R_s\}$ so that $S\cap \mathrm{supp}{f}_j=\mathrm{\varnothing }$. Let $\overrightarrow{x}\in S$. Then $\overrightarrow{x}={\overrightarrow{c}}_s+R\widehat{\mathbf{s}}$ where $\widehat{\mathbf{s}}\in {\mathrm{\Omega }}_0$. Then, substituting the definition of $O_n^m$ from equation (6a) into equation (11a) and using the orthogonality of $Y_n^m(\widehat{\mathbf{s}})$ from equation (5d), we obtain the well known relation (see, for example [14, equations 4-5]):

	${\alpha }_n^m$	$=i^m\sqrt{{\displaystyle \frac{2n+1}{4\pi (n-m)!(n+m)!}}}{R}^{n+1}\cdot$
		${\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\psi }_j({\overrightarrow{c}}_s+R\widehat{\mathbf{s}})Y_n^{m\ast }(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s.$		(12)

Now assume that the degree of approximation used in the FMM is $p>0$. Then, the maximum degree of spherical harmonics under the integral that is relevant to the computations is $2p$. Therefore, if we use a quadrature, such as a composite Gauss-Legendre-Trapezoidal rule, that can integrate spherical harmonics of degree $2p$ exactly, we will be able to compute the coefficients exactly as well. Let $N_q$ be the number of quadrature points, $\{({\widehat{\mathbf{s}}}_k,w_k)\}$, be quadrature points and their corresponding weights on ${\mathrm{\Omega }}_0$, respectively, and let ${\overrightarrow{x}}_k={\overrightarrow{c}}_s+R{\widehat{\mathbf{s}}}_k$, for $1\le k\le N_q$. Then, equation (12) can be approximated as:

	${\tilde{\alpha }}_n^m$	$\approx i^m\sqrt{{\displaystyle \frac{2n+1}{4\pi (n-m)!(n+m)!}}}{R}^{n+1}\cdot$
		${\displaystyle \sum _{k=1}^{N_q}}{w}_k{\psi }_j({\overrightarrow{x}}_k)Y_n^{m\ast }({\widehat{\mathbf{s}}}_k).$		(13)

Therefore, to compute the coefficients ${\alpha }_n^m$, all the library needs are the values of ${\psi }_j(\overrightarrow{x})$, given by equation (10), at a set points on a specified sphere, which the application can readily supply. As remarked earlier, since the observation points ${\overrightarrow{x}}_k$ are chosen to be far away from the support of $f_j$, the integrand in (10a) is non-singular and hence ${\psi }_j({\overrightarrow{x}}_k)$ can be evaluated using simple quadrature rules.

Using the addition theorem [10]:

	${\displaystyle \frac{1}{\parallel \overrightarrow{x}-{\overrightarrow{x}}^{\prime }\parallel }}$	$={\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{\displaystyle \sum _{m=-n}^n}{(-1)}^n{I}_n^{-m}({\overrightarrow{x}}^{\prime }-{\overrightarrow{c}}_s)O_n^m(\overrightarrow{x}-{\overrightarrow{c}}_s),$		(14a)
for $\parallel \overrightarrow{x}-{\overrightarrow{c}}_s\parallel >\parallel {\overrightarrow{x}}^{\prime }-{\overrightarrow{c}}_s\parallel$ in equation (10a) and comparing with equation (11a), it is easy to see that:
	${\alpha }_n^m$	$={(-1)}^n{\displaystyle {\int }_{\mathrm{supp}{f}_j}}{f}_j({\overrightarrow{x}}^{\prime })I_n^{-m}({\overrightarrow{x}}^{\prime }-{\overrightarrow{c}}_s)\mathit{d}{\overrightarrow{x}}^{\prime }.$		(14b)

Comparing equation (14b) with the definition of ${\rho }_n^m$ in equation (11d), we see that:

${\rho }_n^m={(-1)}^n{\alpha }_n^{-m},$

(15)

for the $f_j$ if it were also used as a testing function in a Galerkin MoM. Since:

	${\displaystyle {\int }_{\mathrm{\Omega }}}{g}_m(\overrightarrow{x}){\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}}{f}_n({\overrightarrow{x}}^{\prime })G(\overrightarrow{x},{\overrightarrow{x}}^{\prime })\mathit{d}{\mathrm{\Omega }}^{\prime }\mathit{d}\mathrm{\Omega }$
	$={\displaystyle {\int }_{{\mathrm{\Omega }}^{\prime }}}{f}_n({\overrightarrow{x}}^{\prime }){\displaystyle {\int }_{\mathrm{\Omega }}}{g}_m(\overrightarrow{x})G(\overrightarrow{x},{\overrightarrow{x}}^{\prime })\mathit{d}\mathrm{\Omega }\mathit{d}{\mathrm{\Omega }}^{\prime },$

it follows that we can use exactly the same technique (equations (12) and (13)) for computing ${\rho }_n^m$ of any testing function $g_i$ by evaluating the corresponding Q2M coefficients and using the relation in equation (15).

By this process, we have thus decoupled the application and the library for computing the L2P and Q2M coefficients for the Laplace equation: the application does not have to work with $I_n^m$ functions and the library does not have to know the nature of the basis and testing functions, $(f_j,g_i)$.

A.1. Comparison of performance

Given a multipole degree $p>0$ and using a composite Gauss-Trapezoidal rule with $2{(p+1)}^2$ quadrature points, it is easy to see that the standard approach to computing Q2M (and L2P) coefficients and the potential based approach proposed above have the same asymptotic complexity. However, the proposed method has a larger constant associated with it. We now demonstrate that this additional cost for engineering a simplified interface is negligible in the context of the much larger initialization time (the time for the matrix-vectors are not affected).

Figure 3: Comparing the time for evaluating Q2M coefficients using standard (“Direct Multipole”) and using sampled potentials (“Potential”) for various test cases. We have also plotted the total setup time to provide the larger context.

To do this, we have taken 30 distinct triangular meshes ranging in size from 200 to about 1.2 million. They were generated by Ansys Q3D Extractor for different test cases and do not share a common geometry. For our purposes, it is sufficient to replace the triangles with point sources at their centroids. The canonical matrix in equation (4) is then constructed with a multipole degree of $p=5$ which gives roughly 4 digits of accuracy in the matrix vector products. All the computations were done on a single core of a virtual machine with two Intel(R) Xeon(R) 6238R @2.20GHz processors with 182GB of RAM. In Fig. 3, we compare the time for evaluating the Q2M coefficients using the standard approach and from the potential samples.

It is important to note that the bulk of the initialization/setup time is spent in computing the near matrix and that the near matrix evaluation for concentrated sources gives the least possible time. This is because the commonly used basis functions require more expensive handling of singularities and near singularities. In other words, we are comparing the cost of Q2M evaluations against the best case scenario for the rest of the initialization. It is clear from Fig. 3 that although the proposed method is slightly more expensive, the increase is negligible when compared with the overall setup time. We have also compared the accuracy of the computed matrix vector products and the largest $2-$norm error observed was less than $2\times {10}^{-4}$, which is comparable to the overall precision of the matrix-vector products corresponding to the degree of multipole approximation $p=5$.

B. Evaluation of Q2P and L2P coefficients for the Helmholtz equation

In the case of Helmholtz FMM that uses only multipole expansions involving ${\tilde{O}}_n^m$ and ${\tilde{I}}_n^m$, it is clear that an approach analogous to the one described in section 3.1 can be readily developed. However, most applications of Helmholtz FMM involves electrically large structures where a diagonalized form is employed [2, 10], often combined with multipole expansions to handle multi-scale structures [15]. For such a mixed-form FMM, it is simpler to compute signature functions [10] rather than work with potentials. We now show how signature functions can be used to decouple the FMM library from the applications.

The signature function of the potential ${\psi }_j(\overrightarrow{x})$ of equation (10a), now with the Green’s function for the Helmholtz equation (1), is defined by formula [10]:

	${\tilde{\psi }}_j(\widehat{\mathbf{s}})$	$=\underset{r\to \mathrm{\infty }}{lim}kr{e}^{-ikr}\psi ({\overrightarrow{c}}_s+r\widehat{\mathbf{s}}),$		(16a)
which when approximated by equation (11a) with ${𝒪}_n^m={\tilde{O}}_n^m$, can be shown to be:
	${\tilde{\psi }}_j(\widehat{\mathbf{s}})$	$={\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{\displaystyle \sum _{m=-n}^m}{\alpha }_n^m{\displaystyle \frac{Y_n^m(\widehat{\mathbf{s}})}{i^{n+1}}}.$		(16b)

Therefore, given a signature function $\tilde{\psi }(\widehat{\mathbf{s}})$, the orthogonality of spherical harmonics implies that we can recover the coefficients of the multipole expansion in equation (11a), ${\alpha }_n^m$, using the formula [16]:

${\alpha }_n^m$

$=i^{n+1}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\tilde{\psi }}_j(\widehat{\mathbf{s}})Y_n^{m\ast }(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s.$

(17)

It is well known that [2, 10] the signature function of the Green’s function $G(\overrightarrow{r},{\overrightarrow{r}}^{\prime })=\frac{e^{ik|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}$ for ${\overrightarrow{r}}^{\prime }\in D_s$ is given by:

$$\tilde{G}(\widehat{\mathbf{s}})=k{e}^{ik\widehat{\mathbf{s}}\cdot ({\overrightarrow{c}}_s-{\overrightarrow{r}}^{\prime })}.$$

Therefore, it follows that signature function of ${\psi }_j(\overrightarrow{x})$ of the canonical problem in equation (10a) is given by:

${\tilde{\psi }}_j(\widehat{\mathbf{s}})$

$=k{\displaystyle {\int }_{\mathrm{supp}{f}_j}}{f}_j(\overrightarrow{r^{\prime }})e^{ik\widehat{\mathbf{s}}\cdot ({\overrightarrow{c}}_s-{\overrightarrow{r}}^{\prime })}\mathit{d}{\overrightarrow{r}}^{\prime },$

(18)

and we define the signature of the basis function $f_j$ by ${\tilde{f}}_j(\widehat{\mathbf{s}})={\tilde{\psi }}_j(\widehat{\mathbf{s}})$.

Using the representation [9, 10]:

	${\tilde{I}}_n^m(\overrightarrow{x}-\overrightarrow{c})$	$={\displaystyle \frac{1}{4\pi i^n}}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{e}^{ik(\overrightarrow{x}-\overrightarrow{c})\cdot \widehat{\mathbf{s}}}{Y}_n^m(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,$		(19a)
it is easy to show that the L2P coefficients:
	${\rho }_n^m$	$={\displaystyle \frac{1}{4\pi i^n}}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\tilde{g}}_i(\widehat{\mathbf{s}})Y_n^m(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s$		(19b)
	$\text{where}{\tilde{g}}_i(\widehat{\mathbf{s}})$	$\stackrel{\text{def}}{=}{\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x})e^{ik(\overrightarrow{x}-\overrightarrow{c})\cdot \widehat{\mathbf{s}}}\mathit{d}\overrightarrow{x},$		(19c)
are defined by the signature function of the testing function $g_i$2.

Hence, in a mixed-form broadband FMM that uses both diagonalized and (possibly scaled) multipole expansions, the user only has to provide the signature functions in equation (18) and equation (19c) for setting up the far field calculations in Helmholtz FMM. Furthermore, these functions are easy to define and easy to calculate.

IV. COMPUTATION OF THE DERIVATIVES OF THE FAR FIELD

We claimed in the introduction that it is possible to represent the derivatives of far fields as product of sparse matrices. In this section, we shall demonstrate this claim. The approach relies on the following observation for the Helmholtz equation (A similar argument can be made for the Laplace equation by setting $k=0$ and replacing ${\tilde{I}}_n^m$ with $I_n^m$.)

Consider the potential field ${\psi }_j(\overrightarrow{x})$ due to sources in $D_s$ having an approximation of the form in equation (11b):
	${\psi }_j(\overrightarrow{x})$	$={\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=-n}^n}{\beta }_n^m{\tilde{I}}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r),\overrightarrow{x}\in D_r.$		(20a)

Since:

${\nabla }^2{\displaystyle \frac{\partial {\psi }_j}{\partial u}}+k^2{\displaystyle \frac{\partial {\psi }_j}{\partial u}}$

$={\displaystyle \frac{\partial }{\partial u}}\left\{{\nabla }^2{\psi }_j+k^2{\psi }_j\right\}=0,u=x,y,z,$

(20b)

it follows that $\frac{\partial {\psi }_j}{\partial u}$ also must have an expansion of the form:

${\displaystyle \frac{\partial {\psi }_j(\overrightarrow{x})}{\partial u}}$

$={\displaystyle \sum _{n=0}^{N^{\prime }}}{\displaystyle \sum _{m=-n}^n}{\delta }_n^m{\tilde{I}}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r),\overrightarrow{x}\in D_r,$

(20c)

and that the coefficients $\{{\beta }_n^m\}$ and $\{{\delta }_n^m\}$ must be related through a linear operator $T_u$ so that:

${\psi }_j\mapsto {\displaystyle \frac{\partial {\psi }_j}{\partial u}}$

$⟹\{{\beta }_n^m\}\mapsto \{{\delta }_n^m\}=T_u(\{{\beta }_n^m\}).$

Thus, for non-negative integers $p,q,r\ge 0$, an inductive reasoning shows that:

${\psi }_j\mapsto {\displaystyle \frac{{\partial }^{p+q+r}{\psi }_j}{\partial x^p\partial y^q\partial z^r}}$

$⟹\left\{{\beta }_n^m\right\}\mapsto \left\{{\delta }_n^m\right\}=T_z^r{T}_y^q{T}_x^p\left(\left\{{\beta }_n^m\right\}\right).$

Therefore:

	${\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){\displaystyle \frac{{\partial }^{p+q+r}{\psi }_j}{\partial x^p\partial y^q\partial z^r}}\mathit{d}\overrightarrow{x}$
	$={\displaystyle \sum _{n=0}^{N^{\prime }}}{\displaystyle \sum _{m=-n}^n}{\delta }_n^m\cdot {\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){\tilde{I}}_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r)\mathit{d}\overrightarrow{x}$		(21a)
	$={\displaystyle \sum _{n=0}^{N^{\prime }}}{\displaystyle \sum _{m=-n}^n}{\delta }_n^m{\rho }_n^m,$		(21b)

where ${\rho }_n^m$ are the same receiving coefficients defined in equation (11d). Thus, we see that partial derivatives of arbitrary order manifest as product of linear operators that can be applied to the local expansion coefficients ${\beta }_n^m$ transparent to the application/user.

So far, we have not said anything about the nature of $T_u$. In the next three sections, we shall demonstrate that the finite dimensional representation of $T_u$ is a sparse matrix for the Laplace and low frequency Helmholtz FMM that use explicit multipole expansions and a diagonal matrix for the high frequency Helmholtz FMM that uses signature functions.

A. Derivatives for the Laplace kernel

Let $\overrightarrow{x}=(r,\theta ,\varphi )\in {ℝ}^3$. We begin with the integral representation of the inner function $I_n^m(\overrightarrow{x})$ derived in [10, equation 2.33]:

	$I_n^m(\overrightarrow{x})$	$={\displaystyle \frac{{(-1)}^{n+m}}{n!}}{\displaystyle \frac{1}{2\pi }}\cdot$
		${\displaystyle {\int }_{-\pi }^{\pi }}{\left(z+ix\cos \alpha +iy\sin \alpha \right)}^n{e}^{im\alpha }\mathit{d}\alpha .$		(22)

Without the loss of generality, consider the case of $\frac{\partial I_n^m(\overrightarrow{x})}{\partial x}$. From Appendix 6.1, we have:

${\displaystyle \frac{\partial I_n^m(\overrightarrow{x})}{\partial x}}$

$={\displaystyle \frac{i}{2}}\left(I_{n-1}^{m+1}(\overrightarrow{x})+I_{n-1}^{m-1}(\overrightarrow{x})\right).$

(23)

Substituting equation (23) into equation (20a), we get:

	${\displaystyle \frac{\partial {\psi }_j(\overrightarrow{x})}{\partial x}}$	$={\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=-n}^n}{\beta }_n^m{\displaystyle \frac{\partial I_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r)}{\partial x}}$
		$={\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=-n}^n}{\beta }_n^m\cdot$
		${\displaystyle \frac{i}{2}}\left(I_{n-1}^{m+1}(\overrightarrow{x}-{\overrightarrow{c}}_r)+I_{n-1}^{m-1}(\overrightarrow{x}-{\overrightarrow{c}}_r)\right).$		(24a)

Since $I_n^m=0$ if $n<0$ or $|m|>n$, the above equation can be rewritten as:

${\displaystyle \frac{\partial \psi (\overrightarrow{x})}{\partial x}}$

$={\displaystyle \sum _{n=0}^{N-1}}{\displaystyle \sum _{m=-n}^n}{\displaystyle \frac{i}{2}}({\beta }_{n+1}^{m-1}+{\beta }_{n+1}^{m+1})I_n^m(\overrightarrow{x}-{\overrightarrow{c}}_r).$

(24b)

Identifying:

${\delta }_n^m$

$={\displaystyle \frac{i}{2}}({\beta }_{n+1}^{m-1}+{\beta }_{n+1}^{m+1}),$

and $N^{\prime }=N-1$ in equation (20c), we see that the matrix representation of $T_x$ has at most two non-zero entries per row. It is clear from Appendix 6.1 that a similar reasoning can be made for $T_y$ and $T_z$ as well.

B. Derivatives for the Helmholtz kernel with multipoles

Without the loss of generality, we consider the partial derivative with respect to $x$. Then, it can be shown that $\frac{\partial {\tilde{I}}_n^m}{\partial x}$ can be expressed as a four-term recurrence of the form (see Appendix 6.2 for the derivation and for the explicit expressions):

${\displaystyle \frac{\partial {\tilde{I}}_n^m(\overrightarrow{x})}{\partial x}}$

$={\displaystyle \sum _{\stackrel{\nu =|n-1|}{\nu \ne n}}^{n+1}}{\displaystyle \sum _{p=-1,1}}{T}_{\nu }^{m+p}{\tilde{I}}_{\nu }^{m+p}(\overrightarrow{x}).$

(25a)

A reasoning similar to that made in the previous section shows that each row of the matrix representation of $T_x$ has at most four entries, once again showing that the derivative can be represented as a sparse matrix. Furthermore, a closer look at equation (25a) shows that if equation (20a) has degree $N$, then equation (20c) will have degree $N^{\prime }=N+1$. However, if the degree of the “L2P” coefficients ${\rho }_n^m$ is truncated to $N$, then the summation in equation (21b) will be truncated to degree $N$ and we can set $N^{\prime }=N$.

C. Derivatives for the Helmholtz kernel in diagonalized form

A potential field ${\psi }_j(\overrightarrow{x})$ having an approximation of the form in (20a), can also be represented as [2, 10]:

	${\psi }_j(\overrightarrow{x})$	$={\displaystyle \frac{i}{4\pi }}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\tilde{\psi }}_j(\widehat{\mathbf{s}})\cdot e^{ik(\overrightarrow{x}-{\overrightarrow{c}}_r)\cdot \widehat{\mathbf{s}}}\mathit{d}{\mathrm{\Omega }}_s,$		(26a)
	$\text{where}{\tilde{\psi }}_j(\widehat{\mathbf{s}})$	$\stackrel{\text{def}}{=}{\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=-n}^m}{\beta }_n^m{\displaystyle \frac{Y_n^m(\widehat{\mathbf{s}})}{i^{n+1}}},$		(26b)
so that:
	$P_{ij}$	$={\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){\psi }_j(\overrightarrow{x})\mathit{d}\overrightarrow{x}$
		$={\displaystyle \frac{i}{4\pi }}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\tilde{g}}_i(\widehat{\mathbf{s}})\cdot {\tilde{\psi }}_j(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,$		(26c)
	$\text{where}{\tilde{g}}_i(\widehat{\mathbf{s}})$	$={\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x})e^{ik(\overrightarrow{x}-{\overrightarrow{c}}_r)\cdot \widehat{\mathbf{s}}}\mathit{d}\overrightarrow{x},$		(26d)
and we recognize ${\tilde{g}}_i(\widehat{\mathbf{s}})$ as the signature function defined in equation (19c). Since:
	${\displaystyle \frac{\partial {\psi }_j(\overrightarrow{x})}{\partial x}}$	$={\displaystyle \frac{i}{4\pi }}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}(ik{s}_x)\cdot {\tilde{\psi }}_j(\widehat{\mathbf{s}})\cdot e^{ik(\overrightarrow{x}-{\overrightarrow{c}}_r)\cdot \widehat{\mathbf{s}}}\mathit{d}{\mathrm{\Omega }}_s,$		(26e)
from equation (26a), it follows that:
	${\displaystyle {\int }_{\mathrm{supp}{g}_i}}{g}_i(\overrightarrow{x}){\displaystyle \frac{\partial {\psi }_j(\overrightarrow{x})}{\partial x}}\mathit{d}\overrightarrow{x}$
	$\mathrm{\hspace{1em}\hspace{1em}}={\displaystyle \frac{i}{4\pi }}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{\tilde{g}}_i(\widehat{\mathbf{s}})\cdot (ik{s}_x)\cdot {\tilde{\psi }}_j(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,$			(26f)
demonstrating that the partial differential operator appears as a diagonal modification to the incoming signature function ${\tilde{\psi }}_j(\widehat{\mathbf{s}})$.

Note that this result has been known from the very early days of FMM and has been effectively used for computing with derivatives [17]. The novelty in our approach is that we are using it to transfer the responsibility of computing the derivative to the library instead of requiring the application do it.

V. A PSEUDO-EXAMPLE FOR A LIBRARY INTERFACE

We now discuss how the results from sections 3 and 4 can be used to develop an application independent library of Laplace and Helmholtz FMM. Our objective is not to describe a concrete implementation, but to discuss the essential functionality.

To that end, we divide the library functionality into two parts: the matrix build/initialization functions and matrix apply (“matrix-vector” product) functions. During the initialization, we expect the user to provide what are called “call-back” functions. The library, in turn, repeatedly calls these functions to get the information needed to complete matrix initialization. After the matrix is built, the library provides the functionality to evaluate matrix-vector products, potentially involving different combinations of derivatives.

A. Matrix initialization

Building the matrix consists of three phases:

The tree construction phase

The evaluation of the near interactions

The evaluation of the far field representations (“Q2M” and “L2P” coefficients).

A.1. Constructing the tree

Rationale:

In general, each source/receiver can be assigned a “center” and a spatial extent. The former can be captured as a array of three floating point numbers. Given that the mesh used in the discretization can be extremely non-uniform in practical cases, it is often necessary to take the spatial extent into account to ensure sufficient accuracy. This can lead to non-uniform trees with internal nodes also “hosting” sources and/or receivers. The following listing shows an example of an interface.

A.2. Computing near interactions with multiple matrices

Rationale:

With the ability to compute arbitrary derivatives of the far field, it becomes necessary to handle multiple matrices as well. As a simple example, consider the following summations:

	$f({\overrightarrow{x}}_i)$	$={\displaystyle \sum _{j=1}^N}{a}_{j}G({\overrightarrow{x}}_i,{\overrightarrow{x}}_j),$		(27a)
	$\overrightarrow{F}({\overrightarrow{x}}_i)$	$={\displaystyle \sum _{j=1}^N}{a}_j\nabla G({\overrightarrow{x}}_i,{\overrightarrow{x}}_j).$		(27b)

Although potential and the derivatives of the far field can be computed from a single FMM representation, the four components of the near matrix must be stored simultaneously. In other words, the library needs to allow for the possibility of multiple matrix entries resulting from each pair-wise interaction. If we assume that every pair-wise near interaction results in the same number of matrix entries (although there are exceptions, this usually the case), the following interface would be adequate for most purposes.

In the above example, dim(elements)=4 and a user can return the Green’s function and its gradient as:

$\mathtt{elements}$

$=[G,{\displaystyle \frac{\partial G}{\partial x}},{\displaystyle \frac{\partial G}{\partial y}},{\displaystyle \frac{\partial G}{\partial z}}].$

A.3. Computing Q2M and L2P coefficients

Rationale:

In section 3, we have shown that the required Q2M and L2P coefficients can be evaluated from sampled potential values or signature functions. For the Laplace equation, this can be implemented in the following way. Given a node with sources (or receivers) in it, we can construct a sphere with center, ${\overrightarrow{c}}_s$ at the center of the box and a sufficiently large radius, $R\gg d$ where $d$ is the box length, and introduce a Gauss-Trapezoidal rule as discussed in section 3.1. Then, a function of the following nature can be invoked repeatedly to evaluate the potential, as defined by equation (10a) at each integration point ${\overrightarrow{x}}_k={\overrightarrow{c}}_s+R{\widehat{\mathbf{s}}}_k$ in equation (13).

For the Helmholtz equation, we expect the user to return the signature function, defined by equation (18).

Note that in this case, we need both the center of the box ${\overrightarrow{c}}_s$ and the direction vector ${\widehat{\mathbf{s}}}_k$ passed to the application.

B. Evaluating matrix-vector products

In order to offer maximum flexibility, we split the evaluation of the matrix-vector products (MVP) into three steps:

multipole_pass: performs all the source-to-multipole (Q2M), outer-to-outer (M2M, aggregation), outer-to-inner (M2L) and inner-to-inner (L2L) translations in the standard FMM flow. This is a pre-requisite for calling the far field evaluation step below.

far_axpy: evaluates (the integral of) the potential (or its derivative) for every testing function using the results in Section 4.

near_axpy: evaluates the contribution due the near interaction matrix selected by the user (if there are multiple matrices involved).

The library can offer a functionality similar to the following one to accomplish these tasks:

where the array deriv specifies the orders of the partial derivatives. For example:

$\mathtt{deriv}[1,2,2]⟹{\displaystyle \frac{{\partial }^2}{\partial z^2}}{\displaystyle \frac{{\partial }^2}{\partial y^2}}{\displaystyle \frac{\partial }{\partial x}}.$

The motivation for this is best illustrated by an example. Towards this, consider the evaluation of the two sums in equation (27). A psuedo-code is given below.

Computation of the partial derivatives with respect $y$ and $z$ can be done analogously. Note that the multipole_pass is usually the most dominant part during a matrix-vector evaluation and it is invoked only once.

VI. SUMMARY AND CONCLUSIONS

In this paper, we have proposed an approach to developing a stand-alone, application independent library of Laplace and Helmholtz FMM. Towards this end, we have demonstrated a technique to capture the essential characteristics of the problem needed for setting up the FMM – the basis and testing functions, mesh etc. – using samples of either potentials or the signature functions of a canonical problem. Furthermore, we have shown that it is possible to transfer the responsibility of computing the partial derivatives of the far field to the library, instead of burdening the user to implement complicated integrals. This is accomplished by representing partial derivatives of arbitrary orders as product of certain sparse matrices. We have also outlined a pseudo interface that illustrates how the methods can be used.

Within Ansys Inc., we have developed a flexible library, called Q3Dfastmvpack incorporating these ideas and it is currently being used by multiple Ansys products.

ACKNOWLEDGEMENTS

I would like to thank my colleagues Andy Mathis, Dalian Zheng, Eric Bracken and Indranil Chowdhury for the many discussions. I would also like to thank Werner Thiel for supporting this work.

APPENDICES

A. Partial derivatives of $I_n^m$

Let $h(\alpha ;\overrightarrow{x})=(z+ix\cos \alpha +iy\sin \alpha )$ in equation (22). Consider:

$$\frac{\partial I_n^m(\overrightarrow{x})}{\partial x}=\frac{{(-1)}^{n+m}}{(n-1)!}\frac{i}{2\pi }{\int }_{-\pi }^{\pi }h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{im\alpha }\cos \alpha d\alpha$$

	$={\displaystyle \frac{{(-1)}^{n+m}}{(n-1)!}}{\displaystyle \frac{i}{2\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{im\alpha }{\displaystyle \frac{e^{i\alpha }+e^{-i\alpha }}{2}}\mathit{d}\alpha$
	$={\displaystyle \frac{i}{2}}{\displaystyle \frac{{(-1)}^{n+m}}{(n-1)!}}{\displaystyle \frac{1}{2\pi }}\{{\displaystyle {\int }_{-\pi }^{\pi }}h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{i(m+1)\alpha }d\alpha$
	$+{\displaystyle {\int }_{-\pi }^{\pi }}h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{i(m-1)\alpha }d\alpha \}$
	$={\displaystyle \frac{i}{2}}{\displaystyle \frac{{(-1)}^{n+m}}{(n-1)!}}{\displaystyle \frac{1}{2\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{i(m+1)\alpha }\mathit{d}\alpha$
	$\mathrm{\hspace{1em}}+{\displaystyle \frac{i}{2}}{\displaystyle \frac{{(-1)}^{n+m}}{(n-1)!}}{\displaystyle \frac{1}{2\pi }}{\displaystyle {\int }_{-\pi }^{\pi }}h{(\alpha ;\overrightarrow{x})}^{n-1}{e}^{i(m-1)\alpha }\mathit{d}\alpha$
	$={\displaystyle \frac{i}{2}}\left(I_{n-1}^{m+1}(\overrightarrow{x})+I_{n-1}^{m-1}(\overrightarrow{x})\right),$

where we have used the facts $n-1+m-1=n+m$ and that ${(-1)}^{n+m-2}={(-1)}^{n+m}$. The following results can be similarly derived:

	${\displaystyle \frac{\partial I_n^m(\overrightarrow{x})}{\partial y}}$	$={\displaystyle \frac{1}{2}}\left(I_{n-1}^{m+1}(\overrightarrow{x})-I_{n-1}^{m-1}(\overrightarrow{x})\right),$
	${\displaystyle \frac{\partial I_n^m(\overrightarrow{x})}{\partial z}}$	$=-I_{n-1}^m(\overrightarrow{x}).$

B. Partial derivatives of ${\tilde{I}}_n^m(\overrightarrow{x})$

We have, from equation (19a):

	${\displaystyle \frac{{\partial }^{p+q+r}}{\partial x^p\partial y^q\partial z^r}}{I}_n^m(\overrightarrow{x})$
	$={\displaystyle \frac{{(ik)}^{p+q+r}}{4\pi i^n}}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{s}_x^p{s}_y^q{s}_z^r{e}^{ik\overrightarrow{x}\cdot \widehat{\mathbf{s}}}{Y}_n^m(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,$		(28)

where $\widehat{\mathbf{s}}=(s_x,s_y,s_z)=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta )$. Noting that:

$\sin \theta e^{i\varphi }$

$=(-)\sqrt{{\displaystyle \frac{8\pi }{3}}}{Y}_1^1(\widehat{\mathbf{s}}),$

(29)

from the definition of spherical harmonics $Y_n^m(\theta ,\varphi )$ in equation (5b), and using equation (5c), we get:

	$\sin \theta \cos \varphi$	$=(-)\sqrt{{\displaystyle \frac{2\pi }{3}}}\left[Y_1^1(\widehat{\mathbf{s}})-Y_1^{-1}(\widehat{\mathbf{s}})\right]$		(30a)
	$\sin \theta \sin \varphi$	$=i\sqrt{{\displaystyle \frac{2\pi }{3}}}\left[Y_1^1(\widehat{\mathbf{s}})+Y_1^{-1}(\widehat{\mathbf{s}})\right]$		(30b)
	$\cos \theta$	$=\sqrt{{\displaystyle \frac{4\pi }{3}}}{Y}_1^0(\widehat{\mathbf{s}}).$		(30c)

Let:

${\mathrm{𝒱}}_{\mathrm{\surd }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }$

$={\displaystyle \frac{1}{4\pi i^n}}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{Y}_1^p(\widehat{\mathbf{s}})e^{ik\overrightarrow{x}\cdot \widehat{\mathbf{s}}}{Y}_n^m(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,p=-1,0,1.$

(31)

Then, substituting equation (30) into equation (28) and using equation (31), we get:

	$$\frac{\partial I_n^m(\overrightarrow{x})}{\partial x}=(-)ik\sqrt{\frac{2\pi }{3}}({\mathrm{𝒱}}_{\mathrm{\infty }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }\mathrm{\nwarrow }{\mathrm{𝒱}}_{\mathrm{\nwarrow }\mathrm{\infty }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }\mathrm{\Rightarrow }$$		(32a)
	$$\frac{\partial I_n^m(\overrightarrow{x})}{\partial y}=(-)k\sqrt{\frac{2\pi }{3}}({\mathrm{𝒱}}_{\mathrm{\infty }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }\mathrm{\Downarrow }{\mathrm{𝒱}}_{\mathrm{\nwarrow }\mathrm{\infty }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }\mathrm{\Rightarrow }$$		(32b)
	$$\frac{\partial I_n^m(\overrightarrow{x})}{\partial z}=ik\sqrt{\frac{4\pi }{3}}{\mathrm{𝒱}}_{\mathrm{\prime }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }\mathrm{\swarrow }$$		(32c)

Now, substituting the following well known addition theorem [9, 10]:

$e^{ik\overrightarrow{x}\cdot \widehat{\mathbf{s}}}$

$={\displaystyle \sum _{n^{\prime }=0}^{\mathrm{\infty }}}{\displaystyle \sum _{m^{\prime }=-n^{\prime }}^{n^{\prime }}}4\pi i^{n^{\prime }}{j}_{n^{\prime }}(kx)Y_{n^{\prime }}^{m^{\prime }\ast }(\widehat{\mathbf{x}})Y_{n^{\prime }}^{m^{\prime }}(\widehat{\mathbf{s}}),$

into equation (31), we can express:

	${\mathrm{𝒱}}_{\mathrm{\surd }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }$	$={\displaystyle \frac{1}{4\pi i^n}}{\displaystyle \sum _{\nu =0}^{\mathrm{\infty }}}{\displaystyle \sum _{\mu =-\nu }^{\nu }}4\pi i^{\nu }{(-1)}^{\mu }{j}_{\nu }(k\parallel \overrightarrow{x}\parallel )Y_{\nu }^{-\mu }(\widehat{\mathbf{x}})\cdot$
		${\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{Y}_1^p(\widehat{\mathbf{s}})Y_{\nu }^{\mu }(\widehat{\mathbf{s}})Y_n(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s$		(33a)
	$$=\sum _{\nu =0}^{\mathrm{\infty }}\sum _{\mu =-\nu }^{\nu }{i}^{\nu -n}{(-1)}^{\mu }{j}_{\nu }(k\parallel \overrightarrow{x}\parallel )Y_{\nu }^{-\mu }(\widehat{\mathbf{x}})\mathit{\gamma }\left(\begin{array}{ccc}\hfill \nu \hfill & \hfill 1\hfill & \hfill n\hfill \\ \hfill \mu \hfill & \hfill p\hfill & \hfill m\hfill \end{array}\right),$$			(33b)
where $\mathit{\gamma }$ is the Gaunt coefficient defined as
	$\mathit{\gamma }\left(\begin{array}{ccc}\hfill n_1\hfill & \hfill n_2\hfill & \hfill n_3\hfill \\ \hfill m_1\hfill & \hfill m_2\hfill & \hfill m_3\hfill \end{array}\right)$	$\stackrel{\text{def}}{=}{\displaystyle {\int }_{{\mathrm{\Omega }}_0}}{Y}_{n_1}^{m_1}(\widehat{\mathbf{s}})Y_{n_2}^{m_2}(\widehat{\mathbf{s}})Y_{n_3}^{m_3}(\widehat{\mathbf{s}})\mathit{d}{\mathrm{\Omega }}_s,$		(33c)
for integers $n_i\ge 0$ and $-n_i\le m_i\le n_i$ for $i=1,2,3$. Using the properties of Gaunt coefficients [9], for $\mathit{\gamma }\ne 0$, we have to have $|n-1|\le \nu \le n+1$, $\nu +n+1$ even and $\mu +m+p=0$. Thus, $\nu =|n-1|$ and $\nu =n+1$ and therefore:
	${\mathrm{𝒱}}_{\mathrm{\surd }}\mathrm{\Leftarrow }\overrightarrow{\mathrm{\S }}\mathrm{\Rightarrow }$	$={\displaystyle \sum _{\stackrel{\nu =|n-1|}{\nu \ne n}}^{n+1}}{i}^{\nu -n}{(-1)}^{m+p}{j}_{\nu }(k\parallel \overrightarrow{x}\parallel )Y_{\nu }^{m+p}(\widehat{\mathbf{x}})$
		$\mathit{\gamma }\left(\begin{array}{ccc}\hfill \nu \hfill & \hfill 1\hfill & \hfill n\hfill \\ \hfill -(m+p)\hfill & \hfill p\hfill & \hfill m\hfill \end{array}\right)$		(33d)
		$={(-1)}^{m+p}{\displaystyle \sum _{\stackrel{\nu =|n-1|}{\nu \ne n}}^{n+1}}{i}^{\nu -n}\cdot$
		$\mathit{\gamma }\left(\begin{array}{ccc}\hfill \nu \hfill & \hfill 1\hfill & \hfill n\hfill \\ \hfill -(m+p)\hfill & \hfill p\hfill & \hfill m\hfill \end{array}\right)I_{\nu }^{m+p}(\overrightarrow{x}),$		(33e)
showing that the right hand sides of equations (32a)-(32c) have at most four non-zero terms.

REFERENCES

[1] L. Greengard, “The rapid evaluation of potential fields in particle systems,” Ph.D. thesis, Yale University, 1987.

[2] V. Rokhlin, “Diagonal forms of translation operators for the Helmholtz equation in three dimensions,” Applied and Comp. Harmonic Analysis, vol. 1, pp. 82-93, 1993.

[3] W. C. Chew, J.-M. Jin, E. Michielssen, and J. M. Song, Fast and Efficient Algorithms in Computational Electromagnetics, chap. 9, Boston: Artech House, 2001.

[4] V. Rokhlin, “Rapid solution of integral equations of classical potential theory,” Journal Of Computational Physics, vol. 60, pp. 187-207, 1985.

[5] J. Carrier, L. Greengard, and V. Rokhlin, “A fast adaptive multipole algorithm for particle simulations,” SIAM Journal on Scientific and Statistical Computing, vol. 9, no. 4, pp. 669-686, July 1988.

[6] V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” Journal Of Computational Physics, vol. 86, pp. 414-439,1990.

[7] W.-J. He, X.-W. Huang, M.-L. Yang, and X.-Q. Sheng, “Massively parallel multilevel fast multipole algorithm for extremely large-scale electromagnetic simulations: A review,” Progress In Electromagnetics Research, vol. 173, pp. 37-52,2022.

[8] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions, Cambridge: Cambridge University Press, Cambridge, 2010.

[9] A. Messiah, Quantum Mechanics, New York: Dover Publications, Inc., 2014.

[10] M. A. Epton and B. Dembart, “Multipole translation theory for the three dimensional Laplace and Helmholtz equations,” SIAM J. of Scientific Computing, vol. 16, no. 4, pp. 865-897, July 1995.

[11] C. R. Anderson, “An implementation of the fast multipole method without multipoles,” SIAM Journal of Scientific and Statistical Computing, vol. 13, no. 4, pp. 923-947, July 1992.

[12] T. Eibert, “A diagonalized multilevel fast multipole method with spherical harmonics expansion of the k-space Integrals,” IEEE Transactions on Antennas and Propagation, vol. 53, no. 2, pp. 814-817, Feb. 2005.

[13] B. Dembart and E. Yip, “A 3D Fast multipole method for electromagnetics with multiple levels,” 11th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, vol. 1, pp. 621-628, Mar. 1995.

[14] J. Phillips and J. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 16, no. 10, pp. 1059-1072, 1997.

[15] J.-S. Zhao and W. C. Chew, “Integral equation solution of Maxwell's equations from zero frequency to microwave frequencies,” IEEE Transactions on Antennas and Propagation, vol. 48, no. 10, pp. 1635-1645, Oct. 2000.

[16] B. Dembart and E. Yip, “The accuracy of fast multipole methods for Maxwell’s equations,” IEEE Computational Science and Engineering, vol. 5, no. 3, pp. 48-56, 1998.

[17] J. M. Song and W. C. Chew, “Multilevel fast multipole algorithm for solving combined field integral equations of electromagnetic scattering,” Micro. Opt. Tech. Lett., vol. 10, no. 1, pp. 14-19, Sep. 1995.

FOOTNOTES

¹Note that we are not providing a strict definition for the term “reside.” It is simply a rule that the user provides to determine if a given basis/testing function can be considered as belonging to a given box. For example, in the case of point particles, if the location of the particle falls into a box, it can be considered to be a resident of that box.

²Compare with equation (18) for subtle differences.

BIOGRAPHIES

Sanjay Velamparam I obtained a Ph.D. in Engineering from the Indian Institute of Science, Bangalore, in 1997. My doctoral work was in the then nascent field of fast multipole methods (FMM). Subsequently, working with Dr. Weng Chew at the University of Illinois at Urbana-Champaign, I pioneered the parallelization of Helmholtz FMM on distributed memory computers.

I was with Ansoft Corporation (currently Ansys Inc.) from 2002-2005 working on fast integral equation methods for signal integrity applications. In 2007, I joined Acceleware Corporation, Calgary, and worked on linear algebraic algorithms for GPU computing. In 2010, I joined Apache Design Solutions, which was later acquired by Ansys. I have been working on fast integral equation methods for several Ansys products such as Ansys Q3D Extractor${}^{\text{®}}$ and SIwave™ since 2010.

I enjoy long bicycle rides, especially mountainous routes, running, and generally love the outdoors. I am married and we are fortunate enough to live in the foothills of the majestic Rocky Mountains in Northern Colorado.

ACES JOURNAL, Vol. 39, No. 5, 405–415
doi: 10.13052/2024.ACES.J.390504
© 2024 River Publishers