Discontinuous Galerkin methods for the Laplace-Beltrami operator on point clouds¹¹1Emails: guozhi.dong@csu.edu.cn; hailong.guo@unimelb.edu.au; zqshi@tsinghua.edu.cn

Let $\Gamma_{h}:=\bigcup_{j\in J}\Gamma_{h}^{j}$ represent a subdivision of $\Gamma$, where $h$ denotes a parameter determined by the subdivision’s size, and $J\subset\mathbb{N}$ stands for an index set. Each $\Gamma_{h}^{j}\subset\Gamma$ is homeomorphic to an open subset in $\mathbb{R}^{2}$. Specifically, we assume $\Gamma=\bigcup_{j\in J}\overline{\Gamma_{h}^{j}}$ and $\Gamma_{h}^{j_{1}}\cap\Gamma_{h}^{j_{2}}=\emptyset$ for all index pairs $j_{1}\neq j_{2}\in J$. In practical applications, the explicit representation of the manifold $\Gamma$ might be unavailable. Thus, $\Gamma_{h}$ remains unknown, and only approximations or discrete representations are accessible. Consequently, we consider a collection of patches $\hat{\Gamma}_{h}=\bigcup_{j\in J}\hat{\Gamma}^{j}_{h}$, where each $\hat{\Gamma}^{j}_{h}$ approximates $\Gamma_{h}^{j}$ for every $j\in J$. We assume that every $\hat{\Gamma}^{j}_{h}$ is also homeomorphic to an open subset in $\mathbb{R}^{2}$, and there exists a bijective map between $\hat{\Gamma}^{j}_{h}$ and $\Gamma^{j}_{h}$. Particularly, we consider the following $k$-th order patchwise manifold:

Throughout this article, the letter $C$ or $c$, with or without subscripts, denotes a generic constant which is independent of $h$ and may not be the same at each occurrence. For convenience, we represent $x\leq Cy$ (or $x\geq Cy$) as $x\lesssim y$ (or $x\gtrsim y$). Later in this section, we shall provide the detail about how to construct $\hat{\Gamma}^{k}_{h}$ from point clouds. It is important to note that for patchwise manifolds $\{\hat{\Gamma}^{k}_{h}\}$, global continuity is not always necessary. An example of such manifolds includes local polynomial patches. Algorithm 1 illustrates how to construct them from a given point cloud representation of $\Gamma$. To achieve this, we introduce a parameter domain denoted by $\Omega_{h}=\bigcup_{j\in J}\Omega_{h}^{j}$. We shall abuse a bit of notation in the paper, as later $\Omega_{h}^{j}$ may denote both a planar domain (e.g., a triangle) and also the corresponding hyperplane to which the parameter domain belongs. Subsequently, we establish a patchwise mapping $\pi:\Omega_{h}\to\Gamma_{h}$. One approach to defining this map is by utilizing normal vectors on $\Gamma_{h}$ as in (2.4). On each $\Omega_{h}^{j}$, we construct a local coordinate system where the origin lies within $\Omega_{h}^{j}$, and $\pi$ serves as a mapping from $\mathbb{R}^{2}$ to $\mathbb{R}^{3}$. This mapping $\pi$ induces a metric tensor $g$ of $\Gamma$ through the following relation:

It is important to highlight that $\pi$ isn’t necessarily defined using the global coordinates in $\mathbb{R}^{3}$, but alternatively defined on each local patch $\Omega_{h}^{j}$ using the local coordinates within it. Here, $\partial$ denotes the Jacobian, for instance, $\partial\pi\in\mathbb{R}^{3\times 2}$ signifies the Jacobian of the geometry map evaluated locally based on the local coordinates on $\Omega_{h}^{j}$. Remarkably, (2.2) plays a pivotal role in the error analysis presented in this paper.

A bijective map $\pi^{k}_{h}:\Omega_{h}\to\hat{\Gamma}^{k}_{h}$ can be established analogously to $\pi$. Consequently, we obtain the bijective map $\pi^{k}_{h}\circ\pi^{-1}:\Gamma_{h}\to\hat{\Gamma}^{k}_{h}$. Using the parametrization map $\pi$ (or $\pi^{k}_{h}$), the gradient of a scalar function on manifolds can be computed as follows:

Here, $\bar{v}=v\circ\pi$ (or $\bar{v}^{k}_{h}=v^{k}_{h}\circ\pi^{k}_{h}$), and $\nabla$ represents the gradient operator on the planar parametric domain $\Omega_{h}$. Note that this expression of surface gradient is different from the one used in many references for surface finite element methods like [18, 24]. This novel approach enables the development of a new numerical analysis framework to investigate PDEs on surfaces or (sub-)manifolds more broadly.

One should be aware that the choice of $\Omega_{h}$ is flexible, and the mapping $\pi:\Omega_{h}\to\Gamma_{h}$ is non-unique. Though there are many different choices of $\pi$, the metric tensor $g$ is invariant. Without loss of generality, we consider in this paper $\Omega_{h}$ consisting of triangle faces of some polyhedron close to $\Gamma_{h}$, and the following particular widely used mapping $\pi:\Omega_{h}\to\Gamma_{h}$ in (2.4):

where $\nu(x)=\nu(\pi(x))$ is the unit outward normal vector at the point $\pi(x)$ on $\Gamma$, and $d(x)$ is the sign distance function of $\Gamma$. Note that the expression in (2.4) is coordinate independent. We always consider a patchwise local coordinate system where the origin is on the face of $\Omega_{h}$. Assuming $\Omega_{h}$ to be close enough to $\Gamma$, then the mapping $\pi$ is bijective between $\Omega_{h}^{j}$ and $\Gamma_{h}^{j}\subset\Gamma$. This can be argued as follows: The Jacobian of $\pi$ is $\partial\pi(x)=P_{\tau}-d(x)H$. Here $H$ is the second fundamental form, whose eigenvalues are the principal curvatures $\kappa$ of $\Gamma$. $P_{\tau}$ is the tangential projection and it satisfies $P_{\tau}\partial\pi(x)=\partial\pi(x)=\partial\pi(x)P_{\tau}$. Let $\kappa$ be the scalar curvature of $\Gamma$. It is easy to see that the rank of $\partial\pi$ is the same as the dimension of $\Gamma$ provided that $d\kappa<1$. Thus the pseudo-inverse of $\partial\pi$ is bounded, and we have $\pi:\Omega_{h}\to\Gamma_{h}$ bijective. On the other hand $d\kappa<1$ requires the distance of $\Omega_{h}$ to $\Gamma_{h}$ to be inversely proportional to the curvature.

Eventually, we will develop DG methods for solving PDEs on $\hat{\Gamma}^{k}_{h}$ instead of $\Gamma_{h}$ with a novel error analysis based on local parametrization. Although we have introduced the concept of $k$-th order patchwise manifolds, it remains unclear how to obtain such patchwise manifolds from point clouds and under which conditions to have the approximation property in (2.1). In the rest of this section, we shall address this problem, which is important for designing optimally convergent high-order DG methods on point clouds.

To get the reference mesh, we reconstruct an initial triangular mesh $\Omega_{H}$. For this initial mesh, we only require that it is sufficiently close to the surface such that the map $\pi$ given in (2.4) is well defined. This requirement is very mild. Generating the initial mesh can be easily achieved from a point cloud using surface reconstruction methods such as screened Poisson reconstruction (SPR) [32], Gauss reconstruction [37], and others. Figure 1 illustrates an example of the initial triangular mesh constructed using SPR.

If the initial mesh is too coarse to be used in the computation, we can refine the initial mesh to a satisfactory one. Each triangle in the mesh is subdivided into four triangles by connecting its edge centers until the diameter of the largest triangle reaches the filling distance $h_{s}$ outlined in Definition 2.2.

Note that though the definition of filling distance involves the knowledge of $\Gamma$, we only need a magnitude estimate of it for the convergence analysis. This could be achievable without knowing $\Gamma$; instead, we may use the largest distance of neighbouring points in $\Sigma$, which is close to $2h_{s}$, to have such a bound.

The refined mesh is denoted as $\tilde{\Omega}_{h}$. Subsequently, for each vertex of $\tilde{\Omega}_{h}$, we identify its $k$-nearest neighbors ($k\geq 3$) in the point cloud $\Sigma$ to fit a hyperplane. The vertex is then projected onto the hyperplane to obtain a new vertex. Finally, connecting all the new vertices yields the reference mesh $\Omega_{h}$, maintaining the same connections as those in $\tilde{\Omega}_{h}$. It can be easily verified that $\Omega_{h}$ possesses the property:

This property will be utilized in the error analysis to achieve an optimal convergence rate. Throughout the following discussion, we assume that the filling distance $h_{s}$ in (2.5) is sufficiently small, and the largest diameter of the triangles on $\Omega_{h}$ satisfies $c_{a}h_{s}\leq h\leq c_{b}h_{s}$. Here, $c_{a}$ and $c_{b}$ are constants, controlling the magnitude of $h$ to ensure that the polynomial reconstruction in Algorithm 1 is well-posed, meaning that the polynomial function is uniquely determined from the scattered points there.

Using the previously obtained reference mesh as a foundation, we are able to create local patch approximations by Algorithm 1.

We point out that $\hat{\Gamma}_{h}^{k,j}$ is a piecewise polynomial patch, but not necessarily interpolating $\Gamma_{h}$. Now we provide the details on computing $\psi^{j}_{k}$ in Step (5) Algorithm 1. The idea is inspired from the equation in (2.4). More explicitly, given a point $x\in\Omega_{h}^{j}$ over the parametric domain, we compute the point $\psi\in\hat{\Gamma}_{h}^{k,j}$ such that $\psi(x)-x=-d_{h}^{k}(x)\nu_{h}^{k,j}(x)$, where $\nu_{h}^{k,j}$ denotes the outward unit normal vector of the local polynomial graph, and $d_{h}^{k}$ is the length of the vector $\psi-x$. The point $\psi(x)$ is then analogous to $\pi(x)$. Let $\psi=(\psi_{1},\psi_{2},\psi_{3})$ be a local polynomial graph, i.e., $\psi_{3}=p_{h}^{k,j}(\psi_{1},\psi_{2})$, and we consider the local coordinates, $x=(x_{1},x_{2},0)\in\Omega_{h}^{j}$. Notice that $d_{h}^{k}(x)=\sqrt{(\psi_{1}-x_{1})^{2}+(\psi_{2}-x_{2})^{2}+(\psi_{3})^{2}}$, and $\nu_{h}^{k,j}(x)=\frac{(-\partial_{1}p_{h}^{k,j},\;-\partial_{2}p_{h}^{k,j}\;,1)^{\top}}{\sqrt{(\partial_{1}p_{h}^{k,j})^{2}+(\partial_{2}p_{h}^{k,j})^{2}+1}}$.

This gives us the following nonlinear system

We rewrite the above system in a compact form as $M(\psi)=0$, which can then be solved using Newton’s method. For that, we reduce the above system to a two-variable equation

This system then gives rise to the following iterations for some initial guess, e.g., $(\psi_{1,0},\psi_{2,0})^{\top}=(x_{1},x_{2})^{\top}$:

where $M^{m}_{i}=M_{i}(\psi_{1,m},\psi_{2,m})$, and

which is a symmetric matrix. Note that its inverse matrix can be explicitly calculated

The Newton iterations are terminated once an expected accuracy level is reached (in our examples, we choose it to be $10^{-14}$).

One may notice that a continuous surface could potentially be reconstructed from the above algorithm, where one can add an additional step after Step (5) by averaging the nodal points corresponding to the same edges in the reference mesh. This would unify the nodal points on the common edge, and then use the same polynomial interpolation as in Step (6) of Algorithm 1 to achieve a continuous patchwise reconstruction. The same error bounds will hold for this continuous approximation, though proving it will require a bit more work. Since our analysis can cover both discontinuous and continuous cases, and the former appears to be more general in the error analysis, we will focus on the discontinuous case in this paper, as high-order numerical methods for PDEs on continuous surfaces have already been discussed in the literature, e.g., in [18].

We now demonstrate the approximation property for the nodal points $\{\psi^{i}\}_{i\in\mathbb{N}}$ reconstructed from Algorithm 1. Considering the equation in (2.8), the partial differential operators might seem to reduce the approximation accuracy of the polynomial functions $p_{h}^{j}$ after the Newton iterations. However, in the following, we show that surprisingly this is not the case, as long as the Newton algorithm (2.10) converges sufficiently.

We employ the subsequent vector decomposition in the paper.

Additionally, we introduce the following assumptions concerning $\Gamma_{h}$, $\Omega_{h}$, and $\hat{\Gamma}_{h}^{k}$, based on which we establish that the patches reconstructed from Algorithm 1 constitute a $k$-th order patchwise manifold.

Before proving Proposition 2.7, we present an auxiliary lemma regarding the polynomial approximation of scattering points. This lemma constitutes a standard result available in various references, such as [41]. For the sake of completeness, a proof is included in Appendix A.

The error bounds (2.11) and (2.12) enable us to apply previous results in Lemma 2.11 for the reconstructed patches from Algorithm 1. Note that for continuous interpolated patches, the error in Proposition 2.7 vanishes.

The bounded curvature of $\Gamma$ implies another useful estimate $\lvert\nu_{\Omega_{h}}\rvert\lesssim h$ as demonstrated in the following lemma.

The scaling constant $C$ can be further optimized with respect to the scalar curvature $\kappa$, although it is not pursued here.

To simplify notations, we will not distinguish metric tensors $g$ and $g^{k}_{h}$ with $g\circ\pi$ and $g^{k}_{h}\circ\pi^{k}_{h}$ in the following. In Theorem 2.12, we show that $\bigcup_{j\in J}\hat{\Gamma}_{h}^{k,j}$ given by polynomial interpolating the nodal points $(\psi_{i}^{j})_{i\in I_{j}}$ is a $k$-th order patchwise manifold.