A Lightweight, Geometrically Flexible Fast Algorithm for the Evaluation of Layer and Volume Potentials

Abstract

Over the last two decades, several fast, robust, and high-order accurate methods have been developed for solving the Poisson equation in complicated geometry using potential theory. In this approach, rather than discretizing the partial differential equation itself, one first evaluates a volume integral to account for the source distribution within the domain, followed by solving a boundary integral equation to impose the specified boundary conditions. Here, we present a new fast algorithm which is easy to implement and compatible with virtually any discretization technique, including unstructured domain triangulations, such as those used in standard finite element or finite volume methods. Our approach combines earlier work on potential theory for the heat equation, asymptotic analysis, the nonuniform fast Fourier transform (NUFFT), and the dual-space multilevel kernel-splitting (DMK) framework. It is insensitive to flaws in the triangulation, permitting not just nonconforming elements, but arbitrary aspect ratio triangles, gaps and various other degeneracies. On a single CPU core, the scheme computes the solution at a rate comparable to that of the fast Fourier transform (FFT) in work per gridpoint.

Figures (8)

• Figure 1: (a) A portion of a domain composed of triangular elements. The triangle $T_j$ is a boundary triangle with curved boundary segment $\Gamma_j$ described by at least three points (indicated by blue circles). The green squares in $T_j$ indicate the locations of quadrature nodes for an integration rule over the triangle. For an interior point $\bm{x}$, we denote by $\bm{b}$ the closest point on the boundary. In our asymptotic formulas, we re-express the curve and layer potential densities in terms of the parameter $\xi$, the excursion in the tangent direction $\bm{t}_b$. (b) Given the points defining a curved boundary segment (or chunk) $\Gamma_j$, we reparametrize $\Gamma_j$ in terms of chordal arclength $l$$-$ that is, the distance along the chord connecting the first endpoint $P_1$ with the second endpoint $P_2$. $h(l)$ is the excursion of a boundary node $Q_i$ in the direction of the outward normal to the chord. When we wish to refine/oversample the user-provided chunk, we do this in the chordal parametrization by interpolating $h(l)$. The additional work is negligible and the user interface is unaffected.
• Figure 1: Relative errors as functions of $\delta$ for $\mathcal{S}_L[\sigma]$ (left), $\mathcal{D}_L[\mu]$ (center) and $\mathcal{V}_L[f]$ (right) using asymptotic expansions of order $O\!\left(\delta^{p/2}\right)$ with $p = 4,5$. Solid lines with solid markers are the relative $\ell_{\infty}$ errors, dotted lines with open markers are the relative $\ell_{2}$ errors. Square markers ($\mathsmaller{\blacksquare}$, $\mathsmaller{\square}$) correspond to $p=4$, while round markers ($\bullet$, $\circ$) correspond to $p=5$. The dashed and dot-dash lines, without markers, are reference curves with convergence rate $\delta^{p/2}$ for $p=4, 5$, respectively.
• Figure 2: A volumetric mesh for the domain $\Omega$ in \ref{['ss:exintneu']} (left). There are $1142$ triangles in the discretization, with $h\approx 0.1$ and $208$ nodes along the boundary $\partial\Omega$, marked with black dots. The solution to the interior Neumann problem \ref{['eq:neusolana']}, with $u_{c} = 0$ is plotted on the right.
• Figure 3: For the Neumann problem in \ref{['ss:exintneu']}, the left-hand plot shows the relative $\ell_{\infty}$ error (solid line with solid markers) and the relative $\ell_{2}$ error (dotted line with open markers) using our solver. The dashed line is for reference, with slope corresponding to a convergence of order $O\!\left(\Delta x ^{3}\right)$. The right-hand plot shows the total time for all steps (excluding mesh generation) for solving and evaluating the solution at all volumetric quadrature nodes.
• Figure 4: A volumetric mesh for the domain $D$ (left) in \ref{['ss:exextdir']}, with $36318$ elements and $h\approx 0.1$. There are $228$ nodes along the boundary $\partial\Omega$, marked with black dots. The right-hand figure shows the modulus of the source density \ref{['eq:exextsourcedens']} on a logartihmic scale (base 10).

Theorems & Definitions (20)

• Remark 1.1
• Remark 2.1
• Definition 2.2
• Lemma 2.3
• Proof 1
• Lemma 2.4
• Proof 2
• Lemma 2.5
• Lemma 2.6
• Lemma 2.7