Evaluation of Inner Products of Implicitly-defined Finite Element Functions on Multiply Connected Planar Mesh Cells

Given $v,w\in V_{p}(K)$, note that $\Delta v$ and $\Delta w$ are given polynomials of degree at most $p-2$. Let $P$ and $Q$ be polynomials of degree at most $p$ satisfying

As pointed out in [19], such polynomials $P$ and $Q$ can be explicitly constructed term-by-term by observing that

is a polynomial anti-Laplacian of $x^{\alpha}$ for a multi-index $\alpha$. That is, $\Delta P_{\alpha}(x)=x^{\alpha}$. Note that, in practice, $P_{\alpha}$ is obtained only by manipulation of polynomial coefficients, and poses no computational barrier. The same can be said of other operations involving polynomials, such as gradients, etc.

Since the functions

are harmonic, we have the expansion

For the first two integrals in the final expression, the normal derivatives $\partial\phi/\partial\mathbf{n}$ and $\partial\psi/\partial\mathbf{n}$ may be computed using the Dirichlet-to-Neumann map discussed below. Furthermore, $\nabla P\cdot\nabla Q$ is clearly a polynomial of degree at most $2p-2$. As noted in [3], a straightforward application of the Divergence Theorem shows that

In this fashion, we reduce the volumetric integral $\int_{K}\nabla v\cdot\nabla w~{}dx$ to readily computable boundary integrals.

Consider the problem of determining the normal derivative of a harmonic function $\phi$ given its Dirichlet trace $\phi|_{\partial K}$. Recall that $\widehat{\phi}$ is a harmonic conjugate of a harmonic function $\phi$ whenever $\phi,\widehat{\phi}$ are continuously twice differentiable on $K$ and satisfy the Cauchy-Riemann equations:

Given that $\phi$ is harmonic on a simply connected domain $K$, the existence of a harmonic conjugate of $\phi$ is guaranteed, and $\widehat{\phi}$ is unique up to an additive constant. If $\partial K$ is smooth, for every $x\in\partial K$ it holds that

where $G(x,y)=-(2\pi)^{-1}\ln|x-y|$ is the fundamental solution of the Laplacian in $\mathbb{R}^{2}$. Supposing that the boundary $\partial K$ is traversed counterclockwise, we let $\mathbf{t}$ denote the unit tangent vector and $\mathbf{n}$ denote the outward unit normal vector, so that the normal and tangential derivatives of $\phi$ and $\widehat{\phi}$ are related by

from which the right-hand side in (6) can be computed. Since $\widehat{\phi}$ is unique only up to an additive constant, we impose $\int_{K}\widehat{\phi}~{}ds=0$, which we add to the left-hand side above to obtain

In practice, we solve this integral equation numerically for $\widehat{\phi}$ on $\partial K$ using a Nyström method, where the right-hand side is computed using the tangential derivative of $\phi$, which is readily accessible from its trace $\phi|_{\partial K}$. Having obtained values of the harmonic conjugate $\widehat{\phi}$ on the boundary $\partial K$, we may obtain its tangential derivative $\partial\widehat{\phi}/\partial\mathbf{t}$ via numerical differentiation, which then yields values of the normal derivative $\partial\phi/\partial\mathbf{n}$. Indeed, if $x(t)$ is a sufficiently smooth parameterization of $\partial K$ and we define $G(t)=\widehat{\phi}(x(t))$, then

Since $G(t)$ is periodic, a natural choice to obtain $G^{\prime}(t)$ is to write a Fourier expansion $G(t)=\sum_{k=-\infty}^{\infty}\omega_{k}\,e^{\mathfrak{i}kt}$ and obtain an approximation of $G^{\prime}(t)$ by truncating the series

In practice, a Fast Fourier Transform (FFT) may be used on a discretization of $G(t)$, and an inverse FFT used on the coefficients $\mathfrak{i}k\omega_{k}$ to obtain the approximate values of $G^{\prime}(t)$.

Details of such calculations, including the case where $\partial K$ has corners, are discussed in [22].

Let $v=\phi+P$ and $w=\psi+Q$ be as above. We have the expansion

Notice that the last integral can be computed with (5), whereas the two middle integrals have the form

where $r$ is a polynomial and $\eta$ is harmonic. Using (3), let $R$ be a polynomial such that

then applying Green’s Second Identity, we have

The remaining integral to be computed in (11) is the $L^{2}$ inner product of the two harmonic functions $\phi$ and $\psi$. Toward this end, suppose that $\Phi$ is an anti-Laplacian of the harmonic function $\phi$, that is,

Then using Green’s Second Identity again yields

The problem of determining such a $\Phi$, in particular its trace $\Phi|_{\partial K}$ and normal derivative $\partial\Phi/\partial\mathbf{n}$, is addressed as follows.

Notice that if $\Phi$ is an anti-Laplacian of $\phi$, it holds that $\Phi$ is biharmonic, that is, $\Delta^{2}\Phi=0$. A well-known fact (see, for example, pp. 269 of [12]) is that every biharmonic function is of the form

where $f,g$ are some analytic functions, $\text{Re}\,[z]$ denotes the real part of $z\in\mathbb{C}$, $\overline{z}$ denotes the complex conjugate, and we use the natural identitification of the complex plane with $\mathbb{R}^{2}$ via $x=(x_{1},x_{2})\mapsto z=x_{1}+\mathfrak{i}x_{2}$. Since any anti-Laplacian of $\phi$ will suffice for our purposes, we will take $g=0$ and write

where $\rho=4(\text{Re}\,f)$ is a harmonic function and $\widehat{\rho}=4(\text{Im}\,f)$ is a harmonic conjugate of $\rho$. It follows from the Cauchy-Riemann equations that

and that the gradients of $\rho$ and $\widehat{\rho}$ must take the form

where $\widehat{\phi}$ is a harmonic conjugate of $\phi$.

In order to compute $\rho$ and $\widehat{\rho}$, consider the following. For the sake of illustration, we assume that the boundary $\partial K$ is smooth, but a similar approach works when $\partial K$ is only piecewise smooth, using some minor modifications. Given the traces of $\phi$ and $\widehat{\phi}$ on $\partial K$, we have access to the tangential derivative of $\rho$ via

Given a sufficiently smooth parameterization $x(t):[0,2\pi]\to\partial K$ of the boundary, we define $g:[0,2\pi]\to\mathbb{R}$ by

By the Fundamental Theorem of Calculus, we may obtain an anti-derivative $G$ via

Note that $g$ is $2\pi$-periodic and admits a Fourier expansion

As mentioned above, in practice we truncate this series and compute the Fourier coefficients $\omega_{k}$ using an FFT. Integrating termwise yields

for an arbitrary constant $C$. In light of Remark 2.2, we may pick $C$ arbitrarily; for instance, choose $C=0$. Moreover, since $\rho(x(t))$ is $2\pi$-periodic, we also see that $\omega_{0}=0$. In the computational context, we apply an inverse FFT to the coefficients $-\mathfrak{i}\,\omega_{k}/k$ in order to obtain approximate values of $\rho$ on the boundary. We apply an analogous procedure to obtain values of $\widehat{\rho}$ on the boundary, using the tangential derivative data

and computing an anti-derivative of

In our discussion so far, we have assumed the cell boundary $\partial K$ to be smooth, but here we will briefly address how the calculations described above can be modified when $\partial K$ has one or more corners. Elaboration on these modifications can be found in [22].

Each edge of the mesh cell $K$ is discretized into $2n+1$ points, including the endpoints, so that the boundary is discretized into $N=2n\times(\#~{}\text{edges of}~{}K)$ points, with redundant endpoints being neglected. When an edge $e\subseteq\partial K$ is a $C^{2}$ smooth closed contour, the boundary points are assumed to be sampled according a strongly regular parameterization $x(t)$ of $e$ (i.e. $|x^{\prime}(t)|\geq\delta$ for all $x(t)\in e$ and for some fixed $\delta>0$). In the case where $e$ terminates at a corner, we employ a Kress reparameterization, defined as follows (cf. [20]). Suppose that $x(t)$ is a strongly regular parameterization of $e$ for $t\in[0,2\pi]$, then define $\widetilde{x}(u)=x(\tau(u))$ using

where the Kress parameter $\sigma\geq 2$ is fixed. The Kress reparameterization is not regular, with $\widetilde{x}^{\prime}(u)$ vanishing at the endpoints. Indeed, $\tau^{\prime}(u)$ has roots at $0$ and $2\pi$ of order $\sigma-1$, which leads to heavy sampling of the boundary near corners. This effect is amplified for larger values of $\sigma$.

Recall that whenever $\eta$ has a sufficiently smooth Dirichlet trace, we can compute the weighted tangential derivative

by using, for instance, the FFT-based approach described by (10). Replacing $x(t)$ with a Kress reparameterization leads to difficulty in recovering the values of the tangential derivative from the weighted tangential derivative. The same can be said for the weighted normal derivative

Notice, though, that for the sake of computing the boundary integral

the weighting term $|x^{\prime}(t)|$ appears in the Jacobian anyway, so it is natural to keep the tangential and normal derivatives in their weighted form, inlcuding the case when using a Kress reparameterization.

Note that, for the sake of effectively applying an FFT, we assume that the parameter $t$ is sampled at equispaced nodes $t_{k}=hk$, $h=\pi/n$, $0\leq j\leq 2n+1$, and likewise for the parameter $u$ when using a Kress reparameterization.

Thus far, we have all the necessary tools to compute the goal integrals (1) and (2) in the case where $K$ is simply connected. It is worth reiterating that both of these volumetric integrals have been successfully reduced to contour integrals along the boundary $\partial K$, and there is no need for 2-dimensional quadratures as all necessary computations occur only on $\partial K$.

We have the option, though, of obtaining interior values of $v\in V_{p}(K)$ as follows. Write $v=\phi+P$ as above, and determine a harmonic conjugate $\widehat{\phi}$ of $\phi$. Then $f=\phi+\mathfrak{i}\widehat{\phi}$ is an analytic function, and for any fixed interior point $z=x_{1}+\mathfrak{i}x_{2}\in K$ we have Cauchy’s integral formula

Furthermore, we can obtain interior values of $\nabla\phi$ by observing that

Interior values of higher derivatives, such as the components of the Hessian, can be obtained in similar fashion if so desired.

We conclude this section with a few remarks about computational complexity. Assume the boundary $\partial K$ is parameterized and then discretized using $N$ points. The Nyström system resulting from (8) is dense, though well-conditioned, and simple linear solvers come with a computational cost $\mathcal{O}(N^{3})$. Using more sophisticated methods, such as GMRES, make an improvement, but in general will never be better than $\mathcal{O}(N^{2})$. Although even more sophisticated methods, such as those based on hierarchical matrices [14] or hierarchical semiseparable matrices [31], can reduce the computational complexity even further, for relatively small problems such as those considered here, GMRES is sufficient.

The FFT calls used for numerical differentiation have a computational cost of $\mathcal{O}(N\,\log N)$, and integration along $\partial K$ using (using, say, the trapezoid rule) take $\mathcal{O}(N)$ operations. Operations on polynomials, such as computing anti-Laplacians, can be performed by manipulation of the coefficients and do not meaningfully contribute to the computational cost. So despite the many terms we have encountered in the expansion of the integrals (1) and (2), in practice these expansions are relatively cheap in comparison to the cost of obtaining the trace of the harmonic conjugate. Note that the latter computation need only happen once for each function $v\in V_{p}(K)$ considered.

Additionally, we can use the notion of trigonometric interpolation to reduce computational cost even further, as was explored in [22]. With the boundary discretized into $N$ points, we can solve the the Nyström system obtained from (8) as usual to obtain the harmonic conjugate $\widehat{\phi}$. While performing numerical differentiation with FFT as proposed, we have the Fourier coefficients at our disposal, which allows for rapid interpolation to, say, $M=2^{m}N$ points. We then compute the boundary integrals obtained from expanding (1) and (2) using standard 1D quadratures (e.g. the trapezoid rule, Simpson’s rule, etc.) on the larger collection of $M$ points. The heuristics presented in [22] suggest that similar levels of accuracy are achieved as would be in the case where all $M$ sampled points are used for solving the Nyström system.

In the next two sections, we address how our approach can be modified to accomodate multiply connected mesh cells.

Our present goal is to determine the coefficients $a_{1},\dots,a_{m}$, as in the statement of Theorem 3.1, given the trace of a harmonic function $\phi$. We will see that we simultaneously determine $\widehat{\psi}$ by solving an integral equation similar to (8), and thereby arrive at the Dirichlet-to-Neumann map

Assume for now that the boundary $\partial K$ is $C^{2}$ smooth. The case with corners is handled with Kress reparameterization, as discussed in the previous section. Our current task is to generalize the technique described by (8) in the case where $K$ is multiply connected. An alternative approach to the method we discuss here is presented in [13], which is comparable to our method in terms of cost and accuracy when all boundary edges are smooth, but does not achieve similar levels of accuracy when corners are present.

Let $\widehat{\psi}$ denote a harmonic conjugate of $\psi$ satisfying $\int_{\partial K}\widehat{\psi}~{}ds=0$. Just as in the simply connected case, we have

Making the replacement

and rearranging yields

This integral equation is underdetermined due to the $m$ additional degrees of freedom $a_{1},\dots,a_{m}$ in contrast to (8). To resolve this, we multiply both sides of $\phi=\psi+\sum_{j=1}^{m}a_{j}\lambda_{j}$ by the normal derivative $\partial\lambda_{\ell}/\partial\mathbf{n}$ and integrate over $\partial K$ to obtain

Invoking Green’s Second Identity and the Cauchy-Riemann equations yields

To write this in a form more conducive to computation, observe that the Fundamental Theorem of Calculus for contour integrals implies that

since $\partial K$ consists of $m+1$ closed contours. From the Product Rule, we obtain

Therefore, $\widehat{\psi}$ ought to satisfy

In summary, we have obtained a system of equations (17) and (18) for determining the trace of $\widehat{\psi}$ on $\partial K$. Discretizing the boundary into $N$ points, as we did for the simply connected case, this system of equations yields a square augmented Nyström system in $N+m$ variables, which we may solve with the same techniques used for solving (8). The case when $\partial K$ has corners can be handled using a Kress reparameterization, just as in the simply connected case.

Next, we wish to construct an anti-Laplacian of a harmonic function

as in the statement of Theorem 3.1. It is simple to verify that

is an anti-Laplacian of $\lambda_{j}(x)=\ln|x-\xi_{j}|$, so if $\Psi$ is an anti-Laplacian of $\psi$, then we have that

is an anti-Laplacian of $\phi$. The normal derivative can be computed using

Analogous to the simply connected case, we might seek potentials $\rho,\widehat{\rho}$ of the vector fields

While $\mathbf{F}$ and $\widehat{\mathbf{F}}$ both have vanishing curls, this is not sufficient to guarantee that they are both conservative on a multiply connnected domain—a simple counterexample being

taken on a circular annulus centered at the origin.

An elementary observation from complex analysis is that an analytic function $g=\psi+\mathfrak{i}\widehat{\psi}$ has an antiderivative if and only if $\mathbf{F}=(\psi,-\widehat{\psi})$ and $\widehat{\mathbf{F}}=(\widehat{\psi},\psi)$ are conservative vector fields. Indeed, $G=\rho+\mathfrak{i}\widehat{\rho}$ satisfies $G^{\prime}=g$ for $\nabla\rho=\mathbf{F}$ and $\nabla\widehat{\rho}=\widehat{\mathbf{F}}$. This simple fact inspires us to decompose $g$ as $g_{0}+g_{1}$, where $g_{0}$ has an antiderivative and the real part of $g_{1}$ has an anti-Laplacian which can be computed a priori. The following proposition reveals such a decomposition. The proof is inspired by that given in [4] and is unlikely to surprise a reader familiar with elementary complex analysis, but we include it for the sake of completeness.

Using the above results, we may decompose $g=\psi+\mathfrak{i}\widehat{\psi}$ into one part that has an antiderivative

and another part that is a linear combination of rational functions

where we use $\zeta_{j}=\xi_{j}\cdot(1,\mathfrak{i})$, with $\xi_{j}$ being chosen when Theorem 3.1 is applied. Easy calculations verify that

satisfies $\Delta M_{j}=b_{j}\mu_{j}-c_{j}\widehat{\mu}_{j}$, and whose normal derivative can be directly obtained from