A Meshfree Collocation Scheme for Surface Differential Operators on Point Clouds

We present a meshfree collocation scheme to discretize intrinsic surface differential operators over scalar fields on smooth curved surfaces with given normal vectors and a non-intersecting tubular neighborhood. The method is based on discretization-corrected particle strength exchange (DC-PSE), which generalizes finite difference methods to meshfree point clouds. The proposed Surface DC-PSE method is derived from an embedding theorem, but we analytically reduce the operator kernels along surface normals to obtain a purely intrinsic computational scheme over surface point clouds. We benchmark Surface DC-PSE by discretizing the Laplace–Beltrami operator on a circle and a sphere, and we present convergence results for both explicit and implicit solvers. We then showcase the algorithm on the problem of computing Gauss and mean curvature of an ellipsoid and of the Stanford Bunny by approximating the intrinsic divergence of the normal vector field. Finally, we compare Surface DC-PSE with surface finite elements (SFEM) and diffuse-interface finite elements (DI FEM) in a validation case.


Introduction
Partial Differential Equations (PDEs) on curved surfaces and differentiable manifolds are an important tool in understanding and studying physical phenomena such as surface flows [36,25] and active morphogenesis [21].Analytically solving intrinsic PDEs in curved surfaces, however, quickly becomes impossible for nonlinear PDEs or for surfaces that do not possess a global parameterization.Therefore, numerical methods for solving intrinsic PDEs on curved surfaces are important, and a wide variety of both embedded and embedding-free schemes have been developed to consistently discretize intrinsic differential operators over scalar fields on surfaces.
Embedding-free methods require a (at least local) parameterization of the surface in order to discretize the differential operators via coordinate charts or a local basis of the manifold [37].This includes methods based on local moving frames [6], a concept originally developed in continuous group theory, where the surface geometry is locally represented by intrinsic orthonormal bases.This has been used to solve surface PDEs over meshfree point clouds using Moving Least Squares (MLS) approximations [19].The concept of local moving frames has also been combined with discontinuous Galerkin discretization, e.g., to solve shallow-water equations on arbitrary rotating surfaces [7].Other embedding-free Finite Element Methods (FEM) include Intrinsic Surface FEM (ISFEM), which discretizes differential operators on a triangulation of the surface [2,14], and methods based on Discrete Exterior Calculus (DEC) [24].
Embedding methods discretize the surface problem in an embedding space of codimension 1 and use projections to restrict the differential operators computed in the embedding space to the surface manifold.This includes methods that use explicit tracer points to represent the surface, but interpolate to an embedding mesh to evaluate differential operators [18], diffuse-interface methods based on phase-field representations of the surface [23], embedding FEM such as TraceFEM [26] and diffuse-interface FEM [17,30], narrow-band level-set methods based on orthogonal extension of the surface quantities [9,4], level-set methods based on the closest-point transform [29,20], and volume-of-fluid methods for surface PDE problems [16].
While each of these methods has its specific strengths, embedding methods usually generalize better to complex-shaped or arbitrary surfaces [29].However, they tend to have higher computational cost, because computations are done in the higherdimensional embedding space and additional extension (for level sets), right-hand-side evaluation (for phase fields), or interpolation (for closest-point transforms) steps are required, albeit specific optimizations are available, e.g., for level sets [27].Embeddingfree methods are usually more accurate, because they avoid the interpolation and projection errors arising when the discretization of the embedding space does not trace the surface exactly, but they tend to be more difficult to implement and harder to generalize to complex-shaped or moving/deforming surfaces.
Here, we present a meshfree collocation method for PDEs on smooth and orientable curved surfaces with non-intersecting tubular neighborhood.The method combines elements from embedding and embedding-free approaches.It is algorithmically embedding-free in the sense that surface quantities are represented on tracer points that are contained in the surface.This also discretizes and represents the surface itself as a point cloud.But the method is mathematically related to embedding approaches, since the stencils used to approximate differential operators at the surface points are computed in the embedding space by a reduction operation along the local normal vector, which needs to be known or computed.Intuitively, this projects the discrete operators, rather than projecting the flux vectors as typically done in embedding methods.The resulting method therefore shares properties of moving frame approaches, such as the low dimensionality (and hence low computational cost) and the meshfree character [6,19].It combines these with properties of embedding methods, such as their flexibility in generalizing to complex surfaces [29], and their ability to compute extrinsic differential-geometric quantities.
Our method is based on the Discretization-Corrected Particle Strength Exchange (DC-PSE) collocation scheme for arbitrary (surface) point clouds.DC-PSE is related to Generalized Finite Difference Methods (GFDM) [35] and to MLS [32].Given the local surface normal n, we derive intrinsic discrete operators by first creating an embedding narrow-band and placing collocation points along the normal from each surface point.We then determine the regular DC-PSE operator kernels in the embedding space.These kernels are subsequently reduced under the condition of orthogonal extension ∇f •n = 0 for any (sufficiently) differentiable scalar field f (x) ∈ R to derive intrinsic kernels at the surface points x.This is possible due to the kernel nature of DC-PSE, and it preserves the information from the embedding space in a scheme that only requires computation over surface points.
This paper is organized as follows: Section 2 recollects the DC-PSE method for convenience and introduces the notation.In Section 3, we describe the Surface DC-PSE scheme for numerically consistent discretization of surface differential operators.We present validation and convergence result in Section 4 and conclude in Section 5.

Discretization-Corrected Particle Strength Exchange (DC-PSE)
DC-PSE is a numerical method for discretizing differential operators on irregular distributions of collocation points [32].The method was originally derived as an improvement over the classic Particle Strength Exchange (PSE) [12] scheme, reducing its quadrature error on irregularly distributed collocation points, but mathematically amounts to a generalization of finite differences [32].The PSE/DC-PSE class of collocation methods uses mollification with a symmetric smoothing kernel η ϵ (•) to approximate (sufficiently smooth) continuous functions where f ϵ (x p ) is a regularized approximation of the function f at location x p ∈ Ω of collocation point p.The scalar ϵ is the smoothing length (or the kernel width) of the mollification.Linear differential operators in R d , defined by the multi-index α = (α 1 , . . ., α d ) ∈ Z d with |α| = d i=1 α i are approximated by Taylor series expansion to find a discrete operator at collocation point x p .The order of approximation r depends on the kernel η ϵ used in Eq. (1), and h(x p ) is the average distance between collocation point p and its neighbors within the kernel support.We use the arithmetic mean of the L 1 -distances to compute h, but since all norms are equivalent, the convergence order (but not the actual error magnitude) is independent of the choice of average.The Taylor expansion yields integral constraints (also known as continuous moment conditions), which the kernel η ϵ needs to fulfill in order to reach a certain convergence order r [12].DC-PSE uses different kernels η p ϵ (•, •) for different collocation points p and directly acts on a given quadrature of Eq. (1) with collocation points x q ∈ Ω, resulting in the discrete operator: where N (x p ) are all collocation points in the neighborhood (of a certain radius r c defined by the kernel width) around point x p , as illustrated in Fig. 1a.The positive sign in the parenthesis is used for odd |α|, the negative sign for even |α|.This renders the operator conservative on symmetric collocation point distributions, i.e., when η p ϵ (x p , x q ) = η q ϵ (x q , x p ).In DC-PSE, the kernels η p ϵ are thus not determined from continuous moment conditions, as in PSE, but directly from the discrete moment conditions that result from substituting Eq. ( 4) into the quadrature of Eq. ( 1) [32] for a given set {x q } N q=1 .This adapts the kernels to the specific distribution of collocation points (hence the name "discretization-corrected") and avoids the quadrature error of PSE [12], leading to a scheme that is consistent with order r on almost1 arbitrary collocation point sets.This means that at each collocation point, a potentially different kernel is used for the same differential operator if the neighboring collocation points within the kernel support are distributed differently.Evaluating such a kernel at the locations of the collocation points yields a generalized finite-difference stencil, which reduces to the classic compact finite differences on regular grid arrangements of points [32].
DC-PSE kernels are determined at runtime by solving a small system of linear equations for each collocation point, resulting from the discrete moment conditions in its kernel neighborhood.For this, one can choose the function space such that the kernels are compact and symmetric.A frequent choice are polynomials windowed by truncated exponentials [33] of finite radius r c .The polynomial coefficients a γ are determined for a given α and given collocation points x q ∈ N (x p ), such that the following discrete moment conditions are satisfied: where is the discrete moment of order β of the kernel η α ϵ , and β min is the parity of |α|, because the zeroth moment Z 0 h vanishes for even operators.Under these conditions, DC-PSE is consistent with order r as long as i.e., the kernel width ϵ scales proportionally with the average inter-point distance h around x p [32].

Surface DC-PSE
We generalize DC-PSE to surface differential operators based on the following classic result [29,22]: Let S ⊂ R d be a differentiable manifold that possesses a tubular neighborhood T and is orientable 2 and f : S → R. Define F : T → R, such that the restriction F | S = f , and F is constant along the normal direction n of S, i.e., ∇F • n = 0.Then, on the surface S, where ∇ S f is the intrinsic surface gradient.A similar result is true for the intrinsic divergence operator (∇ S • ) and for tangential vector field that is extended by constant extension to all surfaces displaced along the normal of S [29,22].
Given this result, it is straightforward to see the advantages of a meshfree discretization: it allows for conforming discretization of the surface and for exact constant orthogonal extension by simply copying points along the normal.This creates an embedding narrow-band of exact closest-point function values within the tubular neighborhood T without a need for interpolation.If T is non-intersecting with a radius of at least r c everywhere, the conditions of the result in Eq. ( 9) are satisfied in the constructed embedding.Analogous to the closest-point method [29], one can then discretize differential operators in the embedding space.Due to the additive kernel nature of DC-PSE, the discrete operators in the embedding space can be reduced to only the surface points.
This reduction becomes clear from the formulation of the DC-PSE method.Indeed, we realize that the constant normal extension can be made internal to the operator evaluation by accumulating the kernel coefficients along the normals.To see this, consider the DC-PSE operator in Eq. ( 4) in the embedding space.The neighborhood N for the summation contains both surface points x s and normally extended points x n , as shown in Fig. 1b.Because the f (x n ) are identical copies of the values of the respective surface points, we note that the pre-factors f (x s ) ± f (x p ) in the summation of Eq. ( 4) are the same for all extended normal points and the corresponding surface point x s .Hence, for each given pair of a center point x p and another surface point x s , the interactions with the corresponding normally extended points can be factored out from the kernel summation: defining the surface kernels η S (x p , x s ).These can be evaluated over only the surface points x s = N S (x p ) in the in-surface neighborhood N S (x p ) around the surface point x p , see Fig. 1b, yielding the Surface DC-PSE operator: Importantly, the surface kernels η S (x p , x s ), summed over all orthogonally extended points, can directly be computed when determining the kernel weights and without explicitly creating or storing the normally extended points x n .
Evaluating a Surface DC-PSE operator involves only the neighboring points on the surface and requires no narrow band or normally extended grid, even though the construction of the operators uses an embedding.This leads to a corresponding reduction in computational complexity for operator evaluation, as computations are only performed on a (d − 1)-dimensional surface embedded in d-dimensional space.In comparison, the cost of operator evaluation for embedding methods such as the closestpoint method is O(k(d − 1)), where k > 1 is the narrow-band width, which scales proportionally with the order of convergence.

Surface DC-PSE kernel construction
Surface DC-PSE requires two algorithms that are not part of the standard, flat-space DC-PSE method: an algorithm to create the intrinsic neighborhood of a surface point p, and an algorithm to determine the surface kernels η S at a surface point p.We follow the example of Ref. [5] and use explicit component notation and a concrete example in order to directly relate to implementations in computer code.
Figure 1b illustrates a piece of the tubular neighborhood of radius r c (circle) of a curved surface S embedded in R 2 .The red point p at position x p on the surface is the "center" collocation point at which we derive the discrete Surface DC-PSE operator Q α S f (x p ) for a scalar surface field f .Points in light blue are surface points within the embedding-space neighborhood of radius r c (circle) of x p , and the yellow points are the orthogonal extensions x n .By default, the spatial separation δn between adjacent orthogonal extensions of the surface point p is the arithmetic mean of the distances between p and the other surface points within the kernel support, measured in the embedding space.This favors isotropic-resolution neighborhoods and, thus, low condition numbers of the DC-PSE kernel system matrix.The number of orthogonally extended points should be N n ≈ round(r c /δn) to either side of the surface, which is the default.Only surface points are actually allocated and stored.
In order to determine the DC-PSE kernel η p ϵ in the embedding space, the distances between x p and all collocation points in its embedding-space neighborhood N are required.In the example of Fig. 1b, the neighborhood (circle) includes the surface points N S = {s 1 , s 2 } and the normally extended points {n i } 10 i=1 .The two pale-yellow points are not part of the neighborhood.Algorithm 1 constructs the neighbor set along with the corresponding distances.Surface normals at a given point are indexed by the point index for better readability, i.e., n(x p ) := n p .In the example of the figure, this results in the output where d q is the distance between the collocation points p and q in the embedding space in units of ϵ p := ϵ(x p ).
Using this neighborhood data structure, the embedding-space DC-PSE operator at point p in the example of the figure reads: Since the embedding-space kernels are evaluated at concrete distances, the η p ϵ (d q ) are just scalar numbers.All kernel values that share the same pre-factor are thus summed to the surface kernels η S (x p , x q ), q ∈ N S , obtaining: For even differential operators, i.e., derivatives with even |α|, the third term vanishes identically and can be skipped in the calculations.But this is not the case for oddorder derivatives.With this rearrangement, each evaluation of the operator at a point p only requires three kernel evaluations instead of the 12 that would be required in the embedding case.In addition, the normally extended points never need to be allocated and stored, as all kernel computations can happen on the fly.Algorithm 2 details the procedure for Surface DC-PSE operator construction.For the example from Fig. 1b, this results in the surface DC-PSE kernel values: which can directly be used in Eq. ( 11) to evaluate surface differential operators.

Results
We validate and benchmark the Surface DC-PSE method.First, we verify its convergence in test cases with known analytical solution.Then, we show applications to cases with more general surfaces where no analytical solution is available.Finally, we compare with surface finite-element methods (SFEM) in a validation study.In all cases, orthogonal extension is exact, copying surface points along the known normals.Therefore, ∇F • n in Eq. 9 is always zero by construction.Numerically evaluating this term requires approximating the gradient ∇F , which we confirmed to converge with the order of accuracy of the discretization scheme used.

Laplace-Beltrami operator on a circle and a sphere
We start by verifying convergence for the Laplace-Beltrami operator on the unit circle S 1 .The collocation points are distributed regularly using equi-angular spacing.We use a normal spacing of δn = 3/(N p − 1) to compute the surface operators in Eq. ( 11) in a narrow band with N n = 4 layers on each side of the surface.N p is the number of surface points x s .We choose r c = 4.1δn as the operator support and r = 2 as the desired order of convergence.The Laplace-Beltrami operator is characterized by the multi-index α = (2, 0) + (0, 2).Note that this multi-index is 2-dimensional, despite the circle being one-dimensional, since the operators are constructed in the embedding space, but evaluated intrinsically.
We test the numerical approximation of the surface operator on the function in polar coordinates.The error is computed against the analytical solution The visualization of the numerical solution and the convergence plot of the absolute errors are shown in Fig. 2a,c.We observe second-order convergence to the analytical solution, as expected for r = 2.
We further test the method on the unit sphere S 2 .The collocation points are distributed using the Fibonacci sphere technique [13].We use a normal spacing of δn = 0.8/( 3 N p − 1) to determine the surface operators in Eq. ( 11) in a narrow band with N n = 2 layers on each side of the surface.N p is the number of points on the sphere.We choose r c = 2.9δn as the operator support and r = 2 and r = 4 as the desired orders of convergence.The Laplace-Beltrami operator is characterized by the three-dimensional multi-index α = (2, 0, 0) + (0, 2, 0) + (0, 0, 2).We test the numerical approximation of the surface operator on the scalar spherical harmonic function f (θ, ϕ) = Y lm (16) in spherical coordinates for the mode l = 4, m = 0 (Fig. 3a).The error is computed against the analytical solution and is plotted in Fig. 3c.We also use this test case to benchmark against the Closest Point (CP) method [29] with L 2 and L ∞ errors plotted in Fig. 3c.While Surface DC-PSE is less accurate than the CP method for second-order operators, it outperforms for fourth-order operators, which is likely due to the better condition numbers of the small linear systems to be solved for each point.
Finally, we perform a strong scaling benchmark of the computation time with increasing numbers of CPU cores with both codes implemented in the parallel computing library OpenFPM [15,34] in C++ and run on the same hardware.We plot the parallel efficiency (i.e., the speed-up divided by the number of cores) in Fig. 3b for both the construction and the evaluation of the operators of both methods.We further report the absolute wall-clock times on one and 24 cores in the inset table.These times show that evaluating the Surface DC-PSE operators over all points in the domain is about one order of magnitude faster than evaluating the CP transform [29].In addition, Surface DC-PSE scales better with increasing numbers of parallel CPU cores.This is due to the simpler kernel evaluation of DC-PSE 3 .Eventually, the efficiency for both methods drops, as is expected for strong scaling with constant communication overhead.The time required to construct the Surface DC-PSE operators with N n = 2, however, is about two orders of magnitude larger than that for constructing the CP representation in a narrow band of radius 5 as required by the regression support.However, it still scales better with increasing CPU core count.For Eulerian simulations, where collocation points do not move, the kernels have to be determined only once at the beginning of a simulation, or they can be loaded from files for standard point distributions, potentially leading to an overall lower computational cost of Surface DC-PSE.

Poisson equation on a circle and a sphere
Surface DC-PSE operators can also be used for implicit equations by solving a linear system of equations with a system matrix constructed using the Surface DC-PSE operators.We test this by solving the Poisson equation on the unit circle S 1 : with Dirichlet boundary condition at one point (1, 0) conforming to the analytical solution We use the same Surface DC-PSE operators as in the previous subsection to construct the system of equations, which is then solved using the KSPGMRES solver from PETSc [3].Fig. 2b,d show the solution f and the convergence plot of the absolute error with respect to the analytical solution.
Next, we test the method in three dimensions by solving the Poisson equation on the sphere S 2 : with Dirichlet boundary condition along the great circle parallel to the y − z plane conforming to the analytical solution We solve the resulting linear system with KSPGMRES from PETSc [3] without preconditioning.The convergence plots for orders r = 2 and r = 4 are shown in Fig. 3d.The collocation points on the sphere were constructed using the Fibonacci sphere technique [13].While this generates pseudo-regular point distributions on the sphere, their average spacing does fluctuate a bit, explaining the slight waviness of the curve especially for the fourth-order operators.

Mean and Gauss curvature computation
We further verify Surface DC-PSE by computing the mean curvature H and Gauss curvature K of an ellipsoid with a = 1, b = 0.8, c = 0.75 and parameterization (u, v) We compute both curvatures from the embedded shape tensor ∇ S n, i.e., the extension of the intrinsic shape operator to the embedding space.We numerically compute this tensor in Cartesian coordinates as the 3 × 3 matrix where n = (n 1 , n 2 , n 3 ) are the components of the analytically given normal vectors at the collocation point.All intrinsic derivatives over the scalar fields n 1 , n 2 , n 3 are approximated by Surface DC-PSE operators.Mean curvature H is then computed as the trace of the embedded shape tensor and Gauss curvature K as the product of its non-zero eigenvalues (principal curvatures).These numerical computations are verified against the analytical solutions We approximate the embedded shape tensor ∇ S n using Surface DC-PSE with δn = 3.0/(N p − 1), r c = 2.9δn, N n = 2, and r = 2.The results and the convergence plot of the absolute errors are shown in Fig. 4a,c.As specified by r, we observe secondorder convergence to the analytical solution when decreasing the average spacing h between the points.The relative errors are visualized in Fig. 4b.They concentrate around extremal points of the curvature, as expected.Finally, we apply the same mean-curvature computation to an arbitrary surface with no analytical solution, the Stanford bunny from the Stanford Computer Graphics Laboratory.We use the down-sampled version of the original data set with 2960 points on the surface, obtained from https://www.stlfinder.com/model/stanford-bunny-S4kAUsKI/3091553.The result is visualized in Fig. 4d, showing that the Surface DC-PSE qualitatively works also for non-algebraic surfaces.

Comparison with Surface Finite Element Methods
We validate Surface DC-PSE by comparing it with surface Finite Element Methods (FEM) on a test case specifically developed for benchmarking FEM solutions of surface problems [1].The benchmark problem considers a two-dimensional surface S with an isolated "bump" described by the graph of the function u(x) = αζ(∥x − p∥/r), where α ≥ 0 is the height of the bump, p = (−0.5,0) is the position of its center, and r = 0.25 is the bump radius (see Fig. 5a).The function ζ(d) = exp − 1 1−d 2 for d < 0.975 and 0 otherwise is a cut-off compressed Gaussian.The surface is thus defined as This surface is asymmetric along the x direction, containing regions of both positive and negative Gauss curvature (see Fig. 5b).
The collocation points for Surface DC-PSE are regularly distributed in the flat parts of the surface, while in the region of the bump ([−0.75, −0.25] × [−0.25, 0.25]) we use the Fibonacci sphere technique [13] to place the points.This results in a total of N p = 18 439 points on the surface, of which 16 441 are in the flat parts of the surface.The resulting average point spacing is h = 0.03125.We use a normal spacing of δn = h to generate the surface operators with N n = 3 layers to either side of the surface.We choose r c = 2.9δn as the operator support and r = 2 as the desired order of convergence.The Neumann boundary conditions at the edge of the surface are imposed using the method of images [8] with around 3000 "ghost points" outside the domain.Time integration over t ∈ [0, 1] is done using the explicit fifth-order Dormand-Prince method [10] with a time-step size of δt = 10 −4 .
We qualitatively compare the solution computed by Surface DC-PSE with those obtained by Finite Element Methods using the data reported in Ref. [1] for Surface FEM (SFEM) [11], Intrinsic Surface FEM (ISFEM) [2], trace FEM (TraceFEM) [26], and diffuse interface (DI) FEM [17,30].We report the results obtained by high-resolution SFEM and DI FEM.The solutions obtained using a lower-resolution SFEM, as well as using ISFEM and TraceFEM are visually indistinguishable and reported elsewhere [1].
The comparison is done at two points on the surface above x 0 = p = (−0.5,0) and . There is no analytical solution for this test case, but a high-resolution SFEM solution serves as a reference.This "highRes SFEM" reference solution was obtained using a tentimes finer mesh than Surface DC-PSE, with the finest space resolution (on the bump) set to h ≈ 0.0027, a time step size of δt = 10 −4 , a standard BDF-2 scheme, and a polynomial order of two [1].Results are shown in Fig. 6.
We observe excellent agreement between the Surface DC-PSE solution and the "highRes SFEM" reference solution for α = 0 (dotted lines in Fig. 6).For quantitative comparison, we report the absolute difference in the peak values of the solutions, e peak .At both probe points, e peak ≈ 7 • 10 −4 for the flat case, which is comparable to the other FEM methods [1].
For the curved cases (α ∈ {1, 2}), Surface DC-PSE is still in very good agree-ment with the reference solution and closer to it than DI FEM (i.e., the L 2 difference between highRes SFEM and Surface DC-PSE is smaller in all cases than between high-Res SFEM and DI FEM).For α = 1 (dashed lines in Fig. 6), the absolute difference in the peak values between Surface DC-PSE and highRes SFEM is e peak ≈ 2.6 • 10 −3 above x 0 and e peak ≈ 1.2 • 10 −2 above x 1 .The differences for DI FEM are of the same order.For the highest bump (α = 2, solid lines in Fig. 6), we observe similar e peak values and overall behavior for all methods.Consistently, and as expected, the curves have decreasing peak values for increasing α.

Conclusions
We have presented a meshfree collocation scheme for consistently approximating intrinsic differential operators on smooth curved surfaces represented by point clouds.The present scheme is based on the DC-PSE method for discretizing differential operators on (almost) arbitrarily distributed collocation points in flat spaces [32].We have derived the present surface-intrinsic version by realizing that the kernel evaluations can be factored out across points created by exact constant orthogonal extension, and that the partial sums over the kernels can be precomputed and stored on the surface points only, defining effective surface kernels.This yields a method that is easy to implement and computationally efficient, as it only requires storing points on the surface.In this sense, Surface DC-PSE combines features from embedding methods with features from embedding-free methods.The operators are determined in an embedding formulation, but result in a surface-intrinsic algorithm for operator evaluation.
We have verified the method in different test cases with known analytical solutions.This included evaluating the Laplace-Beltrami operator on the unit circle and the unit sphere, solving Poisson equations on the unit circle and the unit sphere using an implicit solver, and computing mean and Gauss curvature of an ellipsoid via approximation of the embedded shape tensor.In all cases, the Surface DC-PSE method converged as expected.We then applied the method to compute the mean curvature of the Stanford bunny, showing an application to a non-parametric surface.We expect DC-PSE to be more robust (i.e., better conditioned) than Surface MLS in complex geometries [5] and likely computationally more efficient, since Surface MLS uses local moving frames in co-dimension 2, which is different from the present tubular approach in co-dimension 1.Finally, we validated Surface DC-PSE against two different Finite Element Methods (FEM) for surface problems, showing excellent agreement with the reference solution.
Despite its advantages, Surface DC-PSE also has a number of limitations: First, the normal field is required as an input, which can be limiting or introduce additional errors if the normals need to be numerically approximated.Second, for a given point distribution, the numerical error is limited by the curvature of the represented surface and depends on the average spacing between the surface points and the normally extended points (see Figs. 1b and 4b).The required minimum resolution can be determined based on an approximation of the curvature.Third, Surface DC-PSE is only correct for orientable surfaces that possess a non-intersecting tubular neighborhood (i.e., a tubular network) with a radius at least as large as the kernel radius everywhere.This guarantees that the line segments along which the surface points are extended never intersect.
For surfaces that possess any non-intersecting tubular neighborhood, this can always be achieved by choosing the resolution h (locally) sufficiently small, since the tube radius scales with the radius of the DC-PSE kernel.Fourth, we limited our discussion to surfaces in co-dimension 1.While it may be possible to extend Surface DC-PSE to codimension 2 (e.g., curves embedded in R 3 ), the method is likely not computationally efficient in those cases, as it would require constant orthogonal extension in the whole two-dimensional normal space of each surface point.Lastly, constructing the Surface DC-PSE kernels is computationally expensive, as it involves solving a small linear system of equations at each surface point in order to determine the embedding-space DC-PSE kernels.For Eulerian simulations, where the collocation points do not move, the kernels have to be determined once at the beginning of the simulation.However, if points move, e.g. in a Lagrangian simulation or simulations involving deforming surfaces, the kernels need to be recomputed at each time step.While the corresponding cost may be amortized by a gain in accuracy and stability [32], it is still significant.
In future work, we will consider extensions of Surface DC-PSE to Lagrangian problems involving moving and deforming surfaces.This also includes simulations of deformable surfaces, where the surface deformation itself results from intrinsic forcebalance equations [31,21].We will also consider coupling Surface DC-PSE with regular DC-PSE in the surrounding space in order to describe coupled bulk-surface phenomena, as well as adaptive-resolution Surface DC-PSE with resolution and tube radius locally adjusted to the curvature of the surface in order to maintain error equidistribution [28].
In summary, we have extended the meshfree collocation method DC-PSE to scalarvalued problems on curved surfaces, requiring only intrinsic surface points.Like DC-PSE, also Surface DC-PSE computes the operator kernels numerically at runtime and is consistent for any desired order of convergence r.This makes the presented algorithm particularly attractive for higher-order intrinsic operators, such as the fourth-order operators in Fig. 3, and for determining the system matrices of implicit equations on surfaces or implicit time integration schemes.

Figure 1 :
Figure1: (a) Illustration of the DC-PSE method.The collocation points x q (blue) within the symmetric operator support N (x p ) of radius r c around the center point x p (red) are used to approximate the differential operator at x p .The average distance between points in the operator support, h(x p ) defines the accuracy of the approximation.(b) Illustration of the Surface DC-PSE method.The intrinsic differential operator at a point x p (red) on the surface S is evaluated over neighboring surface points x s (blue).The operator kernel is constructed in the tubular neighborhood of radius r c (circle) with 2N n points x n (yellow) replicated along the local surface normal n at distances δn.Point labels are used in the main text for the illustrative example.

Algorithm 1 1 . 4 : 6 : 12 :
Surface DC-PSE: construction of neighborhood of a point p.Input: Point set P on the surface S 2. Cutoff radius for the operator support r c 3. Indices N S of surface points in the neighborhood of p 4. Optional: spacing δn between the normally extended points.Default: average embedding-space distance between surface points 5. Optional: Number of normal copies of each surface point to be used during operator construction N n (symmetric to either sides of the surface).Default: N n = round(r c /δn) Output: List of distances between point p and all surface s i and normal n i points in its neighborhood N (x p ): N dist (x p ) Require: |n si | = |n p | = 1 1: N dist = [ ], k = 0 2: for all s i ∈ N S do 3: N dist .append([]) for i ∈ [−N n , N n ] do 5: d ni = x p − x si − iδn • n si if |d ni | ≤ r c then 7: N dist [k].append(d ni /ϵ p ) ▷ Add to set of the corresponding s i .N dist .append([]) 13: for i ∈ [−N n , N n ]\{0} do ▷ Normal points to p.

14 : d ni = −iδn • n p 15 :Figure 2 :
Figure 2: Visualization and convergence of the Laplace-Beltrami operator on the unit circle.(a) Visualization of ∆ S (sin(θ) + cos(θ)) computed using second-order accurate Surface DC-PSE operators.(b) Visualization of the solution of the Poisson equation ∆ S f = 4 sin(2θ) solved using second-order accurate Surface DC-PSE operators.(c) Convergence plot of the Laplace-Beltrami operator in (a).L ∞ (•) and L 2 (♦) norms of the absolute errors are computed against the analytical solution in Eq. (15) for increasing numbers of points on the circle.(d) Convergence plot of the Poisson equation solution in (b).L ∞ (•) and L 2 (♦) norms of the absolute errors are computed against the analytical solution in Eq.(19) for increasing numbers of points on the circle.

Figure 4 :
Figure 4: Gauss and mean curvature computation using Surface DC-PSE.(a) Visualization of the mean and Gauss curvatures of an ellipsoid, numerically computed as the trace of the embedded shape tensor −0.5∇ S • n = −0.5Tr(∇S n) of the surface normal vector field n and the product of the non-zero eigenvalues of ∇ S n, respectively, using second-order accurate Surface DC-PSE operators.(b) Visualization of the relative errors in the computed curvatures from (a) on 32258 surface points in comparison with the analytical solution in Eq. (24).They range from about 10 −7 to 10 −3 .(c) Convergence plot of the mean and Gauss curvature computations.L ∞ (•, +) and L 2 (♦, ×) absolute errors for mean and Gauss curvatures, respectively, are computed against the analytical solutions in Eq. (24) for decreasing average collocation point spacing h.(d) Visualization of the mean curvature computed on the Stanford bunny with 2960 points using second-order accurate Surface DC-PSE operators.