Assessment of an isogeometric approach with Catmull–Clark subdivision surfaces using the Laplace–Beltrami problems

An isogeometric approach for solving the Laplace–Beltrami equation on a two-dimensional manifold embedded in three-dimensional space using a Galerkin method based on Catmull–Clark subdivision surfaces is presented and assessed. The scalar-valued Laplace–Beltrami equation requires only C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^0$$\end{document} continuity and is adopted to elucidate key features and properties of the isogeometric method using Catmull–Clark subdivision surfaces. Catmull–Clark subdivision bases are used to discretise both the geometry and the physical field. A fitting method generates control meshes to approximate any given geometry with Catmull–Clark subdivision surfaces. The performance of the Catmull–Clark subdivision method is compared to the conventional finite element method. Subdivision surfaces without extraordinary vertices show the optimal convergence rate. However, extraordinary vertices introduce error, which decreases the convergence rate. A comparative study shows the effect of the number and valences of the extraordinary vertices on accuracy and convergence. An adaptive quadrature scheme is shown to reduce the error.


Introduction
Hughes et al. [33] proposed the concept of isogeometric analysis (IGA) in 2005. The early works on IGA [10,18,47] focussed on geometries modelled using Non-Uniform Rational B-Splines (NURBS) as these are widely used in computer aided design (CAD). NURBS can be used to model freeform, two-dimensional curves. However, a NURBS surface is a tensor product surface generated by two NURBS curves, thereby imposing limitations for modelling complex geometries with arbitrary topologies. Complex CAD models are always composed of a number of NURBS patches. These patches are often poorly connected in the design stage. When such models are used for analysis, the unmatched patches must be treated carefully to ensure the geometries are watertight. Furthermore, because NURBS can not be B Zhaowei Liu zhaowei.liu@glasgow.ac.uk 1 Glasgow Computational Engineering Centre, University of Glasgow, Glasgow G12 8LT, UK 2 Chair of Applied Mechanics, Friedrich-Alexander Universität Erlangen-Nürnberg, Paul-Gordan-Str. 3, 91052 Erlangen, Germany locally refined, adaptive mesh refinement method cannot be employed. A number of alternative CAD techniques were developed and adopted in IGA to overcome these limitations, including Hierarchical B-splines [29,52], T-splines [11,49], PHT-splines [23,42], THB-splines [13,30] and LR B-splines [24,34]. Some of these recent techniques are being adopted by the engineering design market. However, the majority are the subject of academic research and not widely used in the CAD community. Moreover, computing the basis functions for analysis using these alternative approaches can be expensive. Catmull and Clark [14] developed a bicubic Bspline patch subdivision algorithm for describing smooth three dimensional objects. The use of Catmull-Clark subdivision surfaces to model complex geometries in the animation and gaming industries dates back to 1978. Catmull-Clark subdivision surfaces can be considered as uniform bi-cubic splines which can be efficiently evaluated using polynomials.
In CAD, distortion of regular parametrizations are inevitable and indeed vital when modelling complex geometries. Allowing 'extraordinary vertices' ensures that Catmull-Clark subdivision surfaces can be used for modelling complex geometries with arbitrary topology. Cirak et al. [17] implemented Loop subdivision surfaces for solving the Kirchhoff-Love shell formulation. This was the first appli-cation of subdivision surfaces to engineering problems. Subdivision surfaces have subsequently been used in electromagnetics [19], shape optimisation [6,7] , acoustics [15,38] and lattice-skin structures [56].
Catmull-Clark subdivision surfaces face a number of challenges when used for analysis. Many of these have been discussed in the literature, however a unified assessment is lacking. This manuscript provides a clear and concise discussion of the challenges and limitations of Catmull-Clark subdivision surfaces.
Engineering designs often require exact geometries including circles, spheres, tori and cones. However, subdivision surfaces can not capture these geometries exactly. Moreover, there are always offsets between the control meshes and the surfaces. Fitting subdivision surfaces [37] aim to overcome this limitation. Although the fitting subdivision surfaces still can not model arbitrary geometries exactly as they are interpolated using cubic splines, they can approximate the given geometries closely through least-square fitting. Another challenge of subdivision surfaces is that they can model smooth closed manifolds easily but require special treatment to model manifolds with boundaries. A common solution is to introduce 'ghost' control vertices to provide bases for interpolating. From the perspective of analysis, the shape functions will span into 'ghost' elements [17]. In addition, the spline basis functions do not possess an interpolating property. Thus it is difficult to directly impose Dirichlet boundary conditions. Meshless methods and extended finite element methods have developed strategies to overcome this problem [28,39]. A common strategy is to modify the weak form of the governing equation. Methods include the Lagrangian multiplier method [5], the penalty method [3] and Nitsche's method [32,43].
Conventional Catmull-Clark subdivision surfaces can not be locally refined. Truncated hierarchical Catmull-Clark subdivision surfaces (THCCS), developed by Wei et al. [54], overcome this limitation. They generalise truncated hierarchical B-splines (THB-splines) to meshes with arbitrary topology. Wei et al. [55] subsequently improved their method using a new basis function insertion scheme and thereby enhanced the efficiency of local refinement. The extraordinary vertices introduce singularities in the parametrisation [41,51]. Catmull-Clark subdivision surfaces have C 2 continuity everywhere except at the surface points related to extraordinary vertices where, as demonstrated by Peters and Reif [45], they possess C 1 continuity. Stam [50] developed a method to evaluate Catmull-Clark subdivision surfaces directly without explicitly subdividing, thus allowing one to evaluate elements containing extraordinary vertices. Although the surface gradients can not be evaluated at the extraordinary vertices, they can be evaluated at nearby quadrature points. Thus, subdivision surfaces can be used as C 1 elements as required, for example, in thin shell theory [17]. Nevertheless, the evaluation of points around extraordinary vertices of Catmull-Clark surfaces introduces error. The conventional evaluation method repeatedly subdivides the element patch until the target point fall into a regular patch allowing a uniform bi-cubic B-spline patch to be mapped the subdivided element patch. The extraordinary vertex also introduces approximation errors because of the singular parameterisations at extraordinary vertices [40,44]. Stam's natural parametrisation only can achieve C 0 continuity at extraordinary vertices. Recently Wawrzinek and Polthier [53] introduced a characteristic subdivision finite element scheme that adopted a characteristic reparameterisation for elements with extraordinary vertices. The evaluated limiting surface is at least C 1 everywhere and the numerical accuracy is improved. Zhang et al. [57] optimised the subdivision scheme to improve its approximation properties when used for thin-shell theory.
Using the finite element method to solve the partial differential equations (PDEs) on surfaces dates back to the seminal work by Dziuk [25], which developed a variational formulation to approximate the solution of the Laplace-Beltrami problems on two dimensional surfaces. This method was extended to solve nonlinear and higher-order equations on surfaces by Dziuk and Elliott [26]. Dziuk and Elliott [27] also provided a thorough review on finite element methods for approximating the solution of PDEs on surfaces. Dedner et al. [22] proposed a discontinuous Galerkin (DG) method for solving a elliptic problem with the Laplace-Beltrami operator on surfaces. Adaptive DG [21] and high-order DG [1] methods were also developed for solving PDEs on surfaces. However, the accuracy of these methods depends on the approximation of the mean curvatures of the surfaces. The geometrical error is dominant when conventional Lagrangian discretisation is used to approximate solutions on complex surfaces. Isogeometric discretisation maintains the exact geometry and overcomes this limitation. Dedè and Quarteroni [20] proposed an isogeometric approach for approximating several surface PDEs involving the Laplace-Beltrami operator on NURBS surfaces. Bartezzaghi et al. [9] solved PDEs with high order Laplace-Beltrami operators on surfaces using NURBS based isogeometric Galerkin method. More accurate results are obtained using an IGA approach over the conventional finite element method. Langer et al. [36] present an isogeometric DG method with non-matching NURBS patches allowing the approximation of PDEs on more complex surfaces.
This work presents a thorough and unified discussion of several major issues related to isogeometric Galerkin formulation based on Catmull-Clark subdivision surfaces. The difficulties associated with imposing Dirichlet boundary conditions, the reduction of the approximation power around extraordinary vertices, and the problem of sufficient numerical integration in the element with extraordinary vertices will be examined and discussed. Previous studies [16,17] on Catmull-Clark subdivision surfaces for analysis introduce ghost degrees of freedoms for constructing basis functions in elements at boundaries. We propose a method which modifies the basis functions at boundaries to ensure they are only associated with given control vertices. No additional ghost degrees of freedom are involved. A penalty method is employed to impose Dirichlet boundary conditions. This does not change the size or symmetry of the system matrix and is straightforward to implement. An adaptive quadrature scheme inspired by [35] is presented to increase the integration accuracy for elements with extraordinary vertices. The proposed method can perform isogeometric analysis on complex geometries using Catmull-Clark subdivision discretisations. A test for approximating Poisson's problem on a square plate is conducted to demonstrate the properties of the method in a simplified setting so as to distill the key features. The approach is also used for solving the Laplace-Beltrami equation which is a benchmark problem for curved manifolds [35,41]. A comparative convergence study is conducted between the Catmull-Clark subdivision method and the conventional finite element method. The effects of the extraordinary vertices and modified bases at boundaries on convergence are examined. Catmull-Clark subdivision surfaces are limiting surfaces generated by successively subdividing given control meshes. They are identical to uniform bi-cubic B-splines. Thus, they have difficulty to represent desired geometries exactly. Here, a least-squares fitting method is used to fit any given geometry using Catmull-Clark subdivision surfaces.
This manuscript first summarises the subdivision algorithm and the evaluation method for Catmull-Clark subdivision surfaces. Then, techniques for using Catmull-Clark for numerical analysis and improving accuracy are presented in Sect. 3. Section 4 presents the Laplace-Beltrami problem and Sect. 5 shows a Galerkin method with Catmull-Clark subdivision surface bases. Section 6 showcases the numerical results.

Catmull-Clark subdivision surfaces
There exist a variety of subdivision schemes, but the basic idea is to use a subdivision scheme to generate a smooth surface through a limiting procedure of repeated refinement steps starting from an initial polygonal grid. The Catmull-Clark algorithm can generate curves and surfaces which are identical to cubic B-splines. The algorithms for curves and surfaces are shown in Appendices A.1 and A.2, respectively. This section will briefly introduce the methods for interpolating and evaluating curves and surfaces using the Catmull-Clark subdivision algorithm.  Figure 1 shows a curve generated using a subdivision algorithm. The interpolated curve is identical to a cubic B-spline curve. The limiting curve can be interpolated using cubic basis splines and associated control points. With a control polygon containing n control points, the curve is naturally divided into n − 1 elements. Each element in the curve is associated with one segment of the control polygon. To interpolate on the target element, four control points including the neighbouring control points are required. For example, if one aims to evaluate the geometry of element 2 in Fig. 1, the four control points P 1 ,P 2 ,P 3 and P 4 are required and the curve point is evaluated as

Curve interpolation and evaluation based on the subdivision algorithm
where ξ ∈ [0, 1] is the parametric coordinate within an element. The basis functions for element 2 are defined by The bases are visualised in Fig. 2a. They are C 2 continuous across element boundaries. Element 1 in Fig. 1 contains the end of the curve, which has an end curve point that coincides with the control point. In order to evaluate this curve, one needs to mirror the point P 2 to P 0 as The curve point can now be evaluated using basis splines with a set of control points shown in Fig. 2b. However, if one adopts a spline discretisation for analysis, this strategy of end element treatment will introduce additional 'ghost-like' degrees of freedom. To avoid this problem, the expression for P 0 (3) is substituted into the interpolating equation yielding Hence only three control points are required to evaluate a curve point and the modified basis functions for interpolating end elements are defined by Figure 2b illustrates the modified basis functions. It achieves the same basis functions as the cubic B-Spline with p + 1 multiple knots at the two end points. The new basis functions do not possess the Kronecker delta property but do have the interpolating property at the boundary. The performance of modified bases in analysis will be discussed in Sect. 6.1. The global basis functions for interpolating the curve in Fig. 1 are shown in Fig. 2c. It is worth noting that this subdivision curve is a cubic B-spline curve and represents a special case of Lane-Riesenfeld subdivision it can not model conical shapes exactly. This property is significantly different to NURBS and motivates Sect. 3.1 on geometry fitting.

Interpolating and evaluating Catmull-Clark subdivision surfaces
One defines the number of elements connected with the vertex as the valence. A regular vertex in a Catmull-Clark surface mesh has a valence of 4. A vertex with a valence not equal to 4 is called an extraordinary vertex. This allows subdivision surfaces to handle arbitrary topologies. In their seminal paper [14], Catmull and Clark proposed a way to modify the weight distributions for extraordinary vertices in order to describe complex geometries. With this simple solution, Catmull-Clark surfaces can use a single mesh to present surfaces of arbitrary geometries while other splinebased CAD tools, such as NURBS surfaces, need to link multiple patches. The limiting surface of the Catmull-Clark subdivision algorithm has C 2 continuity over the surface except at the extraordinary vertices where they have C 1 continuity as proven by Peters and Reif [45]. This section will illustrate the methods of interpolating and evaluating both regular element and element with an extraordinary vertex in Catmull-Clark subdivision surfaces. Figure 3a shows a subdivision surface element (dashed) which does not contain an extraordinary vertex. In order to evaluate a point in this Catmull-Clark element, an element patch must be formed. The patch consists of the element itself and the elements which share vertices with it. A regular element patch has 9 elements with 16 control vertices. The surface point can be evaluated using the 16 basis functions associated with these control points as

Element in a regular patch
where ξ := (ξ, η) is the parametric coordinate of a Catmull-Clark subdivision surface element. A Catmull-Clark surface is obtained as the tensor product of two Catmull-Clark curves. The basis functions are defined by where N (ξ ) or N (η) are the basis functions defined in Eq. (2) and presented in Fig. 3a.
• is the modulus operator and % denotes the remainder operator which gives the remainder of the integer division. Figure 3b shows the element patch of a subdivision surface element (shaded) that has an edge on the physical boundary. This type of element has only 5 neighbour elements so that it belongs to an element patch which has 12 control vertices. To evaluate this element, a common solution is to generate a set of 'ghost' vertices outside the domain to form a full element patch [17]. However, this method involves additional degrees of freedom in numerical analysis. Instead, the curve basis functions in Eq. (5) are adapted to deal with the element on the boundary. The same strategy is used for elements which have two edges on the physical boundary as shown in Fig. 3c.

Element in a patch with an extraordinary vertex
Extraordinary vertices are a key advantage of Catmull-Clark subdivision surfaces which allows them to model complex geometries with arbitrary topologies. However, it increases the difficulty of evaluating the surfaces. Figure 4a shows a Catmull-Clark subdivision element which contains one extraordinary vertex.
In order to evaluate this element, one needs to re-numerate the control points as shown in Fig. 4a. After applying one level of subdivision, new control points are generated and this element is subdivided into four sub-elements, as shown in Fig. 4b. The sub-elements Ω 1 , Ω 2 and Ω 3 are now in a regular patch. However, the last sub-element (grey) still has an extraordinary vertex. If the target point to be evaluated is in this region, we must repeatedly subdivide the element until the point falls into a sub-element with a regular patch. Then, the point can be evaluated within the sub-element with the new set of control points P n,k , where n is the number of subdivision required and k = 1, 2, 3 is the sub-element index shown in Fig. 4b. The new control point set is computed as where D k is a selection operator to pick control points for the sub-elements. A andĀ are two types of subdivision operators. P 0 is the initial set of control points. The detailed approach is given in [50] and also can be found in Appendix A.3. P n,k contains 16 control points. Then, a surface point in the element with an extraordinary vertex can be computed as whereξ is the parametric coordinates of the evaluated point in the sub-element, which can be mapped from ξ as Equation (9) can thus be rewritten as whereN is the Catmull-Clark subdivision surfaces basis function. DefineN as a set of 2κ + 8 basis functions in an element with an extraordinary vertex and N is a set of 16 regular basis functions defined in Eq. (7).N can be calculated in a vector form aŝ The derivatives of the Catmull-Clark subdivision surfaces basis functions for elements containing extraordinary vertices are expressed as and can be computed by where ∂ξ ∂ξ can be considered as a mapping matrix defined by

Remark 1
The calculation of the basis functionsN at a physical point x involves two mappings. The first is from the physical domain to the parametric domain of an element with an irregular patch, x → ξ . Because the irregular patch does not have the tensor-product nature, n levels of subdivisions are required and the point is mapped to the parametric domain of a sub-element, ξ →ξ . This second mapping is defined in Eq. (10). The value of n approaches positive infinity when ξ approaches the extraordinary vertex which has the parametric coordinate (0, 0). Hence the diagonal terms in the mapping matrix (15) tend to positive infinity as n → ∞. This results in the basis functionsN not being differentiable at ξ = 0. This problem is termed singular configuration in [35], and singular parameterisation in [41,51].

Techniques for analysis and improving accuracy
This section presents three techniques which are essential for using Catmull-Clark subdivision surfaces in numerical analysis. A geometry fitting method using Catmull-Clark surfaces is introduced in Sect. 3.1. Section 3.2 illustrates an adaptive quadrature scheme for integrating element with an extraordinary vertex to improve accuracy. Section 3.3 introduces the penalty method for applying essential boundary conditions.

Geometry fitting
Catmull-Clark subdivision surfaces are CAD tools which can construct limiting surfaces from control polygons and meshes. However, in a number of engineering problems, the geometry is given as an industry design and a limit surface that is a "best approximation" of this desired geometry required. Litke et al. [37] introduced a method for fitting a Catmull-Clark subdivision surface to a given shape. They employed, both a least-squares fitting method and a quasi-interpolation method to determine a set of control points for a given surface. The least-square fitting method is used here. One first chooses a set of sample points S = {s 1 , s 2 , . . . , s n s } ∈ Γ , where Γ is the geometry, n s is the number of sample points. Each sample point should be evaluated using Catmull-Clark subdivision bases with control points as where n b = 2κ + 8 is the number of local basis functions.
Then the set of sample points can be evaluated as where P = {P 1 , P 2 , . . . , P n c } is a set of control points with n c control points. L is an evaluation operator of Catmull-Clark curves or surfaces. Set ξ = (0, 0) to ensure the sample points correspond to the control vertices and n s ≡ n c , then L is a square matrix. The control points can be calculated as If more sampling points n s are chosen than the required number of control points n c , then L is invertible, a least-squares method is used to obtain a set of control pointsP that minimises S − LP 2 aŝ shows that 6 sample points are chosen from the given curve and one assembles the evaluation operator for these sampling points. The control points can be obtained by solving (18). Using these control points, the limit curve can be interpolated. Since 6 sample points is not sufficient to capture the given curve, the limit curve is significantly different to the given curve. Figures 5b and c show the curve fitting with 11 and 21 sample points, respectively. Increasing the number of samples points, the limit curve converges to the given curve.

Adaptive quadrature rule for element with an extraordinary vertex
In numerical analysis, a Gauss quadrature rule is applied to integrate over Catmull-Clark subdivision elements. A one dimensional quadrature rule with n q Gauss points can exactly evaluate the integrals for polynomials of degree up to 2n q −1.
The polynomial degree of a cubic B-spline function is 3.
Because the basis functions of a Catmull-Clark subdivision element in regular element patch are generated as the tensor product of two cubic splines, 2×2 Gauss points can be used in this case. However, if a Catmull-Clark subdivision element has an extraordinary vertex, the basis functions are generated by Eq. (12). In this case, basis functions are not polynomials and the derivatives of the basis functions suffer from the singular parametrisation problem, see Remark 1. Thus, the standard Gauss quadrature can not be used to evaluate the element integral. Inspired by [35], an adaptive quadrature rule, well suited to Catmull-Clark subdivision surfaces is adopted by integration at a number of levels of subdivided elements.
With n d levels of subdivisions, the element is subdivided into 3n d + 1 sub-elements as shown in Fig. 4d. The sub-elements can be evaluated using cubic B-splines with new control vertices except for the ones having an extraordinary vertex. Thus the Gauss quadrature rule can be used to evaluate the integrals in 3n d sub-elements. With a number of subdivisions, the integration error can be reduced. In this work, n d = 7 is chosen in order to obtain sufficiently accurate values of the integrals.

Penalty method for applying boundary condition
The basis functions do not have the Kronecker delta and interpolating properties, so boundary conditions can not be directly applied using conventional methods. The method used here is a penalty method which uses a penalty parameter and boundary mass matrix to apply the boundary conditions approximately. It preserves the symmetry of the system matrix and does not increase its size. However, the penalty parameter should be carefully selected. If fine meshes with more degrees of freedom are adopted, a larger penalty parameter must be chosen. The Dirichlet boundary condition is defined as An L 2 projection is used for applying the Dirichlet boundary condition, where for test function v, one obtains Using the cubic B-spline functions in Eq.
(2) to discretise u and v and the same strategy for formulating the system matrix, one introduces a boundary mass matrix as where n b e is the number of boundary elements, and The right hand side vector for applying the boundary conditions is thus Then the discrete system of equations arising from (21) is We note that the elements for applying boundary conditions are the discretisation of the surface boundary which are one dimensional cubic B-spline curves and only one-dimensional Gauss quadrature rule is used for integration. However, one uses the global degrees of freedom indices to assemble M b and f b , so that they have the same size as the system matrix and global right-hand side vector, respectively. Assume the system of equations is expressed as Ku = f, where K is the system matrix, u is the global coefficients vector to be solved for, and f is global right-hand side vector. Then, we scale M b and f b using a penalty factor β and combine them with the systems of equations as The Dirichlet boundary condition (20) is here weakly applied to the system of equations. A relatively large penalty factor β = 10 8 is selected for all numerical examples. It is sufficiently large to ensure good satisfaction of the constraint but not too large so as to significantly impact the conditioning of the system.

Laplace-Beltrami problem
The governing partial differential equation which we want to solve to illustrate fundamental features of subdivision surfaces is given by where Γ is a two dimensional manifold (with outward unit normal vector n) in three dimensional space R 3 and Δ Γ (•) is the Laplace-Beltrami operator (also called surface Laplacian operator). The Dirichlet boundary condition is expressed in (20). We will use a manufactured solution to compute against the approximate solution. The Laplace-Beltrami operator is defined by where ∇ Γ (•) is the surface gradient operator defined by Hence the surface gradient of a scalar function v can be calculated as the spatial gradient subtracted by its normal part as where ∇(•) is the spatial gradient operator. Hence the surface Laplacian of v is given by where ∇ 2 v is the Hessian matrix of v, and ∇n is the gradient of the normal vector, which is arranged in a matrix as We define the total curvature at a surface point x ∈ Γ as the surface divergence of the normal, that is c(x) := ∇ Γ · n. For a given manufactured solution u m , the right hand side of Eq. (27) can thus be computed as

Galerkin formulation
The weak formulation of problem (27) is where v is an admissible test function. The weak formulation is partitioned into n e number of elements, as Discretising v, ∇u and ∇v using the Catmull-Clark basis functions N given in Eq. (7) produces where J is the surface Jacobian for the manifold, given in a matrix form as For details on the computation of J −1 see [46] and for a discussion of superficial tensors such as J in the context of Laplace-Beltrami equation, see [31]. If the element contains an extraordinary vertex, the shape functions N A are replaced byN A in Eq. (12). The surface gradient of the shape functions is computed as and J = ∂x Integrating the discrete problem using Gauss quadrature, the system of Eq. 34 becomes where A is the assembly operator and n q is the number of quadrature points in each element, w i is the weight for i th quadrature point, n e is the number of elements and n b is the number of basis functions of the element. The basis functions N e are replaced byN e if the element e contains an extraordinary vertex. In this case, the basis functions are not differentiable and their derivatives approach positive infinity when points are close to the extraordinary vertex (see Remark 1). Thus |J| approaches positive infinity at extraordinary vertices. Errors result if quadrature is adopted to integrate the contributions from element containing extraordinary vertices.
The discrete system of equations to solve is thus given by

Numerical results
A 'patch test' [58] on a two-dimensional plate is first presented to assess the consistency and stability of the proposed formulation in a simplified setting. Then, the Laplace-Beltrami equation is solved on both cylindrical and hemispherical surfaces. Convergence studies are conducted. The influence of extraordinary vertices is also investigated.

'Patch test'
The 'patch test' is performed on a two dimensional flat plate where the Laplace-Beltrami operator reduces to the Laplace operator. The problem proposed in Sect. 4 reduces to the Poisson problem expressed given by This partial differential equation is solved on the square plate shown in Fig. 6a with the essential boundary conditions The essential boundary conditions are imposed using the penalty method. Natural homogeneous boundary conditions are applied on the remaining two edges of the plate. Four different manufactured functions for f are used. The functions, analytical solutions for u and their gradients ∂u ∂ x 2 are given in Table 1. We investigate both a regular and an irregular mesh. The regular mesh is a 4 × 4 element patch without extraordinary vertices as shown in Fig. 6b. In all of the tests, a geometry error is absent.
For Test 1, the right hand side f = 0 so that ∂u ∂ x 2 = 2. Solving the equation using the proposed Catmull-Clark subdivision method, the numerical result u h is exactly 2 everywhere as shown in Fig. 7b. Thus passes the consis- π sin(π x 2 ) 1 π sin(π x 2 ) + 2x 2 cos(π x 2 ) + 2 Test 1 has no right-hand side term, thus the analytical solution u is linear and its gradient is a constant. The analytical solutions for Tests 2 and 3 are quadratic and cubic, respectively, and their gradients are linear and quadratic, respectively. Test 4 has a sine function as the right-hand side term which gives a cosine function as the gradient of the analytical solution tency test and the eigenvalue of the system matrix are all positive and non-zero after application of the essential boundary conditions. The gradient ∂u ∂ x 2 for Test 2 and 3 are linear and quadratic respectively. Recall that when interpolating functions in elements with edges on physical boundaries, the basis functions are modified, see Eqs. (3) and (4). In other words, the gradients of the function u are expected to be constant at boundaries. Figure 7a, c and e show the numerical results for these tests. The results are smooth and capture the analytical solutions well. Figure 7d and f compare the numerical results of ∂u ∂ x 2 to the analytical solution for Tests 2 and 3. The Catmull-Clark subdivision method is also compared to linear and quadratic Lagrangian finite element methods. There is a substantial error in both boundary regions in Test 2 for Catmull-Clark subdivision method. This is because the method imposes the gradient to be constant at both boundaries. The numerical result of the Catmull-Clark subdivision method in Test 3 has a substantial error in the region close to the top boundary (x 2 = 2) but captures the gradient in the region close to the bottom boundary (x 2 = 0) well because the analytic solution for the gradient in the bottom boundary region is near-constant. These errors at the boundaries will pollute the numerical result in the interior of the domain, which will reduce the convergence rate. The gradients approximated by the linear and quadratic Lagrangian finite elements are piecewise constant and piecewise linear, respectively. The results of the Catmull-Clark subdivision methods for these two tests lies between the linear and quadratic Lagrangian elements. The gradient ∂u ∂ x 2 in Test 4 is a cosine function which is non-polynomial and it behaves as a constant in both boundary regions shown in Fig. 7h. The Lagrangian elements only possess C 0 continuity across elements and their gradients hence have jumps between elements. The Catmull-Clark subdivision elements capture the gradients of the given function better as they are C 1 smooth. Figure 8 shows the plots of normalised global L 2 and H 1 errors against the element size. The normalised global L 2 error is defined by where • L 2 is the L 2 norm defined as • L 2 = Γ | • | 2 dΓ . The normalised global H 1 error is computed as where • H 1 is the H 1 norm defined as We set the element size of the coarsest mesh as 1. Then, the normalised element size for the refined meshes are 1 2 , 1 4 , . . .. The convergence rate of Tests 2 and 3 are sub-optimal at 2.5 (L 2 error) and 1.5 (H 1 error). The optimal convergence rate for cubic elements should be p + 1 = 4 (L 2 error) and p = 3 (H 1 error), where p is the polynomial degree of the basis functions. The numerical result captures the analytical solution well and the convergence rate for Test 4 is optimal. The same convergence study is now repeated starting from a mesh containing extraordinary vertices as shown in Fig. 6c. Figure 8a and b show the plots of normalised element sizes against the L 2 and H 1 errors, respectively. The same convergence rates are obtained for Tests 2 and 3. However, the convergence rate of Test 4 is also reduced to 2.5 (L 2 error) and 1.5 (H 1 error). Figure 8c and d show the plots of normalised element sizes against L 2 and H 1 errors, respectively, for the mesh with an extraordinary vertex.
The Catmull-Clark subdivision method can pass the patch test when the function gradient is a constant but has difficulties to capture the gradients in boundary regions when they do not behave like a constant. When the gradient behaves like a constant in the boundary regions, the optimal convergence rate can be obtained. If this is not the case, a reduction of the convergence rate is observed. The presence of the extraordinary vertex in the patch also reduces the convergence rate. It is also important to note that the Catmull-Clark subdivision elements have advantages in describing non-polynomial functions since their basis functions are cubic and C 2 continuous.

Comparison with NURBS and Lagrangian elements
We now compare the convergence rate associated with Catmull-Clark elements against conventional Lagrangian elements and NURBS. Bézier extraction [12] is adopted to decompose a NURBS surfaces into C 0 Bézier elements to provide an element structure for the isogeometric Galerkin method. This is a widely-used method for isogeometric analysis using T-splines [48]. As the Lagrangian and Bézier elements can fully pass the 'patch test', they both have no approximation error for Tests 1, 2 and 3. Figure 9 compares their behaviour in approximating non-polynomial solution in Test 4. Mesh 1 is used for all methods. All methods exhibit an optimal convergence rate. Since no geometry error is involved in the 'patch test', the Bézier element provides the same performance as the Lagrangian element without the advantages of exact geometry representation. The Catmull-Clark element is slightly more accurate than other two methods for this specific test.

Cylindrical surface example
The first numerical example considered is a cylindrical surface. The analysis domain of the problem is the cylindrical Fig. 10 The geometry is given in (a) and (b) is the control mesh which constructs the best approximating Catmull-Clark subdivision surface of the given geometry. The control mesh is generated using least-squares fitting. c Shows the numerical result u h on the cylindrical surface   Fig. 10a. Surfaces fitting methods are used to construct the control mesh, see Sect. 3.1. The first level control mesh is shown in Fig. 10b. This has no extraordinary vertices. The Laplace-Beltrami problem on this manifold domain is solved using the Galerkin formulation presented in Sect. 5. Essential boundary conditions are applied on ∂Γ . The right-hand side function f is computed using the definition in Eq. (33). Figure 10c shows the numerical result u h which matches the manufactured analytical solution (47) very well. A convergence study is now conducted for this geometry. The refined control meshes are constructed using the least-squares fitting method described in Sect. 3.1. Figure 11 compares the convergence rates between Catmull-Clark subdivision surfaces with two different order Lagrangian elements. In this example, the shortcoming caused by extraordinary vertices and boundary gradients are not present, and the Catmull-Clark subdivision surfaces have the same convergence rate p + 1 as cubic Lagrangian elements.

Hemispherical surface example
The second geometry investigated is a hemispherical surface with radius equal to 1 as shown in Fig. 12a. We use the same strategy to fit the Catmull-Clark subdivision surfaces to the hemispherical surface. The control mesh shown in Fig. 12b is generated to discretise the surface into a number of Catmull-Clark elements. The control mesh has four extraordinary vertices. Figure 12e shows the solution u h .

Convergence study with an isogeometric approach
In engineering, designers usually do not know the geometry of the product in advance. The geometry information is purely from the CAD model. Catmull-Clark subdivision surfaces, as a design tool, provide the geometry which is the Fig. 12 a is a hemispherical surface. b is the control mesh for constructing subdivision surfaces to fit the hemispherical surface. c is 1-Level refined mesh for the hemispherical surface. d is 2-level refined mesh for the hemispherical surface. e shows the numerical result u h on this surface design of the engineering product. In this case, engineers do not need to approximate the given geometry with Catmull-Clark elements. They can directly adopt the discretisation from the CAD model for analysis. For example, we adopt the control mesh shown in Fig. 12b as the initial control mesh. It can be used to generate a limit surface approximating a hemisphere, as shown in Fig. 12a, with Catmull-Clark subdivision bases. It is important to note the limit surface is not an exact hemisphere since it is evaluated using cubic basis spline functions. However, this surface is the domain of our problem and it will stay exact the same during the entire analysis (isogeometric) and h-refinement with subdivision algorithm will not change the geometry.
The same problem is solved on the subdivision surfaces. A convergence study is done with another two levels of subdivision control mesh as shown in Fig. 12c and d. Note, refinement does not change the number of extraordinary vertices. The two new meshes still have four extraordinary vertices. The two control meshes can be used to evaluate the same limit surface shown in Fig. 12a. The Catmull-Clark  subdivision surfaces are compared with quadratic and cubic Lagrangian elements. Generally, Catmull-Clark subdivision elements can achieve higher accuracy per degree of freedom than Lagrangian elements. From the initial to the second level of mesh refinement, the Catmull-Clark subdivision elements have a similar convergence rate to cubic Lagrangian elements. After that, the convergence rate is equivalent to quadratic Lagrangian elements as shown in Fig. 13. Figure 14a shows the sparsity pattern of the system matrix K for the Catmull-Clark subdivision discretisation. The size of the matrix is the same as the system matrix assembled using a linear Lagrange discretisation. However, because the Catmull-Clark subdivision discretisation uses cubic basis functions with non-local support and there are 16 shape functions in a subdivision element with no extraordinary vertex, the number of non-zero entries in columns and rows is more than the linear Lagrange discretisation (i.e. the sparsity is decreased and the bandwidth increased). Thus, the system matrix of a Catmull-Clark subdivision discretisation has the same size but is denser than the linear Lagrange discretisation shown in 14b. Figure 14c is the sparsity patterns of cubic Lagrange discretisations. p-refinement increase the number of degrees of freedom as well as the number of non-zero entries in rows and columns. Thus there is no significant change in the density of the system matrices. The Catmull-Clark subdivision discretisation has the same number of non-zero entries in each row or column as the cubic Lagrangian discretisation but has a much smaller size.

Quadrature error
The presence of extraordinary vertices leads to difficulties in integration as described in Sect. 3.2. Figure 15a, c and e show the point-wise errors at surface points for three levels of mesh refinement using the standard Gauss quadrature rule. The number of extraordinary vertices remains 4 after refinement.
For the analysis using the initial mesh, the error in the regions around extraordinary vertices have similar magnitudes to the other regions. However, after a level of refinement, the error in the other regions is reduced more than the area around the four extraordinary vertices. After the second refinement, the error is concentrated in the areas around the four extraordinary vertices. Figure 15b, d and f plot the point-wise errors on the same mesh analysed with the adaptive quadrature rule shown in Sect. 3.2. The errors around extraordinary vertices are now decreased.

Approximation error
The presence of extraordinary vertices introduces approximation errors. Then we investigate the effect of the number and valence of extraordinary vertices on numerical accuracy. Figure 16a, b and c are three control meshes with different numbers of extraordinary vertices. Figure 16a shows a control mesh without an extraordinary vertex. Figure 16b shows a control mesh with four extraordinary vertices, including two vertices with a valence of 3 and two vertices with a valence of 5. The control mesh in Fig. 16c has seven extraordinary ver-tices, including four vertices with a valence of 4, two vertices with a valence of 5 and one vertex with a valence of 6. It is important to note the three different control meshes construct different but similar geometries. The Laplace-Beltrami problem is solved using the Galerkin formulation with the same right-hand side function f computed in (33). Both standard and adaptive Gauss quadrature rules are used for all cases. Figure 17a, c and e show the solution of u on the surfaces constructed using the three meshes. Because of the similarity of the geometries and solutions, the three cases are used to inves-  Fig. 17b, d and f. Meshes with extraordinary vertices have larger maximum point-wise errors close to the extraordinary vertices, while the mesh without extraordinary vertices has increased uniform point-wise error. Figure 18 shows the convergence rates for the three cases. Meshes without extraordinary vertices can achieve the optimal p + 1 convergence rate and p = 3. In general, the more extraordinary vertices a mesh contains, the more error results. The extraordinary vertices increase the global errors in the results and reduce the convergence rate. Since the global errors also include quadrature errors, the adaptive quadrature rule serves to reduce the quadrature errors. With the adaptive quadrature rule, the convergence rates are improved for the 4 and 7 extraordinary vertices cases but the results still agree with our assumption that increasing the number and valence of extraordinary vertices will produce higher error. Table 2 compares the computational cost for assembling the system matrix for the standard and adaptive quadrature rules. Because the number of extraordinary vertices remains constant after subdivision, the difference in computational time between the standard and adaptive quadrature schemes diminishes.

Complex geometry
This final example considers the ability of the Catmull-Clark method to provide high-order discretisations of complex geometry. The model considered is that of a racing car from CAD and imported into Autodesk Maya [4] for removal of extraneous geometry and the generation of the surface mesh shown in Fig. 19a. Modelling such geometry using NURBS surfaces would require a number of patches to be spliced together. A model based on Catmull-Clark subdivision surface can directly evaluate the smooth limit surface in Fig. 19b where n e = 9152 for this example. Figure 19c indicates the domain where essential (Dirichlet) and natural boundary conditions are applied. The essential boundary Γ d is composed of two parts as where Domains defined for applying boundary conditions.   x 2 (f ) Pointwise error on Γ n .
The natural boundary conditions is applied to the rest of the domain Γ n = Γ \Γ d . The numerical result matches the analytical solution well as shown in Fig. 19d. Figure 19e shows the results on Γ n and a maximum point-wise error 2.8% is observed in Fig. 19f.

Conclusions
A thorough study of the isogeometric Galerkin method with Catmull-Clark subdivision surfaces has been presented. The same bases have been used for both geometry and the Galerkin discretisation. The method has been used to solve the Laplace-Beltrami equation on curved two-dimensional manifold embedded in three dimensional space using the Catmull-Clark subdivision surfaces. An approach to fit given geometries using Catmull-Clark subdivision scheme has been outlined. A method to model open boundary geometries without involving 'ghost' control vertices, but involving errors in function gradients close to boundary regions, has also been described. The penalty method has been adopted to impose the Dirichlet boundary conditions. The optimal convergence rate of p + 1 has been obtained when using a cylindrical control mesh without extraordinary vertices.
A reduction of convergence rates has been observed when the function gradients at the boundaries do not behave like constant, or control meshes contain extraordinary vertices. The adaptive quadrature scheme significantly improves the accuracy. The effect of the number and valence of the extraordinary vertices in convergence rates has been investigated and an adaptive quadrature rule implemented. This successfully improved the convergence rates for the proposed method. The convergence rate of the proposed method is not worse than 2.5 (L 2 error) and 1.5 (H 1 error).
In future work, this method will be investigated with problems requiring C 1 continuity such as the deformations of thin shells.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-

A.1 Lane-Riesenfeld subdivision algorithm for curves
The Lane-Riesenfeld algorithm successively refines a curve starting from an initial control polygon. After a number of subdivisions, the curve is limited to a B-spline. Figure 20 illustrates a special case of this subdivision algorithm. The control point P i 2 j in the i th level of refinement is computed from the upper level control points as Point P i 2 j is the mid-point of P i−1 j -P i−1 j+1 , and is called an 'edge point'. The control point P i 2 j+1 is computed as To compute this point, one needs to connect the mid-points of P i 2 j -P i−1 j+1 and P i−1 j+1 -P i 2 j+2 . The point P i 2 j+1 is the midpoint of the connecting line. This type of point is called 'vertex point'. Each 'vertex point' is associated with an upper level control point. Figure 21 shows two levels of refinements using the Lane-Riesenfeld algorithm and the limiting result which is a cubic B-spline curve.

A.2 Catmull-Clark subdivision algorithm for surfaces
The application of the subdivision algorithm to surfaces follows in a similar manner to curves. One face in the original Equipped with these formulae, the new control points on the ith level of refinement P i can be computed as: S is a subdivision operator -a matrix consisting of a set of weights. Each weight is associated with a control point in P i−1 . The weight distributions for different types of control points are shown in Fig. 23. The weight distributions for extraordinary point are shown in Fig. 24. After successive levels of refinements, a smooth B-spline surfaces is obtained.

A.3 Computing control point set for sub-elements
We denote the control points of an irregular patch in Fig. 4a as a set P. The initial control points of the patch are expressed as P 0 = P 0 0 , P 0 1 , . . . , P 0 2κ+6 , P 0 2κ+7 .
The subdivision step is represented as where A is the subdivision operator given by The terms S, S 11 , S 12 , S 21 and S 22 are defined in [50] and S is given in Eq. 61. To evaluate the sub-element Ω 1 , Ω 2 and Ω 3 in Fig. 4b, one needs to pick 2κ + 8 control points out of the new 2κ + 17 control point patch. A selection operator D + κ for sub-element Ω k and k = 1, 2, 3 is used to select the necessary control points from P 1 , that is As shown in Fig. 4c, after successive subdivisions, the nonevaluable element can be limited to a negligible region.
(68) The sub-element index k is determined as