Aspects of Isogeometric Analysis with Catmull-Clark Subdivision Surfaces

09/23/2019 ∙ by Zhaowei Liu, et al. ∙ University of Glasgow 0

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 studied. 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.



There are no comments yet.


page 13

page 17

page 23

page 25

page 26

page 27

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

hughes2005isogeometric proposed the concept of isogeometric analysis (IGA) in 2005. The early works on IGA [cottrell2009isogeometric, reali2006isogeometric, bazilevs2008isogeometric]

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 any arbitrary two-dimensional curve. 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 connected poorly 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 NUBRS can not be 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 

[forsey1988hierarchical, vuong2011hierarchical], T-splines [sederberg2003t, bazilevscalo2010], PHT-splines [deng2008polynomial, nguyen2011rotation] and LR B-splines [dokken2013polynomial, johannessen2014isogeometric]. However, theses recent techniques are currently not widely used in the CAD community. Catmull and Clark [catmull1978recursively] developed a bicubic B-spline 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 are derived as tensor-products of two Catmull-Clark subdivision curves. 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. However, extraordinary vertices introduce singularities in the parametrisation [takahashi2012application, nguyen2014comparative].

Peters:2008aa proved that Catmull-Clark subdivision surfaces have continuity everywhere except for the surface points related to extraordinary vertices, which possess continuity. stam1998exact developed a method to evaluate Catmull-Clark subdivision surfaces directly without explicitly subdividing, thus providing a method which can 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 elements as required, for example, in thin shell theory. cirakortiz2000 implemented Loop subdivision surfaces for solving the Kirchhoff-Love shell formulation. This was the first application of subdivision surfaces to engineering problems. Subdivision surfaces have subsequently been used in electromagnetics [Dault2015], shape optimisation [bandara2016shape] and acoustics [liu2018isogeometric].

Catmull-Clark subdivision surfaces face a number of challenges when used for analysis. 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 [litke2001fitting]

aim to overcome this limitation. Although the fitting subdivision surfaces still can not model exact geometries since they are interpolated using cubic splines, they can approximate the given geometries closely through least-square fitting. Another challenge of the subdivision surfaces is that they can model smooth closed manifolds easily but need special treatment to model manifolds with boundaries. The common method is to introduce ‘ghost’ control vertices to provide bases for interpolating. From the perspective of analysis, the shape functions will span into ‘ghost’ elements 

[cirakortiz2000]. 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 [fernandez2004imposing, moes2006imposing]. A common strategy is to modify the weak form of the governing equation. Strategies include the Lagrangian multiplier method [babuvska1973finite], the penalty method [arnold1982interior] and Nitsche’s method [nitsche1971variationsprinzip, hansbo2002unfitted].

The primary aim of this paper is to present a thorough analysis of Catmull-Clark subdivision surfaces and an accompanying isogeometric Galerkin formulation. Catmull-Clark subdivision surfaces are limiting surfaces generated by successively subdividing given control meshes. Thus, they have difficulties to represent desired geometries exactly. Here, a least-squares fitting method is used to fit any given geometry using Catmull-Clark subdivision surfaces. Previous studies [cirakortiz2000, Cirak:2001aa]

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 truncates 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 

[juttler2016numerical] is presented to increase the integration accuracy for elements with extraordinary vertices. The proposed method integrates all the aforementioned features and 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 [juttler2016numerical, nguyen2014comparative]. 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 truncated bases on convergence are examined.

This manuscript first summarises the subdivision algorithm and the evaluation method for Catmull-Clark curves and surfaces. Then, techniques for using Catmull-Clark for numerical analysis and improving accuracy are presented in Section 3. Section 4 presents the Laplace-Beltrami problem and Section 5 shows a Galerkin method with Catmull-Clark subdivision surface bases. Section 6 showcases the numerical results.

2 Catmull-Clark subdivision curves and surfaces

Subdivision surfaces is a mature geometry discretisation technology with history dating back to the late 1970s [catmull1978recursively]. 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 Appendix A.1 and A.2, respectively. This section will briefly introduce the methods for interpolating and evaluating Catmull-Clark subdivision curves and surfaces.

2.1 Interpolating and evaluating Catmull-Clark subdivision curves

Figure 1 shows a curve generated using Catmull-Clark subdivision algorithm. The limiting curve can be interpolated using cubic basis splines and associated control points. With a control polygon containing control points, the curve is naturally divided into elements.

Figure 1: A Catmull-Clark curve is interpolated using basis splines and its control polygon.

Each element in the Catmull-Clark subdivision 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 in Figure 1, the four control points ,, and are required and the curve point is evaluated as


where is the parametric coordinate within an element. The basis functions for element are defined by


The bases are visualised in Figure 1(a). They are continuous across element boundaries. Element in Figure 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 to as


The curve point can now be evaluated using basis splines with a set of control points shown in Figure 1(b). However, if one adopts Catmull-Clark subdivision discretisation for analysis, this strategy of end element treatment will introduce additional ‘ghost-like’ degrees of freedom. To avoid this problem, the expression for  (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 1(b) illustrates the modified basis functions. The new basis functions do not possess the Kronecker delta property but do have the interpolating property at the boundary.

Figure 2: (a) Basis splines for interpolating element as a Catmull-Clark curve. (b) The construction of mirroring ghost point to maintain the location of the end point. Basis splines are reconstructed for interpolating an end-element in a Catmull-Clark curve. (c) Global basis functions for interpolating a curve.

The global basis functions for interpolating the curve in Figure 1 are shown in Figure 1(c). It is worth noting that the Catmull-Clark subdivision algorithm can not model arbitrary shapes exactly. It can only model a geometry using cubic B-splines. This property is significantly different to NURBS and motivates section 3.1 on geometry fitting.

2.2 Interpolating and evaluating Catmull-Clark subdivision surfaces

Catmull-Clark subdivision surfaces are derived as tensor-products of two Catmull-Clark subdivision curves but allowing for the presence of extraordinary vertices. The number of elements connected with the vertex is defined 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 [catmull1978recursively], 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 spline-based CAD tools, such as NURBS surfaces, need to link multiple patches. peters2008subdivision have proven that the limiting surface of the Catmull-Clark subdivision algorithm has continuity over the surface except at the extraordinary vertices. The corresponding surface location of the extraordinary vertices possesses continuity. This section will illustrate the methods of interpolating and evaluating both regular element and element with an extraordinary vertex in Catmull-Clark subdivision surfaces.

Element in a regular patch

Figure 2(a) 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 elements with control vertices. The surface point can be evaluated using the basis functions associated with these control points as


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 or are the basis functions defined in Equation (2) and presented in Figure 2(a). is the modulus operator and denotes the remainder operator which gives the remainder of the integer division.

Figure 2(b) shows the element patch of a subdivision surface element (shaded) that has an edge on the physical boundary. This type of element has only neighbour elements so that it belongs to an element patch which has 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 [cirakortiz2000]. However, this method involves additional degrees of freedom in numerical analysis. Instead, the curve basis functions in Equation (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 Figure 2(c).

Figure 3: Element patches for evaluating a Catmull-Clark subdivision element. (a) A regular element. (b) An element with a face on the boundary. (c) An element with two faces on the boundary.

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 3(a) 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 Figure 3(a). After applying one level of subdivision, new control points are generated and this element is subdivided into four sub-elements, as shown in Figure 3(b). The sub-elements , and 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 , where is the number of subdivision required and is the sub-element index shown in Figure 3(b). The new control point set is computed as


where is a selection operator to pick control points for the sub-elements. and are two types of subdivision operators. is the initial set of control points. The detailed approach is given in [stam1998exact] and also can be found in Appendix A.3. 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

Figure 4: (a) An irregular patch for a Catmull-Clark subdivision element with an extraordinary vertex. (b) One level of subdivision of the element patch divides the element into four sub-elements; three of them are in regular patches. (c) Successive subdivisions of the element until the evaluated point falls into a sub-element with a regular patch. (d) Adaptive Gauss quadrature scheme for the element with an extraordinary vertex.

Equation (9) can thus be rewritten as


where is the Catmull-Clark subdivision surfaces basis function. Define as a set of basis functions in an element with an extraordinary vertex and is a set of regular basis functions defined in Equation (7).

can be calculated in a vector form as


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 functions at a physical point involves two mappings. The first is from the physical domain to the parametric domain of an element with an irregular patch, . Because the irregular patch does not have the tensor-product nature, levels of subdivisions are required and the point is mapped to the parametric domain of a sub-element, . This second mapping is defined in Equation (10). The value of approaches positive infinity when approaches the extraordinary vertex which has the parametric coordinate . Hence the diagonal terms in the mapping matrix (15) tend to positive infinity as . This results in the basis functions not being differentiable at . This problem is termed singular configuration in [juttler2016numerical], and singular parameterisation in [takacs2012h2, nguyen2014comparative].

3 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 Section 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.

3.1 Geometry fitting

The Catmull-Clark subdivision curves and surfaces are CAD tools which can construct limiting curves and 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 as an approximation of this desired geometry is required. litke2001fitting introduced a method for fitting a Catmull-Clark subdivision surface to a given shape. In their work, both a least-squares fitting method and a quasi-interpolation method are adopted to determine a set of control points for a given curve or surface. The least-square fitting method is used here. One first chooses a set of sample points , where is the geometry, is the number of sample points. Each sample point should be evaluated using Catmull-Clark subdivision bases with control points as


where is the number of local basis functions. Then the set of sample points can be evaluated as


where is a set of control points with control points. is an evaluation operator of Catmull-Clark curves or surfaces. Set to ensure the sample points correspond to the control vertices and , then is a square matrix. The control points can be calculated as


If more sampling points are chosen than the required number of control points , then is invertible, a least-squares method is used to obtain a set of control points that minimises as


Figure 5 shows an example that fits a geometry using Catmull-Clark subdivision curves. The given curve is defined as . Figure 4(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 4(b) and 4(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.

(a) Fitting curve with sampling points.
(b) Fitting curve with sampling points.
(c) Fitting curve with sampling points.
Figure 5: The process of constructing Catmull-Clark subdivision surfaces to approximate a given curve.

3.2 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 Gauss points can exactly evaluate the integrals for polynomials of degree up to . The polynomial degree of a cubic B-spline function is . Because the basis functions of a Catmull-Clark subdivision element in regular element patch are generated as the tensor product of two cubic splines, 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 Equation (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 [juttler2016numerical], an adaptive quadrature rule, well suited to Catmull-Clark subdivision surfaces is adopted by integration at a number of levels of subdivided elements. With levels of subdivisions, the element is subdivided into sub-elements as shown in Figure 3(d). 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 sub-elements. With a number of subdivisions, the integration error can be reduced. In this work, is chosen in order to obtain sufficiently accurate values of the integrals.

3.3 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


We first use a test function to form a weak formulation for the Dirichlet boundary condition as


Using the Catmull-Clark subdivision curve basis functions in Equation (2) to discretise and and the same strategy for formulating the system matrix, one introduces a boundary mass matrix as


where 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 Catmull-Clark subdivision curves and only one-dimensional Gauss quadrature rule is used for integration. However, one uses the global degrees of freedom indices to assemble and , 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 , where is the system matrix, is the global coefficients vector to be solved for, and is global right-hand side vector. Then, we scale and 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.

4 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 ) in three dimensional space 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 can be calculated as the spatial gradient subtracted by its normal part as


where is the spatial gradient operator. Hence the surface Laplacian of is given by


where is the Hessian matrix of , and is the gradient of the normal vector, which is arranged in a matrix as


We define the total curvature at a surface point as the surface divergence of the normal, that is . For a given manufactured solution , the right hand side of Equation (27) can thus be computed as


5 Galerkin formulation

The weak formulation of problem (27) is


where is an admissible test function. The weak formulation is partitioned into number of elements, as


Discretising , and using the Catmull-Clark basis functions given in Equation (7) produces


where is the surface Jacobian for the manifold, given in a matrix form as


For details on the computation of see [peterson2005mapped] and for a discussion of superficial tensors such as in the context of Laplace-Beltrami equation, see [gurtin2002interface]. If the element contains an extraordinary vertex, the shape functions are replaced by in Equation (12). The surface gradient of the shape functions is computed as




Integrating the discrete problem using Gauss quadrature, the system of equations 34 becomes


where is the assembly operator and


is the number of quadrature points in each element, is the weight for quadrature point, is the number of elements and is the number of basis functions of the element. The basis functions are replaced by if the element 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 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


6 Numerical results

A ‘patch test’ [zienkiewicz1997finite] 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.

6.1 ‘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 Section 4 reduces to the Poisson problem expressed given by


This partial differential equation is solved on the square plate shown in Figure 5(a) with the essential boundary conditions


Four different manufactured functions for are used. The functions, analytical solutions for and their gradients are given in Table 1. We investigate both a regular and an irregular mesh. The regular mesh is a element patch without extraordinary vertices as shown in Figure 5(b). The boundary conditions are imposed using the penalty method.

(a) Plate problem setup.
(b) A regular mesh (Mesh 1).
(c) A mesh with extraordinary vertices (Mesh 2).
Figure 6: Schematic of the patch test on a plate.
Test No.
Table 1: Four test case functions for the plate problem. Test 1 has no right-hand side term, thus the analytical solution 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.

For Test 1, the right hand side so that . Solving the equation using the proposed Catmull-Clark subdivision method, the numerical result is exactly everywhere as shown in Figure 6(b)

. Thus passes the consistency test and the eigenvalue of the system matrix are all positive and non-zero after application of the essential boundary conditions. This also proves the stability of the method. The gradient

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 truncated, see Equations (3) and (4). In other words, the gradients of the function are expected to be constant at boundaries. Figures 6(a)6(c) and 6(e) show the numerical results for these tests. The results are smooth and capture the analytical solutions well. Figures 6(d) and 6(f) compare the numerical results of to the analytical solution for Test 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 () but captures the gradient in the region close to the bottom boundary () 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 in Test 4 is a cosine function which is non-polynomial and it behaves as a constant in both boundary regions shown in Figure 6(h). The Lagrangian elements only possess 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 smooth.

Figure 8 shows the plots of normalised global and errors against the element size. The normalised global error is defined by


where is the norm defined as . The normalised global error is computed as


where is the norm defined as . We set the element size of the coarsest mesh as . Then, the normalised element size for the refined meshes are . The convergence rate of Test 2 and 3 are sub-optimal at ( error) and ( error). The optimal convergence rate for cubic elements should be ( error) and ( error), where 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 Figure 5(c). Figures 7(a) and 7(b) show the plots of normalised element sizes against the and errors, respectively. The same convergence rates are obtained for Tests 2 and 3. However, the convergence rate of Test 4 is also reduced to ( error) and ( error). Figure 7(c) and 7(d) show the plots of normalised element sizes against and errors, respectively, for the mesh with an extraordinary vertex.

(a) Solution for Test 1.
(b) Gradient for Test 1.
(c) Solution for Test 2.
(d) Gradient for Test 2.
(e) Solution for Test 3.
(f) Gradient for Test 3.
(g) Solution for Test 4.
(h) Gradient for Test 4.
Figure 7: Solution and the gradient plotted along the line for the plate test. The numerical results are compared to the analytical solutions for tests 2, 3 and 4.
(a) Normalised global error (Mesh 1).
(b) Normalised global error (Mesh 1).
(c) Normalised global error (Mesh 2).
(d) Normalised global error (Mesh 2).
Figure 8: Convergence study for Tests 2, 3 and 4 using the regular mesh (Mesh 1) and the mesh with an extraordinary vertex (Mesh 2).

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 continuous.

6.2 The Laplace-Beltrami equation

The following sections will solve the Laplace-Beltrami equation (27) on different two dimensional Catmull-Clark subdivision manifolds. An analytical solution of function is manufactured as


where is a point on the two dimensional manifold in three dimensional space. Applying the Laplacian operator on gives


Then, the Hessian matrix is computed as


The second term in (33) includes the normal vector and its gradient which can not be computed analytically. In the present work, an projection is used to compute the coefficients of normal vectors associated with control points, denoted by , in order to numerically interpolate the normal vector derivatives at any surface points, thus


6.2.1 Cylindrical surface example

The first numerical example considered is a cylindrical surface. The analysis domain of the problem is the cylindrical surface shown in Figure 8(a). Surfaces fitting methods are used to construct the control mesh, see Section 3.1. The first level control mesh is shown in Figure 8(b). This has no extraordinary vertices. The Laplace-Beltrami problem on this manifold domain is solved using the Galerkin formulation presented in Section 5. The right-hand side function is computed using the definition in Equation (33). Figure 8(c) shows the numerical result which matches the manufactured analytical solution (47) very well.

(a) Geometry
(b) Control mesh
(c) Numerical result.
Figure 9: The geometry is given in Figure (a) and Figure (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. Figure (c) shows the numerical result on the cylindrical surface.

A convergence study is now conducted for this geometry. The refined control meshes are constructed using the least-squares fitting method described in Section 3.1. Figure 10 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 as cubic Lagrangian elements.

(a) Global error against number of degrees of freedom.
(b) Global error against normalised element size.
Figure 10: Convergence study for cylindrical example: comparison of the Catmull-Clark elements to the quadratic and cubic Lagrangian elements.

6.2.2 Hemispherical surface example

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

Figure 11: (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 on this surface.
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 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 Figure 10(b) as the initial control mesh. It can be used to generate a limit surface approximating a hemisphere, as shown in Figure 10(a), 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 -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 Figure 10(c) and 10(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 Figure 10(a). 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.

(a) against number of degrees of freedom.
(b) against normalised element sizes.
Figure 12: Convergence study: comparison of the Catmull-Clark elements with the quadratic and cubic Lagrangian elements.
Sparsity patterns
(a) Catmull-Clark subdivision elements (1280 elements, 1313 DoFs).
(b) Linear Lagrangian elements (1280 elements, 1313 DoFs).
(c) Cubic Lagrangian elements (1280 elements, 11617 DoFs).
Figure 13: Comparison of sparsity patterns between the Catmull-Clark elements and the Lagrangian elements

Figure 12(a) shows the sparsity pattern of the system matrix 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 12(b). Figure 12(c) is the sparsity patterns of cubic Lagrange discretisations. -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 rows and columns as the cubic Lagrangian discretisation but has a much smaller size.

6.3 Investigation of extraordinary vertices

The presence of extraordinary vertices leads to difficulties in integration as described in Section 3.2. Figures 13(a)13(c) and 13(e) show the point-wise errors at surface points for three levels of mesh refinement using the standard Gauss quadrature rule. 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. Figures 13(b)13(d) and 13(f) plot the point-wise errors on the same mesh analysed with the adaptive quadrature rule shown in Section 3.2. The errors around extraordinary vertices are now decreased.

(a) Initial mesh with standard Gauss quadrature.
(b) Initial mesh with adaptive Gauss quadrature.
(c) level subdivision mesh with standard Gauss quadrature.
(d) level subdivision mesh with adaptive Gauss quadrature.
(e) level subdivision mesh with standard Gauss quadrature.
(f) level subdivision mesh with adaptive Gauss quadrature.
Figure 14: Point-wise error plots over spherical surfaces.

The effects of the number and valence of extraordinary vertices are investigated. Figures 14(a)14(b) and 14(c) are three control meshes with different numbers of extraordinary vertices. Figure 14(a) shows a control mesh without an extraordinary vertex. Figure 14(b) 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 Figure 14(c) has seven extraordinary vertices, 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 computed in (33). Both standard and adaptive Gauss quadrature rules are used for all cases. Figures 15(a), 15(c) and 15(e) show the solution of on the surfaces constructed using the three meshes. Because of the similarity of the geometries and solutions, the three cases are used to investigate the effects of the number of extraordinary vertices in numerical results. The point-wise errors on the three surfaces are shown in Figures 15(b),  15(d) and 15(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 17 shows the convergence rates for the three cases. Meshes without extraordinary vertices can achieve the optimal convergence rate. 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 are dominated by the integration errors on the elements with extraordinary vertices, they can be reduced by using an adaptive quadrature rule. With the adaptive quadrature rule, the convergence rates are improved for the 4 and 7 extraordinary vertices cases.

(a) Regular mesh has no extraordinary vertex.
(b) Mesh has four extraordinary vertices (two valence 5 vertices and two valence 3 vertices).
(c) Mesh has seven extraordinary vertices (one valence 6 vertex, two valence 5 vertices and four valence 3 vertices).
Figure 15: Three control meshes with different number of extraordinary vertices.
(a) Numerical result on the cylindrical surface with no extraordinary vertex.
(b) Point-wise error on the cylindrical surface with no extraordinary vertex.
(c) Numerical result on the cylindrical surface with 4 extraordinary vertices.
(d) Point-wise error on the cylindrical surface with 4 extraordinary vertex.
(e) Numerical result on the cylindrical surface with 7 extraordinary vertices.
(f) Point-wise error on the cylindrical surface with 7 extraordinary vertex.
Figure 16: Numerical results and point-wise errors on the cylindrical surfaces constructed by thee different meshes. The white grids are used to indicate the locations of extraordinary vertices.
Figure 17: Convergence study: comparison of the Catmull-Clark subdivisions with different number of extraordinary vertices. The adaptive Gauss quadrature achieves better convergence rates over the standard Gauss quadrature rule.

7 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 new 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 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 ( error) and ( error).

In future work, this method will be used to compute solutions for problems requiring continuity such as the deformations of thin shells.


Appendix A Appendix

a.1 Catmull-Clark subdivision algorithm for curves

The Catmull-Clark algorithm successively refines a curve starting from an initial control polygon. After a number of subdivisions, the curve is limited to a cubic B-spline.

Figure 18: Computing new control points through the Catmull-Clark algorithm.

Figure 18 illustrates the subdivision algorithm. The control point in the level of refinement is computed from the upper level control points as


Point is the mid- point of , and is called an ‘edge point’. The control point is computed as


To compute this point, one needs to connect the mid-points of and . The point is the mid-point of the connecting line. This type of point is called ‘vertex point’. Each ‘vertex point’ is associated with an upper level control point. Figure 19 shows two levels of refinements using the Catmull-Clark algorithm and the limiting result which is a Catmull-Clark subdivision curve.

Figure 19: Schematics of the Catmull-Clark Algorithm

a.2 Catmull-Clark subdivision algorithm for surfaces

The application of the Catmull-Clark algorithm to surfaces follows in a similar manner to curves. One face in the original control mesh is split into four new faces. For a closed surface, the numbers of faces and control vertices are doubled. Figure 20 shows an example of generating a new control mesh through the Catmull-Clark algorithm.

Figure 20: Catmull-Clark subdivision algorithm for surfaces.

Similar to the one-dimensional Catmull-Clark curve, the new refined control points can be classified into three types: face points, edge points and vertex points. The face points in the

refinement are computed as


where and are indices of control points for orthogonal directions. The face point is the central point of the original face. The edge point is computed as


or likewise


The ‘vertex point’ is computed as


Equipped with these formulae, the new control points on the level of refinement can be computed as:


is a subdivision operator – a matrix consisting of a set of weights. Each weight is associated with a control point in . The weight distributions for different types of control points are shown in Figure 21. The weight distributions for extraordinary point are shown in Figure 22. After successive levels of refinements, a smooth B-spline surfaces is obtained.

(a) Edge point case 1.
(b) Edge point case 2.
(c) Vertex point.
(d) Face point.
Figure 21: The weight distribution for computing different types of new control points.
Figure 22: Weight distributions for computing an extraordinary point with valence .

a.3 Computing control point set for sub-elements

We denote the control points of an irregular patch in Figure 3(a) as a set . The initial control points of the patch are expressed as


Through one level of subdivision we generate new control points ( is the valence), as shown in Figure 3(b), denoted by


The subdivision step is represented as


where is the subdivision operator given by


The terms , , , and are defined in [stam1998exact] and is given in Equation 57. To evaluate the sub-element , and in Figure 3(b), one needs to pick control points out of the new control point patch. A selection operator for sub-element and is used to select the necessary control points from , that is


Then a surface point can be evaluated with the cubic spline basis functions as


As shown in Figure 3(c), after successive subdivisions, the non-evaluable element can be limited to a negligible region.

Assume the target point has parametric coordinates . One first determines how many subdivisions are required for this point by:


The sub-element index is determined as


The surface point is located in the regular sub-element after the refinement. The patch for this element is picked with the selection operator as


The enlarged set contains control vertices, which is generated from the subdivision of as


The set has control vertices. It is successively refined from the initial set as


where is a square matrix operator which subdivides the patch for computing the new patch for the irregular element, and is defined by