Solving Partial Differential Equations on Evolving Surfaces via the Constrained Least-Squares and Grid-Based Particle Method

This section briefly summarizes the Grid Based Particle Method (GBPM) proposed and developed in a series of studies [18, 17, 16]. We refer the interested readers to the reference and thereafter. The overall algorithm of the GBPM is given as follows, and a graphical illustration is given in Figure 1.

In the GBPM, we represent the interface by meshless particles associated with an underlying Eulerian mesh. Each sampling particle on the interface is chosen to be the closest point from each underlying grid point in a small neighborhood of the interface. This one-to-one correspondence gives each particle an Eulerian reference during the evolution.

In the first step, we define an initial computational tube for active grid points and use their corresponding closest points as the sampling particles for the interface. A grid point is called active if its distance to the interface is smaller than a given tube radius, and we label the set containing all active grids $\Gamma$. We associate the corresponding closest point on the interface to each active grid point. This particle is called the footpoint associated with this active grid point. Furthermore, we can also compute and store certain Lagrangian information of the interface at the footpoints, including normal, curvature, and parametrization, which will be helpful in various applications.

This representation is illustrated in Figure 1(a) using a circular manifold as an example. We plot the underlying mesh in a solid line, all active grids using small circles, and their associated footpoints on the manifold using squares. To each grid point near the interface (blue circles), we associate a footpoint on the interface (red squares). The relationship of each of these pairs is shown by a solid line link. To track the motion of the interface, we move all the sampling particles according to a given motion law. This motion law can be very general. Suppose the interface is moved under an external velocity field. Since we have a collection of particles on the interface, we move these points just like all other particle-based methods, which is simple and computationally efficient. We can solve a set of ordinary differential equations using a high-order scheme, which gives a very accurate interface location. For more complicated motions, the velocity may depend on the geometry of the interface. This can be done through local interface reconstruction.

It should be noted that a footpoint after motion may not be the closest point on the interface to its associated active grid point anymore. For example, Figure 1 (b) shows the location of all particles on the interface after the constant motion $\textbf{u}=(1,1)^{T}$ with a small time step. As we can see, these particles on the interface are no longer the closest point from these active grid points to the interface. More importantly, the motion may cause those original footpoints to become unevenly distributed along the interface. This may introduce stiffness when particles are moving toward each other and a large error when particles are moving apart. To maintain a quasi-uniform distribution of particles, we need to resample the interface by recomputing the footpoints and updating the set of active grid points ($\Gamma$) during the evolution. During this resampling process, we locally reconstruct the interface, which involves communications among different particles on the interface. This local reconstruction also provides geometric and Lagrangian information at the recomputed footpoints on the interface. The key step in the method is a least squares approximation of the interface using polynomials at each particle in a local coordinate system, $\{(\textbf{n}^{\prime})^{\perp},\textbf{n}^{\prime}\}$, with the associated particle as the origin. Using this local reconstruction, we find the closest point from this active grid point to the local approximation of the interface. This gives the new footpoint location. Further, we also compute and update any necessary geometric and Lagrangian information, such as the normal vector, curvature, and possibly an updated interface parametrization at this new footpoint.

We approximate the surface locally using a polynomial from each grid point $\mathbf{p}$. For example, in the three dimensions, we could use a second-order polynomial given by

$$h(x,y)=a_{0,0}+a_{1,0}x+a_{0,1}y+a_{2,0}x^{2}+a_{1,1}xy+a_{0,2}y^{2}\,,$$

where $(x,y,h(x,y))$ are coordinates in the local coordinate system, $\{(\textbf{n}^{\prime})^{\perp},\textbf{n}^{\prime}\}$. The original idea in [18] fits the polynomial using a collection of $(n+1)$-footpoints, $(x_{i},y_{i},h(x_{i},y_{i}))$ with $i=0,1,\cdots,n$, in the neighborhood of $\mathbf{p}$. With these points, we determine the least squares solution of the linear system:

$$\begin{pmatrix}1&x_{0}&y_{0}&x_{0}^{2}&x_{0}y_{0}&y_{0}^{2}\\ 1&x_{1}&y_{1}&x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}\\ 1&x_{2}&y_{2}&x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 1&x_{n}&y_{n}&x_{n}^{2}&x_{n}y_{n}&y_{n}^{2}\end{pmatrix}\begin{pmatrix}a_{0,0}\\ a_{1,0}\\ a_{0,1}\\ a_{2,0}\\ a_{1,1}\\ a_{0,2}\end{pmatrix}=\begin{pmatrix}h(x_{0},y_{0})\\ h(x_{1},y_{1})\\ h(x_{2},y_{2})\\ \vdots\\ h(x_{n},y_{n})\end{pmatrix}\,.$$

This set of footpoints in the local reconstruction must be chosen carefully. As discussed in the original GBPM formulation, this set of points has to satisfy the following conditions. The first one is that these points have to be distinct. In the previous formulation, we have conditioned that these footpoints must be separated by at least a certain distance on the local tangent plane. This constraint avoids putting an extensive amount of weight on these sampling points. The second constraint is that we would like to sample the same segment on the interface. When the curvature of the interface is too large compared to the underlying mesh size, one might collect footpoints from different interface segments, which will result in a wrong local reconstruction of the surface. In the original work, we proposed checking if the associated normal vectors from these footpoints were pointing in a similar direction. In particular, for each of these footpoints, $\mathbf{x}_{i}$ has an associated vector $\mathbf{n}^{\prime}_{i}$ corresponding to the surface normal before the motion. We have checked if $\mathbf{n}^{\prime}_{i}\cdot\mathbf{n}^{\prime}_{j}>\cos\theta_{\max}$ for some angle $\theta_{\max}$.

The above strategies are in some sense imposing some hard constraints in the selection phase. We either use the neighboring sampling point or eliminate it when choosing the subset of footpoints in the local reconstruction step. This paper proposes relaxing this by imposing some soft constraints. The first condition on the well-separation can be partly handled by the moving least squares (MLS) method [15], which introduces a weighting to the footpoints based on the distance between points on the tangent plane. We will not concentrate on this part. To avoid information from different segments on the surface, instead of simply getting rid of those points pointing at an angle larger than $\theta_{\max}$, we propose to add a weighting in the least squares further to emphasize the contribution from footpoints points with similar normal directions. Mathematically, we let $\mathbf{x}_{0}=(x_{0},y_{0},h(x_{0},y_{0}))$ be the closest footpoint to the grid point $\mathbf{p}$ and denote those normal vectors at the footpoints before motion by $\mathbf{n}^{\prime}_{i},i=0,1,\cdots,n$. We now introduce a weight $c_{i}$ to each collected footpoints given by

$$c_{i}=\left\{\begin{array}[]{cc}\mathbf{n}^{\prime}_{0}\cdot\mathbf{n}^{\prime}_{i}&\mbox{ if $\mathbf{n}^{\prime}_{0}\cdot\mathbf{n}^{\prime}_{i}>\cos\theta_{\max}$,}\\ 0&\mbox{ otherwise.}\end{array}\right.$$

(3)

Once we have this set of weights, the $i$-th equation in the least squares system will be modified to

$$c_{i}\begin{pmatrix}1&x_{i}&y_{i}&x_{i}^{2}&x_{i}y_{i}&y_{i}^{2}\end{pmatrix}\begin{pmatrix}a_{0,0}\\ a_{1,0}\\ a_{0,1}\\ a_{2,0}\\ a_{1,1}\\ a_{0,2}\end{pmatrix}=c_{i}\,h(x_{i},y_{i})\,.$$

Let $\Gamma$ be a regular surface in $\mathbb{R}^{2}$ parametrized by $\mathbf{x}$. Mathematically, the surface gradient of a smooth function $\xi$ is given by

$$\nabla_{\Gamma}\xi=\frac{\xi_{u}G-\xi_{v}F}{EG-F^{2}}\mathbf{x}_{u}+\frac{\xi_{v}E-\xi_{u}F}{EG-F^{2}}\mathbf{x}_{v},$$

where $E,F$ and $G$ are the coefficients of the first fundamental form of the regular surface $\Gamma$, the quantities $\xi_{u}=\frac{\partial}{\partial u}\xi(\mathbf{x}(u,v))$ and $\xi_{v}=\frac{\partial}{\partial v}\xi(\mathbf{x}(u,v))$. For any local vector field $X=A\mathbf{x}_{u}+B\mathbf{x}_{v}$, we have the surface divergence defined as

$$\nabla_{\Gamma}\cdot X=\frac{1}{\sqrt{EG-F^{2}}}\left[\frac{\partial}{\partial u}\left(A\sqrt{EG-F^{2}}\right)+\frac{\partial}{\partial v}\left(B\sqrt{EG-F^{2}}\right)\right]\,.$$

The surface Laplacian (Laplace-Beltrami) operator $\Delta_{\Gamma}\xi$ is defined as $\Delta_{\Gamma}\xi=\nabla_{\Gamma}\cdot(\nabla_{\Gamma}\xi)$. In the local representation, it can be expressed as

	$\displaystyle\Delta_{\Gamma}\xi$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{EG-F^{2}}}\left[\frac{\partial}{\partial u}\left(\frac{G}{\sqrt{EG-F^{2}}}\xi_{u}\right)-\frac{\partial}{\partial u}\left(\frac{F}{\sqrt{EG-F^{2}}}\xi_{v}\right)\right.$
			$\displaystyle+\left.\frac{\partial}{\partial v}\left(\frac{E}{\sqrt{EG-F^{2}}}\xi_{v}\right)-\frac{\partial}{\partial v}\left(\frac{F}{\sqrt{EG-F^{2}}}\xi_{u}\right)\right]\,.$

Therefore, to design a numerical method for computing any differential operator on the surface, we have to approximate the first fundamental form of the surface and also provide a way to approximate the derivatives in the local coordinate systems.

The first issue is to design a good local representation of surface $\Gamma$. One natural representation is the local tangential projection as shown in Figure 2. We assume that a set of sampling points represents the interface. These sampling points can be obtained from various approaches. For example, it can be the GBPM representation we summarized above, a typical point cloud dataset, or an interface represented by triangulation. The representation itself does not play a crucial role in the following discussions. Assume further that for any individual sampling point $\mathbf{y}^{*}$, we can collect a set of $M$-neighboring sampling points denoted by $\{\mathbf{y}_{j}:j=1,\cdots,M\}$. It is certainly easier to have the interface represented by a triangulation since we have the connectivities of all vertices (the sampling points) on the interface. For the GBPM representation, we can collect these neighboring points through the underlying mesh. Some KNN algorithms can also be applied to point cloud datasets.

For the sampling point $\mathbf{y}^{*}$, we let $\mathbf{n}$ be the corresponding (approximated) normal vector of the surface. With $\mathbf{n}$ and its tangent plane, we introduce a new local coordinate system $(x,y,z)$ with the $z$-direction representing the component in the $\mathbf{n}$-direction and $\mathbf{y}^{*}$ becomes $(\mathbf{x}_{0},0)$ in the new coordinate system. In this local coordinate system, a function can locally approximate the surface, and we denote it as $h(x_{j},y_{j})=z_{j}$. Once we have this local approximation $h$, we can approximate the derivatives $h_{x},h_{y},h_{xx},h_{xy}$ and $h_{yy}$ for the first fundamental form in the surface differential operators.

Numerically, this leads to one of the following problems. The first one is an interpolation problem: Given $(x_{j},y_{j},z_{j})$ for $j=1,\cdots,n$, determine an interpolant $h(x,y)$ which passes through all data points such that $(x_{j},y_{j},h(x_{j},y_{j}))$. One idea is to use the polynomial interpolation. For example, when the neighbors are carefully chosen such that $n=6$, we can use a second-order polynomial to interpolate these data points. This approach better estimates the local error if the interpolation uses the polynomial basis. However, we generally have no control over the distribution of $(x_{j},y_{j})$ from the projection of the sampling points onto the local tangent plane. In particular, some of these projection points could be very close to each other or could collide on the tangent plane. Designing a good subset of points might be challenging, which could lead to a good interpolation problem. Therefore, we do not consider this class of approach further.

The second is the least squares approximation problem: Given $(x_{j},y_{j},z_{j})$ for $j=1,\cdots,n$, determine a function $h(x,y)$ that best explains the data points. Mathematically, we obtain the set of parameters $\alpha_{i}$ to minimize the mismatch between the data points and $\sum_{\psi_{i}\in\mathcal{B}}\alpha_{i}\psi_{i}(\mathbf{x})$ where $\mathcal{B}$ is the space for the basis, i.e. we solve

$$\min_{\alpha_{i}}\sum_{j=0}^{M}\left[h(\mathbf{x}_{j})-\sum_{\psi_{i}\in\mathcal{B}}\alpha_{i}\psi_{i}(\mathbf{x}_{j})\right]^{2}\,.$$

The choice of the basis is clearly not unique. For example, the original GBPM developed in [18, 17, 16] chose the standard polynomial basis for $\mathcal{B}$. Of course, one could also replace the standard least squares approximation with a weighted least squares approach so that those points far away from $\mathbf{x}_{0}$ do not contribute equally like the point $\mathbf{x}_{0}$ itself.

In a recent work [37, 38], we have developed a general strategy to impose a hard constraint in the least squares reconstruction by enforcing the least squares reconstruction passing through $(\mathbf{x}_{0},h(\mathbf{x}_{0}))$ exactly. We call this a CLS approach. Under this constraint, the least squares method can be modified to fit the function value $h(\mathbf{x})-h(\mathbf{x}_{0})$ with functions in the form of $\sum_{\psi_{i}\in\mathcal{B}}\beta_{i}[\psi_{i}(\mathbf{x})-\psi_{i}(\mathbf{x}_{0})]$. This only requires a simple modification of the standard least squares approximation method. We will give out more information below. Let $\mathcal{L}$ be some linear differential operator. We can also follow a similar idea to obtain an approximation $\mathcal{L}(h(\mathbf{x}_{0}))$ based on the CLS, i.e.

$$\mathcal{L}(h(\mathbf{x}_{0}))=\sum_{j=1}^{M}w_{j}(h(\mathbf{x}_{j})-h(\mathbf{x}_{0}))=\sum_{j=0}^{M}w_{j}h(\mathbf{x}_{j})\,,$$

with the coefficient $w_{0}=-\sum_{j=1}^{M}w_{j}$.

The choice of the basis is not unique. In this paper, we present two different sets of bases and give the technical implementation details. The first is the simple polynomial basis, which leads to an approximation result similar to the original GBPM method for PDEs but with an extra hard constraint at the target location during reconstruction. The second basis follows our previous development in [37, 38], which uses the radial basis functions. A different basis will give slightly different results. We are not trying to conclude which of the following two basis sets will yield better numerical results in all applications. Still, we will only show some examples to explain better the framework that we are proposing.

To end this section, we summarize our overall approach to solving PDEs on evolving surfaces. The original GBPM handles both the interface sampling and the interface motion. At the same time, the proposed CLS methods help to approximate the surface, the local reconstructions of the function value $u$ on the interface, and the coefficients in the first fundamental form for the surface differential operator.

Our first example is the motion of a sphere under a single vortex flow [10, 12, 25, 18]. This is a challenging test case that can test the stability of our GBPM with the modified local reconstruction. The velocity field is given by

	$\displaystyle v_{1}(x,y,z)$	$\displaystyle=$	$\displaystyle 2\sin(\pi x)^{2}\sin(2\pi y)\sin(2\pi z)\cos\left(\frac{\pi t}{T}\right)$
	$\displaystyle v_{2}(x,y,z)$	$\displaystyle=$	$\displaystyle-\sin(2\pi x)\sin(\pi y)^{2}\sin(2\pi z)\cos\left(\frac{\pi t}{T}\right)$
	$\displaystyle v_{3}(x,y,z)$	$\displaystyle=$	$\displaystyle-\sin(2\pi x)\sin(2\pi y)\sin(\pi z)^{2}\cos\left(\frac{\pi t}{T}\right)\,,$

for some given final time $T$. The initial shape is a sphere centered at $(0.35,0.35,0.35)$ with radius $r_{0}=0.15$. In Figure 4, we compute the solution at $t=0,0.375,0.75$ and $1.5$ with $T=1.5$ using $\Delta x=0.01$. The mean distance from the final footpoints to the point $(0.35,0.35,0.35)$ at the final time is $0.1506$ with the standard deviation $5.114\times 10^{-4}$. The largest absolute error in the location is $2.862\times 10^{-3}$. There is a total of $8530$ initially at $t=0$. The number gradually increases to the maximum value $13801$ when $t=0.75$, finally dropping to $8615$ at the final time. Figure 5 shows the detailed evolution of footpoints and the error in the sphere’s radius $r-r_{0}$ at the final time.

To see the effect of the new weighting (3) introduced in the local reconstruction, we compare two numerical solutions in Figure 6. In this example, we repeat the similar setup but choose a slightly smaller final time with $T=1$. Figure 6 shows the changes in the number of footpoints and the solution at the final time $T=1$. For the original scheme without any weighting in the local reconstruction, we found that the mean distance from the final footpoints to the center at the final time is $0.1557$ with the standard deviation $5.832\times 10^{-3}$. The modified local reconstruction gives a more accurate solution with the corresponding mean distance $0.1502$ with the standard deviation $2.247\times 10^{-4}$.

In this section, we give two numerical examples to show that the proposed CLS method works well in solving PDEs on surfaces sampled by the GBPM representations. The first example is the Cahn-Hilliard equation on the unit sphere. The surface Cahn-Hilliard equation has the form

$\displaystyle u_{t}=\frac{\nu}{Pe}\Delta_{\Gamma}\left[g^{\prime}(u)\right]-\frac{\nu\,Cn^{2}}{Pe}\Delta_{\Gamma}^{2}u\,.$

(6)

where $Cn$ is the Cahn-Hilliard number, $Pe$ is the Peclet number and $\nu$ is a diffusion parameter. In this example, we solve the equation on a unit sphere with the parameters in the equation given by $Pe=\nu=1$, $Cn=0.06$, and the double well potential function $g(u)=u^{2}(1-u)^{2}$. Our numerical solution at different times is shown in Figure 7. The initial condition is a uniform random perturbation with a magnitude of 0.01, with an average concentration of $u$ of 0.5.

In the second example, we simulate the Turing patterns on a static surface. Consider the following coupled system of equation (2) given by

	$\displaystyle\frac{Du}{Dt}+u\nabla_{\Gamma}\cdot\mathbf{v}-\mathcal{D}_{1}\Delta_{\Gamma}u$	$\displaystyle=$	$\displaystyle f_{1}(u,w)$
	$\displaystyle\frac{Dw}{Dt}+w\nabla_{\Gamma}\cdot\mathbf{v}-\mathcal{D}_{2}\Delta_{\Gamma}w$	$\displaystyle=$	$\displaystyle f_{2}(u,w)\,,$		(7)

with

$\displaystyle f_{1}(u,w)=\gamma(a-u+u^{2}w)\,\mbox{ and }\,f_{2}(u,w)=\gamma(b-u^{2}w)\,.$

(8)

The parameters $\gamma$, $a$, and $b$ are all positive constants. In this example, we follow [1] and use the following parameters for the kinetics function in (8): $a=0.1$, $b=0.9$, $\mathcal{D}_{1}=1$ and $\mathcal{D}_{2}=10$. All three examples use the same initial condition given by $(u,w)=(1,0.9)$ with a uniform random perturbation of magnitude 0.01. The parameter $\gamma$ is given in different scales. The grid size in the GBPM representation is $\Delta x=0.1$, and the stopping time is $t=10$.

Figure 8 shows our numerical solutions corresponding to various $\gamma$ at the fixed final time $t=10$. Even though the exact solution to these examples is not available, these results are similar to those obtained in various articles, including [3, 24, 1]. Numerically, the CLS method is flexible and matches the GBPM well in solving PDEs on a static surface.

This three-dimensional example is taken from [7, 8, 28] to test the accuracy of a PDE solver on moving surfaces. The time-dependent oscillating ellipsoid has the following exact solution

$$\Gamma(t)=\left\{\mathbf{x}=(x,y,z)\in\mathbb{R}^{3}\left|\frac{x^{2}}{a(t)}+y^{2}+z^{2}=1\right.\right\}\,,$$

where $a(t)=1+\sin 2t$. The surface $\Gamma$ can be explicitly parametrized as

$$\Gamma(t,\theta,\phi):=\left(\sqrt{a(t)}\cos\theta\cos\phi,\sin\theta\cos\phi,\sin\phi\right)\,.$$

The associated velocity field is given by $\mathbf{v}(\mathbf{x},t)=\frac{a^{\prime}(t)}{a(t)}(x,0,0)^{T}$. Therefore, the equation involves both the advection along the surface and the deformation by the motion in the normal direction. We choose the exact continuous solution $u(\mathbf{x},t)=e^{-6t}xy$ so that the corresponding source term on the right-hand side of the equation is given by

$\displaystyle f(\mathbf{x},t)=\left\{-6+\frac{a^{\prime}(t)}{a(t)}\left(1-\frac{x_{1}^{2}}{2N}\right)+\frac{1+5a(t)+2a^{2}(t)}{N}-\frac{1-a(t)}{N^{2}}\left[x_{1}^{2}+a^{3}(t)\left(x_{2}^{2}+x_{3}^{2}\right)\right]\right\}u(\mathbf{x},t)$

with $N=x^{2}+a^{2}(t)(y^{2}+z^{2})$. In this simulation, we choose the mesh size $\Delta x=0.1$ and a fixed time step size $\Delta t=0.1\,(\Delta x)^{2}$. We plot our solution at different time $t=\pi/8,\pi/4,3\pi/8$ and $\pi/2$ in Figure 9. Figure 10 gives the $L^{\infty}$-error in the numerical solution computed using the CLS with different mesh sizes at different times. Table 1 shows the errors in the computed solution at different times $t=0.01,0.02,0.03$ and $0.04$ computed on a different mesh. We see that our numerical approach is roughly second-order accurate.

Concerning the computational efficiency, the CPU time for one marching step on the mesh with $\Delta x=0.1$ is $3.1$ seconds, including $1.5$ seconds spent on the CLS reconstruction and $1.6$ seconds on the GBPM computations. This computational time is recorded under a 4-core laptop in MATLAB with a parallel modulo implemented.

In this case, the interface velocity is strongly coupled with the solution of the PDE. This example is taken from [28], where we evolve an initial torus

$$\left(1-\sqrt{x^{2}+y^{2}}\right)^{2}+z^{2}=0.3^{2}$$

according to the given velocity $\mathbf{v}=(0.1\kappa+5u)\mathbf{n}$ where $\kappa$ is the mean curvature, $\mathbf{n}$ is the unit normal vector, and $u$ is the solution to the PDE

$\displaystyle u_{t}+V\frac{\partial u}{\partial\mathbf{n}}-\kappa Vu+\Delta_{\Gamma}u=0\text{\; on\; }\Gamma\,.$

The initial value of $u$ is given by $u(\mathbf{x},t)=1+20xyz$. Figure 11 shows our computed solution of $u$ on the evolving surface at different times from $0.02$ to $0.08$ in an increment of $0.02$. Even though the exact solution to this problem is unavailable, these numerical solutions are similar to those presented in [28].