High-order Finite Element–Integral Equation Coupling on Embedded Meshes

To fix notation and for the benefit of the reader unfamiliar with boundary integral equation methods, we briefly summarize the approach taken by this family of methods when solving boundary value problems for linear, homogeneous, constant-coefficient partial differential equations.

Let $\Gamma\subset\mathbb{R}^{d}$ ($d=2,3$) be a smooth, rectifiable curve. For a scalar linear partial differential operator $\mathcal{L}$ with associated free-space Green’s function $G(x,x_{0})$, the single-layer and double-layer potential operators on a density function $\gamma$ are defined as

Here $\nabla_{x_{0}}$ denotes the gradient with respect to the variable of integration and $\hat{n}(x_{0})$ is the outward-facing normal vector. In addition, the normal derivatives of the layer potentials are denoted

respectively. For the Laplacian $\triangle$ in two dimensions, $G(x,x_{0})=-{(2\pi)}^{-1}\log(x-x_{0})$. As $\mathcal{D}\gamma(x)$ and $\mathcal{S}^{\prime}\gamma(x)$ are discontinuous across the boundary $\partial\Omega$, we use the notation $\bar{\mathcal{D}}\gamma(x)$ and $\bar{\mathcal{S}}^{\prime}\gamma(x)$ to denote the principal value of these operators for target points $x\in\partial\Omega$. Consider, as a specific example, the exterior Neumann problem in two dimensions for the Laplace equation:

where $\lim_{y\to x_{+}}$ denotes a limit approaching the boundary from the exterior of $\Omega$. By choosing the integration surface $\Gamma$ in the layer potential as $\partial\Omega$ and representing the solution $u$ in terms of a single layer potential $u(x)=\mathcal{S}\gamma(x)$ with an unknown density function $\gamma$, we obtain that the Laplace PDE and the far-field boundary condition are satisfied by $u$. The remaining Neumann boundary condition becomes, by way of the well-known jump relations for layer potentials (see [27]) the integral equation of the second kind

The boundary $\Gamma$ and density $\gamma$ may then be discretized and, using the action of $\mathcal{\bar{S}}^{\prime}$ and, e.g. an iterative solver, solved for the unknown density $\gamma$. Once $\gamma$ is known, the representation of $u$ in terms of the single-layer potential (2a) may be evaluated anywhere in $\mathbb{R}^{d}\setminus\Omega$ to obtain the sought solution $u$ of the boundary value problem.

In composing a representation of the solution to the interface problem out of single- and double-layer potentials, the objective is to obtain an integral equation of the second kind — i.e., of the form $(I+\mathcal{A})\gamma=g$, where $\mathcal{A}$ is a compact integral operator. The single- and double-layer operators, as operators on the boundary, are indeed compact. This typically results in benign conditioning that is independent of mesh size and allows the application of a Nyström discretization [27].

First, we solve the FE problem, which, in the form of (4), is not coupled to the IE part. The weak form of (4) requires finding a

where $\mathcal{F}$ and $\mathcal{M}$ are defined for any $v\in H^{1}_{0}(\hat{\Omega})$ as

For the discrete FE solution of the continuous weak problem (5), we consider a Galerkin formulation with $u_{1}$, $v\in V^{h}\subset H^{1}_{0}(\hat{\Omega})$.

We represent $u_{2}=\mathcal{D}\gamma$ in terms of the double layer potential of an unknown density $\gamma$.

The resulting integral equation problem in (4) is

This operator is known to have a trivial nullspace and thus we are guaranteed existence and uniqueness of a density solution $\gamma\in C^{0}(\partial\Omega)$ for $g-u_{1}\rvert_{\partial\Omega}\in C^{0}(\partial\Omega)$ by the Fredholm alternative for a sufficiently smooth curve $\partial\Omega$ [27].

For concreteness, we next discuss the Nyström discretization for this problem type, omitting analogous detail for subsequent problem types. We fix a family of composite quadrature rules with weights $w_{h,i}$ and nodes $\xi_{h,i}\subset\partial\Omega$ parametrized by the element size $h$ so that

for a kernel $K^{\mathcal{A}}$ associated with an integral operator $\mathcal{A}$ and densities $\gamma\in C^{q}(\partial\Omega)$. (We use this notation throughout, e.g. we use $K^{\mathcal{D}}$ to denote the kernel of the double layer potential $\mathcal{D}$.) We let

for general layer potential operators $\mathcal{A}$. Using a conventional Nyström discretization, the unknown discretized density

satisfies the linear system given by

Once $\gamma^{h}$ is known, the solution $u_{2}$ can be computed as

We note that, although the density is numerically represented only in terms of pointwise degrees of freedom, $\gamma^{h}$ can be extended to a function defined everywhere on $\partial\Omega$ by making use of the fact that it solves an integral equation of the second kind (6), yielding

while, on account of (8), agreeing with the prior definition of $\gamma^{h}$.

In this section we establish a decoupling estimate that allows us to express the error in the overall solution in terms of the errors achieved by the IE and FE solutions to their associated sub-problems given in (4). For notational convenience, we introduce $r_{1}=u_{1}\rvert_{\partial\Omega}$ for the restriction of $u_{1}$ to the boundary of $\Omega$ and $r_{1}^{h}=u_{1}^{h}\rvert_{\partial\Omega}$ for the restriction of the approximate solution $u_{1}^{h}$ in the finite element subspace $V^{h}\subset H^{1}_{0}(\hat{\Omega})$.

The purpose of Lemma 1 is to reduce the error encountered in the coupled problem to a sum of errors of boundary value problems each solved by a single, uncoupled method, so that standard FEM and IE error analysis techniques apply to each part.

Analogous bounds can be derived in Sobolev spaces, specifically the $H^{1/2}(\partial\Omega)$ norm on the boundary and the $H^{1}(\Omega)$, which in turn can be related to $L^{2}$ the norms included in the results of the numerical experiments described below.

The method we develop in this paper is a combined FE-IE solver that improves on the approach in [22] in a number of ways. First, we make use of a representation of the solution that gives rise to an integral equation of the second kind, leading to improved conditioning and the applicability of the Nyström method. Second, we use high-order accurate quadrature for the evaluation of the layer potentials, leading to improved accuracy. The earlier method [22] uses a coordinate transformation to remove the singularity and employs adaptive quadrature for points near the singularity. We make use of quadrature by expansion (QBX) [32] using the pytential [33] library. QBX evaluates layer potentials on and off the source surface by exploiting the smoothness of the potential to be evaluated. It forms local/Taylor expansions of the potential off the source surface using the fact that the coefficient integrals are non-singular. Compared with the classical singularity removal method based on polar coordinates, QBX is more general in terms of the kernels it can handle, and it unifies on- and off-surface evaluation of layer potentials. It is also naturally amenable to acceleration via fast algorithms [34]. The finite element terms are evaluated in standard $Q^{n}$ spaces using the deal.II library [35, 36].

We consider three different right-hand sides with varying smoothness: a constant function, a $C^{0}$ piecewise bilinear function, and a piecewise constant function.

These cases are selected to expose different levels of regularity in the problem. A classical $C^{2}$ solution is expected from $f_{\text{c}}$, while $f_{\text{bl}}$ and $f_{\text{pw}}$ admit solutions only in $H^{3}$ and $H^{2}$, respectively.

All problems are defined on the domain $\Omega=\{x:|x|_{2}\leq 0.5\}$, a circle of radius $0.5$. In addition, the domain is embedded in a square domain $\hat{\Omega}=[-0.6,-0.6]\times[0.6,0.6]$, as illustrated in Figure 6 along with the solution obtained when using the right-hand side $f_{\text{c}}$.

Table 1 reports the self-convergence error in the finite element and integral equation portions of the solution for each test case compared to a fine-grid solution whose parameters are given. We see that the method exhibits the expected order of accuracy given the smoothness of the data. In particular, the method is high-order even near the embedded boundary, in contrast with standard immersed boundary methods. The degree of the FE polynomials space $p$ and the truncation order $p_{\text{QBX}}$ in Quadrature by Expansion are chosen so as to yield equivalent orders of accuracy in the solution–for instance $p=1$ and $p_{\text{QBX}}=2$ yield second-order accurate approximations [37, 38]. In the lower-smoothness test cases, we note a marked difference between the error in the $\infty$- and the 2-norm of the error, both shown. Also, for linear basis functions, the finite element convergence rate in $\|\cdot\|_{\infty;\Omega}$ is suboptimal. This matches analytical expectations as known error estimates of the error in this norm for this basis include a factor of $\log(1/h)$ [39, 40].

We solve (3) as before by representing

and posing a system of BVPs for $u_{1}$ and $u_{2}$:

with a given constant $A$. In two dimensions, (18) includes a far-field boundary condition for $u_{2}$ that differs from the standard far-field BC for the exterior Dirichlet problem, $u(x)=O(1)$ as $x\to\infty$ [27]. There are two reasons for this modification. First, the BVP (3) allows solutions containing a logarithmic singularity within $\Omega$. Without permitting logarithmic blow-up at infinity, such solutions would be ruled out by the splitting (18). Neither $u_{1}$ (nonsingular throughout $\Omega$) nor $u_{2}$ would be able to represent them. Second, allowing nonzero additive constants in $u_{2}$ would lead to a non-uniqueness, since for any given constant $C$, $(u_{1}^{\text{new}},u_{2}^{\text{new}})=(u_{1}+C,u_{2}-C)$ would likewise be an admissible solution.

Next, we support the existence of a solution of the coupled BVPs (18) with the stricter-than-conventional decay condition $o(1)$, we let $A=0$ without loss of generality. We remind the reader that any solution of the exterior Dirichlet problem in two dimensions may be represented as $\mathcal{D}\gamma+C$, for some constant $C$ [27, Thm. 6.25]. Since $(\mathcal{D}\gamma)(x)=O(|x|^{-1})$ ($x\to\infty$) [27, (6.19)], the only loss from our more restrictive decay condition is a constant, which, as discussed above, may be contributed by $u_{1}$.

From (18), we see that the two subproblems are now fully coupled. We cast the subproblem for $u_{1}$ in variational form in anticipation of FEM discretization and the subproblem for $u_{2}$ in terms of layer potentials. To arrive at the coupled system, we first define the operator $\mathcal{R}$ as the restriction to the boundary $\partial\Omega$ and $\hat{\mathcal{R}}$ as the restriction to the boundary $\partial\hat{\Omega}$.

We write $u_{2}$ in terms of an unknown density $\gamma$ using a layer potential operator $\mathcal{A}$ such that $u_{2}=\mathcal{A}\gamma$, while ensuring that the resulting integral equation is of the second kind:

for some constant $C$.

Next, we decompose $u_{1}$ as

where $\tilde{u}_{1}\in H^{1}_{0}(\hat{\Omega})$ is zero on the boundary $\partial\hat{\Omega}$ and $\hat{u}_{1}\in H^{1}(\hat{\Omega})$ is used to enforce the boundary conditions. $\hat{u}_{1}$ is defined by a lifting operator $\mathcal{E}:H^{1/2}(\partial\hat{\Omega})\rightarrow H^{1}(\hat{\Omega})$ that selects a specific $\hat{u}_{1}$ in the volume from its boundary restriction $\hat{r}_{1}\equiv\hat{\mathcal{R}}\hat{u}_{1}$. (The precise choice of the lifting operator within these guidelines has no influence on the obtained solution $\hat{u}_{1}$.)

The coupled problem is then to find $\gamma\in C(\partial\Omega)$, $\tilde{u}_{1}\in H^{1}_{0}(\hat{\Omega})$, and $\hat{r}_{1}\in H^{1/2}(\partial\hat{\Omega})$ such that

where $\mathcal{A}\gamma$ for $u_{2}$ is used in (18).

Next, we isolate the density equation in (21) using a Schur complement, which results in

with the solution operator $\mathcal{U}:L^{2}(\hat{\Omega})\times H^{1/2}(\partial\hat{\Omega})\rightarrow H^{1}(\hat{\Omega})$, where $\mathcal{U}\,[\zeta;\;\hat{\rho}]$ is defined as the function $\mu=\tilde{\mu}+\mathcal{E}\hat{\rho}$, and where $\tilde{\mu}\in H^{1}_{0}(\hat{\Omega})$ satisfies

This allows us to express the IE solution to the coupled problem (21) in terms of input data $f$, $g$, and $\hat{g}$, along with the action of $\mathcal{U}$. Once $\gamma$ is known, the FE solution is found as $u_{1}=\mathcal{U}\,[f;\;\hat{g}-\hat{\mathcal{R}}\mathcal{A}\gamma]$.

The form in (22) identifies two remaining issues. The first is the choice of $C$, which is determined by selecting a layer potential representation $\mathcal{A}$ of $u_{2}$. The conventional choice, a double layer potential, is not suitable because the exterior limit of the double layer operator $\mathcal{D}$ has a nullspace spanned by the constants. A common way of addressing this issue involves adding a layer potential with a lower-order singularity [27]; however, this is inadequate for our coupled FE-IE system (for $d=2$), as we explain below. Instead, we choose $\mathcal{A}=\mathcal{D}+\mathcal{S}$, the sum of the double and single layer potentials, each with the same density. This choice, by the exterior jump relations for the single and double layer potentials [27] establishes $C=1/2$.

Second, uniqueness and existence for (21) hinges on joint compactness of the composition of operators $\mathcal{R}\mathcal{U}\,[0;\;\hat{\mathcal{R}}\mathcal{A}]$ using the next lemma.

Using slightly different machinery, a discrete version of Lemma 2 is available at least in $\mathbb{R}^{2}$. To this end, let $\hat{\Omega}\subset\mathbb{R}^{2}$ be convex and polygonal and define a finite element subspace $V^{h}\subset H^{1}(\hat{\Omega})$ of continuous polynomials of degree $\geq 1$ on a quasi-uniform triangulation of $\hat{\Omega}$ (in the sense of [40]). Also define $V^{h}_{0}:=H^{1}_{0}(\hat{\Omega})\cap V^{h}$. Further define the discrete lifting operator $\mathcal{E}^{h}:H^{1/2}(\partial\hat{\Omega})\to V^{h}$ and the discrete solution operator $\mathcal{U}^{h}:V^{h}\times H^{1/2}(\partial\hat{\Omega})\to V^{h}$ where $\mathcal{U}^{h}\,[\zeta;\;\hat{\rho}]$ is defined as the function $\tilde{\mu}^{h}+\mathcal{E}^{h}\hat{\rho}$, where $\tilde{\mu}^{h}\in V^{h}_{0}$ satisfies

(Again, the precise choice of the discrete lifting operator within these guidelines has no influence on the obtained solution.)