Tuned Hybrid Non-Uniform Subdivision Surfaces with Optimal Convergence Rates
Xiaodong Wei Xin Li Yongjie Jessica Zhang Thomas J. R. Hughes
math.NA cs.NA


language=Matlab,
    basicstyle=,
    breaklines=true,
    morekeywords=matlab2tikz,
    keywordstyle=,
    morekeywords=[2]1, keywordstyle=[2],
    identifierstyle=,
    stringstyle=,
    commentstyle=,
    showstringspaces=false,
    numbers=left,
    numberstyle=,
    numbersep=9pt, 
    emph=[1]for,end,break,emphstyle=[1], 




International Journal for Numerical Methods in Engineering



2]Xiaodong Wei
1]Xin Licor1
3]Yongjie Jessica Zhang
4]Thomas J.R. Hughes
[cor1]Corresponding author, lixustc@ustc.edu.cn, tel: +86-551-63607202
[1]School of Mathematical Sciences, University of Science and Technology of China, Hefei, Anhui, China
[2]Institute of Mathematics, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
[3]Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA
[4]Oden Institute, The University of Texas at Austin, Austin, TX 78712, United States



This paper presents an enhanced version of our previous work, hybrid non-uniform subdivision surfaces <cit.>, to achieve optimal convergence rates in isogeometric analysis. We introduce a parameter λ (1/4<λ<1) to control the rate of shrinkage of irregular regions, so the method is called tuned hybrid non-uniform subdivision (tHNUS). Our previous work corresponds to the case when λ=1/2. While introducing λ in hybrid subdivision significantly complicates the theoretical proof of G^1 continuity around extraordinary vertices, reducing λ can recover the optimal convergence rates when tuned hybrid subdivision functions are used as a basis in isogeometric analysis. From the geometric point of view, the tHNUS retains comparable shape quality as <cit.> under non-uniform parameterization. Its basis functions are refinable and the geometric mapping stays invariant during refinement. Moreover, we prove that a tuned hybrid subdivision surface is globally G^1-continuous. From the analysis point of view, tHNUS basis functions form a non-negative partition of unity, are globally linearly independent, and their spline spaces are nested. We numerically demonstrate that tHNUS basis functions can achieve optimal convergence rates for the Poisson's problem with non-uniform parameterization around extraordinary vertices.



Non-Uniform Subdivision Extraordinary Vertex Optimal Convergence Rates Isogeometric Analysis








§ INTRODUCTION


Isogeometric analysis (IGA) has emerged as a powerful technology to unify geometric modeling and numerical
simulation <cit.>, which employs the same basis functions used in computer-aided
design (CAD) and simulations. IGA has grown into a large family of numerical methods incorporating various spline techniques, such as NURBS (Non-Uniform Rational B-Splines) <cit.>, hierarchical B-splines <cit.>, T-splines <cit.>, polynomial splines over T-meshes <cit.>, and locally refinable B-splines <cit.>.

The study of extraordinary vertices[An interior vertex in a quadrilateral mesh is called an extraordinary vertex if it is shared by other than four faces.] has been one of the most active research directions in IGA because they are inevitable in complex watertight geometric representations. Along this direction, simultaneously fulfilling the requirements from both design and analysis is a significant challenge. Numerous methods have been developed over the past few years, but among them, only a few constructions can achieve optimal convergence rates in IGA, such as geometrically smooth multi-patch construction <cit.>, degenerated Bézier construction <cit.>, manifold-based construction <cit.>, and blended C^0 construction for unstructured hexahedral meshes <cit.>. A common simplification in all these constructions is to adopt uniform parameterization around extraordinary vertices, i.e., the surrounding knot intervals are assumed to be the same. While the support of non-uniform parameterization is a necessary step forward to be compatible with the current industry standard in CAD, i.e., NURBS, the related study on the above-mentioned constructions has not been reported in the literature.

On the other hand, subdivision methods, as a generalization of splines, provide a flexible means to deal with extraordinary vertices, where an infinite series of spline patches are smoothly joined around extraordinary vertices. The combination of flexibility and global smoothness makes them not only the standard in the computer animation industry but also a promising candidate for IGA. Indeed, some of the subdivision methods have been studied in the context of IGA, such as the use of Loop subdivision in thin-shell analysis <cit.> and the development of Catmull-Clark solids <cit.>. However, several challenging problems need to be carefully investigated before we can fully leverage the power of subdivision methods, such as developing efficient quadrature rules to integrate infinite piecewise polynomials around extraordinary vertices <cit.>, supporting non-uniform parameterizations to be compatible with NURBS <cit.>, and recovering optimal convergence rates <cit.>. This paper intends to address both non-uniform parameterization and optimal convergence behavior at the same time.

The present work is a follow-up of our preceding work on hybrid non-uniform subdivision (HNUS) <cit.>, which generalizes bicubic NURBS to arbitrary topology with proved G^1 continuity around extraordinary vertices. HNUS features high quality in geometric modeling under non-uniform parameterization. When applied to IGA, HNUS basis functions are not optimal but lead to improved convergence rates compared to Catmull-Clark subdivision.

Motivated by the idea of tuned Catmull-Clark subdivision <cit.> under uniform parameterization, we introduce a parameter λ∈ (1/4,1) in HNUS to control the shrinkage rate in irregular regions such that we can recover optimal convergence. The enhanced version of HNUS is therefore called tuned hybrid non-uniform subdivision (tHNUS). In fact, the parameter λ is the subdominant eigenvalue (the 2nd and 3rd eigenvalues which are equal) of the tHNUS subdivision matrix, that plays a crucial role in surface continuity <cit.> as well as the convergence performance <cit.>. Note that tHNUS coincides with the original HNUS when λ=1/2. From the geometric point of view, tHNUS retains comparable shape quality as HNUS. Its basis functions are refinable and the geometric mapping stays invariant during refinement. Moreover, we prove that the tHNUS surface is globally G^1-continuous. From the analysis point of view, tHNUS basis functions form a non-negative partition of unity, are globally linearly independent, and their spline spaces are nested. Moreover, we numerically demonstrate that tHNUS can achieve optimal convergence rates in the Poisson's problem by reducing λ, regardless of whether parameterization around extraordinary vertices is uniform or not. As an interesting side product, we also show that simply applying the standard Gauss quadrature rule to every element (close to or far away from extraordinary vertices) in tHNUS does not influence simulation accuracy or convergence.

The reminder of the paper is organized as follows. Section <ref> presents the subdivision rules of tHNUS. The proof of G^1 continuity for tHNUS surfaces is given in Section <ref>. The tHNUS basis functions are derived and their properties are discussed in Section <ref>. In Section <ref>, we present numerical tests of both geometric modeling and IGA. Section <ref> concludes the paper and discusses the future work.



§ TUNED HYBRID NON-UNIFORM SUBDIVISION SURFACES



Our discussion assumes that the input control mesh is a regular manifold mesh where all the faces are quadrilaterals. If initially a mesh has polygonal faces, we apply a single NURSS (Non-Uniform Recursive Subdivision Surface) refinement <cit.> to obtain an all-quadrilateral mesh. A non-negative scalar, which is called the knot interval, is assigned to each edge of the control mesh. We further assume that in each face, the knot intervals on the opposite edges coincide. A non-uniform parameterization is obtained by assigning different knot intervals to different edges as long as the assumption for knot intervals holds.

The tHNUS consists of two sets of rules: the topological rules to manipulate mesh connectivity, and the geometric rules to update the coordinates of involved control points. Each set of rules can be further divided into the first level and the subsequent levels. All the rules of tHNUS coincide with those of the original HNUS <cit.> except for the geometric rule corresponding to the subsequent  levels. We will concisely cover all the rules in the following to keep the explanation self-contained. One may refer to <cit.> for more details.

We start with the topological rules of tHNUS, which consist of rules for the first level and the subsequent levels, as illustrated in Figure <ref>. The rule corresponding to the first level converts the input quadrilateral mesh to its hybrid counterpart. Each extraordinary vertex is replaced by a polygonal face, whereas each spoke edge[A spoke edge is an edge touching a certain extraordinary vertex.] is replaced by a quadrilateral face. To make the resulting mesh conforming, additional vertices and edges are further replaced by certain faces; see Figure <ref>(a). Note that all the edges of a newly added polygonal face have a zero knot interval. Under the assumption of knot intervals, this means that all the newly added faces have a zero (parametric) measure. In regular regions, introducing zero-knot-interval edges leads to a reduction in continuity of basis functions from C^2 to C^1.



The topological rule for the subsequent levels is about how to split the initial hybrid mesh as in Figure <ref>(a). Overall, every edge with a nonzero knot interval is split equally into two. As a result, a quadrilateral face is split into four or two subfaces, depending on the number of nonzero-knot-interval edges it has. All the polygonal faces[We refer to non-quadrilaterals as polygonal faces.] stay unchanged (topologically).

We next introduce the geometric rules of tHNUS, which again are divided into the first level and the subsequent levels. At the first level, the rule to update regular vertices is the same as NURBS refinement, whereas the rule to compute polygon vertices is derived such that the limit surface of tHNUS has the same limit point and tangent plane as that of the non-uniform subdivision via eigen-polyhedron <cit.>. We take the eigen-polyhedron-based subdivision as the reference because it shows demonstrated shape quality under non-uniform parameterization. However, the computation is rather complicated and there is no explicit formula available. Alternatively, a simple explicit rule was provided in <cit.>, where each polygon vertex is computed as a convex combination of neighboring vertices. However, this explicit rule does not guarantee shape quality. Note that the geometric rule for the first level plays a crucial role in determining shape quality, but it has nothing to do with the proof of surface G^1 continuity or the convergence performance in IGA.



Now we provide the geometric rule of tHNUS for the subsequent levels, which differs from HNUS in that there is an additional tuning parameter λ in the formula to update polygon vertices. Referring to Figure <ref> and given knot intervals a_i, d_i, the points P_i^0, 0, P_i^1, 0, P_i^1, 1 and P_i^0, 1 in the refined mesh are defined as


    {P_i^0, 0   = (1 - λ)C + λ P_i^0, 0 + 2λα_i( -n P_i^0, 0 + ∑_j = 0^n-1( 1 + 2cos(2(j - i)π/n))P_j^0, 0), 
    P_i^1, 1   = d_id_i+1P_i^1, 1 + d_i(d_i+1+2a_i+1)P_i^1, 0 +
    d_i+1(d_i+2a_i)P_i^0, 1 + (d_i+2a_i)(d_i+1+2a_i+1)P_i^0, 0/4(d_i+a_i)(d_i+1+a_i+1), 
    P_i^1, 0   = d_id_i+1P_i-1^0, 1 + d_i(d_i+1+2d_i-1)P_i^1, 0 +
    d_i+1(d_i+2a_i)P_i-1^0, 0 + (d_i+2a_i)(d_i+1+2d_i-1)P_i^0, 0/4(d_i-1+d_i+1)(d_i+a_i), 
    P_i-1^0, 1   = d_id_i-1P_i^1, 0 + d_i(2d_i+1+d_i-1)P_i-1^0, 1 +
    d_i-1(d_i+2a_i)P_i^0, 0 + (d_i+2a_i)(2d_i+1+d_i-1)P_i-1^0, 0/4(d_i-1+d_i+1)(d_i+a_i),
    .


where α_j = 1/nd_j-1d_j+2/(d_j-1+d_j+1)(d_j + d_j+2), then

    C = ∑_i=0^n-1(d_iP_i+1^0, 0 + d_i+2P_i^0, 0)(d_i-1+d_i+3)/∑_j=0^n-1(d_j+d_j+2)(d_j-1+d_j+3)≐∑_i=0^n-1β_i P_i^0, 0.

The remaining points are computed by the NURBS mid-knot insertion. For example,

    P_i^2, 0   = a_iP_i^1, 0/2(d_i + a_i) + 1/4(d_i+1 + 2d_i-1)P_i^1, 0 + d_i+1P_i-1^0, 1/d_i+1+d_i-1 + d_iP_i^3, 0/2(d_i + a_i), 
    P_i^2, 1   = a_iP_i^1, 1/2(d_i + a_i) + 1/4d_i+1 P_i^1, 1 + (d_i+1 + 2a_i+1)P_i^1, 0/d_i+1 + a_i+1 + d_iP_i^3, 1/2(d_i + a_i), 
    P_i^2, 2   =  P_i^1, 1/4 + a_ia_i+1P_i^1, 1 + d_ia_i+1P_i^3, 1 + a_id_i+1P_i^1, 3 + d_id_i+1P_i^3, 3/4(d_i+a_i)(d_i+1+a_i+1)
       + a_i+1/4(d_i+1+a_i+1)M_1 + d_i+1/4(d_i+1+a_i+1)M_3 + d_i/4(d_i+a_i)M_2 + a_i/4(d_i+a_i)M_4,

where

    M_1   = d_i+1 P_i^1, 1 + (d_i+1 + 2a_i+1)P_i^1, 0/2(d_i+1 + a_i+1),
       M_2   = P_i^1, 1 + P_i^2, 1/2, 
    
    M_4   = P_i^1, 1 + (d_i + 2a_i)P_i^0, 1/2(d_i + a_i),
       M_3   = P_i^1, 1 + P_i^1, 2/2.



All the computations of these points are the same as those in HNUS except for polygon vertices P_i^0, 0, where the tuning parameter λ is introduced to control the size of shrinkage in the updated polygon: the smaller λ is, the more the polygon shrinks. As a result, isoparametric lines become more concentrated around extraordinary vertices. We will see examples in Section <ref>. When λ=1/2, the computation of P_i^0, 0 coincides with that in HNUS, and thus tHNUS is equivalent to HNUS in this particular case. However, the generalization via introducing λ is not as straightforward as it appears. The key insight is that λ turns out to be the subdominant eigenvalues (i.e., the 2nd and 3rd eigenvalues) of the subdivision matrix in tHNUS. As has been reported in <cit.>, convergence behavior in a subdivision scheme is mostly influenced by the subdominant eigenvalues. Therefore, tuning λ is equivalent to “controlling" convergence. We will have more detailed discussion about how λ improves convergence with specific examples in Section <ref>.


Tuned subdivision is a well studied subject aiming to optimized subdivision stencils (i.e., coefficients in the subdivision matrix) to improve certain properties of a subdivision scheme, for example, to minimize curvature variations to achieve a better surface fairness <cit.>. Recently, it has been explored in the context of IGA to improve accuracy <cit.> as well as convergence <cit.>. In particular, the tuned Catmull-Clark subdivision <cit.> is the first work in IGA that is able to use a subdivision scheme to achieve optimal convergence rates (in the L^2-norm error by solving the Poisson's equation). However, the optimization framework proposed in <cit.> only works for uniform parameterization and cannot be extended to non-uniform subdivision schemes because the subdivision stencils in a uniform subdivision scheme like Catmull-Clark only depend on the valence of a given extraordinary vertex, and the optimization can be focused on a finite number of stencils of interest. Thus, optimization only needs to be done once and the optimized stencils can be stored for future use. On the other hand, the subdivision stencils in a non-uniform subdivision depend on not only the valence of an extraordinary vertex, but also the surrounding knot intervals, leading to infinite possible cases of stencils. Therefore, it is not feasible to apply optimization to non-uniform subdivision because otherwise it would be very time-consuming and also problem specific.


We introduce λ explicitly to the formula of P_i^0, 0. Note that λ is a single parameter for all situations. It is independent of the valence of extraordinary vertices and the choice of knot intervals. However, this is only possible when bounded curvature is not of primary interest. tHNUS generally does not have bounded curvature. Nonetheless, bounded curvature under non-uniform parameterization remains an open problem and may not be available at all.




§ PROOF OF CONTINUITY



In order to prove tHNUS surfaces to be G^1 continuous, we need to prove that the spectrum of the subdivision matrix satisfies certain constraints and the associated characteristic map is regular and injective. Note that introducing λ to HNUS indeed significantly complicates the proof of G^1 continuity. Referring to Figure <ref> for the notations, the subdivision rule can be written into the following equations since the neighbor knot intervals a_i equals to d_i with enough subdivision levels.




    {P_j^0, 0   = (1 - λ)C + λ P_j^0, 0 + 2λα_j[ -n P_j^0, 0 + ∑_i = 0^n-1 ( 1 + 2cos(2(j-i)π/n))P_i^0, 0], 
    P_j^1, 0   = 3(2d_j-1+d_j+1)/8(d_j-1+d_j+1)P_j^0, 0 + 3d_j+1/8(d_j-1+d_j+1)P_j-1^0, 0 +
    (2d_j-1+ d_j+1)/8(d_j-1+d_j+1)P_j^1, 0 + d_j+1/8(d_j-1+d_j+1)P_j-1^0, 1 , 
    P_j-1^0, 1   = 3d_j-1/8(d_j-1+d_j+1)P_j^0, 0 + 3(d_j-1+2d_j+1)/8(d_j-1+d_j+1)P_j-1^0, 0 +
    d_j-1/8(d_j-1+d_j+1)P_j^1, 0 + (d_j-1+2d_j+1)/8(d_j-1+d_j+1)P_j-1^0, 1 , 
    P_j^1, 1   = 9/16P_j^0, 0 + 3/16P_j^0, 0 + 3/16P_j^0, 0 + 1/16P_j^1, 1.
    .


We arrange them in a matrix form M = S_n M, i.e.,

    (
    [   P_0^0, 0;          ⋮; P_n-1^0, 0;   P_0^1, 0;          ⋮; P_n-1^0, 1;   P_0^1, 1;          ⋮; P_n-1^1, 1 ])
    = (
        [                                                      ;         Q_n                 0                 0      ;                                                      ;                     E_0     …     0                  ;           ⋇           ⋮     ⋱     0           0      ;                       0     0 E_n-1                  ;                                      1/16     0     0;           ⋇                 ⋇           ⋮     ⋱      ;                                         0     0  1/16;       ])
    (
      [   P_0^0, 0;          ⋮; P_n-1^0, 0;  P_0^01, 0;          ⋮; P_n-1^0, 1;   P_0^1, 1;          ⋮; P_n-1^1, 1 ]) .

Denote Q_n = (Q_i, j), where i, j ∈{0,1,…, n-1}, and then we have

    Q_i, j = {[ (1-λ)β_j + 2( 1 + 2cos(2(j-i)π/n))λα_i,                                   j ≠ i;               λ + (1-λ)β_i- 2(n-3)λα_i,                                   j = i ].

and

    E_j = (
                 [ 2d_j-1+ d_j+1/8(d_j-1+d_j+1)         d_j+1/8(d_j-1+d_j+1);         d_j-1/8(d_j-1+d_j+1)  d_j-1+2d_j+1/8(d_j-1+d_j+1);                              ]).





Given an extraordinary vertex of any valence and an arbitrary choice of positive knot intervals, the eigenvalues of Q_n satisfy

    λ_1 = 1 > λ_2 = λ_3 = λ > |λ_k|, k = 4, 5, …, n.


We use the discrete Fourier transform to compute the eigenvalues of Q_n. Let p_k and p_k (k = 0, …, n-1) be the Fourier vectors corresponding to P_j and P_j, respectively, i.e.,

    p_k = 1/n∑_j = 0^n-1P_j^0, 0ω^jk,    p_k = 1/n∑_j = 0^n-1P_j^0, 0ω^jk, 
       P_k^0, 0 = ∑_j = 0^n-1p_jω^jk,    P_k^0, 0 = ∑_j = 0^n-1p_jω^jk,

where ω = e^2π/n and ω = e^-2π/n. Now the subdivision rule can be formulated in terms of the Fourier vectors,

    ∑_k=0^n-1p_kω^jk = ∑_k=0^n-1(∑_j=0^n-1(1-λ)β_jω^jk)p_k
    + p_0 + λω^jp_1 + λω^j(n-1)p_n-1 + 2λ(1/2 - nα_j)∑_k=2^n-2p_kω^jk.

Using the inverse discrete Fourier transform, we obtain

    ([   p_0;   p_1;     ⋮; p_n-1 ]) = (
                [          1   (1-λ)β_1          ⋯ (1-λ)β_n-1;          0          λ          ⋇          0;          0          0      B_n-3          0;          0          0          ⋇          λ;            ])
    ([   p_0;   p_1;     ⋮; p_n-1 ]) ,

where

    B_n - 3 = λ I - 2λ(
                                 [         ∑_j=0^n-1α_j      ∑_j=0^n-1α_jω^j                    … ∑_j=0^n-1α_jω^(n-4)j; ∑_j=0^n-1α_jω^(n-1)j         ∑_j=0^n-1α_j                    ⋯ ∑_j=0^n-1α_jω^(n-5)j;                    ⋮                    ⋮                    ⋱                    ⋮;     ∑_j=0^n-1α_jω^4j     ∑_j=0^n-1α_jω^5j                    ⋯         ∑_j=0^n-1α_j;                      ]) =: λ I - 2λ G_n.

Similar to <cit.>, both G_n and I - G_n are positive definite. Denote λ_B, i the eigenvalues of B_n-3. As G_n is positive definite, 1/2I - 1/2λ B_n-3 (=G_n) is also a positive definite matrix, which means that λ_B, i < λ. On the other hand, I - G_n is a positive definite matrix as well because μ_k < 1, and equivalently, I - (1/2I - 1/2λB_n-3) is positive definite, which means that λ_B, i > -λ. Therefore, we complete the proof.




Given an extraordinary vertex of any valence and an arbitrary choice of positive knot intervals, if λ > 1/4, then the eigenvalues of S_n satisfy

    λ_1 = 1 > λ_2 = λ_3 = λ > |λ_k|, k = 4, 5, …, 4n.



The eigenvalues of S_n consist of those of Q_n, E_i and 1/16I_n, where I_n is an n× n identity matrix. As proved in Lemma <ref>, the first three eigenvalues of Q_n are 1, λ, λ, and the remaining ones are less than λ. 1/16I_n has n equal eigenvalues 1/16 (< λ). It is also straightforward to verify that the eigenvalues of the 2×2 matrix E_i are 1/4 (< λ) and 1/8 (< λ). Therefore, we conclude that the eigenvalues of S_n are

    λ_1 = 1 > λ_2 = λ_3 = λ > |λ_k|, k = 4, 5, …, 4n.



The next step is to compute the characteristic map and prove that it is regular and injective. We first prove the following lemma.


Let P_i = (cos(2iπ/n), sin(2iπ/n)) ∈ℝ^2 (i = 0, …, n-1), C_P = ∑_i=0^n-1β_iP_i,
and P be an n×2 vector containing all P_i, i.e., P = [P_0, P_1, …, P_n-1]^T. Then we have

    S_n(P - C_P) = λ(P - C_P).



Denote P = S_nP, and we can obtain


    P_j - C_P=   λ( P_j - C_P) + 2λα_j[ -n (cos(2jπ/n), sin(2jπ/n)) + ∑_i = 0^n-1( 1 + 2cos(2(j-i)π/n))(cos(2iπ/n), sin(2iπ/n))] 
    
    =   λ( P_j - C_P) + 2λα_j[ -n (cos(2jπ/n), sin(2jπ/n)) + ∑_i = 0^n-12cos(2(j-i)π/n)(cos(2iπ/n), sin(2iπ/n))] 
    
    =   λ( P_j - C_P) + 2λα_j[ -n (cos(2jπ/n), sin(2jπ/n)) + 
       ∑_i = 0^n-1(cos(2jπ/n) +
    cos(2(j-2i)π/n), sin(2jπ/n) - sin(2(j-2i)π/n))] 
    
    =   λ( P_j - C_P) .


Since the above equation holds for any 0 ≤ j ≤ n-1, we conclude

    S_n(P - C_P) = λ(P - C_P) .





The characteristic map of tHNUS is regular and injective for any valence extraordinary vertices and any positive knot intervals if λ∈ (1/4, 1).


To prove that the characteristic map is regular and injective, we need a 4×4 grid of control points. We first compute the coordinates of this control grid that is used to define the characteristic map. The key idea is based on the fact that applying subdivision to the control grid of a characteristic map is equivalent to scaling the control grid by λ.

Referring to Figure <ref>, we have control points P_i^j, k, where 0 ≤ j, k ≤ 3 (0 ≤ i ≤ n-1). According to Lemma <ref>, if
we let P_i^0, 0 = (cos(2iπ/n), sin(2iπ/n)) ∈ℝ^2, C = ∑_i=0^n-1β_iP_i^0, 0, then we have
S_n[P_0^0, 0-C, …, P_n-1^0, 0-C]^T = λ[P_0^0, 0-C, …, P_n-1^0, 0-C]^T.

Further let E_i = d_i/d_i+d_i+2P_i+1^0,0+d_i+2/d_i+d_i+2P_i^0,0, p = P_i^0,0 - C, v = E_i-1 - C and w = E_i - C.
By definition, we have

    1/4(d_i+1+2d_i-1/2d_i+1+2d_i-1(P_i^1,0 - P_i^0,0)+
    d_i+1/2d_i+1+2d_i-1(P_i-1^0,1 - P_i-1^0,0)) + 1/2(E_i-1-C)    = λ(P_i^1,0 - P_i^0,0), 
    1/4(d_i-1/2d_i+1+2d_i-1(P_i^1,0 - P_i^0,0)+
    d_i-1+2d_i+1/2d_i+1+2d_i-1(P_i-1^0,1 - P_i-1^0,0)) + 1/2(E_i-1-C)    = λ(P_i-1^0,1 - P_i-1^0,0).

Solving the linear systems, we obtain

    P_i^1, 0 - P_i^0, 0 = P_i-1^0, 1 - P_i-1^0, 0 = 4(1-λ)/4λ-1v + 4(1-2λ)/8λ-1(p - v) .

Similarly, we compute P_i^2, 0, P_i^3, 0, P_i-1^0, 3, P_i-1^0, 3 as follows,

    P_i^2, 0 - P_i^1, 0   = 18(1-λ)/(8λ-1)(4λ-1)v + 18(1-2λ)/(16λ-1)(8λ-1)(p-v) ,
    
    P_i^3, 0 - P_i^2, 0   = 6(1-λ)(1+λ)/(8λ-1)λ(4λ-1)v + 3(1-4λ^2)/(16λ-1)λ(8λ-1)(p-v) .

We can also compute the remaining control points P_i^j, k (1 ≤ j, k ≤ 3), whose coefficients are complex expressions in λ. The detailed expressions are given in the Appendix.

With all these control points, we can now extract the Bézier control points for patches P_1, P_2 and P_3; see Figure <ref>(b). For example, in the patch P_2, let B_2^j, k (j, k = 0, …, 3) be the 4×4 Bézier control points. We denote S_2^j, k = B_2^j+1, k - B_2^j, k and T_2^j, k = B_2^j, k+1 - B_2^j, k. All S_2^j, k and T_2^j, k can be written as linear combinations of p, v and w, where the coefficients are again complex expressions in λ; see Appendix. We further plot some of these coefficients as functions of λ∈ (1/4, 1); see Figures <ref> and <ref>. We observe that S_2^j, k are convex combinations of vectors p, v and -w, while T_2^j, k are convex combinations of p, -v and w. Moreover, C is a convex combination of the points E_i from Equation (<ref>), so the patch P_2 is regular and injective. As a result, all the control points P_i^j, k (0 ≤ j, k ≤ 3) lie in the region bounded by two rays CE_i-1 and CE_i, which means that any two different patches must not intersect with one another. Similar results can also be achieved for patches P_1 and P_3. Therefore, the characteristic map of tHNUS is regular and injective for any λ∈ (1/4, 1), any valence extraordinary vertices and any positive knot intervals.








Given an arbitrary 2-manifold control mesh with any choice of positive knot intervals and any λ∈ (1/4, 1), the corresponding tHNUS limit surface is globally G^1-continuous.


The theorem is a direct result of Lemma <ref>, Lemma <ref> and Lemma <ref>.





§ HYBRID SUBDIVISION BASIS FUNCTIONS



In this section, we introduce basis functions of hybrid non-uniform subdivision. The derivation of such subdivision functions essentially follows Stam's method for Catmull-Clark subdivision <cit.>. However, there are two major differences. First, Catmull-Clark basis functions are associated with the input quadrilateral control mesh, whereas tHNUS basis functions are associated with the hybrid control mesh; see Figure <ref>(a). Second, Catmull-Clark subdivision features uniform knot intervals everywhere, leading to a subdivision matrix that only depends on the valence of a particular extraordinary vertex. In contrast, tHNUS (or HNUS) supports general non-uniform knot intervals, so the subdivision matrix depends not only on the valence of the extraordinary vertex, but also on the surrounding knot intervals.



 §.§ Definition of basis functions



We now introduce how tHNUS basis functions are defined on a hybrid control mesh. We start with distinguishing different types of faces. Recall that there exists both quadrilateral and polygonal faces in the hybrid mesh, and each edge in a polygonal face is assigned with a zero knot interval by construction. The knot intervals of other edges inherit from the input quadrilateral mesh and are constrained by the assumption that opposite edges in a quadrilateral face have the same knot interval. Moreover, note that edges perpendicular to the boundary also have zero knot intervals to make use of open knot vectors. An example of the knot interval configuration is shown in Figure <ref>(b), where the hybrid mesh is obtained from the input mesh in Figure <ref>(a).



We identify faces of zero-measure and nonzero-measure according to their parametric areas, which are computed using knot intervals. Zero-measure faces are not used in geometric representation and have no contribution to analysis. Note that all the polygonal faces and boundary faces have a zero-measure. The nonzero-measure faces, on the other hand, are divided into regular and irregular faces. An irregular face is a nonzero-measure face that shares a vertex with a certain polygonal face; all the other nonzero-measure faces are regular; see Figure <ref>(c). The tHNUS basis functions defined on a regular element[We use face and element interchangeably, but “face" emphasizes mesh topology whereas “element" is IGA-oriented.] are simply B-spline basis functions. In what follows, we restrict our attention to those defined on an irregular element. For simplicity of explanation, we assume that there is only one polygonal face in the 1-ring neighborhood of an irregular element. The 1-ring neighborhood of a face is a collection of faces that share vertices with this face, and recursively, the n-ring (n≥ 2) neighborhood consists of faces in the (n-1)-ring neighborhood as well as the faces sharing vertices with the (n-1)-ring neighborhood.


In a hybrid control mesh, each interior vertex is shared by four faces (or edges) and thus it has a regular valence of four. However, it does not mean that mesh irregularities are removed by converting the input quadrilateral mesh to its hybrid counterpart. In fact, irregularities are now manifested in the polygonal faces, which will be detailed in the following.






Given an irregular element Ω, let N denote the number of vertices in its adjacent polygonal face which is equivalent to the valence of the corresponding extraordinary vertex in the input quadrilateral mesh. There are K:=N+12 basis functions defined on Ω, associated with a local mesh around the polygonal face; see Figure <ref>(a). We denote

    𝐁_0(u,v) = [B_0,1(u,v), B_0,2(u,v), …, B_0,K(u,v)]^T

and

    𝐏_0 = [P_0,1, P_0,2, …, P_0,K]^T

the basis functions and the corresponding control vertices, respectively. Their indices are ordered according to Figure <ref>(a). The surface patch, i.e., the geometric mapping restricted to Ω is then represented by

    s(u,v) = 𝐏_0^T𝐁_0(u,v),    (u,v)∈Ω.

Note that Ω naturally has a parametric domain [0,d_1]× [0,d_2] that is determined by the corresponding knot intervals. We rescale it to Ω=[0,1]^2 to unify the treatment of irregular elements. The influence of the rescaling will be discussed in Remark <ref>.





Our focus is to derive 𝐁_0, which relies on subdivision of the corresponding control mesh 𝐏_0. Applying subdivision once yields Level-1 control vertices, denoted by

    𝐏_1    = [P_1,1, P_1,2, …, P_1,K]^T = 𝐒_1 𝐏_0, 
    𝐏̅_1    = [P_1,1, P_1,2,  …, P_1,K, P_1,K+1,…, P_1,M]^T = 𝐒̅_1 𝐏_0,

where M:=K+9=N+21. The subdivision matrices 𝐒_1 and 𝐒̅_1 have the dimension of K× K and M× K, respectively. Clearly, 𝐒̅_1 yields additional 9 vertices compared to 𝐒_1. The entries of 𝐒_1 and 𝐒̅_1 come from the tHNUS geometric rules as well as the mid-knot insertion of B-splines, see Equations (<ref>, <ref>) and (<ref>, <ref>), respectively. Among the four subelements at Level 1, three of them Ω_k^1 (k=1,2,3) are regular and correspond to regular C^1 B-spline patches. In other words, the surface patch restricted to Ω_k^1 is given by

    s(u,v) = 𝐏_1,k^T𝐍_1,k(u,v),    (u,v)∈Ω_k^1 ⊂Ω,

where 𝐏_1,k is a subvector of 𝐏̅_1 and 𝐍_1,k(u,v) is the vector of B-splines defined on Ω_k^1 (k=1,2,3); see Figure <ref>(c–e). Both 𝐏_1,k and 𝐍_1,k have a dimension of 16 due to the bicubic degree setting. 𝐏_1,k can be obtained with the help of a permutation matrix 𝐓_k, i.e., 𝐏_1,k=𝐓_k 𝐏̅_1 = 𝐓_k 𝐒̅_1 𝐏_0. Therefore, Equation (<ref>) becomes

    s(u,v) = (𝐓_k 𝐒̅_1 𝐏_0)^T 𝐍_1,k(u,v) = 𝐏_0^T (𝐓_k 𝐒̅_1)^T𝐍_1,k(u,v),    (u,v)∈Ω_k^1.

For Equations (<ref>, <ref>) to be equivalent under arbitrary choice of 𝐏^0, we need

    𝐁_0(u,v) = (𝐓_k 𝐒̅_1)^T𝐍_1,k(u,v),    (u,v)∈Ω_k^1 .

In other words, we have found the definition of 𝐁_0 on three quarters (Ω_1^1, Ω_2^1 and Ω_3^1) of Ω.



Now we are left to find the definition of 𝐁_0 on the remaining quarter [0,1/2]^2, and we proceed with the same idea explained above. As a result, the domain Ω is partitioned into an infinite series of tiles,

    Ω = ⋃_n=1^∞⋃_k=1^3 Ω_k^n,

where

    Ω_1^n    = [1/2^n, 1/2^n-1] ×[0,1/2^n], 
    Ω_2^n    = [1/2^n, 1/2^n-1] ×[1/2^n,1/2^n-1], 
    Ω_3^n    = [0,1/2^n] ×[1/2^n,1/2^n-1].

Analogous to Equation (<ref>), we have

    𝐏_n    =𝐒_n 𝐏_n-1 = 𝐒_n 𝐒_n-1⋯𝐒_1 𝐏_0, 
    𝐏̅_n    =𝐒̅_n 𝐏_n-1 = 𝐒̅_n 𝐒_n-1⋯𝐒_1 𝐏_0.

Again, 𝐒_n and 𝐒̅_n (n≥ 1) have the dimension of K× K and M× K, respectively. Note that 𝐒_2=𝐒_3=⋯ =𝐒_n (n≥ 2) and 𝐒̅_3 = 𝐒̅_4 = ⋯ = 𝐒̅_n (n≥ 3) because the ratios of knot intervals around an irregular subelement become fixed as the subdivision level increases. Therefore, following the same argument in deriving Equation (<ref>), we have the general expressions for 𝐁_0,

    𝐁_0(u,v) =
    
    (𝐓_k 𝐒̅_1)^T𝐍_1,k(u,v)     n=1 
    ( 𝐓_k 𝐒̅_2𝐒_1 )^T𝐍_2,k(u,v)     n=2 
    ( 𝐓_k 𝐒̅_3(𝐒_2)^n-2𝐒_1 )^T𝐍_n,k(u,v)     n≥ 3
    ,

where (u,v)∈Ω_k^n, and 𝐍_n,k(u,v) (n≥ 1) is the vector of B-splines defined on Ω_k^n.


Rescaling Ω to [0,1]^2 only affects the knot vectors (or equivalently, the vectors of knot intervals) of B-splines 𝐍_n,k(u,v) (n≥ 1). For example, when n=1 without scaling, the vector of knot intervals in the u direction is

    {d_3/2,d_3/2,0,d_1/2,d_1/2,a_1/2,a_1/2,b_1/2},

whereas after scaling with respect to d_1, it becomes

    U_1 = {d_3/2 d_1,d_3/2 d_1,0,1/2,1/2,a_1/2 d_1,a_1/2 d_1,b_1/2 d_1},

and similarly, we have the vector of knot intervals in the v direction,

    V_1 = {d_N/2 d_2,d_N/2 d_2,0,1/2,1/2,a_2/2 d_2,a_2/2 d_2,b_2/2 d_2}.

Moreover, when n=2, we have

    U_2 = {d_3/2^2 d_1,d_3/2^2 d_1,0,1/2^2,1/2^2,1/2^2,1/2^2,a_1/2^2 d_1},V_2  = {d_N/2^2 d_2,d_N/2^2 d_2,0,1/2^2,1/2^2,1/2^2,1/2^2,a_2/2^2 d_2},

and when n≥ 3,

    U_n = {d_3/2^n d_1,d_3/2^n d_1,0,1/2^n,1/2^n,1/2^n,1/2^n,1/2^n}, V_n  = {d_N/2^n d_2,d_N/2^n d_2,0,1/2^n,1/2^n,1/2^n,1/2^n,1/2^n}.

𝐍_n,k(u,v) are defined using U_n and V_n.



In fact, it is practically useful to rescale each tile Ω_k^n to [0,1]^2 by

    [           k=1  ξ = 2^n u-1,    η = 2^n v,;           k=2  ξ = 2^n u-1, η = 2^n v -1,;           k=3    ξ = 2^n u, η = 2^n v -1.;               ]

Correspondingly, the vectors of knot intervals are rescaled to

    [                                                       n=1 Ξ_1 = {d_3/d_1, d_3/d_1, 0,1,1, a_1/d_1,a_1/d_1,b_1/d_1}, Θ_1 = {d_N/d_2, d_N/d_2, 0,1,1, a_2/d_2,a_2/d_2,b_2/d_2},;                                                       n=2            Ξ_2 = {d_3/d_1, d_3/d_1, 0,1,1, 1, 1,a_1/d_1},             Θ_2 = {d_N/d_2, d_N/d_2, 0,1,1, 1,1,a_2/d_2},;                                                      n≥ 3                Ξ_3 = {d_3/d_1, d_3/d_1, 0,1,1, 1, 1,1  },                  Θ_3 = {d_N/d_2, d_N/d_2, 0,1,1, 1,1,1 },;                                                           ]

where the rescaled knot intervals are now independent of the subdivision level n when n ≥ 3. In summary, the basis functions of interest are defined as

    𝐁_0(u,v) =
    
    (𝐓_k 𝐒̅_1)^T𝐛_1,k(ξ(u),η(v))     n=1 
    ( 𝐓_k 𝐒̅_2𝐒_1 )^T𝐛_2,k(ξ(u),η(v))     n=2 
    ( 𝐓_k 𝐒̅_3(𝐒_2)^n-2𝐒_1 )^T𝐛_3,k(ξ(u),η(v))     n≥ 3
    ,

where (ξ,η)∈ [0,1]^2 and 𝐛_l,k(ξ,η) are B-splines defined using Ξ_l and Θ_l (l= 1,2,3). Note that when n=1,2, we have 𝐍_n,k(u,v) = 𝐛_n,k(ξ(u),η(v)), and when n≥ 3, 𝐍_n,k(u,v) = 𝐛_3,k(ξ(u),η(v)).








Evaluation of 𝐁_0(u,v) at (0,0) in an irregular element needs the computation of lim_n→∞ (𝐒_2)^n-2. Following <cit.>, we need to eigen-decompose 𝐒_2 such that 𝐒_2 = 𝐕Λ𝐕^-1, where Λ is a diagonal matrix containing the eigenvalues of 𝐒_2 and 𝐕 is an invertible matrix with columns being the corresponding eigenvectors. Accordingly, we have lim_n→∞ (𝐒_2)^n-2 =lim_n→∞𝐕Λ^n-2𝐕^-1. Recall that all the eigenvalues are smaller than 1 (and greater than 0) except the first one λ_1=1, so lim_n→∞Λ^n-2 is a matrix whose entries are all zero except the first-row-first-column entry, which is one. On the other hand, the derivatives of 𝐁_0(u,v) are not bounded around (0,0). We can see this by applying the chain rule, for example,

    𝐁_0(u,v)u = 2^n ( 𝐓_k 𝐒̅_3(𝐒_2)^n-2𝐒_1 )^T𝐛_3,k(ξ,η)ξ,

where the factor comes from dξ/d u = 2^n. A differentiable version of 𝐁_0 (with respect to certain parameters) can be obtained via characteristic-map-based reparameterization <cit.>. However, in this paper, we are interested in applying tHNUS basis functions in the context of IGA, so we only need derivatives at quadrature points that are away from (0,0). Moreover, what we eventually need is derivatives with respect to the physical coordinates, for example,

    𝐁_0 ∘ s^-1(x,y)x = 𝐁_0uux + 𝐁_0vvx,

where s^-1 is the inverse mapping of s(u,v). The troublesome factor 2^n, which may cause overflow when n becomes too large, is canceled out with that from ∂ u/∂ x (which is 2^-n) and does not cause any numerical issues. The same argument applies to higher order derivatives.



In <cit.>, the eigen structure (Λ,𝐕) is precomputed for different valence numbers and stored in a file for repeated use. However, the same scheme cannot be applied to tHNUS because the subdivision matrix 𝐒_2 depends on not only the valence number but also the surrounding knot intervals, leading to infinite possible cases of 𝐒_2. Therefore, the eigen structure of 𝐒_2 needs to be found in real time for every irregular element. Alternatively, we can directly perform matrix multiplications to compute (𝐒_2)^n-2, especially when the valence number is small and basis functions need to be computed at points other than (0,0). This is indeed the case in IGA where evaluation is needed at quadrature points. In practice, we adopt a near-machine-precision tolerance (e.g., 10^-13) to prevent a potential overflow issue.












In the previous discussion, 𝐁_0 is derived under the assumption that there is only one polygonal face next to an irregular element, which, however, is not a necessary condition. When an irregular element has multiple adjacent polygonal faces, we treat it as a macro element and pseudo-subdivide it once. Each of the resulting four subelements only has one polygonal face, where basis functions are defined according to our previous discussion. In other words, basis functions are well defined on each quarter of the original macro element. This extension follows the same idea proposed in <cit.>, which extends Stam's derivation <cit.> to arbitrary unstructured quadrilateral meshes.








 §.§ Quadrature


To apply the standard Gauss quadrature rule, we need to guarantee that the involved basis functions are polynomials (rather than piecewise polynomials) on each integration cell. However, the functions in 𝐁_0 are piecewise smooth polynomials defined on an infinite series of subdomains, i.e., {Ω_k^n}_n=1^∞ (k=1,2,3). The straightforward way is to apply the Gauss quadrature rule on each cell Ω_k^n up to a certain fine level, which was adopted in several subdivision-based isogeometric methods <cit.>. We call such a quadrature the full quadrature scheme. In our patch test, we observe that the solution achieves machine precision (∼ 10^-16) when the 4-point rule is used and the level n is set to be 10. As a result, a total number of 496 quadrature points are needed for a single irregular element. In contrast, only 16 Gauss quadrature points are used for a regular element.

Alternatively, we can “brutally" apply the Gauss quadrature rule to the entire irregular element. In other words, only 16 Gauss quadrature points are placed on an irregular element. We observe that such a reduced quadrature scheme does not influence convergence. In fact, it does not introduce noticeable numerical error in terms of the L^2- or H^1-norm error compared to the full quadrature. We will numerically compare the two schemes in the next section.



 §.§ Properties


Now we briefly discuss several properties of tHNUS basis functions, including non-negative partition of unity, refinability (equivalent to nested spline spaces), and global linear independence. The non-negative partition of unity of tHNUS basis functions follows from the fact that all the entries in the subdivision matrix are non-negative and each row sum of the subdivision matrix is one. Refinability states that each basis function of a given mesh can be represented as a linear combination of those defined on a refined mesh. In fact, we can see this property in the derivation of 𝐁_0, where basis functions are always expressed as linear combinations of functions in the refined meshes.

Finally, the global linear independence implies linear independence on the entire domain, and it can be easily shown under the mild assumption that each irregular element has at least one regular element as its direct neighbor. Under this assumption, every basis function has support on a certain regular element, where it is simply a B-spline. As B-splines are linearly independent on such an element, we can conclude that all the basis functions are linearly independent on the entire domain by going through all the regular elements. The proof on general meshes becomes more involving because we need to resolve different configurations of polygonal faces, or equivalently, configurations of extraordinary vertices in the input mesh. A complex configuration usually occurs when the mesh is very coarse such that many extraordinary vertices may be next to one another. When this is the case, we can perform global refinement to guarantee linear independence.







§ NUMERICAL EXAMPLES



In this section, we present several numerical examples using tHNUS surfaces in both geometric modeling and IGA.





 §.§ Geometric modeling with tHNUS surfaces

We show some tHNUS limit surface examples and compare them with the existing non-uniform subdivision schemes.
We first show the graphs of blending functions for the extraordinary points (EPs) with different valences, such as valence-5 EP in
Figure <ref>, valence-6 EP in Figure <ref> and valence-7 EP in Figure <ref>.
As stated in <cit.>, the approaches in <cit.>, <cit.> and <cit.>
produce limit surfaces with very similar quality in all the examples. Therefore, we only show the limit surface
comparisons in one example as shown in Figure <ref>. All the rest of the examples only show the limit surface of
the new tHNUS with different λ. We can observe that all different λ can produce better shape quality
than those approaches in <cit.>, <cit.> and <cit.>, but the small λ
produces worse shape quality surround the EPs, see Figures <ref> and <ref> for the
details.












 §.§ IGA applications using tHNUS basis functions


In this section, we test the performance of tHNUS basis functions in the context of IGA. We solve the Poisson's equation with several unstructured quadrilateral meshes as the input. We start with convergence tests on a unit square, whose input control mesh has two EPs, one of valence 3 and the other of valence 5; see Figure <ref>(a). These tests are aimed at studying: (1) the role of λ in convergence, (2) the feasibility of using reduced quadrature, and (3) the influence of non-uniform parameterizations on convergence.



First, we study the influence of the tuning parameter λ on convergence behavior, where we choose λ to be 0.65, 0.5, and 0.26. Recall that tHNUS is equivalent to the original HNUS when λ=0.5. We adopt uniform parameterization (i.e., same knot intervals) around EPs as well as full quadrature in this study. With the manufactured solution u(x,y)=sin(π x) sin(π y), we summarize the convergence plots in Figure <ref>. We observe that a smaller λ delivers a better convergence behavior, and particularly, optimal convergence rates are achieved when λ=0.26. The tuned Catmull-Clark subdivision (with uniform parameterization) was studied in <cit.>, where optimal convergence rates in the L^2-norm were observed in the Poisson's problem when λ=0.39. It indicates that the tuning parameter in tHNUS plays a less sensitive role than that in <cit.> because tHNUS requires a smaller λ to recover optimal convergence. The reason may be that λ brings more vertices to move further towards each EP than in tHNUS. As a result, the tuned Catmull-Clark subdivision has a faster shrinkage in irregular regions. We will provide insights about why reducing λ recovers optimal convergence later when we study the meshes with high-valence EPs.







Second, we compare two quadrature schemes in irregular elements, full quadrature versus the reduced quadrature, under uniform parameterization with λ=0.26 and λ=0.5. We observe in Figure <ref>(a) that there is no noticeable difference in terms of both L^2- and H^1-norm errors. In other words, both quadrature schemes deliver the same level of accuracy when λ=0.26. In contrast, quadrature plays an important role when λ=0.5, where the full quadrature yields nearly one-order higher convergence rates than the reduced quadrature. This indicates that when λ=0.26, the corresponding basis functions (piecewise polynomials) in irregular elements can be better approximated by polynomials than those using λ=0.5, and 16 quadrature points seem to suffice to retain accuracy. However, further study is needed to fully understand the mechanism behind.







Third, we study several different non-uniform parameterizations for convergence test, which can be obtained by assigning different knot intervals to the edges in the input control mesh. Semi-uniform knot intervals are usually adopted in the literature, where all the edges are assigned a unit knot interval except for those perpendicular to the boundary, which are assigned a zero knot interval. To have non-uniform parameterization around EPs, we modify the semi-uniform setting in two ways: (1) the knot interval (denoted by d) of highlighted spoke edges takes values d∈{1,2,5,10}; and (2) every spoke edge is assigned a different knot interval; see Figure <ref>(a, b). In both cases, we observe in Figure <ref>(c, d) that tHNUS basis functions can achieve optimal convergence rates with λ=0.26. We also observe that the convergence plots corresponding to a larger d slightly shift up, meaning that larger difference in knot intervals yields larger approximation error. In other words, the “distortion" in parameterization influences accuracy rather than convergence.





Now, we consider meshes with high-valence EPs (valence 6, 7 and 8), where each spoke edge is assigned a different knot interval. We again observe optimal convergence rates with λ=0.26; see Figure <ref>. Moreover, let us have a close look at how λ influences parameterization around EPs. In particular, we compare isoparametric lines using two different λ's (0.5 versus 0.26) around a valence-6 EP, where the same input control mesh in Figure <ref>(a) is used in both cases. We find in Figure <ref> that isoparametric lines are overlaid with one another in most regions, and with a smaller λ, the isoparametric lines (blue curves) are more bent towards the extraordinary surface point s(0,0). Equivalently speaking, a smaller λ yields smaller refined irregular elements in the physical domain. Therefore, the mesh around s(0,0) becomes denser than that using a larger λ, and as a result, the asymptotic approximation error controlled by s(0,0) can be reduced using such a denser mesh. Ideally, optimal convergence rates can be achieved by reducing λ, which indeed is the case in all our numerical tests when λ=0.26.


Although the globally smooth tHNUS basis functions can be applied to solve 4th-order partial differential equations (PDEs), our preliminary tests only show suboptimal convergence in solving the biharmonic equation, where obtained convergence rates in terms of L^2-, H^1- and H^2-norm errors are around 2, 2, and 1, respectively. This is consistent with the result reported in <cit.>, where a thin-shell problem was solved and the reported convergence rates in L^2- and energy norm errors are 2 and 1, respectively. In other words, reducing λ alone is not sufficient for high-order PDEs. We conjecture that to recover the optimal convergence in this case, we may need more degrees of freedom around EPs following similar ideas in <cit.>. However, this would further complicate the current subdivision framework, so we postpone related results in a follow-up work.




§ CONCLUSIONS AND FUTURE WORK



We have presented a tuned version of hybrid non-uniform subdivision, tHNUS, by introducing a parameter λ∈ (1/4, 1), which is also the second and third eigenvalues of the subdivision matrix. The tHNUS surface is proved to be G^1-continuous for any positive knot intervals and extraordinary vertices of any valence. The tHNUS surface has satisfactory shape quality for any λ under non-uniform parameterization. However, the highest shape quality is achieved when λ=0.5. In other words, the original HNUS generally performs better in geometric modeling than tHNUS. On the other hand, tHNUS basis functions can achieve optimal convergence rates when λ is reduced to 0.26, regardless of which quadrature scheme is used and whether parameterizations are uniform or non-uniform around EPs.

In the future, we can extend tHNUS in the following three issues. First, converting an input quadrilateral mesh to its hybrid counterpart can be restricted locally to irregular regions without introducing zero-measure faces throughout the entire mesh. This can be done by allowing T-junctions <cit.> in the hybrid mesh, but support of T-junctions in a hybrid mesh requires a much more sophisticated data structure to accommodate both polygonal faces and quadrilaterals with T-junctions. Second, tHNUS can be adapted to hierarchical splines <cit.> due to its refinability property, where the initial level corresponds to the initial hybrid mesh. The construction of hierarchical tHNUS essentially follows those proposed for truncated hierarchical Catmull-Clark subdivision surfaces <cit.>, but the differences lie in dealing with hybrid meshes and non-uniform knot intervals. Third, improving tHNUS to achieve optimal convergence rates in solving high-order PDEs is another challenging but very interesting direction to pursue. Currently, we can only show optimal convergence in solving the 2nd-order PDE. In the case of high-order PDEs, our preliminary tests suggest that it is not sufficient to tune λ alone and additional treatment is needed. We plan to our investigation by adding more control points around extraordinary vertices.



§ ACKNOWLEDGEMENTS



X. Wei was supported in part by ERC AdG project CHANGE n. 694515. X. Li was supported by the NSF of China (No.61872328), NKBRPC (2011CB302400), SRF for ROCS SE, and the Youth Innovation Promotion Association CAS. Y. Zhang was supported in part by the NSF grant CBET-1804929. T.J.R. Hughes was partially supported by the Office of Naval Research, USA (Grant Nos. N00014-17-1-2119 and N00014-13-1-0500).




§ APPENDIX


The control points P_i^j,k (1≤ j,k≤ 3) of the characteristic map are listed as follows.


    P_i^1,1   = -(λ-2) ((32 λ^3-100 λ^2+31 λ-2) p-12 λ^2 (v+w))/(λ-4) λ (4 λ-1) (8 λ-1)
    
    P_i^1,2   = 6 λ^2 ((256 λ^3+136 λ^2-2251 λ-232) v+2 (128 λ^3-184 λ^2-149 λ+10) w)+(-4096 λ^6+18688 λ^5-14576 λ^4-22112 λ^3+8275 λ^2-608 λ+20) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1)
    
    P_i^1,3   = 6 λ^2 (2 (320 λ^3-919 λ^2-800 λ+7) v+(160 λ^3-314 λ^2-13 λ+2) w)+(-2560 λ^6+15904 λ^5-29348 λ^4+3535 λ^3+928 λ^2-62 λ+2) p/(λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1)
    
    P_i^2,1   = 6 λ^2 ((256 λ^3-368 λ^2-298 λ+20) v+(256 λ^3+136 λ^2-2251 λ-232) w)+(-4096 λ^6+18688 λ^5-14576 λ^4-22112 λ^3+8275 λ^2-608 λ+20) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1)
    
    P_i^2,2   = 6 (4096 λ^4+33024 λ^3-96320 λ^2-11271 λ-1160) λ^2 (v+w)+(-65536 λ^7-194560 λ^6+2362752 λ^5-4183824 λ^4+1165584 λ^3-71505 λ^2+976 λ+100) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    P_i^2,3   = 6 λ^2 ((74240 λ^4-139104 λ^3-145684 λ^2-8118 λ-7) v+(58880 λ^4-106224 λ^3-32914 λ^2-2208 λ-127) w)/2 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-942080 λ^7+4816384 λ^6-6465888 λ^5+394032 λ^4+345858 λ^3-33192 λ^2+1313 λ-1) p/2 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    P_i^3,1   = 6 λ^2 ((160 λ^3-314 λ^2-13 λ+2) v+2 (320 λ^3-919 λ^2-800 λ+7) w)+(-2560 λ^6+15904 λ^5-29348 λ^4+3535 λ^3+928 λ^2-62 λ+2) p/(λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1)
    
    P_i^3,2   = 6 λ^2 ((58880 λ^4-106224 λ^3-32914 λ^2-2208 λ-127) v+(74240 λ^4-139104 λ^3-145684 λ^2-8118 λ-7) w)/2 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-942080 λ^7+4816384 λ^6-6465888 λ^5+394032 λ^4+345858 λ^3-33192 λ^2+1313 λ-1) p/2 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    P_i^3,3   = 6 (25600 λ^4-30144 λ^3-61166 λ^2-5685 λ-236) λ^2 (v+w)+(-409600 λ^7+1923584 λ^6-1925280 λ^5-1010388 λ^4+470631 λ^3-36036 λ^2+1060 λ+16) p/(λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)



The expressions of all S_2^j, k (0≤ j,k≤ 3) are listed as follows.


    S_2^0,0   = 3 λ(28672 λ^4-44160 λ^3-104888 λ^2-4242 λ-59) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)+3 λ(-20480 λ^4+125568 λ^3+2296 λ^2-2886 λ-107) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-589824 λ^7+2557952 λ^6-1417856 λ^5-723088 λ^4+318940 λ^3-25264 λ^2+853 λ+3) p/4 (λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^1,0   =-6 λ(14336 λ^4-53888 λ^3+58744 λ^2+2009 λ+30) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-12 λ(5120 λ^4-11808 λ^3+1904 λ^2-234 λ+41) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (81920 λ^7-467968 λ^6+851904 λ^5-540792 λ^4+116238 λ^3-7683 λ^2+211 λ+1) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^2,0   =-(163840 λ^5-360448 λ^4+20512 λ^3+2732 λ^2+571 λ+47) w/4 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(458752 λ^6-1579008 λ^5+1185408 λ^4+775632 λ^3-25056 λ^2+567 λ-4) v/8 (λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(2621440 λ^8-14696448 λ^7+25499648 λ^6-13590144 λ^5+1352256 λ^4+248436 λ^3-23683 λ^2+897 λ-5) p/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^0,1   =18 λ(2048 λ^4-1856 λ^3-7028 λ^2-669 λ-13) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)+18 λ(-2048 λ^4+12288 λ^3+1904 λ^2-556 λ-17) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       -3 (32768 λ^7-129024 λ^6+47872 λ^5+108496 λ^4-38716 λ^3+3018 λ^2-95 λ-1) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^1,1   = -6 λ(12288 λ^4-38912 λ^3+43904 λ^2+4127 λ+76) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-12 λ(6144 λ^4-14080 λ^3+1036 λ^2-299 λ+80) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (98304 λ^7-553984 λ^6+985152 λ^5-637296 λ^4+138348 λ^3-9009 λ^2+216 λ+4) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^2,1   = -3 (32768 λ^5-93184 λ^4+70336 λ^3+61032 λ^2-210 λ+7) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-3 (32768 λ^5-70656 λ^4-4768 λ^3+888 λ^2+236 λ+15) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(524288 λ^8-2883584 λ^7+4827136 λ^6-2557248 λ^5+206544 λ^4+64080 λ^3-6056 λ^2+233 λ-1) p/4 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^0,2   =6 λ(4096 λ^4-1024 λ^3-19348 λ^2-1877 λ-117) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)+6 λ(-4096 λ^4+21888 λ^3+20552 λ^2-3786 λ-97) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-65536 λ^7+215040 λ^6+122880 λ^5-518496 λ^4+160368 λ^3-11706 λ^2+295 λ+9) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^1,2   =-6 λ(8192 λ^4-28928 λ^3+34888 λ^2+4205 λ+228) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-48 λ(1024 λ^4-2272 λ^3-868 λ^2-65 λ+39) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (65536 λ^7-393216 λ^6+768768 λ^5-550584 λ^4+120330 λ^3-6531 λ^2+32 λ+12) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^2,2   =-(32768 λ^5-100352 λ^4+78176 λ^3+82984 λ^2+1862 λ+7) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(32768 λ^5-63488 λ^4-49120 λ^3+2596 λ^2+845 λ+43) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(524288 λ^8-2998272 λ^7+5300224 λ^6-2834112 λ^5-22440 λ^4+158118 λ^3-13784 λ^2+513 λ-1) p/6 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^0,3   =(393216 λ^7+233472 λ^6-2491392 λ^5-285864 λ^4-21918 λ^3-258 λ^2) v/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)+(-393216 λ^7+1769472 λ^6+3983616 λ^5-415440 λ^4-81576 λ^3+1449 λ^2-12 λ) w/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-1048576 λ^8+2555904 λ^7+6412288 λ^6-14197248 λ^5+3735840 λ^4-144900 λ^3-5573 λ^2+717 λ-1) p/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^1,3   =-(65536 λ^5-262144 λ^4+334720 λ^3+52988 λ^2+3718 λ+41) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(65536 λ^5-139264 λ^4-137984 λ^3-23368 λ^2+8098 λ-7) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(1048576 λ^8-6782976 λ^7+14569472 λ^6-11210496 λ^5+2223408 λ^4-32748 λ^3-8632 λ^2+615 λ+1) p/12 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    S_2^2,3   =-(65536 λ^5-219136 λ^4+173056 λ^3+229220 λ^2+9745 λ+131) v/3 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(65536 λ^5-108544 λ^4-208832 λ^3-11944 λ^2+6694 λ+239) w/3 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(1048576 λ^8-6291456 λ^7+11761664 λ^6-6099072 λ^5-1021344 λ^4+664680 λ^3-54730 λ^2+1881 λ+7) p/18 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)



The expressions of all T_2^j, k (0≤ j,k≤ 3) are listed as follows.


    T_2^0,0   = -3 λ(20480 λ^4-125568 λ^3-2296 λ^2+2886 λ+107) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-3 λ(-28672 λ^4+44160 λ^3+104888 λ^2+4242 λ+59) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-589824 λ^7+2557952 λ^6-1417856 λ^5-723088 λ^4+318940 λ^3-25264 λ^2+853 λ+3) p/4 (λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^1,0   =-18 λ(2048 λ^4-12288 λ^3-1904 λ^2+556 λ+17) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-18 λ(-2048 λ^4+1856 λ^3+7028 λ^2+669 λ+13) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       -3 (32768 λ^7-129024 λ^6+47872 λ^5+108496 λ^4-38716 λ^3+3018 λ^2-95 λ-1) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^2,0   =-6 λ(4096 λ^4-21888 λ^3-20552 λ^2+3786 λ+97) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-6 λ(-4096 λ^4+1024 λ^3+19348 λ^2+1877 λ+117) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-65536 λ^7+215040 λ^6+122880 λ^5-518496 λ^4+160368 λ^3-11706 λ^2+295 λ+9) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^3, 0   =-(-65536 λ^5-38912 λ^4+415232 λ^3+47644 λ^2+3653 λ+43) w/4 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(131072 λ^6-589824 λ^5-1327872 λ^4+138480 λ^3+27192 λ^2-483 λ+4) v/8 (λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(-1048576 λ^8+2555904 λ^7+6412288 λ^6-14197248 λ^5+3735840 λ^4-144900 λ^3-5573 λ^2+717 λ-1) p/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^0, 1   = -12 λ(5120 λ^4-11808 λ^3+1904 λ^2-234 λ+41) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-6 λ(14336 λ^4-53888 λ^3+58744 λ^2+2009 λ+30) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (81920 λ^7-467968 λ^6+851904 λ^5-540792 λ^4+116238 λ^3-7683 λ^2+211 λ+1) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^1,1   =-12 λ(6144 λ^4-14080 λ^3+1036 λ^2-299 λ+80) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-6 λ(12288 λ^4-38912 λ^3+43904 λ^2+4127 λ+76) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (98304 λ^7-553984 λ^6+985152 λ^5-637296 λ^4+138348 λ^3-9009 λ^2+216 λ+4) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^2,1   =-48 λ(1024 λ^4-2272 λ^3-868 λ^2-65 λ+39) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-6 λ(8192 λ^4-28928 λ^3+34888 λ^2+4205 λ+228) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +2 (65536 λ^7-393216 λ^6+768768 λ^5-550584 λ^4+120330 λ^3-6531 λ^2+32 λ+12) p/(λ-4) λ (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^3,1   = -(65536 λ^5-139264 λ^4-137984 λ^3-23368 λ^2+8098 λ-7) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(65536 λ^5-262144 λ^4+334720 λ^3+52988 λ^2+3718 λ+41) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(1048576 λ^8-6782976 λ^7+14569472 λ^6-11210496 λ^5+2223408 λ^4-32748 λ^3-8632 λ^2+615 λ+1) p/12 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^0,2   = (-983040 λ^7+2162688 λ^6-123072 λ^5-16392 λ^4-3426 λ^3-282 λ^2) v/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)+(-1376256 λ^7+4737024 λ^6-3556224 λ^5-2326896 λ^4+75168 λ^3-1701 λ^2+12 λ) w/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(2621440 λ^8-14696448 λ^7+25499648 λ^6-13590144 λ^5+1352256 λ^4+248436 λ^3-23683 λ^2+897 λ-5) p/24 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^1,2   =-3 (32768 λ^5-70656 λ^4-4768 λ^3+888 λ^2+236 λ+15) v/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-3 (32768 λ^5-93184 λ^4+70336 λ^3+61032 λ^2-210 λ+7) w/2 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(524288 λ^8-2883584 λ^7+4827136 λ^6-2557248 λ^5+206544 λ^4+64080 λ^3-6056 λ^2+233 λ-1) p/4 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^2,2   = -(32768 λ^5-63488 λ^4-49120 λ^3+2596 λ^2+845 λ+43) v/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(32768 λ^5-100352 λ^4+78176 λ^3+82984 λ^2+1862 λ+7) w/(λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(524288 λ^8-2998272 λ^7+5300224 λ^6-2834112 λ^5-22440 λ^4+158118 λ^3-13784 λ^2+513 λ-1) p/6 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
    
    T_2^3,2   =-(65536 λ^5-108544 λ^4-208832 λ^3-11944 λ^2+6694 λ+239) v/3 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)-(65536 λ^5-219136 λ^4+173056 λ^3+229220 λ^2+9745 λ+131) w/3 (λ-4) (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)
       +(1048576 λ^8-6291456 λ^7+11761664 λ^6-6099072 λ^5-1021344 λ^4+664680 λ^3-54730 λ^2+1881 λ+7) p/18 (λ-4) λ^2 (4 λ-1) (8 λ-1) (16 λ-1) (32 λ-1) (64 λ-1)



The following mathematica code has been used to compute the above control points.
mathematical/hybrid_subd.m




plain

