Isogeometric Analysis for Surface PDEs with Extended Loop Subdivision

TL;DR

This paper introduces an isogeometric analysis method for surface PDEs using extended Loop subdivision, achieving high accuracy and efficiency on complex surfaces, and outperforming traditional finite element methods.

Contribution

The paper develops a novel isogeometric analysis approach based on extended Loop subdivision, providing exact geometry representation and improved accuracy for surface PDEs.

Findings

Second-order accuracy in L2-norm demonstrated

Outperforms standard linear finite element methods

Effective on both open and closed surfaces

Abstract

We investigate the isogeometric analysis for surface PDEs based on the extended Loop subdivision approach. The basis functions consisting of quartic box-splines corresponding to each subdivided control mesh are utilized to represent the geometry exactly, and construct the solution space for dependent variables as well, which is consistent with the concept of isogeometric analysis. The subdivision process is equivalent to the $h$-refinement of NURBS-based isogeometric analysis. The performance of the proposed method is evaluated by solving various surface PDEs, such as surface Laplace-Beltrami harmonic/biharmonic/triharmonic equations, which are defined on different limit surfaces of the extended Loop subdivision for different initial control meshes. Numerical experiments demonstrate that the proposed method has desirable performance in terms of the accuracy, convergence and computational cost for solving the above surface PDEs defined on both open and closed surfaces. The proposed approach is proved to be second-order accuracy in the sense of $L^2$-norm by theoretical and/or numerical results, which is also outperformed over the standard linear finite element by several numerical comparisons.

Figures

Equations

$$T_{\mathbf{x}}\mathcal{S}=\mathrm{span}\{{\mathbf{x}}_u,{\mathbf{x}}_v\}\mathrm{and}\mathbf{n}=\frac{{\mathbf{x}}_u\times {\mathbf{x}}_v}{\Vert {\mathbf{x}}_u\times {\mathbf{x}}_v\Vert },$$

$$I=⟨\mathrm{d}\mathbf{x},\mathrm{d}\mathbf{x}⟩=g_{11}\mathrm{d}u\mathrm{d}v+2g_{12}\mathrm{d}u\mathrm{d}v+g_{22}\mathrm{d}u\mathrm{d}v,$$

$$II=-⟨\mathrm{d}\mathbf{x},\mathrm{d}\mathbf{n}⟩=b_{11}\mathrm{d}u\mathrm{d}u+2b_{12}\mathrm{d}u\mathrm{d}v+b_{22}\mathrm{d}v\mathrm{d}v$$

$$g=\det [g_{\alpha \beta }],[g^{\alpha \beta }]=[g_{\alpha \beta }{]}^{-1}\mathrm{and}b=\det [b_{\alpha \beta }]\text{for}\alpha ,\beta =1,2.$$

$${\int }_{\mathcal{S}}\mathrm{d}A={\iint }_{\Omega }|{\mathbf{x}}_u\wedge {\mathbf{x}}_v|\mathrm{d}u\mathrm{d}v$$

$$\mathcal{W}({\mathbf{x}}_u)=-{\mathbf{n}}_u\mathrm{and}\mathcal{W}({\mathbf{x}}_v)=-{\mathbf{n}}_v.$$

$$H:=\frac{k_1+k_2}{2}=\frac{\mathrm{tr}(S)}{2}=\frac{b_{11}{g}_{22}-2b_{12}{g}_{12}+b_{22}{g}_{11}}{2g}.$$

$$f_u=\frac{\partial f(\mathbf{x}(u,v))}{\partial u}\mathrm{and}{f}_v=\frac{\partial f(\mathbf{x}(u,v))}{\partial v}.$$

$${\nabla }_{\mathcal{S}}f(\mathbf{x})=[{\mathbf{x}}_u,{\mathbf{x}}_v][g^{\alpha \beta }][f_u,f_v{]}^{\mathrm{T}}\in {\mathbb{R}}^3.$$

$$\nabla f(\mathbf{x})=[f_u,f_v{]}^{\mathrm{T}}=[{\mathbf{x}}_u,{\mathbf{x}}_v{]}^{\mathrm{T}}{\nabla }_{\mathcal{S}}f.$$

$${\mathrm{div}}_{\mathcal{S}}(\mathbf{v})=\frac{1}{\sqrt{g}}\left[\frac{\partial }{\partial u},\frac{\partial }{\partial v}\right]\left[\sqrt{g}[g^{\alpha \beta }][{\mathbf{x}}_u,{\mathbf{x}}_v{]}^{\mathrm{T}}\mathbf{v}\right].$$

$${\Delta }_{\mathcal{S}}f={\mathrm{div}}_{\mathcal{S}}({\nabla }_{\mathcal{S}}f).$$

$${\Delta }_{\mathcal{S}}f=\frac{1}{g}(g_{22}{f}_{11}+g_{11}{f}_{22}-2g_{12}{f}_{12}).$$

$${\mathbf{x}}_{\beta }^{-1}\circ {\mathbf{x}}_{\alpha }:{\Omega }_{\alpha \beta }={\mathbf{x}}_{\alpha }^{-1}({\Omega }_{\alpha }\cap {\Omega }_{\beta })\subset {\Omega }_{\alpha }⟶{\Omega }_{\beta \alpha }={\mathbf{x}}_{\beta }^{-1}({\Omega }_{\alpha }\cap {\Omega }_{\beta })\subset {\Omega }_{\beta }$$

$${\int }_{\Omega }⟨\mathbf{v},{\nabla }_{s}f⟩+f{\mathrm{div}}_{\mathcal{S}}(\mathbf{v})\mathrm{d}A={\int }_{\partial \Omega }f⟨\mathbf{v},{\mathbf{n}}_c⟩\mathrm{d}s.$$

$$C_k(\mathcal{S})=\left\{f\in C^{\infty }(\mathcal{S}):{\int }_{\mathcal{S}}|{\nabla }^{j}f{|}^2\mathrm{d}A\le \infty \text{for}j=0,\dots ,k\right\}.$$

$$\Vert f{\Vert }_{H^k(\mathcal{S})}:={\left(\sum _{j=0}^k{\int }_{\mathcal{S}}|{\nabla }^{j}f{|}^2\mathrm{d}A\right)}^{\frac{1}{2}}.$$

$$C_k(\mathcal{S})=C^{\infty }(\mathcal{S})\subset C^k(\mathcal{S})\subset H^k(\mathcal{S}).$$

$${\mathbf{p}}_i^{k+1}=\frac{1}{8}(3{\mathbf{p}}_0^k+3{\mathbf{p}}_i^k+{\mathbf{p}}_{i-1}^k+{\mathbf{p}}_{i+1}^k),i=1,\dots ,n,$$

$${\mathbf{p}}_0^{k+1}=(1-n\alpha ){\mathbf{p}}_0^k+\alpha \sum _{j=1}^n{\mathbf{p}}_j^k,$$

$$\alpha =\frac{1}{n}\left[\frac{5}{8}-{\left(\frac{3}{8}+\frac{1}{4}\cos \frac{2\pi }{n}\right)}^2\right]$$

$${\mathbf{p}}_0^{\infty }=(1-nl){\mathbf{p}}_0^0+l\sum _{j=1}^n{\mathbf{p}}_j^0,l=\frac{1}{n+3/(8\alpha )},$$

$$\mathbf{p}(\xi ,\eta )=\sum _{i=1}^{12}{B}_i(\xi ,\eta ){\mathbf{p}}_i^k,$$

$${\mathbf{x}}^{k+1}=S{\mathbf{x}}^k=S^k{\mathbf{x}}^1,k=1,2,\dots .$$

$${\mathbf{x}}^1={\mathbf{e}}_0{\mathbf{a}}_0^{\mathrm{T}}+{\mathbf{e}}_1{\mathbf{a}}_1^{\mathrm{T}}+{\mathbf{e}}_2{\mathbf{a}}_2^{\mathrm{T}}+\cdots +{\mathbf{e}}_n{\mathbf{a}}_n^{\mathrm{T}},{\mathbf{a}}_j\in {\mathbb{R}}^3,$$

$${\mathbf{x}}^{k+1}={\lambda }_0^k{\mathbf{e}}_0{\mathbf{a}}_0^{\mathrm{T}}+{\lambda }_1^k{\mathbf{e}}_1{\mathbf{a}}_1^{\mathrm{T}}+{\lambda }_2^k{\mathbf{e}}_2{\mathbf{a}}_2^{\mathrm{T}}+\cdots +{\lambda }_n^n{\mathbf{e}}_n{\mathbf{a}}_n^{\mathrm{T}},$$

$$(1-n_i{l}_i)\mathrm{u}({\mathbf{p}}_i^k)+l_i\sum _{j=1}^{n_i}\mathrm{u}({\mathbf{p}}_j^k)=u({\mathbf{p}}_i^{\infty }),i=1,\dots ,m,$$

$$\Vert u-u_I{\Vert }_{H^l(\mathcal{S})}\le c{h}^{2-l}\Vert u{\Vert }_{H^2(\mathcal{S})},l=0,1.$$

$$\left\{\begin{array}[]{l}-\Delta_{\cal S}u=f,\quad\text{in}~{}{\cal S},\\[5.0pt] u=0,\quad\text{on}~{}\partial{\cal S},\end{array}\right.$$

$$\left\{\begin{array}[]{l}{\rm Find}\ u\in H^{1}_{0}({\cal S})\ {\rm such\ that}\\[5.0pt] \displaystyle\int_{\cal S}\nabla_{\cal S}u\cdot\nabla_{\cal S}\varphi\ {\rm d}A=\int_{\cal S}f\varphi\ {\rm d}A,\quad\forall\varphi\in H^{1}_{0}({\cal S}).\end{array}\right.$$

Full text

Keywords:

Isogeometric Analysis, Extended Loop Subdivision, Surface PDEs

1 Introduction

The isogeometric analysis (IGA) was introduced by Hughes et al. [1, 2] to replace the traditional finite elements by Non-Uniform Rational B-Splines (NURBS) [3, 4] or T-splines [5, 6, 7, 8]. The concept of IGA shows great potential in developing the seamless integration between computer-aided design (CAD) and finite element method (FEM), which possesses higher accuracy than traditional FEM. The $h$-refinement, $p$-refinement and even $k$-refinement can be easily implemented by knot insertion and/or order elevation to improve the simulation accuracy without changing the geometry, which vanish the necessity of subsequent communication with CAD system. From the analytical point of view, hierarchical bases are preferred. Several constructions of such spaces have been proposed, the most common of which are (Truncated) Hierarchical B-Splines (THB) [9], Locally Refined (LR)-Splines [10, 11] and Polynomial/Rational Splines over Hierarchical T-Meshes (PHT/RHT splines) [12, 13, 14, 15, 16]. A review including an efficient IGA code can be found in [17]. A computation framework is reused in IGA on a set of three-dimensional models with similar semantic features in [18].

On the other hand, surface PDEs have been attained extensive attentions in a variety of subjects, such as structural mechanics, biology, fluid dynamics, and image processing, etc. [19, 20, 21, 22, 23, 24, 25, 26, 27]. A comprehensive review about this topic can be founded in [32]. Although the FEMs were developed to solve such kind of problems (see, for instance, [28, 29, 30, 31]), which cannot possess the ability to represent the geometry exactly [2]. To overcome this deficiency, the NURBS-based IGA was proposed to solve the high-order surface PDEs in [32]. In this article, we will further study the subdivision-based IGA, and use this method to solve surface PDEs.

Subdivision is a powerful technique in surface modeling, which provides a simple and efficient method to generate smooth surfaces with arbitrary topology structures. Moreover, it is capable of recovering sharp features of surfaces with creases and corners. Subdivision surfaces play a key role in computer graphics and numerical analysis. A class of piecewise-smooth surface representations was introduced based on subdivision to reconstruct smooth surface from scattered data [33]. Thin-shell finite element analysis was used to describe both the geometry and associated displacement fields [34]. The limit representation of Loop subdivision for triangular meshes was combined with the diffusion model to arrive at a discretized version of the diffusion problem [35]. Mixed finite element methods based on subdivision technology were used to construct high-order smooth surfaces with specified boundary conditions [36]. Truncated hierarchical Catmull-Clark subdivision was developed recently to support local refinement and generalize truncated hierarchical B-splines to arbitrary topology [37].

As a compatible technique of NURBS, subdivision surfaces are capable of the refinability of B-spline techniques. There recently have been a few works on the application of the subdivision-based IGA approach. Volumetric IGA based on Catmull-Clark solids was investigated in [38]. For complex physical domain, Powell-Sabin splines were used as IGA tools for advection-diffusion-reaction problems [39]. The bivariate splines in the rational Bernstein-Bézier form over the triangulation were applied in [40]. A reproducing kernel triangular B-spline-based FEM was proposed to solve PDEs [41]. Collocated isogeometric boundary element methods and unstructured analysis-suitable T-spline surfaces were coupled for linear elastostatic problems in [42]. A new generalized surface and IGA elements having the vertices of the irregular quad mesh through complementing bicubic splines and biquartic splines near irregularities in the mesh layout was presented in [43]. A framework for geometric design and IGA on unstructured quadrilateral meshes was proposed in [44]. A new type of Hermite bases for bicubic spline defined over a rectangular mesh with arbitrary topology was investigated in the framework of IGA [45]. Second-order PDEs on surfaces with the IGA was considered in [46].

Contributions

In the present paper, we introduce an alternative numerical method, namely, the isogeometric analysis based on the extended Loop subdivision approach (IGA-Loop), for solving the PDEs on subdivision surfaces. The finite elements consisting of quartic box-splines corresponding to each subdivided control mesh are utilized to represent the geometry of interest, and construct the solution space for dependent variables as well. The subdivision process is equivalent to the $h$-refinement of NURBS-based isogeometric analysis. The basis functions induced by the extended Loop subdivision possess the ability to represent the geometry exactly which is consistent with the concept of isogeometric analysis. The performance of the proposed method is evaluated by solving various surface PDEs, such as surface Laplace-Beltrami harmonic/biharmonic/triharmonic equations, which are defined on various limit surfaces of the extended Loop subdivision for different initial control meshes. Numerical experiments demonstrate that the proposed method has desirable performance in terms of the accuracy, convergence and computational cost for solving the above classical surface PDEs defined on both open and closed surfaces. The proposed approach is proved to be second-order accuracy in the sense of $L^{2}$-norm by theoretical and numerical results for the surface Laplace-Beltrami harmonic equation, which is outperformed over the standard linear finite element by several numerical comparisons. Furthermore, the $L^{2}$-norm error with the second-order convergence rate can be observed for the surface Laplace-Beltrami biharmonic/triharmonic equations.

Outline

Let us outline the content of the paper. The required mathematical preliminaries are introduced in Section 2. The isogeometric analysis based on the extended Loop subdivision is investigated in Section 3. We apply the proposed method to solve second/fourth/sixth-order surface PDEs in Section 4. The Section 5 presents several numerical experiments with comparison to the linear finite element method. Finally, the Section 6 concludes the paper.

2 Mathematical Preliminaries

Here we introduce some mathematical preliminaries which are requisite for treating the PDEs on surface.

2.1 Surface Parameterization

Let ${\cal S}=\{{\bm{x}}(u,v)\in\mathbb{R}^{3}:[u,v]^{{\rm T}}\in\Omega\subset\mathbb{R}^{2}\}$ be a piece of regular smooth parametric surface. The tangent plane and unit normal vector at ${\bm{x}}$ of ${\cal S}$ are defined as

[TABLE]

respectively, where ${\bm{x}}_{u}$ and ${\bm{x}}_{v}$ are the coordinate tangent vectors. The first fundamental form of $\cal S$ is

[TABLE]

where the coefficients $g_{\alpha\beta}=\langle{\bm{x}}_{u^{\alpha}},{\bm{x}}_{u^{\beta}}\rangle$, $\alpha,\beta=1,2$ with $(u^{1},u^{2})=(u,v)$. The second fundamental form of $\cal S$ is

[TABLE]

with the coefficients $b_{\alpha\beta}=\langle{\bm{x}}_{u^{\alpha}u^{\beta}},{\bm{n}}\rangle=-\langle{\bm{x}}_{u^{\alpha}},{\bm{n}}_{u^{\beta}}\rangle$. Set

[TABLE]

Area

The positive value

[TABLE]

is called the area of surface $\cal S$.

Mean Curvature

First we introduce Weingarten transformation ${\cal W}:T_{\bm{x}}{\cal S}\rightarrow T_{\bm{x}}{\cal S}$ which is a linear mapping, and satisfies

[TABLE]

The linear mapping can be represented by a $2\times 2$ matrix $S=[b_{\alpha\beta}][g^{\alpha\beta}].$ The eigenvalues $k_{1}$ and $k_{2}$ of $S$ are the principal curvatures of $\cal S$. Their arithmetic average is the mean curvature:

[TABLE]

2.2 Surface Differential Operators

Let $C^{1}({\cal S})$ denote the set of function $f:{\cal S}\rightarrow\mathbb{R}$ with continuous derivatives at any ${\bm{x}}\in{\cal S}$. Analogously, we can define the sets $C^{2}({\cal S})$, $C^{3}({\cal S}),$ $\ldots$. For the sufficiently differentiable $f({\bm{x}})=f({\bm{x}}(u,v))$ defined on surface $\cal S$, we denote

[TABLE]

For simplicity, we set $f_{1}:=f_{u}$, $f_{2}:=f_{v}$, $f_{11}:=f_{uu}$, $f_{12}:=f_{uv}$ and $f_{22}:=f_{vv}$.

Tangential Gradient Operator

Let $f\in C^{1}({\cal S})$, then the tangential gradient operator $\nabla_{\cal S}$ acting on $f$ at ${\bm{x}}\in{\cal S}$ is defined as

[TABLE]

By the definition of $g^{\alpha\beta}$ we have

[TABLE]

Obviously, $\nabla_{\cal S}f\in T_{\bm{x}}{\cal S}$ and hence $\langle\nabla_{\cal S}f({\bm{x}}),{\bm{n}}\rangle=0$.

Divergence Operator

Let ${\bm{v}}$ be a $C^{1}$ smooth vector field on surface $\cal S$. Then the divergence operator ${\rm div}_{\cal S}$ acting on ${\bm{v}}$ is defined as

[TABLE]

Laplace-Beltrami Operator

Let $f\in C^{2}({\cal S})$, then $\nabla_{\cal S}f$ is a smooth vector field on surface $\cal S$. The Laplace-Beltrami operator (LBO) $\Delta_{s}$ acting on $f$ is defined as

[TABLE]

With the definitions of $\Delta_{\cal S}$ and ${\rm div}_{\cal S}$, we derive

[TABLE]

2.3 Global Properties of Surface

Here we introduce the concept of the global surface, and derive some global properties, for example Green formula. There will be quite useful for the numerical analysis of surface PDEs.

Definition 2.1

The ${\cal S}=\bigcup\limits_{\alpha\in\Lambda}{\cal S}_{\alpha}$ is called a global orientable surface in $\mathbb{R}^{3}$ if it satisfies:

• (1)

Each ${\cal S}_{\alpha}=\{{\bm{x}}_{\alpha}(u,v)\in\mathbb{R}^{3}:[u,v]^{\rm T}\in{\Omega}_{\alpha}\subset\mathbb{R}^{2}\}$ is a regular surface in $\mathbb{R}^{3}$.

• (2)

For any $\alpha,\beta\in\Lambda$, if ${\cal S}_{\alpha}\cap{\cal S}_{\beta}\neq\varnothing$, ${\cal S}_{\alpha}\cap{\cal S}_{\beta}$ is still a piece of surface and the composite mapping

[TABLE]

is differentiable.

• (3)

For every ${\cal S}_{\alpha}$, there is an orientation ${\bm{n}}_{\alpha}$. Moreover ${\bm{n}}_{\alpha}={\bm{n}}_{\beta}$ for ${\bm{x}}\in{\cal S}_{\alpha}\cap{\cal S}_{\beta}$.

Theorem 2.1** (Green formula for LBO)**

Let ${\cal S}$ be an orientable surface, and $\Omega$ be a subregion of ${\cal S}$ with a piecewise smooth boundary $\partial{\Omega}$. Let ${\bm{n}}_{c}$ be the outward unit normal along the boundary $\partial{\Omega}$. Then for a given smooth vector field ${\bm{v}}\in C^{1}(\Omega)$, we have

[TABLE]

2.4 Sobolev Spaces on Surface

The FEM is applied to solve PDEs on surface, therefore the theory of Sobolev spaces on surface is requisite. Assume that ${\cal S}$ is a sufficiently smooth surface. For a given constant $k$ and a function $f\in C^{\infty}({\cal S})$, denote $\nabla^{k}f$ the $k$-th order covariant derivative of function $f$, with the convention $\nabla^{0}f=f$. Let

[TABLE]

We have the following definition of Sobolev space $H^{k}({\cal S})$.

Definition 2.2

Let ${\cal S}$ be a compact surface with at least $k$-th order smoothness. Sobolev space $H^{k}({\cal S})$ is the completion of $C_{k}({\cal S})$ in the sense of norm

[TABLE]

For the compact surface $\cal S$, we have

[TABLE]

3 Isogeometric Analysis Based on Extended Loop Subdivision

Isogeometric analysis is a recently developed computational approach that the solution space for dependent variables is represented in terms of the same functions which represent the geometry. The isogeometric procedures can be developed based on NURBS, A-patches or subdivision surfaces [2]. In this work, we investigate isogeometric analysis for surface PDEs based on subdivision surfaces.

A subdivision surface is a method of representing a smooth surface, which can be generated through a progressively iterative subdivision process starting from an initial control mesh by a prescribed subdivision scheme. The subdivision schemes used to generate smooth surfaces can be separated into two categories: interpolatory and approximatory. In the first category, the vertex positions of the initial mesh are fixed, and only the positions of new added vertices need to be computed in each subdivision step (see [47, 48]). However, for the second category, both the old and new vertex positions are required to update during each refinement (see [49, 50]). In general, the surfaces generated by approximating schemes have better quality than those generated by interpolating ones. In this article, the extended Loop subdivision is considered, which falls into the second category.

3.1 Extended Loop Subdivision Scheme

Let $\Omega^{0}_{h}$ be the input closed triangular mesh, which denotes as the initial control mesh of the conventional Loop subdivision. In each refinement step, each triangle is subdivided into four sub-triangles in which all the new vertices are generated by the weighted average of the old vertices, and all the old vertex positions on the refined mesh $\Omega^{k+1}_{h}$ are computed by the weighted average of the vertex positions on the mesh $\Omega^{k}_{h}$. Let ${\bm{p}}^{k}_{0}$ be a vertex with one-ring neighbors ${\bm{p}}_{j}^{k}$ ($j=1,\dots,n$) on $\Omega^{k}_{h}$, where $n$ is the valence of ${\bm{p}}^{k}_{0}$. The positions of the new generated vertices ${\bm{p}}^{k+1}_{i}$ on $\Omega^{k+1}_{h}$ corresponding to the edges of the previous mesh are computed by

[TABLE]

where the subscript index $i$ is regarded as modulo $n$. The old vertex gets new position according to the rule

[TABLE]

where

[TABLE]

proposed by Loop in [50] is named as Loop’s weight [52].

The conventional Loop subdivision scheme is suitable only for subdividing closed triangular control meshes without boundaries. Actually, for a lot of geometric modeling problems, the surfaces are usually constructed in a piecewise manner with fixed boundaries. In such a case, Loop subdivision scheme cannot be implemented near the boundaries of the control mesh. To overcome this deficiency, the extended Loop subdivision scheme was proposed by Biermann et al. in [51], which can treat triangular control meshes with boundaries. Note that the subdivision rule for boundaries is the same as that of cubic B-spline, and the control vertices on boundaries are treated as the control vertices of cubic B-spline curves with spaced knots. The extended Loop subdivision scheme is described as follows.

3.1.1 Vertex Refinement

The vertices can be separated into three types: corner vertex, boundary vertex and interior vertex. For different type of vertex, the refined strategy is given as:

• (1)

Corner vertex: The corner vertices are to be interpolated, meaning which are fixed.

• (2)

Boundary vertex: Let ${\bm{x}}_{i}$ be a boundary vertex, and let ${\bm{x}}_{l}$ and ${\bm{x}}_{r}$ be its two neighbor vertices on the boundary, then ${\bm{x}}_{i}$ is updated by $\frac{1}{8}{\bm{x}}_{l}+\frac{3}{4}{\bm{x}}_{i}+\frac{1}{8}{\bm{x}}_{r}$.

• (3)

Interior vertex: Use the conventional Loop scheme as (3.1) and (3.2).

3.1.2 Edge Refinement

Overall, we divide the edges into boundary, sub-boundary and interior edges. Boundary edges lie on the boundaries, which are the features of control mesh in general. Sub-boundary edges are not boundary edges but adjacent to the boundary vertices. The rest are the interior edges. For different type of edge, the edge refinement is performed as:

• (1)

Boundary edge: Let $[{\bm{x}}_{l}{\bm{x}}_{r}]$ be a boundary edge. The newly added vertex on this edge is the average of ${\bm{x}}_{l}$ and ${\bm{x}}_{r}$, i.e., $\frac{1}{2}({\bm{x}}_{l}+{\bm{x}}_{r}$).

• (2)

Sub-boundary edge: Let $[{\bm{x}}_{i}{\bm{x}}_{j}]$ be a sub-boundary edge with ${\bm{x}}_{i}$ being a boundary vertex, ${\bm{x}}_{l}$ and ${\bm{x}}_{r}$ being the two wing neighbor vertices of $[{\bm{x}}_{i}{\bm{x}}_{j}]$, then the newly added vertex on this edge is defined as $(\frac{3}{4}-\gamma){\bm{x}}_{i}+\gamma{\bm{x}}_{j}+\frac{1}{8}{\bm{x}}_{l}+\frac{1}{8}{\bm{x}}_{r}$ where $\gamma=\frac{1}{2}-\frac{1}{4}\cos\theta_{k}$ with $\theta_{k}=\frac{\pi}{k}$ for a boundary point, $\theta_{k}=\frac{\alpha}{k}$ for a convex corner point, and $\theta_{k}=\ \frac{2\pi-\alpha}{k}$ for a concave corner point.

• (3)

Interior edges: Use the conventional Loop scheme as (3.1) and (3.2).

3.2 Limit Form of the Extended Loop Subdivision

The explicit representation for the limit position of each vertex can be derived, which is stated as the following lemma.

Lemma 3.1** ([50])**

Let ${\bm{p}}^{0}_{0}$ be a control vertex of valence $n$ on the mesh $\Omega^{0}_{h}$, and ${\bm{p}}^{0}_{i}$, $i=1,\ldots,n$, be its one-ring neighbor control vertices. The vertices sequence $\{{\bm{p}}^{k}_{0}\}$ converges to a unique point

[TABLE]

as the subdivision step $k\rightarrow\infty$, where the $\alpha$ is defined by (3.3).

The limit surface of the extended Loop subdivision is $C^{2}$ everywhere except at the irregular vertices (i.e., having valence other than six) where is $C^{1}$ (see [53]). By Lemma (3.1), we can evaluate the position of the limit surface at any finite subdivision level and at any vertex by weighted averaging the vertex and its neighbors. However, it is not feasible to compute any point on the limit surface. Fortunately, there exists the explicit expression for the regular Loop subdivision surface patch. If all vertices of the triangle have valence six and none of its two-ring neighbor vertices is a boundary vertex, the resulting surface patch is called regular. The regular patch can be exactly described by a quartic box-spline, more precisely, which is formulated by the linear combination of 12 basis functions:

[TABLE]

where $(1-\xi-\eta,\xi,\eta)$ are the barycentric coordinates of the shaded unit reference triangle and ${\bm{p}}^{k}_{i}$ are the indexed vertices of the control mesh in Figure 1 ($a$). The analytic expressions of the basis functions $B_{i}$ are given as those in [52]. Each triangle of the control mesh can be regarded as a parametric domain, which corresponds to one triangular patch of the limit surface.

If a triangle is irregular, i.e., at least one of its vertices has a valence other than six or one of its two-ring neighbor vertices is a boundary vertex, the resulting surface patch cannot be represented by a quartic box-spline. For evaluating irregular patches, we use the fast scheme proposed by Stam in [56], in which the mesh needs to be subdivided repeatedly until the parameter values of interest are interior to a regular patch (See Figure 1 ($b$)). Hence, the basis functions given in (3.5) can be used to represent the geometry exactly.

3.3 Finite Element Space for IGA-Loop

Let ${\cal S}$ be a considered smooth surface, which is generated by the limit of the extended Loop subdivision scheme. And we denote $\Omega_{h}^{k}$ the $k$-th refined control mesh of the discretized representation ${\cal S}_{h}$, where ${\bm{p}}^{k}_{i}$ denotes its $i$-th control vertex and $h$ is the length of maximal edge. Each control vertex ${\bm{p}}^{k}_{i}$ on $\Omega_{h}^{k}$ can be evaluated recursively, which is described in Section 3.1. Each triangle patch of $\Omega_{h}^{k}$ can be locally represented by a linear combination of explicit box-spline tiles as (3.5). The boundaries of ${\cal S}_{h}$ are represented as the cubic B-spline curves which are preserved as the subdivision proceeds. The subdivision step is equivalent to the $h$-refinement of the NURBS-based isogeometric analysis.

Let us define the basis functions of the finite element space for IGA-Loop. For each control vertex ${\bm{p}}^{k}_{i}$ on the control mesh $\Omega^{k}_{h}$ of the surface ${\cal S}_{h}$, including boundary control vertices, we associate it with a basis function $\phi_{i}$, where $\phi_{i}$ is defined as the limit of the extended Loop subdivision scheme applying to the zero control values everywhere except at ${\bm{p}}^{k}_{i}$ where it is one. Note that the basis functions have a compact support within two-ring neighbors. The mesh $\Omega^{k}_{h}$, formulated as piecewise patches, is served as the parametric domain of the basis functions. Actually the set of quartic box-splines corresponding to each subdivided computable control patch are utilized to represent the geometry of interest and the solution space for dependent variables.

The control mesh $\Omega^{k}_{h}$, as a piecewise linear surface, is served as the definition domain of the basis function $\phi_{i}$. The mapping from $\Omega^{k}_{h}$ to $\phi_{i}$ is defined by a dual subdivision process. More precisely, when the extended Loop subdivision scheme is applied to the control function values recursively, the linear subdivision scheme (each triangle is partitioned into four equal-sized sub-triangles) is applied to the control mesh correspondingly. The limit of the former is $\phi_{i}$ and that of later is $\Omega^{k}_{h}$ itself.

The basis functions share some properties with the well-known B-spline basis. These properties are important in the proposed method. Now let us describe them as follows.

• (1)

Positivity. The weights of the extended subdivision rules are positive. Hence the basis function $\phi_{i}$ is nonnegative everywhere and positive around ${\bm{p}}^{k}_{i}$.

• (2)

Locality. It is known that the limit value at a control vertex is a linear combination of the one-ring neighbor values. Hence, the limit value is zero at a control vertex if the control values on the one-ring neighbor control vertices are zero. Therefore, the support of the basis function is within the two-ring neighborhood.

• (3)

Partition of Unity. Since all the subdivision rules have the properties that the weights are summed to one. Therefore, if we choose all the control values as one. The control values after one subdivision step are still one. This implies that $\sum\phi_{i}=1$. This property is called partition of unity.

• (4)

Interpolatory Properties at the Boundary. The extended subdivision rules on the boundary do not involve the interior control vertices. Hence the basis functions for the interior control vertices are zero at the boundary.

• (5)

Tangential Property. Let ${\bm{p}}^{k}_{i}$ be a control vertex, with non of its one-ring neighbor control vertices are boundary control vertices. Then $\nabla_{s}\phi_{i}$ vanishes on the boundary. This fact can be observed by considering the eigen-decomposition of the control vertices. Let ${\bm{x}}^{k}\in\mathbb{R}^{(n+1)\times 3}$ be a vector consisting of one-ring neighbor control vertices of ${\bm{p}}_{i}^{k}$ at the subdivision level $k$, $S\in\mathbb{R}^{(n+1)\times(n+1)}$ be the local subdivision matrix that converts ${\bm{x}}^{k}$ to ${\bm{x}}^{k+1}$, i.e.,

[TABLE]

Here $n$ stands for the valence of ${\bm{p}}_{i}^{k}$. Suppose ${\bm{x}}^{1}$ is decomposed into

[TABLE]

where ${\bm{e}}_{0}$, ${\bm{e}}_{1}$, $\ldots$, ${\bm{e}}_{n}$ are the eigenvectors of $S$. Here we assume that these eigenvectors are arranged in the order of non-increasing eigenvalues $\lambda_{j}$. Then

[TABLE]

where $\lambda_{0}=1$, $\lambda_{1}=\lambda_{2}<1$. It is well-known that the limit position at the center is ${\bm{a}}_{0}$. The tangent direction at this vertex are ${\bm{a}}_{1}$ and ${\bm{a}}_{2}$, and ${\bm{a}}_{j}$ is given by ${\bm{a}}_{j}^{\rm T}=\tilde{\bm{e}}_{j}^{\rm T}{\bm{x}}^{1}$. $\tilde{\bm{e}}_{j}$ are the left eigenvectors of $S$ with normalized condition $\tilde{\bm{e}}_{j}^{\rm T}{\bm{e}}_{j}=1$. The analysis above is valid for control function values. The fact that $\nabla_{\cal S}\phi_{i}$ vanishes on the boundary implies that the tangent vector of the subdivision surface on the boundary determined by the boundary and sub-boundary control vertices. This is similar to Bézier and B-spline surfaces.

• (6)

Linear independency. The functions $\phi_{i}$, $i=0,\ldots,m$, are defined by the extended Loop subdivision scheme. This fact can be derived from the unique solvability of the following interpolation problem:

Lemma 3.2** ([54])**

Let ${\bm{p}}^{\infty}_{i}$ be the limit position of the control vertex ${\bm{p}}^{k}_{i}$ with its 1-ring neighbor vertices ${\bm{p}}^{k}_{j},j=1,2,\ldots,n_{i}$, $u({\bm{p}}^{\infty}_{i})$ be the $i$-th interpolation function value, and ${\rm u}({\bm{p}}^{k}_{i})$ be the $i$-th control function value, where $i=1,\ldots,m$. The system

[TABLE]

is always solvable uniquely. Here the $l_{i}$ is defined as (3.4).

• (7)

Interpolant Error. With the basis of Lemma 3.2, we derive the following interpolant error estimation:

Lemma 3.3** ([55])**

Assume that $\Omega_{h}^{k}$ is the $k$-th refined control mesh for the smooth surface ${\cal S}$ generated by the extended Loop subdivision, where $h$ stands for the length of the maximal edge. For $u\in H^{2}({\cal S})$, there exists an interpolation function $u_{I}\in H^{2}({\cal S}_{h})$ such that

[TABLE]

4 Applications to Surface PDEs

We consider the approximations for the Laplace-Beltrami harmonic, biharmonic and triharmonic problems on different surfaces by the proposed IGA-Loop strategy. The numerical examples we provide in Section 5 include both open and closed surfaces, more concretely, which are the Laplace-Beltrami harmonic problem on a quarter of cylinder and one-eighth of sphere, the Laplace-Beltrami biharmonic problem on a cylinder, and the Laplace-Beltrami triharmonic problem on a unit sphere. The numerical error will be analyzed which shows the consistent convergence rate with (4.11). For comparison, we also consider the approximations of these three problems by standard linear element discretization. Hence we adopt the mixed formulations of the high-order Laplace-Beltrami biharmonic and triharmonic problems.

4.1 Surface PDEs

Laplace-Beltrami Harmonic Problem

The Laplace-Beltrami harmonic equation with zero boundary condition can be written as

[TABLE]

where the boundary $\partial{\cal S}$ of surface $\cal S$ is Lipschitz continuous. Here $f:{\cal S}\rightarrow\mathbb{R}$ is a given sufficiently regular function, and $\Delta_{\cal S}$ is the Laplace-Beltrami operator. The weak form of (4.1) is given as follows:

[TABLE]

Laplace-Beltrami Biharmonic Problem

The Laplace-Beltrami biharmonic equation with homogeneous boundary conditions reads as

[TABLE]

where ${\bm{n}}$ is the outward directed unit vector normal to the boundary $\partial{\cal S}$. Let $\Delta_{\cal S}u=-v$, then the mixed weak formulation of (4.3) reads as

[TABLE]

where $v:{\cal S}\rightarrow\mathbb{R}$ are the auxiliary unknowns.

Laplace-Beltrami Triharmonic Problem

The Laplace-Beltrami triharmonic equation with homogeneous boundary conditions is formulated as

[TABLE]

Let $\Delta_{\cal S}u=-v$ and $\Delta_{\cal S}v=-w$, then the mixed formulation of (4.5) reads as

[TABLE]

where $v,w:{\cal S}\rightarrow\mathbb{R}$ are the auxiliary unknowns.

4.2 IGA-Loop Discrertization of Surface PDEs

The finite element space defined by the limit form of the extended Loop subdivision is employed to represent the geometry and also represent the solution space for the discretization of the surface PDEs. The input triangular mesh serves as the initial control mesh of the extended Loop subdivision. Such function is $C^{2}$-continuous everywhere except at the extraordinary vertices with $C^{1}$-continuity. Using globally $C^{1}$-continuous basis functions induced by the extended Loop subdivision yields our IGA finite dimensional space that is the subspace of the Sobolev space $H^{2}({\cal S})$. For each control vertex ${\bm{x}}_{i}$ of the limit surface ${\cal S}$, we associate it with a basis function $\phi_{i}$, where $\phi_{i}$ is defined by the limit of the extended Loop subdivision for the zero control values everywhere except at ${\bm{x}}_{i}$ , where it is one. Hence the support of $\phi_{i}$ is a piecewise function, and covers the two-ring neighborhood of the vertex ${\bm{x}}_{i}$.

The function can be locally parameterized on the unit triangle defined by $T=\{(\xi,\eta)\in\mathbb{R}^{2}:\xi\geq 0,\eta\geq 0,\xi+\eta\leq 1\}$ where $(1-\xi-\eta,\xi,\eta)$ are the barycentric coordinates of the unit triangle. Using this parameterization, the discretized representation of the surface domain is ${\cal S}_{h}=\bigcup_{i=1}^{k}{\mathcal{T}}_{i},\ \mathring{\mathcal{T}}_{i}\cap\mathring{\mathcal{T}}_{j}=\varnothing\ {\rm for}\ i\neq j$, where $\mathring{\mathcal{T}}_{i}$ is the interior of triangular patch $\mathcal{T}_{i}$. Each triangular patch can be parameterized locally as

[TABLE]

Denote ${\bm{e}}_{j},j=1,\ldots,m_{i}$ be the two-ring neighborhood patches around the control vertex ${\bm{x}}_{i}$. Then if ${\bm{e}}_{j}$ is regular, the explicit box-spline expression as in (3.5) exists for $\phi_{i}$ on ${\bm{e}}_{j}$. If ${\bm{e}}_{j}$ is irregular, local subdivision, as described in Section 3.2, is needed around ${\bm{e}}_{j}$ until the parameter values of interest are interior to a regular patch. The parameterization has no overlap. Each control vertex ${\bm{x}}_{i}\in{\cal S}$ has its unique parameter coordinates. Using the set of the basis functions $\{\phi_{i}\}$, the limit surface ${\cal S}$ of the extended Loop subdivision is expressed as

[TABLE]

Each control triangular surface patch of the surface ${\cal S}$ is defined locally by only a few related basis functions, since the supports of the basis are compact. The surface boundaries are represented as the cubic B-spline curves which are preserved as the subdivision proceeds. With the parameterization, the differential operators on the surface as described in Section 2 can be computed directly, and the computation of the function integration on the surface is replaced by

[TABLE]

4.2.1 Numerical Discretization

In what follows, we discretize the weak formulations of the Laplace-Beltrami based problems. Let ${\bm{x}}_{j},j=1,\ldots,n_{0}$ be the interior control vertices of the surface $\cal S$, and ${\bm{x}}_{j},j=n_{0}+1,\ldots,n$ be its boundary control vertices. Recalling (4.2), (4.4) and (4.6), by means of the basis functions $\phi_{j}\in H^{2}({\cal S}_{h})$ introduced by the limit form of the extended Loop subdivision, we have the following discrete description for the function $u^{h}$ to be determined

[TABLE]

where $u^{h}_{j}\ (j=n_{0}+1,\ldots,n)$ is the known boundary conditions. The other two auxiliary functions $v^{h}$ and $w^{h}$ are represented as

[TABLE]

For simplicity, denote

[TABLE]

the solution vector, and the other two auxiliary solution vectors

[TABLE]

The coefficient matrices and the right-hand side terms are represented respectively as

[TABLE]

whose elements are

[TABLE]

Take the test functions $\varphi=\phi_{i},i=1,\ldots,n_{0}$, for the Laplace-Beltrami harmonic problem (4.2) on a quarter of cylinder and one-eighth of sphere in Test Suite 1. Considering the zero boundary condition, we obtain the linear system as

[TABLE]

Choose the two sets of the test functions $\varphi=\phi_{i},i=1,\ldots,n_{0}$, and $\psi=\phi_{i},i=1,\ldots,n$, for the Laplace-Beltrami biharmonic problem (4.4) on a cylinder in Test Suite 2. Considering the homogeneous boundary condition, we achieve the following linear system

[TABLE]

Taking the three sets of the test functions $\varphi=\psi=\zeta=\phi_{i},i=1,\ldots,n$, for the Laplace-Beltrami triharmonic problem (4.6) on a unit sphere in Test Suite 3, the linear system reads as

[TABLE]

4.2.2 Precompute the Basis Functions

The framework of our IGA-Loop is the same standard process as the classical FEM. The related basis functions and their derivatives need to be precomputed for each control patch of ${\cal S}_{h}$ before solving the linear systems. However, those computations are not intuitive in comparison with standard linear elements since the required two-ring neighbors around each patch have arbitrary topological structure. And additional geometric data are reflected in the subdivision schemes around boundaries. We classify the control mesh into interior patches, sub-boundary patches, and boundary patches. The patches containing boundary vertices are named as boundary patches. The patches adjacent to boundary patches are called sub-boundary patches. And all of the other patches are called interior pathes.

We use the following algorithms to treat those patches. For interior patches, the Stam’s algorithm is applied to this case (see [56]). The sub-boundary patches can be subdivided into four interior sub-patches by once, then the algorithm for interior case can be used. In addition, the boundary patches can be subdivided repeatedly till their sub-patches belong to the sub-boundary case, then the above methods can be applied to evaluate them.

Note that the Stam’s fast evaluation scheme is always suitable for interior patches with only one extraordinary vertex. Therefore, it is necessary to first subdivide once each patch of the initial mesh. The evaluation of basis functions over their support elements uses general Gaussian integration, which just needs a few subdivision steps to bring Gaussian quadrature knots into a box-spline patch. The integrations for computing the matrix elements are computed by a 12-point Gaussian quadrature rule. That is, each triangle is subdivided into four sub-triangles and a 3-point Gaussian quadrature rule is employed on each of the sub-triangles. The 3-point Gaussian quadrature rule has error bound $O(h^{3})$, where $h$ is the maximal edge length. We adopt the adaptive numerical method developed in our former work [55].

4.2.3 Error Estimate

The interpolant error (3.6) obtained in our former work can be applied to the case of second-order surface PDEs with zero boundary condition (4.1). By means of the general analysis method of the classical finite elements, we can obtain the $H^{1}$-norm estimation of the exact solution $u$ and the approximate solution $u^{h}$ resulting from IGA-Loop, and the $H^{0}$-norm estimation by means of the Aubin-Nitsche duality argument. The error estimation is stated as the following theorem.

Theorem 4.1

Consider the Laplace-Beltrami harmonic problem defined on the limit form ${\cal S}$ of the extended Loop subdivision, endowed with the zero boundary condition. Let $u\in H^{2}({\cal S})$ be the exact solution of the problem, and $u^{h}$ be the approximate solution obtained with IGA-Loop strategy, the following error estimates hold:

[TABLE]

where $c$ is a constant independent of $h$.

Furthermore, we also analyze the numerical errors for the cases of the fourth-order and the sixth-order surface PDEs, even defined on closed surface like the sphere, which show the consistent convergence rates.

5 Numerical Experiments

In this section, we present several numerical experiments to show the performance of our IGA-Loop method by solving the three surface PDEs on different surfaces. It should be indicated that the numerical computations are conducted on the limit surface generated by the extended Loop subdivision, therefore the integration evaluation of the Gauss-Legendre knots is performed on the triangulation of the limit surface.

We compare the proposed method (IGA-Loop) with the linear finite element method (FEM-Linear) in terms of accuracy, convergence and computational cost. That is because FEM-Linear is the most often used finite elements in industry, which has the same number of control vertices and error estimate as those of IGA-Loop. The following test suites seek to assess the performance of the proposed IGA-Loop against different surface PDEs.

5.1 Test Suite 1: Surface Laplace-Beltrami Harmonic Equation

Here we consider the tests for the Laplace-Beltrami harmonic equation. The first example is

[TABLE]

with the exact solution

[TABLE]

on a quarter of the cylinder ${\cal S}_{1}=\{(x,y,z):x^{2}+y^{2}=1\ \&\ x\geq 0\ \&\ y\geq 0\ \&\ 0\leq z\leq 2\}$ with calculated $f$.

We solve the above problem based on three different control meshes from coarse to refined, as shown in the top row of Figure 1. The total numbers of vertices/patches are 221/384, 825/1536, and 3185/6144 respectively. They have the same limit surface ${\cal S}_{1}$, namely, one quarter of a cylinder. The computational accuracy is indicated by the error distribution $u-u^{h}$. As shown in Figure 1, the accuracy of IGA-Loop is higher than that of FEM-Linear. In addition, we also compare the convergence rate of IGA-Loop against successive refinement. The error is progressively decreasing as the control mesh becoming finer and finer as shown in Figure 2. Apart from the visual comparison in Figure 1, the quantitative comparison of the error between IGA-Loop and FEM-Linear is performed in Figure 2. As we can see that IGA-Loop has higher accuracy.

Consider the second example as the same Laplace-Beltrami harmonic equation as the first one, but with different exact solution

[TABLE]

on a one-eighth of the sphere ${\cal S}_{2}=\{(x,y,z):x^{2}+y^{2}+z^{2}=1\ \&\ x\geq 0\ \&\ y\geq 0\ \&\ z\geq 0\}$ for suitable $f$.

As before, we compute the above surface PDE based on three different control meshes from coarse to refined, as shown in the top row of Figure 3, which has the same limit surface ${\cal S}_{2}$. The total numbers of vertices/patches are 107/176, 389/704, and 1481/2816 respectively. For numerical comparison, FEM-Linear is used to solve the same problem. As shown in Figure 3, the accuracy of IGA-Loop is higher than that of FEM-Linear by the profile of error distribution. We further compare the convergence rate of the proposed method against successive refinement. The error is gradually decreasing as the control mesh becoming finer and finer as shown in Figure 4. The quantitative comparison of the error between the proposed method and FEM-Linear is given in Figure 4. As we can see that IGA-Loop has higher accuracy.

In Figure 2 and 4, the errors $e=\|u-u^{h}\|_{L^{2}({\cal S})}$ are computed with the $L^{2}$-norm of the surface operator, where both methods have the convergence rate 2. From the comparison of the approximation errors involved in FEM-Linear and IGA-Loop discretizations, for which the errors versus the number of subdivision times are plotted, we observe the former is approximately 1.5 times more than the latter. The FEM-Linear discretization requires the number of degree of freedoms more than the IGA-Loop method, which means IGA-Loop is potentially more efficient. The explanation of the phenomenon is that the approximation error is introduced by the discretization of the surface geometries using the linear finite elements, instead, exact representation of the surfaces can be achieved by means of IGA-Loop.

We performed the above two numerical tests to demonstrate the mathematical results in Theorem 4.1. For the Laplace-Beltrami biharmonic problem defined on a cylinder and the Laplace-Beltrami triharmonic problem defined on a unit sphere, the error estimate (4.11) in norm $L^{2}({\cal S})$ of Theorem 4.1 holds, which will be studied in the following contents.

5.2 Test Suite 2: Surface Laplace-Beltrami Biharmonic Equation

In what follows, we solve a fourth-order example of the Laplace-Beltrami biharmonic equation as

[TABLE]

where the exact solution

[TABLE]

with $\theta={\rm atan}(y/x)$ on the surface ${\cal S}_{3}=\{(\theta,z):0\leq\theta\leq 2\pi\ \&\ 0\leq z\leq 1\}$ for suitable $f$.

As shown in the top row of Figure 5, three different control meshes from coarse to refined are given. The total numbers of vertices/patches are 432/768, 1632/3072, and 6336/12288 respectively. They have the same limit surface ${\cal S}_{3}$ as the extended Loop subdivision proceeds. We solve the Laplace-Beltrami biharmonic equation based on the proposed IGA-Loop. The computational accuracy is shown by the error distribution $u-u^{h}$. As shown in Figure 5, the accuracy of IGA-Loop is higher than that of FEM-Linear.

We further compare the convergence rate of IGA-Loop against successive refinement. The error is gradually decreasing as the control mesh becoming finer and finer as shown in Figure 6. The quantitative comparison of the error in norm $L^{2}({\cal S})$ versus the number of subdivision times, obtained with IGA-Loop and FEM-Linear is given in Figure 6. As shown that IGA-Loop has higher accuracy. In Figure 6, the errors of norm $L^{2}({\cal S})$ under three times refinement level is plotted. We observe that the convergence rate is 2 for the error in the sense of norm $L^{2}({\cal S})$, which is consistent with the error estimate (4.11). As shown in Figure 6, the approximation errors obtained by the FEM-Linear discretization is approximately 1.6 times more than that of our IGA-Loop discretization. It means that the IGA-Loop approximation only requires a smaller number of degree of freedoms than FEM-Linear to achieve the same accuracy.

5.3 Test Suite 3: Surface Laplace-Beltrami Triharmonic Equation

Furthermore, a sixth-order Laplace-Beltrami triharmonic equation is solved as following

[TABLE]

where for calculated $f$, the exact solution is

[TABLE]

on the closed sphere ${\cal S}_{4}$ with the radius $r=1$.

As shown in the top row of Figure 7, three different control meshes from coarse to refined are provided. The total numbers of vertices/patches are 706/1408, 2818/5632, and 11266/22528, respectively. They have the same limit surface ${\cal S}_{4}$ by the extended Loop subdivision. We solve the Laplace-Beltrami triharmonic equation based on the proposed IGA-Loop. The computational accuracy is shown by the error distribution $u-u^{h}$. As shown in Figure 7, the accuracy of IGA-Loop is higher than that of FEM-Linear.

Different from the above two tests, the test model for the sixth-order problem is on a closed surface. We further observe the convergence rate of IGA-Loop against successively finer meshes of the unit sphere in Figure 8. The approximate error (4.11) in norm $L^{2}({\cal S})$ holds and the same convergence rate 2 is obtained. Furthermore, we can observe the approximation error obtained with the FEM-Linear discretization is approximately 1.8 times more than that of the IGA-Loop discretization. Therefore, IGA-Loop obtains the same level of accuracy more efficiently than FEM-Linear.

5.4 Computational Cost

Actually, the solver of our IGA-Loop can be merged into the framework of FEM-Linear, however the support of each basis function of IGA-Loop is two-ring neighbors described in Section 3. In this section, for the four test examples of Figures 1, 3, 5 and 7, we list the corresponding time costs of computing the basis functions in Tables 1 and 2. To show that the extended Loop subdivision scheme does not require structured meshes and can support the same meshes with any topological structure as the standard finite elements, the valences of the control vertices in these surfaces are in the range of 3 to 12 depicted in Figure 1, 3, 5 and 7.

Table 1 corresponds to the examples of the Laplace-Beltrami harmonic problem where the number of control vertices and triangular patches is shown in the first and the third columns, the second and the fourth columns list the time cost (in seconds) of computing the basis functions and their derivatives because they can be pre-computed and saved in a data structure. Once refinement makes the number of triangular patches on the refined meshes increases four times and their sizes approximately decrease by half. The computational cost of IGA-Loop with the same control meshes is larger than FEM-Linear. One reason is the computation of the derivatives of the linear basis functions is unnecessary, as it happens instead with IGA-Loop. The other reason is the computation around the boundary patches is more complex for our IGA-Loop than FEM-Linear because the extended Loop subdivision scheme embodies the angle information different of every patch.

The data for the examples of the Laplace-Beltrami biharmonic problem and the Laplace-Beltrami triharmonic problem are listed in Table 2. Similar results to the Laplace-Beltrami harmonic problem can be observed. Here we notice the data of the Laplace-Beltrami triharmonic problem which is listed in the third and the fourth columns. The example is a closed sphere, the computation cost of IGA-Loop is almost equal to FEM-Linear. We will explain the phenomenon. Firstly the computation of the basis functions around boundary patches for IGA-Loop is removed, and the rest is only the case of computing on the interior patches. Secondly, interior patches share the same set of basis functions which depend only on the valence list of their control vertices, so we can merge the patches according to their valence list, and then use Stam’s fast evaluation to treat them. Finally, as the mesh refinement proceeds, the added patches are ordinary whose valences of the three control vertices are six.

6 Conclusion

The isogeometric analysis for surface PDEs based on the extended Loop subdivision is presented. The set of quartic box-splines corresponding to each subdivided control mesh are utilized to represent the geometry of interest, and to construct the solution space for dependent variables as well. The finite elements induced by the extended Loop subdivision possess the ability to represent the geometry exactly which fully agrees with the concept of isogeometric analysis.

We have evaluated the performance of the proposed IGA-Loop in dealing with the various surface PDEs, including the Laplace-Beltrami harmonic/biharmonic/triharmonic equations on different limit surfaces. As shown in the visual and quantitative results, the proposed method can yield desirable performance in terms of the accuracy, convergence and computational cost for solving the above classical surface PDEs. As demonstrated with both the theoretical and numerical results, the proposed IGA-Loop approach is proved to be second-order accuracy in the sense of $L^{2}$ norm of the surface. Through various comparisons, the performance of IGA-Loop is also outperformed over the standard linear finite elements. The solution of these problems will help to establish a computational framework for the geometric PDEs solved by IGA-Loop in our future work.

Acknowledgments

Qing Pan is supported by National Natural Science Foundation of China (NSFC) (No.11671130), Hunan Provincial Natural Science Foundation of China (No.2018JJ2248, 2019JJ50395), and the Open Project Program of the National Laboratory of Pattern Recognition of China (NLPR). Chong Chen is supported by the Key Research Project of Beijing National Natural Science Foundation (No.Z180002). Gang Xu was supported by NSFC (No.61761136010, 61772163) and Zhejiang Provincial Natural Science Foundation (No.LR16F020003).