The Laplace-Beltrami problem along surfaces in three dimensions appears in many fields, including plasma physics [beltrami3], fluid mechanics [rahimian2010petascale, veerapaneni_2011], electromagnetics [gendeb1, gendeb2, superconductor1, nedelec], surface diffusion, pattern formation [kim2020pattern, plaza2004effect], and computational geometry [kromer2018highly]. Suppose that is an at least twice differentiable surface without boundary embedded in . Let denote the intrinsic surface gradient, and let it’s adjoint on the space of functions on the boundary, denoted by , be the the intrinsic surface divergence. Given a mean zero function defined on , the Laplace Beltrami problem is to find a mean zero function satisfying , where is the Laplace-Beltrami operator defined as .
In many of the applications discussed above, the Laplace-Beltrami problem naturally arises for the computation of the Hodge decomposition of tangential vector fields. Any smooth tangential vector field on the surface can be written as
where , are smooth functions on , is a harmonic field satisfying , and is the unit normal to the surface . The curl free component and the divergence free component can be computed via the solution of the following Laplace Beltrami problems:
A related problem is that of determining a basis for harmonic vector fields on non-contractible surfaces, which plays a crucial role for several integral representations in computational electromagnetics. On a genus object, the space of harmonic vector fields forms a dimensional vector space. One can compute elements of this basis via the solution of the two Laplace-Beltrami problems above for a randomly generated smooth vector field
, which with high probability tends to have a non-trivial projection onto the harmonic vector field subspace.
Numerical methods for the Laplace-Beltrami problem can be typically classified into two categories. The first category of methods rely on the direct discretization of the differential operator on surface. Examples include finite element methods[bansch2005finite, bonito2020finite, burman2017cut, demlow2007adaptive], virtual element methods [beirao2014hitchhiker, frittelli2018virtual], differencing methods [wang2018modified], and level set methods [bertalmio2001variational, chen2015closest, greer2006fourth, macdonald2008level, macdonald2010implicit]. The level set methods differ from the rest of the direct discretizations of the PDE as they rely on using an embedded finite difference grid in the volume in to discretize an extension of the surface. The other category of methods is to reformulate the Laplace Beltrami problem as an integral equation on the surface. On a flat surface, the Green’s function for the Laplace-Beltrami equation is , which can be used as a parametrix for the solution on general surfaces as well [kropinski_2014_fast]. This idea has also been used for the solution of the related Yukawa-Beltrami problem on general surfaces including the surfaces with boundary in [kropinski_2016_integral], and for the solution of Laplace-Beltrami problem on non-smooth surfaces [goodwill2021numerical]. An alternate integral equation based approach is to use the Green’s function of the free-space Laplace equation in three dimensions as either left or right preconditioners for the Laplace Beltrami problem [oneil2018].
In this work, we present a fast multipole method accelerated iterative solver for the solution of the Laplace Beltrami problems on complicated three dimensional surfaces based on the integral formulation presented in [oneil2018], and using the locally corrected quadrature approach of [greengard2020fmm]. In some respects, this work can be considered part II of the earlier work [oneil2018]. The primary motivation of choosing this approach is two-fold: first that the integral equation tends to be as well-conditioned as the underlying physical problem, and secondly, since the representation uses layer potentials related to the free space Green’s function in three dimensions, the numerical solver can be coupled to existing optimized fast multipole method (FMM) libraries [greengard-1997, greengard2020fmm]. We demonstrate the well-conditioned nature of the integral equation through a battery of numerical experiments, each of which incorporates a high-order representation of the surface and high-order quadratures. We also show the importance of numerical second kindness when discretizing the integral equation using locally corrected quadrature methods that rely on the accurate computation of the principal value part of the layer potentials in the representation. In particular, while there are several equivalent mathematical formulations for the same integral equation on surface, we show that the formulation which explicitly handles the identity term of the integral equation tends to be the most stable numerical implementation. Lastly, with the acceleration provided by fast multipole methods, our integral formulation can be used to compute a basis for harmonic vector fields on relatively high-genus geometries.
The rest of the paper is organized in follows. In section 2, we describe the Laplace-Beltrami problem. In section 3, we review the integral equation formulation presented in [oneil2018], and in section 4, we review a high order FMM-accelerated iterative solver. In section 5, we discuss an application of the solver to the computation of a basis for harmonic vector fields on non-contractible surfaces. In section 6, we present several numerical examples to demonstrate the efficiency of our approach, including a discussion on the importance of numerical second kindness. Finally, in section 7, we conclude by discussing avenues for future research.
2 The Laplace-Beltrami problem
The Laplace-Beltrami operator, also known as the surface Laplacian, is the generalization of the Laplace operator to general (smooth) manifolds. In this work, we will focus on smooth surfaces embedded in three dimensions. Suppose is an oriented surface with at least two continuous derivatives, and let denote a local parametrization of
. The first fundamental form, or equivalently the metric tensor at, is given by
where and are abbreviations for and respectively and are the tangents to the surface at . The unit normal at can written as
The unit vectors in the directions of , , and form a local coordinate system (not necessarily orthogonal) at each point .
Next, suppose that is a smooth tangential vector field defined along and is expressed in terms of the local tangent vectors , as
The surface divergence of , denoted by is then [nedelec, frankel] given by
where is the determinant of the metric tensor . The surface gradient operator , formally defined as the adjoint of the surface divergence operator on the space of functions with respect to the induced metric on , is given by
where is a scalar function defined on and , , denote the components of the inverse of the metric tensor, denoted by .
The Laplace-Beltrami operator is then defined as , and is explicitly given by the expression
The standard Laplace-Beltrami problem is to solve the following PDE for along ,
where is a given function. Conditions on the regularity required in , and the subsequent regularity obtained in have been studied in [rosenberg, jost, epstein-2010]. We will assume that all functions and surfaces are sufficiently smooth so as to obtain high order accuracy in the subsequently numerical examples.
As is clear, the Laplace-Beltrami operator maps constant functions along to 0. Thus, the equation (2.7) is rank-one deficient [imbertgerard_2017, sifuentes_2015]. The following lemma[rosenberg, oneil2018] summarizes a classical result regarding the well-posedness of the above Laplace-Beltrami problem.
Let be a smooth, closed, and orientable boundary. Along , the Laplace-Beltrami operator is uniquely invertible as a map from , the space of mean-zero functions defined on ,
Indeed, if the function defined along has mean zero, then there exists a unique mean-zero solution, , to the Laplace-Beltrami problem .
3 Integral formulation
In this section, we review an equivalent well-conditioned second kind integral equation for the solution of the Laplace-Beltrami problem. This formulation relies on preconditioning the Laplace-Beltrami operator on both the left and right using the Laplace single layer potential. The resulting Fredholm operator can be expressed in terms of standard layer potentials for the free space Laplace equation in three dimensions.
Before describing this formulation, we first define the relevant layer potentials used in the formulation. The Green’s function for Laplace’s equation in three dimensions given by
The single and double layer potentials are then given by:
When the target point lies on , the integral defining the double layer should be interpreted in a principal value sense. For smooth surfaces , both and are compact operators as maps from . Related to the single and double layer potentials are the restrictions of their normal derivatives to given by
The integral defining is to be interpreted in a principal value sense, and that for is to be interpreted in a finite-part sense[kress_2014, nedelec]. The operator is also a compact operator on , whereas is a hypersingular operator.
Finally, let denote the restriction of the second normal derivative of to given by:
where, as for the integral defining , the integral should be interpreted in a finite-part sense [kress_2014, nedelec, oneil2018]. It can be shown that the difference operator is a compact operator on [nedelec], as the dominant singularities in the integrand cancel each other out.
In the following theorem, we now review the equivalent second-kind boundary integral equation presented in [oneil2018], for the solution of the Laplace-Beltrami equation (2.7).
Suppose that is a mean-zero function in . Then there exists a unique solution in which satisfies
and . Here is the mean curvature along , and is the operator given by
Furthermore, is the unique mean zero solution to .
A proof of this result can be found in [oneil2018].
4 An FMM-accelerated solver
In this section, we briefly describe the details of our numerical solver for integral equation formulations of the Laplace-Beltrami problem. A full detailed discussion would merely be a reproduction of much of the content of [greengard2020fmm], which presents a FMM-accelerated library of solvers for various boundary integral equations. Our numerical solver is based on the software packages fmm3dbie[FMM3DBIE] and FMM3D [FMM3D]; subsequently, examples for the Laplace-Beltrami problem will be included in the package fmm3dbie. We now give a very brief overview of the solver, and make a few remarks which are specific to the Laplace-Beltrami problem.
For simplicity, consider a second kind integral equation of the form
where is either or one it’s directional derivatives. Suppose that the surface is represented via a disjoint union of patches , and that each patch parametrized by a non-degenerate chart , where is the standard simplex given by
An order discretization of is one where the components of , and are represented by a polynomial expansion of total degree less than on . One way to obtain such a representation is to sample the chart and the derivative information at order Vioreanu-Rokhlin nodes [vioreanu_2014]
which are stable interpolation nodes for orthogonal polynomials[bremer_2012, koornwinder_1975] on .
Triangulations of surfaces obtained from CAD or standard meshing packages often tend to be low order and introduce several artificial edges and corners on the discretized surface. In our examples, with the aim of showing high-order convergence, we use the surface smoothing algorithm discussed in [vico2020surface] to convert these low-order triangulations of complex geometries to arbitrarily high order triangulations of a nearby related surface. In some of our examples, analytic parameterizations of the domains are known and are used instead of the algorithm of [vico2020surface]. The algorithm was only applied to complex geometries, such as the multi-holed object in Figure 3.
Let denote the union of the samples of the charts at the Vioreanu-Rokhlin nodes of order . A Nyström discretization of equation 4.1 is given by
where is an approximation for the solution . In [greengard2020fmm], the layer potential for is approximated using a combination of generalized Gaussian quadrature for the contribution of to the layer potential, adaptive integration for patches which are close to and target independent oversampled quadratures for all the remaining patches. Let denote the oversampled quadrature nodes and weights for integrating smooth functions on , then the Nyström discretization in equation 4.2 takes the form
where and are the interpolated value of at and respectively, obtained via local polynomial interpolation of the samples , and are the quadrature nodes and weights (which can be precomputed and stored) used in the adaptive integration procedure to accurately compute the contribution of the layer potential due to self panel and the near panels at the target location . There are patches which are in the near field of each target, and since the far field part of the layer potential is computed using a target independent quadrature rule, the corresponding contribution can be computed using standard fast multipole methods [greengard-1997, greengard2020fmm]. If equation 4.3 is well-conditioned, then we can obtain the solution in time using an iterative algorithm such as GMRES [saad-1986].
In the context of solving the Laplace-Beltrami problem, there are two details that aren’t already a part of the fmm3dbie [FMM3DBIE] package: a) the contribution of needs to be handled as a single operator since both nor are hypersingular operators, but their sum is a compact operator of order ; and b) the mean curvature needs to be computed by spectrally differentiating and .
5 Harmonic vector fields
For a surface with unit normal , a harmonic vector field is one which satisfies
On a genus surface, it is well known that the space of harmonic vector fields is dimensional [dai_2014, epstein-2013, epstein-2010]. Moreover, it follows from the definition that if is a harmonic vector field, then is also a harmonic vector field which is linearly independent of . Furthermore, harmonic vector fields are completely characterized by their integrals along and cycles of the geometry. Any tangential vector field along a smooth surface admits a Hodge decomposition [epstein-2010, epstein-2013]:
where and are scalar functions defined on the surface and is a harmonic vector field. Given a smooth tangential vector field , one can use this decomposition to compute linearly independent harmonic vector fields , by solving the following Laplace-Beltrami equations and computing :
A convenient way of choosing such smooth vector fields is to first define where each component is a random polynomial, and then set . For computational efficiency, we use as a low a polynomial degree as possible to define while making sure that the integrals over and cycles are linearly independent. Moreover, suppose that is contained in , then each component is defined to be a scaled tensor product Legendre polynomial [olver2010] , where is the standard Legendre polynomials of degree on . The extra factor of
in the choice of Legendre polynomials is to avoid requiring additional degrees of freedom to representon the surface due to large gradients of the Legendre polynomials near the end points and .
With high probability, different random polynomial vector fields result in different projections on the different elements of the dimensional space of vector fields. Recall that, if is a harmonic vector field then is also a harmonic vector field which is linearly independent of . Thus, in practice one can construct a basis for the harmonic vector field subspace by computing for different vector fields , and the remaining vector fields of the basis are defined as .
6 Numerical examples
In this section, we provide several numerical examples demonstrating the accuracy and computational efficiency of our solver for the Laplace-Beltrami problem, as well as an application of the solver for computing harmonic vector fields.
6.1 Numerical second-kindness
We first turn out attention to the differences in numerical discretizations of the following three mathematically equal integral representations for the Laplace-Beltrami PDE:
For convenience of exposition, we have dropped the rank-1 operator as it is not relevant to the discussion in this section. While the above operators (Eq. (6.1a), (6.1b), (6.1c)) are all equal as maps on smooth functions on , their performance when discretized numerically differs significantly on patch based discretizations of the surface.
The first approach in equation 6.1a of composing with is extremely unstable on patch based discretizations where smooth functions are represented as independent polynomial expansion on each patch. The discretization does not enforce any smoothness constraint on functions across the patch interfaces which results in a spurious null space proportional to the number of degrees of freedom used to discretize the problem. In particular, the surface divergence is computed after projecting the values of onto a truncated Koornwinder polynomial expansion. On each patch, the surface divergence operator has a null space proportional to the order of the discretization. To circumvent this issue, one could use discretizations commonly used in finite element methods which preserve smoothness of functions across patches. Equation 6.1b is an equivalent formulation which avoids this issue.
Even though equation 6.1b avoids the spurious null-space issue, both equations 6.1b and 6.1a tend to perform poorer when using iterative solvers than equation 6.1c. One of the key reasons for this difference is that any numerical discretization of integral operators is inherently compact since the operator is approximated to a finite rank operator. This results in the spectrum of the discretized integral operator to cluster close to the origin causing numerical conditioning issues. Explicitly adding the identity term, as is the case in equation 6.1c, results in the spectrum of the discretized integral operator to cluster at . This greatly enhances the convergence of iterative methods such as GMRES [saad-1986].
We demonstrate this difference in behavior, by inverting a known solution to the Laplace Beltrami problem on the unit sphere which is assumed to be discretized using number of patches of order . We suppose that the data is given by
where are the spherical harmonics [vico_2014] of degree and order . We discretize the layer potentials in representations in equations 6.1c, 6.1b and 6.1a based on the procedure discussed in section 4 and solve the linear system using GMRES, where at each iteration the matrix vector application is accelerated using an FMM. We verify the correctness of our discretization by computing the error in
where we have used the fact that on the unit sphere [vico_2014]. In table 1, we report which is the error in applying . For all three representations, we note that the error converges to at the expected rate of .
To study the relative performance of these representation in iterative solvers, we first compare the spectral properties of the three discretized linear systems. We form the dense matrix forfigure 1. As expected, we see that equation 6.1a has a large spurious null space and that the other two discretizations are full rank. Moreover, the singular values of equation 6.1b cluster at , while the singular values of equation 6.1c cluster at 0.25. We then solve the discretized linear systems for the given function using GMRES. To the right in figure 1
, we plot the estimated residual as a function of iteration number for the three representations. The GMRES residual forequations 6.1b and 6.1a stagnates at the level of discretization error, while equation 6.1c is more robust and the residual converges to machine precision independent of discretization error. The analytical solution for the given data is given by
In table 1, we report the error in the computed solution , denoted by . For equations 6.1c and 6.1b, the error converges to at the expected rate of (even though GMRES stagnates at discretization error for equation 6.1b), while for equation 6.1a, the spurious null space completely pollutes the solution.
6.2 Order of convergence
Next, we demonstrate the accuracy and convergence of our solver following the representation in equation 6.1c. To this end, we use the solver to compute a numerical solution of Laplace-Beltrami equation on a torus with major radius and minor radius and compare it with the known analytical solution. Suppose that the torus is parameterized as , where . A simple calculation shows that for this geometry,
Suppose that the torus is discretized with patches of order resulting in discretization points. We compute the numerical solution to Laplace-Beltrami problem on torus . In table 2, we tabulate the numerical results by presenting the error in the computed numerical solution (denoted by ), the computed rate of convergence and the mean of the computed solution (denoted by ). We observe expected order convergence to the analytical solution where the observed order is appropriately dependent on the discretization order of the patches.
6.3 An alternative integral equation
An alternative second-kind integral equation can be derived for the Laplace-Beltrami problem merely by right preconditioning by . In short, if we let , then after a suitable application of various Calderón identities [contopanagos-2002, oneil2018], the Laplace-Beltrami problem is replaced by
This formulation is also second-kind, but whose compact part differs ever so slightly from that of equation 6.1c. In table 3, we compare the two integral representations (equation 6.1c and equation 6.6) with respect to number of GMRES iterations required, and the accuracy in computing on the torus. Since the compact part of both the representations don’t differ by much, the numerical performance of both these representations is practically indistinguishable.
|error with||# iterations||error with||# iterations|
|100||5||1500||6.3 E-5||9||6.5 E-5||10|
|400||5||6000||2.0 E-6||8||2.0 E-6||8|
|1600||5||24000||6.5 E-8||8||6.4 E-8||8|
As can be seen from the results above, there isn’t any notable difference in solving the Laplace-Beltrami PDE by representing the solution and preconditioning on the left with another application of , or by representing , both in terms of the accuracy achieved and the performance of the iterative solvers. However, in the context of solving Maxwell’s equations or static currents in type I superconductors using the Debye source representation, we find that using is more faithful to the spectral properties of the solution [gendeb1, gendeb2, superconductor1]. As a map from mean zero functions to itself, if , then , which implies that the map from is a pseudo-differential operator of order . However, if we represent , then the map from is a pseudo-differential operator of order , while for the representation , the map from is indeed a pseudo differential operator of order . This significantly affects the spectral properties of some of the operators appearing in the generalized-Debye representation.
For example, one of the operators in the Generalized-Debye representation is
where is a surface charge, and is the single-layer potential for the Helmholtz equation [papas] with wave number given by
Note that is a pseudo-differential operator of order . If we represent the solution of using , then we need to discretize the following operator instead of ,
It turns out that is a second kind Fredholm operator of the form for some compact operator as opposed to which is compact operator. To the best of our knowledge, explicit Calderon identities are not known for this operator. However, one can deduce the second-kindness of the operator and the strength of the identity term by analyzing the action of the operator on spherical harmonics when the boundary is the unit sphere. In particular, if is the unit sphere, then
where and are the spherical Bessel and Hankel functions of order respectively. This result can be derived using the results in [vico_2014] for example. From the asymptotic expressions of spherical Bessel and Hankel functions, it follows that as . Not handling the identity term explicitly can affect the performance of quadrature methods (such as Generalized Gaussian quadrature which is the method used in fmm3dbie) as these methods rely on discretizing just the principal value parts of operators. No such issue arises when is replaced with .
6.4 Harmonic fields
Now we present the numerical results for the computation of the harmonic vector fields on a twisted torus geometry whose boundary is parametrized by with
where the non-zero coefficients are , and . For this example, we define , with and and then define . In figure 2, we plot the tangential vector field , it’s curl-free component, it’s divergence-free component and the harmonic vector field . The second linearly independent harmonic vector field can be obtained by computing . The surface was discretized using patches of order . In table 4, we tabulate the results of a convergence study. We estimating the error in the computation of the harmonic vector fields by computing the norms of and computed using spectral differentiation. For this example , so the above error estimate serves as a proxy for a relative error estimate as well. We observe the surface divergence and surface curl converge to zero as indicated in table 4.
We also compute the 20 linearly independent harmonic vector fields on a genus 10 surface. For this geometry, we choose a random surface vector field whose components are random tensor product Legendre polynomials as discussed in section 5. In figure 3, we plot the surface along with these 20 linearly independent harmonic vector fields computed using our solver. The geometry was discretized using patches of order . As before, we verify the computations by calculating the norms of the surface divergence of and , where is the computed harmonic field and observing that they converge to zero. In table 5, we tabulate convergence results for a harmonic vector field on the genus 10 surface computed using our solver.
7 Conclusions and future work
In this paper, we have presented a high-order and fast multipole accelerated iterative solver for the numerical solution of the Laplace-Beltrami equation on complicated surfaces in three dimensions. The Laplace Beltrami PDE is converted to a second kind integral equation by representing the solution or , and using appropriate Calderon identities. The resulting integral equation which requires the application of various Laplace layer potentials is then discretized using a high-order method with locally corrected quadratures for computing the weakly-singular layer potentials, and accelerated using standard fast multipole methods.
While the integral equation can be written in several equivalent analytic forms, we demonstrate the necessity of maintaining numerical second-kindness by explicitly isolating the identity term of the Fredholm equation in order to avoid stagnation of iterative solvers like GMRES. In some of these equivalent formulations, the operator ends up having a large dimensional spurious null space (proportional to the number of patches used). This issue is related to the difficulty of directly discretizing the differential form of the Laplace-Beltrami equation when using patch based discretizations of surfaces and not enforcing smoothness across patches. However, both of these issues can be mitigated by explicitly isolating the identity part of the Fredholm equation.
Finally, we also presented numerical examples demonstrating the computation of harmonic vector fields on surfaces using the Laplace-Beltrami solver. These vector fields are both of mathematical interest and are also required for solving problems in type I superconductivity and electromagnetism on topologically non-trivial geometries.
There are still several open questions that remain to be addressed. These include coupling this high order discretization to fast direct solvers; subsequently using these solvers to develop iterative solvers for the solution of Maxwell’s equations using the generalized Debye formulation; and extending these ideas for the solution of surface diffusion problems that arise in pattern formation and cell biology. These are all ongoing areas of research.