Subdivision surfaces with isogeometric analysis adapted refinement weights

by   Qiaoling Zhang, et al.
University of Cambridge

Subdivision surfaces provide an elegant isogeometric analysis framework for geometric design and analysis of partial differential equations defined on surfaces. They are already a standard in high-end computer animation and graphics and are becoming available in a number of geometric modelling systems for engineering design. The subdivision refinement rules are usually adapted from knot insertion rules for splines. The quadrilateral Catmull-Clark scheme considered in this work is equivalent to cubic B-splines away from extraordinary, or irregular, vertices with other than four adjacent elements. Around extraordinary vertices the surface consists of a nested sequence of smooth spline patches which join C^1 continuously at the point itself. As known from geometric design literature, the subdivision weights can be optimised so that the surface quality is improved by minimising short-wavelength surface oscillations around extraordinary vertices. We use the related techniques to determine weights that minimise finite element discretisation errors as measured in the thin-shell energy norm. The optimisation problem is formulated over a characteristic domain and the errors in approximating cup- and saddle-like quadratic shapes obtained from eigenanalysis of the subdivision matrix are minimised. In finite element analysis the optimised subdivision weights for either cup- or saddle-like shapes are chosen depending on the shape of the solution field around an extraordinary vertex. As our computations confirm, the optimised subdivision weights yield a reduction of 50% and more in discretisation errors in the energy and L_2 norms. Although, as to be expected, the convergence rates are the same as for the classical Catmull-Clark weights, the convergence constants are improved.



There are no comments yet.


page 5

page 6

page 9

page 10

page 11


A convergent evolving finite element algorithm for Willmore flow of closed surfaces

A proof of convergence is given for a novel evolving surface finite elem...

Convergence of Dziuk's semidiscrete finite element method for mean curvature flow of closed surfaces with high-order finite elements

Dziuk's surface finite element method for mean curvature flow has had si...

Tuned Hybrid Non-Uniform Subdivision Surfaces with Optimal Convergence Rates

This paper presents an enhanced version of our previous work, hybrid non...

Tangential Errors of Tensor Surface Finite Elements

We discretize a tangential tensor field equation using a surface-finite ...

Reissner-Mindlin shell theory based on tangential differential calculus

The linear Reissner-Mindlin shell theory is reformulated in the frame of...

Discrete Total Variation of the Normal Vector Field as Shape Prior with Applications in Geometric Inverse Problems

An analogue of the total variation prior for the normal vector field alo...

Aspects of Isogeometric Analysis with Catmull-Clark Subdivision Surfaces

An isogeometric approach for solving the Laplace-Beltrami equation on a ...
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

Isogeometric analysis aims to provide a seamless engineering design-analysis workflow by using a single common representation for geometric modelling and analysis. This is usually achieved by representing geometry and discretising analysis models with the same kind of basis functions Hughes et al. (2005). The prevailing feature-based CAD modelling systems rely on trimmed NURBS and boundary representations (B-Reps). The resulting non-watertight geometries consisting of several trimmed patches pose unique challenges to finite element analysis. As a generalisation of splines, subdivision surfaces can provide watertight representations for geometries with arbitrary topology. After their early success in computer animation and graphics they are now supported in many CAD systems, including Catia, PTC Creo and Autodesk Fusion 360. Before the advent of isogeometric analysis, it had already been realised that subdivision surfaces provide also ideal basis functions for finite element analysis, in particular, of thin-shells Cirak et al. (2000); Reif and Schröder (2001); Cirak and Ortiz (2001); Cirak et al. (2002), see also more recent work Cirak and Long (2011); Bandara and Cirak (2018).

Subdivision schemes for generating smooth surfaces were first described in the late 1970s as an extension of low degree B-splines to control meshes with non-tensor-product connectivity 

Doo and Sabin (1978); Catmull and Clark (1978)

. In subdivision a geometry is described with a control mesh and a limiting process of repeated refinement. For parts of the mesh containing only regular vertices, with each adjacent to four quadrilateral faces, the refinement rules are adapted from knot insertion rules for B-splines. For the remaining parts with extraordinary vertices the refinement rules are chosen such that they yield in the limit a smooth surface. Subdivision refinement is a linear mapping of coordinates of the coarse control mesh to the coordinates of the refined mesh with a subdivision matrix. Hence, the local limit surface properties can be inferred from the eigenstructure of the subdivision matrix after a discrete Fourier transform 

Doo and Sabin (1978); Reif (1995). The

continuity of the surface and its curvature behaviour at the extraordinary vertex depend on eigenvalues and the ordering, i.e. Fourier indices, of the corresponding eigenvectors. In turn, both depend on the coefficients of the subdivision matrix that encodes the specific refinement rules applied.

As known, around extraordinary vertices short-wavelength surface oscillations, i.e. ripples, may occur irrespective of continuity and boundedness of curvature Peters and Reif (2004); Karčiauskas et al. (2004). There have been many attempts to improve the fairness of subdivision surfaces, that is, to minimise curvature variations, by carefully tuning the refinement rules, earlier works include Halstead et al. (1993); Kobbelt (1996). More recently, in Augsdörfer et al. Augsdörfer et al. (2006) the refinement rules for Catmull-Clark and other quadrilateral schemes have been optimised such that the variation of the Gaussian curvature is minimised while ensuring bounded curvatures. Different from the direct search method used in Augsdörfer et al. (2006), the refinement rules can also be obtained from a nonlinear constrained optimisation problem. Barthe et al. Barthe and Kobbelt (2004) apply such a procedure to triangular Loop and -subdivision schemes with a multi-objective cost function comprised of terms penalising divergence of curvatures and aiming local quadratic precision. In Ginkel et al. Ginkel and Umlauf (2008) a fairness increasing cost function containing the third derivatives of the surface in combination with continuity and bounded curvature constraints is optimised.

In the present paper, we optimise the subdivision refinement rules so that their approximation properties are improved when used in finite element analysis of thin-shells. Thin-shells are prevalent in many engineering applications, most prominently in aerospace, automotive and structural engineering, and are equivalent to thin-plates when their unstressed geometry is planar Ciarlet (2005)

. The thin-shell energy functional, and weak form, depend on the second order derivatives of the stressed surface. Consequently, it is crucial to reduce any short-wavelength oscillations in the subdivision surface. As the included examples demonstrate, meshes with extraordinary vertices usually lead to lower convergence rates than meshes with tensor-product connectivity. For obtaining the improved isogeometric analysis adapted refinement rules we postulate a constrained optimisation problem with a cost function measuring the errors in approximating cup- and saddle-like quadratic shapes. Three of the weights in the Catmull-Clark subdivision scheme around an extraordinary vertex are chosen as degrees of freedom for optimisation. As constraints the

continuity of the surface is strictly enforced and bounded curvatures are enforced as long as non-negative real weights are feasible. The eigenstructure of the subdivision matrix is extensively used in formulating the optimisation problem as usual in previous related work (Peters and Reif, 2008, Chapter 4,5) and (Sabin, 2010, Chapter 15). We compute the eigenvalues and eigenvectors numerically after applying a discrete Fourier transform that exploits the local circular symmetry around the extraordinary vertex. The local parameterisation of the subdivision surface required for evaluating the finite element integrals and the cost function is obtained with the algorithm proposed by Stam Stam (1998). Two sets of optimised weights for cup- and saddle-like shapes are obtained. The weights for finite element analysis are chosen depending on the dominant shape of the solution field around an extraordinary vertex.

For completeness, we note that subdivision is not the only approach for creating smooth surfaces on arbitrary connectivity control meshes. Over the years numerous and smooth constructions with have been proposed, too many to name here. The search for sufficiently flexible smooth surface representations, especially with and , is still open. It is worth mentioning that none of the existing constructions is widely used in commercial CAD systems. This may well be because their implementation is too complicated. The application of basis functions resulting from smooth constructions for isogeometric analysis is currently a very active area of research. For instance, the utility of constructions with NURBS has recently been explored in Nguyen et al. (2016); Collin et al. (2016); Kapl et al. (2017). Alternatively, constructions relying on manifold-based surface constructions Majeed and Cirak (2017); Ying and Zorin (2004); Grimm and Hughes (1995) and constructions relying on singular parameterisations have also been investigated Toshniwal et al. (2017); Nguyen and Peters (2016); Reif (1998). Some of these schemes are able to provide optimal convergence rates.

The outline of this paper is as follows. In Section 2 the Catmull-Clark subdivision is introduced, with a review of the relevant theory on eigenanalysis of the subdivision matrix. Specifically, the necessary conditions for smoothness and boundedness of the curvature are motivated, and the local parameterisation of subdivision surfaces using the characteristic map is introduced. These are all classical results and concepts which are mostly unknown in isogeometric analysis. In Section 3 the proposed constrained optimisation problem and its numerical solution are discussed. Two sets of subdivision weights are derived that minimise the thin-plate energy norm errors in approximating locally cup- and saddle-like shapes. Subsequently, it is shown how a finite element solution can be locally decomposed into cup- and saddle-like components. Depending on this decomposition and the following choice of optimal weights, a second more accurate finite element analysis can be performed. In Section 4 the proposed approach is applied to transversally loaded thin-plate problems using meshes with extraordinary vertices and the convergence of the errors in and energy norms is reported.

2 Catmull-Clark subdivision surfaces

2.1 Refinement weights and the subdivision matrix

Catmull-Clark subdivision is a generalisation of cubic tensor-product B-splines to unstructured meshes Catmull and Clark (1978). On non-tensor-product meshes the number of faces connected to a vertex, i.e. valence , can be different from four. The vertices with are referred to as extraordinary or star vertices. During subdivision refinement each face of the control mesh is split into four faces and the coordinates of the old and new control vertices are computed with the subdivision weights given in Figure 1. The weights in each of the three diagrams have to be normalised so that they add up to one. The unnormalised weights assigned to the extraordinary vertex (empty circle) are denoted by , and respectively. For and bivariate cubic B-splines the three weights take the values , and . The new vertices introduced by the subdivision process are all regular (with ) and the total number of irregular vertices in the mesh remains constant. That is, the irregular vertices are more and more surrounded by regular vertices.

Figure 1: Subdivision weights for the Catmull-Clark scheme with the empty circle denoting the extraordinary vertex. The weights in each of the three diagrams have to be normalised so that they add up to one. For Catmull-Clark scheme the three weights take the values , and , where is the valence.

In order to study the smoothness behaviour of subdivision surfaces near an extraordinary vertex, it is sufficient to consider only the vertices in its immediate vicinity. A -neighbourhood of a vertex is formed by the union of faces that contain the vertex. The -neighbourhood is defined recursively as the union of all -neighbourhoods of the -neighbourhood vertices. It is assumed that the considered -neighbourhood has only one single extraordinary vertex located at its centre. The -neighbourhood control vertices at the refinement level are mapped to control vertices with the subdivision matrix ,


The square subdivision matrix can be readily derived from the weights indicated in Figure 1. The control point coordinates at level are arranged in this form


with each row containing the coordinates of one control point with the index .

2.2 Eigendecomposition of the subdivision matrix

For the tuning approach to be introduced in Section 3, it is necessary to consider the -neighbourhood around an extraordinary vertex, see Figure 2.

Figure 2: Three-rings of faces around an extraordinary vertex with valence and the numbering of the vertices. In the index pair the first is the segment number and the second is the vertex number.

The -neighbourhood consists of segments with each segment containing vertices, excluding the extraordinary vertex with index . Hence, there are vertices so that the subdivision matrix has the dimensions . In establishing the subdivision matrix it is assumed that the index pair is converted to a scalar index as . The eigenvalues and eigenvectors of are closely related to the smoothness and other properties of the subdivision surface. The eigendecomposition of the asymmetric subdivision matrix reads




where are the eigenvalues and , are the right and left eigenvectors respectively. Throughout the paper it is assumed that the eigenvalues are sorted in descending order with largest being . The subdivision matrix has a cyclical structure due to the cyclic symmetry of the weights given in Figure 1. We assume that the vertices in the -neighbourhood are enumerated according to Figure 2.

Because of the cyclical structure, the eigendecomposition of can be best computed with a discrete Fourier transform (DFT). As pioneered in Doo and Sabin (1978), DFT is crucial in identifying the different geometric shapes described by the different eigenvectors. For transforming the following extended DFT matrix is considered:


with the complex number where

, the identity matrix 

of size

and the zero vector 

of size . The first row and column of  have been introduced for the extraordinary vertex. In obtaining (5) the standard relations and with the complex conjugate are used. The inverse transform is obtained by replacing with its complex conjugate . The subdivision matrix is Fourier transformed according to


leading to a block diagonal matrix


where the blocks are of size for and of size for . Due to its block-diagonal structure the eigendecomposition of the transformed matrix


can be more readily determined. Namely, it is sufficient to consider the eigenvalue problems for each of the blocks separately, i.e.,


where the eigenvalues within each block are also sorted in descending order, i.e, is the largest eigenvalue in the block . The eigenvalues and eigenvectors and are the union of all the eigenvalues and the block-wise eigenvectors. In obtaining the eigenvectors and , each of size , the corresponding block-wise vectors and

are suitably padded with zeros. The vectors

and are of size for and of size for . Moreover, the subdivision matrix and its Fourier transform have the same eigenvalues and their eigenvectors are related by


Each block corresponds to a specific rotational frequency . As pointed out, the eigenvectors and can have non-zero entries only in the components corresponding to a specific and . Hence, the transformation of and according to (10) yields always a column of each of which corresponds to a specific rotational frequency111The first columns and rows of and are assigned to the extraordinary vertex and do not represent harmonics.. To this end, recall the Euler identity


Hence, for a fixed angular frequency the vectors and will assign each control vertex with a fixed index , c.f. Figure 2, a value that oscillates with the angular frequency while circumnavigating the extraordinary vertex by incrementing .

Furthermore, for geometric interpretation of the eigendecomposition it is helpful to realise that most of the eigenvalues have the multiplicity of two. That is, the eigenvalues and are identical in the blocks . The corresponding eigenvectors and have the same eigenfrequency because the columns and of the DFT matrix  are the complex conjugates of each other.

2.3 Limit analysis and smoothness

The eigenvalues and eigenvectors and of the Fourier transformed subdivision matrix have to satisfy certain conditions for a subdivision scheme leading to a smooth well-defined surface, see Doo and Sabin (1978); Peters and Reif (2008). To understand this, consider the projection of control mesh vertex coordinates at subdivision level

into the eigenspace of the subdivision matrix, using the orthogonality of left and right eigenvectors (



where each of the scalar products yield a row vector with components. Subdividing the -neighbourhood in the eigenspace, while considering the eigendecomposition (4), gives


Hence, the repeated subdivision of the -neighbourhood can be simply achieved with


From this equation it is evident that the properties of a subdivision surface are widely governed by the eigenstructure of the subdivision matrix. The subdivision matrix

is a stochastic matrix, i.e. only positive entries and each row adds up to 1, so that its largest eigenvalue is

and the components of the corresponding eigenvector are all equal to . In the limit all control vertices converge to . The first term in (14) can be eliminated by translating the initial control vertex coordinates by . Without loss of generality, in the following we assume that the coordinate system for -neighbourhood has been chosen so that the first term in (14) is zero, that is,


For a (symmetric) -continuous subdivision surface the subdominant eigenvalues and have to satisfy the following relationship:


In addition, the corresponding eigenvectors , , and have to come from the eigendecomposition of the blocks and  Doo and Sabin (1978); Reif (1995). As discussed in Section 2.2, owing to the symmetry properties of the Fourier transformation, and have the same eigenvalues, and the eigenvectors and are the complex conjugates of and . All the four eigenvectors , , and are usually complex and have the angular frequency . A set of real eigenvectors each of size representing vertex values can be obtained as the linear combination of the complex ones, e.g., with and where . To avoid a proliferation of symbols we will use the same symbols for the so-computed real and complex eigenvectors.

A necessary condition for the -continuity of a subdivision surface (with no artificial flat spots) is that the subsubdominant eigenvalues satisfy


and the corresponding eigenvectors come from the eigendecomposition of the blocks , and  Peters and Reif (2008). Remember that is naturally satisfied due to the duplicity of eigenvalues from blocks and when , as mentioned in Section 2.2.

2.4 Characteristic map

As first proposed in Reif Reif (1995) the characteristic map provides a means for parameterisation of the surface generated by a subdivision scheme. Parameterisation of the subdivision surface, at least a local one, is essential in order to associate the so far discrete representation based on control vertices with a continuous differentiable representation. This is, for instance, required for finite element analysis using subdivision surfaces. The characteristic map is defined using the two real right eigenvectors and corresponding to the subdominant eigenvalue . The characteristic control mesh shown in Figure 3 representing the 3-neighbourhood around an extraordinary vertex has the coordinates


where the third out-of-plane coordinate is chosen as . As discussed in Section 2.2, recall that the two eigenvectors and have the angular frequency and have been chosen so that they are orthogonal in the plane spanned by the corresponding two complex eigenvectors. This can be done without loss of generality because the subdivision construction is invariant under affine transformations. Hence, the coordinates of control vertices oscillate with in the horizontal direction and with in the vertical direction leading to the shown characteristic control mesh.

Figure 3: Characteristic control mesh of Catmull-Clark scheme for valence (left) and (right).

As suggested in Reif Reif (1995), the planar surface described by the characteristic control mesh can be used for the parameterisation of subdivision surfaces. To this end, first consider the subdivision refinement of the characteristic mesh. According to (15) and the orthogonality of left and right eigenvectors (4), the subdivision refinement of the characteristic mesh simply yields a scaled version of the same mesh:


Hence, the refined control mesh is simply obtained by scaling the control mesh by . Repeated subdivision yields repeated scaling of the control mesh. This combined with the fact that the Catmull-Clark scheme leads to bivariate cubic B-splines in patches with only ordinary vertices is used for parameterising the subdivision surface. During subdivision refinement each patch is split into four patches. In particular, in the patches adjacent to the extraordinary vertex three of the created patches have only regular vertices and can be parameterised with bivariate cubic B-splines.

With repeated refinement more and more of the subdivision surface can be parameterised with cubic B-splines. A practical algorithm for efficient implementation of this parameterisation has been introduced in Stam Stam (1998). Without going into details we define the bijective characteristic map


with and , which maps a set of square domains with representing the faces in the control mesh into the characteristic domain . The smooth parameterisation provided by the characteristic map is illustrated in Figure 4.

Figure 4: Characteristic map from a unit square to the characteristic domain.

With the subdivision basis functions , consisting of cubic B-splines and obtained according to Stam (1998), the characteristic map can be written as


For brevity, in the following we omit the face index in the basis function and write .

3 Optimisation of subdivision weights

We aim to modify the subdivision weights and of the Catmull-Clark scheme to improve its approximation properties when used in finite element analysis. As known in CAD, not all parameters yield visually appealing surfaces even when they give -continuous surfaces, for a quantitative analysis see Peters and Reif (2004). Small surface oscillations, i.e. ripples, appear when the represented surface is not planar.

3.1 Preliminaries

First, we consider the representation of a polynomial scalar field over the characteristic domain . It is assumed that the scalar field is given in the form


where and the functions on the second line are introduced for notational convenience. The chosen functions , , , , and can represent all quadratics and their choice will be discussed further below. The approximation of over can be studied by comparing with it the limit surface resulted from the control points


where the third coordinate is a vector formed by the scalar function evaluated at the vertex locations , row by row. The linear functions  and can be exactly represented so that we are mainly concerned about the quadratic terms , and .

The specific form of the quadratic functions in (22) is motivated by the eigenstructure of the subdivision matrix , see Section 2.2. Specifically, the control point values , and and the corresponding control point values , and have matching angular frequencies over the 3-neighbourhood of the extraordinary vertex222In order for the phase to match, the indexing of the vertices has to begin along the edge aligned with the -axis.. It is straightforward to confirm the orthogonality relations


According to (12), the projection of the control vertex coordinates into the eigenspace of the subdivision matrix , while neglecting the terms with higher orders than quadratic, yields


With the eigendecomposition (4) the subdivision refinement of this control mesh gives


That is, the subdivision refinement of the first two components yields the characteristic domain and the third component yields the graph of the surface approximating , see Figure 5. The corresponding limit surface has, according to (21), the following form:


To compare the shapes and quantitatively, we introduce the thin-plate energy norm


with the Poisson ratio .

Figure 5: Quadratic shapes over a characteristic control mesh with valence . Note that  has the same shape like , but only rotated by in the -plane.

Moreover, the necessary conditions for -continuity given in (17), repeated here for convenience,

can now be related to the curvature of the three quadratic limit surfaces resulted from repeated refinement of with . In order for the limit surfaces to have finite curvature at the extraordinary vertex, when the first two control vertex components scale with the third has to scale with .

3.2 Constrained optimisation

The constrained optimisation problem for determining the subdivision weights , and that minimise the error in approximating quadratic surfaces is formulated as

subject to: (29b)

with and the constraints representing the necessary -continuity conditions (16) and (17). As mentioned in Section 2.3, owing to the symmetries of the DFT, the constraints (29b) and (29e) are automatically satisfied. Hence, the constraints reduce to two independent equations for the three unknowns. To reduce the constrained optimisation problem into an unconstrained one, it is convenient to first solve the nonlinear system of equations


That is, to determine the dependence of the weights , and  on the variable . To solve (30) we use in our implementation the Python library SciPy, to be more specific, the quasi-Newton method with a BFGS update with a suitable cost function. However,  and  can become complex for some  values Augsdörfer et al. (2006). For Catmull-Clark, it is smaller  values which result in complex weights. For instance, there is no real solution for   and  for in case of valence . Instead of excluding  values leading to complex weights we relax the second constraint (30b) by considering the modified constraint equations


In implementations where the boundedness of curvature must be satisfied, one can constrain the value in a valid range or consider more degrees of freedom for optimisation Augsdörfer et al. (2006); Donatelli et al. (2016). In our numerical experiments, we found that considering the modified constraint equations leads to smaller energy norm errors in comparison to constraining the range of possible values. After solving (30) or  (31) and determining , , and the constrained optimisation problem (29a) can now be restated as an unconstrained problem


which is a one-dimensional optimisation problem that can be solved by direct search.

3.3 Optimised weights for valence vertices

As an example for obtaining optimised weights, we consider the valence vertex case. The proposed optimisation follows the same procedure regardless of valence. It is sufficient to consider only the approximation of the quadratic functions and with cup-like and saddle-like geometries, respectively. The function has the same saddle-like geometry as , only rotated by  in the -plane. During optimisation the thin-plate energy norms in (32) are evaluated in the 2-neighbourhood of the extraordinary vertex, which is the same as the support size of the basis functions. This explains why we consider -neighbourhood around an extraordinary vertex, see Figure 2, because the evaluation in the second-ring elements needs the third-ring control vertices.

Figures (a)a and (b)b show the relative energy norm errors in approximating cup- and saddle-like geometries, respectively, when the subdominant eigenvalue and number of Gauss integration points are varied. It can be seen that while has a significant influence on the error the number of integration points appears to be irrelevant. Figure 7 shows the relative energy norm error both in cup- and saddle-like geometries when integration points are used. In comparison to Catmull-Clark weights, also indicated in Figure 7, for the optimised subdivision weights lead to a reduction of errors in both cup- and saddle-like geometries. Moreover, the most optimal value for the cup-like geometry is and for the saddle-like geometry is , see Table 1 for the values of the optimised weights. According to Peters and Reif Peters and Reif (2004), the obtained subdivision surfaces, same as original Catmull-Clark scheme, are -continuous at extraordinary vertices and everywhere else, because the eigenvalues satisfy the required relations and the characteristic map is regular and injective.

Cup 13.4575 0.999938 0.999938 0.550
Saddle 13.9851 0.824885 0.824885 0.585
Original Catmull and Clark (1978) 15 1 1 0.550
Table 1: Optimised weights in Catmull-Clark subdivision scheme for valence vertices.
Figure 6: Relative energy norm error in dependence of the sub-dominant eigenvalue and number of integration points.
Figure 7: Relative energy norm error for cup- and saddle-like geometries in dependence of the sub-dominant eigenvalue and for the Catmull-Clark(CC) scheme. The number of integration points used for all data values is . The most optimal value for the cup-like geometry is , while for the saddle-like geometry it is . See Table 1 for the values of the optimised weights.

3.4 Application-dependent choice of refinement weights

When subdivision surfaces are used for finite element analysis, the solution field has quite often a mixture of both cup- and saddle-like components. And the solution field at a specific extraordinary vertex is only known after the finite element analysis. Therefore, in a first step we use the optimal weights for the cup-like geometry to obtain an initial finite element solution. Afterwards, for each extraordinary vertex with valence , a local shape decomposition is performed to determine whether the local solution is cup or saddle dominated. If the cup component dominates, the optimal weights for the cup-like geometry are chosen. If instead the saddle component dominates, the optimal weights for the saddle-like geometry are chosen. After the optimal weights for each extraordinary vertex are chosen, a second finite element analysis is performed to obtain the final solution with smaller discretisation errors.

For thin-plate and thin-shell finite element problems the decomposition of a solution into cup- and saddle-like components may be accomplished as described in the following. Suppose that the local finite element solution in a 3-neighbourhood of an extraordinary vertex is denoted as and has the dimensions . In a coordinate system centred at the limit position of the extraordinary vertex, according to  (15), we can write


The corresponding limit surface has at the extraordinary vertex the normal vector , defined by


where the vectors and represent the two, usually non-orthogonal, tangent vectors. Multiplying (33) with the normal vector gives by eliminating its first two terms


The vector represents the out-of-plane coordinates of the control points in a coordinate system aligned with the tangent plane at the extraordinary vertex. The corresponding limit surface over the characteristic domain has the following representation:


and can be approximated with quadratic functions , and , see (22), such that


where denotes the refinement level required to evaluate at the point using the Stam Stam (1998) algorithm. The refinement level dependent factors always vanish when, as required for curvature continuity, for . After computing the energy densities, i.e. the integrand in (28), for each of the three quadratic components their ratio can be determined. Moreover, the two last terms with saddle-like geometries are energetically equivalent so that their components can be combined. This gives the following ratio between cup- and saddle-like energies:


In numerical computations the ratio is used to decide which set of subdivision refinement weights to use.

4 Examples

We consider the finite element analysis of thin plates to demonstrate the benefits of the optimised subdivision weights over Catmull-Clark weights. The plates are square shaped, simply supported and subjected to either uniform or sinusoidal distributed transversal loads, see Table 2. The thin plate energy functional depends on the second derivatives of the displacement field, c.f. (28). Hence, the accurate approximation of the quadratic terms in the solution field is crucial. For details of finite element implementation we refer to Cirak et al. Cirak et al. (2000); Cirak and Long (2011). The analytical solutions of all the computed problems are known and can be found in Timoshenko et al. (Timoshenko and Woinowsky-Krieger, 1959, Chapter 5).

Length ,
Young’s modulus
Poisson’s ratio
Uniform loading
Sinusoidal loading
Table 2: Geometry, material and loading of the computed thin plates.

Two different unstructured control meshes shown in Figure 8 are used in the numerical computations. Both meshes have extraordinary vertices with valences and . We do not use optimised subdivision weights for valence vertices with subdominant eigenvalues . During subdivision refinement their 1-neighbourhoods shrink faster than the 1-neighbourhoods of other vertices with , see also Figure 3. Hence, it can be expected that the benefits of optimising the weights of vertices with will be negligible. For each valence vertex the optimal subdivision weights are chosen independently based on the dominant component of the quadratic shape at the vertex.

Figure 8: Initial unstructured coarse control meshes with level .

We demonstrate the benefits of the optimised subdivision weights over the Catmull-Clark weights by plotting the convergence of and energy norm errors. The successively refined meshes are obtained by subdivision using the Catmull-Clark weights. In all examples Gauss quadrature points are used to evaluate the finite element integrals, which appears to be sufficiently accurate as shown in Figure (a)a and Figure (b)b. See also Jüttler et al. (2016) for a systematic study on numerical integration of subdivision surfaces.

4.1 Uniform loading, symmetric unstructured mesh

As the first example, we compute the deformation of a simply supported square plate subjected to uniform transversal loading. The plate is discretised with the symmetric unstructured mesh shown in Figure (a)a. Since the displacement field is usually not known prior to a finite element analysis, we use the optimal weights for cup-like shapes to solve the plate bending problem on a level control mesh. Afterwards, for extraordinary vertices with valence we decompose the local displacement field energetically to determine whether it is cup- or saddle-dominated and choose the optimal subdivision weights accordingly. For the considered uniform loading the shape decomposition shows that saddle dominates at all valence vertices with cup-saddle ratio . Therefore, we choose the optimal weights for saddle-dominated shapes for all vertices and study the convergence of the finite element solution using meshes from levels to .

Figure 9 shows the level control mesh and the deformed plate. The convergence of and energy norm errors are plotted in Figure 10. The optimised refinement weights reduce the norm error by more than and the energy norm error by more than in comparison to Camtull-Clark subdivision weights.

Figure 9: Control mesh with fourfold symmetry and the deformed plate under uniform loading. The shown control mesh is obtained by subdividing the symmetric coarse control mesh in Figure (a)a twice using Catmull-Clark weights.
Figure 10: Uniform loading with symmetric unstructured mesh. Saddle dominates () at all valence vertices. Optimisation reduces the  norm error by more than  and energy norm error by more than . See Table 1 for the values of the optimised weights corresponding to .

4.2 Sinusoidal loading, symmetric unstructured mesh

Next, we compute the deformation of a simply supported square plate discretised with the unstructured mesh shown in Figure (a)a and subjected to sinusoidal loading. Compared with the first example, the only difference is that sinusoidal loading is applied instead of a uniform loading. After the first finite element analysis with the optimal weights for the cup-like geometry the local shape decomposition shows that the cup component dominates at all valence vertices with cup-saddle ratio . Therefore, we choose the optimal weights for cup-dominated shapes for all vertices and study the convergence of the finite element solution using meshes from levels to .

Figure 11 shows the level control mesh and the deformed plate. The convergence of and energy norm errors are plotted in Figure 12. The optimised refinement weights reduce the norm error by more than and the energy norm error by more than in comparison to Camtull-Clark weights.

Figure 11: The deformed plate under sinusoidal loading on a level symmetric control mesh.
Figure 12: Sinusoidal loading with symmetric unstructured mesh. Cup dominates () at all valence vertices. Optimisation reduces the error by more than and energy norm error by more than . See Table 1 for the values of the optimised weights corresponding to .

4.3 Sinusoidal loading, asymmetric unstructured mesh

In this last example, we compute the deformation of a simply supported square plate discretised with the asymmetric unstructured mesh shown in Figure (b)b and subjected to sinusoidal loading. As in the first two examples, we obtain the first finite element solution using the optimal weights for cup-like shapes on a level control mesh. Subsequently, for each extraordinary vertex with valence the local displacement field is decomposed to determine whether it is cup or saddle dominated and the optimal weights are chosen accordingly. The local shape decomposition shows that the local solution at vertex is cup dominated with cup-saddle ratio and at vertex