The Weighted Basis for PHT-Splines

1  Introduction

Polynomial splines over hierarchical T-meshes (PHT-splines for short) [1] are a kind of global C1 continuous polynomial splines that can be locally refined. PHT-splines possess a very efficient local refinement algorithm and inherit many good properties of T-splines such as adaptivity and locality. PHT-splines are now widely applied in adaptive geometric modeling [2–4], adaptive finite element [5], iso-geometric analysis [6–13] and topology optimization [14].

With further in-depth application, the improvement of PHT-splines has also been underway. One of the limitations of PHT-splines is the restrictions from hierarchical T-meshes, which require the underlying meshes to start with tensor product meshes and be refined by inserting crosses. The generalization of PHT-splines on general T-meshes [15] and modified hierarchical T-meshes (allowing split-in-half in mesh refinement) [16] is precisely aimed at overcoming the restrictions from meshes. Another improvement regard to PHT-splines is the construction of basis functions. The original basis functions [1] have a decay phenomenon (the function values decay as the level increases) when the level difference is big in the underlying support. The decay of basis functions when applied in iso-geometric analysis leads to ill-conditioned stiffness matrices, and thus unstable numerical solutions. A new set of basis functions is defined by the B-splines that defined on the 2 × 2 tensor product meshes associated with basis vertices [2] to overcome the decay phenomenon that occurs in basis functions, but at the cost of losing the partition of unity.

The partition of unity is a default ingredient in the construction of approximating functions in finite element methods [17–19]. Moreover, in the field of topology optimization [14], the density variable takes value in the range of [0,1]. The basis functions that form a partition of unity are more suitable for designing density variables that meet the range condition. In this paper, we aim to provide basis functions that form a partition of unity and also have no decay phenomenon, facilitating more effective applications for PHT-splines.

The truncation mechanism is a common way of constructing basis functions for hierarchical spline spaces that form a partition of unity, such as the original PHT-splines [1], THB-splines (the truncated hierarchical B-splines) [20] and truncated T-splines [21]. But as mentioned above, the truncation mechanism causes the decay in basis functions. In this paper, we will show that the basis functions that form the partition of unity and have no decay can be constructed without the truncation mechanism. The weighted bases that satisfy the required properties are defined by the bases constructed in [2] multiplied by positive weights. The weight for each basis function can be computed explicitly and efficiently based on the geometric information of the basis function. Moreover, when the level difference is not greater than 2, the weights can also be computed by the PHT subdivision scheme.

In [22], the authors also proposed a method to define new basis functions to overcome the decay problem and guarantee the partition of unity. However, their method only considers specific cases that the original truncated basis functions have rectangular supports, in which case they are replaced by the associated tensor product B-splines. Instead, our method is more versatile and can be applied to arbitrary T-meshes easily.

The rest of this paper is organized as follows. In Section 2, some basics of PHT-splines are reviewed. In Section 3, we introduce the weighted bases in detail. In Section 4, the subdivision scheme of PHT-splines is introduced as an alternative algorithm for calculating weights. In Section 5, some hierarchical T-meshes are demonstrated to verify the effectiveness of the proposed method. Finally, we conclude the paper with a summary in Section 6.

2  PHT-Splines

In this section, we give a brief review of the definition of PHT-splines. The readers are recommended to refer to [1] for more details.

2.1 Definition of PHT-Splines

Given a T-mesh T, ℱ denotes all the cells in T and Ω denotes the region occupied by ℱ. The polynomial spline space over T is defined as

𝒮(m,n,α,β,T):={s(x,y)∈Cα,β(Ω)|s(x,y)|ϕ∈Pmnfor any ϕ∈ℱ},(1)

where Pmn is the space of all the polynomials with bi-degree (m,n), and Cα,β(Ω) is the space consisting of all the bivariate functions which are continuous in Ω with order α in the x-direction and with order β in the y-direction. In particular, in the case of m≥2α+1,n≥2β+1, 𝒮(m,n,α,β,T) is well studied and its dimension is given in [23].

PHT is defined by 𝒮(3,3,1,1,T) with T being a hierarchical T-mesh. Here, a hierarchical T-mesh is a special type of T-mesh that has a natural hierarchical structure and defined in a recursive fashion. One generally starts from a tensor product mesh (level 0). From level k to level k+1, a cell at level k is subdivided into four sub-cells which are cells at level k+1.

The dimension formula of PHT has a concise expression with

𝒮(3,3,1,1,T)=4(Vb+V+),(2)

where Vb and V+ represent the number of boundary vertices and interior crossing vertices in T respectively. As defined in [1], a boundary vertex or an interior crossing vertex is called a basis vertex. This dimension formula implies that each basis vertex should be associated with four basis functions. Thus, the construction of PHT is to construct four bicubic C1,1 continuous functions for each basis vertex so that the functions span the PHT and retain some good properties.

2.2 Basis Construction

The original PHT basis functions are constructed by truncation mechanism in [1] level by level in order to make basis functions vanish at other basis vertices. Fig. 1 shows two typical hierarchical meshes, where the basis functions associated the marked vertices have undesired decay as level increases. In Fig. 1a, the maximum values decay rapidly as the level increases. In Fig. 1b, the maximum value does not decrease, but the basis function decays sharply along the refinement direction as level increases. This decay phenomenon of basis functions is first noted in [2]. The decay phenomenon makes the stiffness matrix have large condition number and thus produces unstable solutions.

Figure 1: The decay phenomenon of PHT-spline basis functions associated with the vertex marked by a yellow circle

A new basis functions are proposed in [2] to amend the decay phenomenon. In the following, we give a brief review of the construction.

Definition 1 Suppose v is a basis vertex in T, the support mesh of v, denoted by Tv, is defined as the minimal 2×2 tensor product mesh, where v is the central vertex and all edges are subordinate to T.

Note that if vis a boundary vertex, the support mesh is actually a 2×1, 1×2 or 1×1 tensor product mesh. For any given hierarchical T-mesh, the support mesh of each basis vertex exists and is unique. The support mesh can be found efficiently.

Fig. 2 illustrates the support meshes of six basis vertices including three boundary vertices and three interior crossing vertices, where the six basis vertices are marked by solid circles and the corresponding support meshes are marked by bold lines.

Figure 2: Support meshes of six basis vertices

The basis functions associated with a basis vertex are then defined by the four B-spline functions that defined on the support mesh. Specifically, for a basis vertex v=(si,tj), suppose the support mesh is defined by (si0,si,si1)×(tj0,tj,tj1), then the four bases are defined as follows:

b0(s,t)=N3[ si0,si0,si,si,si1 ][ s ]×N3[ tj0,tj0,tj,tj,tj1 ][ t ],b1(s,t)=N3[ si0,si,si,si1,si1 ][ s ]×N3[ tj0,tj0,tj,tj,tj1 ][ t ],b2(s,t)=N3[ si0,si0,si,si,si1 ][ s ]×N3[ tj0,tj,tj,tj1,tj1 ][ t ],b3(s,t)=N3[ si0,si,si,si1,si1 ][ s ]×N3[ tj0,tj,tj,tj1,tj1 ][ t ].(3)

These four basis functions are linearly independent and have the same support [si0,si1]×[tj0,tj1]. This way of definition has been confirmed in [2] without decay phenomenon.

3  Weighted Bases for PHT-Splines

3.1 Geometric Information

For any basis function b(s,t), we define a geometric information (the function value, the first order partial derivatives and the mixed partial derivative) operator by

Lb(s,t)=(b(s,t),bs(s,t),bt(s,t),bst(s,t)).(4)

For a basis vertex v=(si,tj), the associated four PHT-splines basis functions are defined by (3). We collect the geometric information of b0,b1,b2,b3 at the basis vertex v into the matrix G,

G=(Lb0(v)Lb1(v)Lb2(v)Lb3(v))=(λμ−3αμ−3βλ9αβ(1−λ)μ3αμ−3β(1−λ)−9αβλ(1−μ)−3α(1−μ)3βλ−9αβ(1−λ)(1−μ)3α(1−μ)3β(1−λ)9αβ),(5)

where

α=1Δs1+Δs2,β=1Δt1+Δt2,λ=αΔs2,μ=βΔt2,

Δs1=si−si0,Δs2=si1−si,Δt1=tj−tj0,Δt2=tj1−tj.(6)

The sum of the geometric information of the four basis functions b0,b1,b2,b3 at the associated basis vertex is equal to (1,0,0,0), because the sum of the four basis functions is identically equal to 1. This can also be verified by the column sums in G. As non-vanishing basis functions at basis vertices are not limited to the associated four basis functions, so the basis functions defined in [2] do not form a partition of unity generally.

Suppose that the maximal subdivision level of T is N and Tk denotes the hierarchical T-mesh at level k,k=0,⋅⋅⋅,N. Then we have T=TN. For a basis vertex v, if v∈Tk but v∉Tk−1, then it is called a basis vertex at level k,k=1,⋅⋅⋅,N. Specially, every vertex in T0 is a basis vertex at level 0. The ith basis vertex at level k is denoted by vik and the four basis functions associated with vik are denoted as b4i+jk, j=0,1,2,3.

According to the definitions of support meshes and hierarchical T-meshes, the basis function proposed in [2] has the following property:

•   A basis vertex at level l can not be in the interior of the support of any basis function at level k>l, that is

whenl<k,b4i+jk(vgl)=0,forj=0,1,2,3,i=1,⋯,nk,g=1,⋯,nl,(7)

where nk and nl are the number of basis vertices on level k and l, respectively.

•   The basis functions b4i+jk vanish at the basis vertices on level k, except for the vertices associated with them, that is

whenl=k,b4i+jk(vgl)=0,forg≠i,j=0,1,2,3,k=0,⋯,N.(8)

3.2 Computing Weights

The proposed weighted bases are defined by multiplying the basis functions constructed in [2] by positive weights. This ensures that the weighted bases not only have no decay, but also form the positive partition of unity. In the following, we are going to compute weights w4i+jl>0 for each basis function such that

f(s,t):=∑l=0N∑i=1nl∑j=03w4i+jlb4i+jl(s,t)≡1,(s,t)∈[0,1]×[0,1].(9)

This problem can be understood as the PHT fitting problem of a constant function. Since a PHT-spline space spans the bicubic polynomial space P33, the weights w4i+jl satisfying (9) exist. As described in [1], the control points of the PHT-spline surface, here we call them weights, can be determined based on the geometric information at the basis vertices.

The weights are computed level by level. First, we evaluate the PHT-spline surface at the basis vertices on level 0, that is

∑l=0N∑i=1nl∑j=03w4i+jlb4i+jl(vg0)=∑j=03w4g+j0b4g+j0(vg0)=1,g=1,⋯,n0,(10)

since the basis functions on level l>0 vanish on any basis vertex on level 0. Furthermore, B-spline basis functions b4g+j0,j=0,1,2,3 form a partition of unity, thus we obtain w4g+j0=1, for j=0,1,2,3 and g=1,⋯,n0.

For a basis vertex vgk on level k>0, the non-vanishing basis functions at vgk are composed of two parts: the non-vanishing basis functions at vgk on level l<k and the four basis functions associated with vgk. The indices of basis vertices associated with the first part are denoted by K,

K={j|bjl(vgk)≠0,l<k},

and the weighted sum of the basis functions in K is denoted by h(s,t),

h(s,t)=∑l<k,j∈Kwjlbjl(s,t).(11)

Because the geometric information operator is linear, one has

(Lf)(vgk)=L(h(s,t)+∑j=03w4g+jkb4g+jk(s,t))|vgk=Lh(vgk)+∑j=03w4g+jkLb4g+jk(vgk).(12)

Thus, the weights of the four bases b4g+jk are determined by the following linear equations:

GT(w4gkw4g+1kw4g+2kw4g+3k)=(1−h(vgk)−∂h∂s(vgk)−∂h∂t(vgk)−∂2h∂s∂t(vgk)),G=(Lb4gk(vgk)Lb4g+1k(vgk)Lb4g+2k(vgk)Lb4g+3k(vgk)).(13)

The matrix G is invertible, so the solution of (13) exists and is unique. Specifically, the weights are explicitly expressed by

(w4gkw4g+1kw4g+2kw4g+3k)=(1111λ−13αλ3αλ−13αλ3αμ−13βμ−13βμ3βμ3β−λ+μ−λμ−19αβ−λ−λμ9αβ−μ−λμ9αβλμ9αβ)T(1−h(vgk)−∂h∂s(vgk)−∂h∂t(vgk)−∂2h∂s∂t(vgk)).(14)

The non-vanishing basis functions at vgk can be found easily owing to the simple construction. The weights can be computed explicitly based on (14) without the need to solve a linear system of equations.

Fig. 3a shows a typical hierarchical mesh where the diagonal elements are refined. The support meshes of some marked basis vertices are shaded in orange color. We discuss the weights of the basis functions associated with colored vertices.

Figure 3: (a) The weights associated v00,v01,v02 (marked by red circles) are all equal to [1, 1, 1, 1]. The weights for vi2 are equal to [716,1316,1316,1516]. There are two non-vanishing basis functions (associated with the vertices marked by green squares) at vi2. (b) One element on level 1 is refined. The weights associated vi2 are equal to [0, 0, 0, 0]. (c) The weights associated vi2 are equal to [0.5625, 0.1875, 0.1875, 0.0625]

•   The weights of v00, v01 and v02 (marked by red circles) are all equal to [1, 1, 1, 1], since all basis functions are zeros at these three vertices, except for the basis functions associated with themselves, which means h(s,t)≡0.

•   The weights for vi11 are computed as follows: (1) the right term h(s,t)=∑j=03w4i0+j0b4i0+j0 since only the support of vi00 contains vi11; (2) the elements in G are constructed with Δs1=Δs2=0.125,Δt1=Δt2=0.125,λ=μ=12,α=β=4, based on (6). To sum up,

(w4i11w4i1+11w4i1+21w4i1+31)=(11−12412411−124124−124−1241576−1576124124−15761576).(15)

Then, we obtain [w4i11,w4i1+11,w4i1+21,w4i1+31]=[716,1316,1316,1516].

•   The non-vanishing basis functions at vi2 (marked by red circles) are denoted by vi00 and vi11 (marked by green squares). The weighted sum is expressed as

h(s,t)=∑j=03(w4i0+j0b4i0+j0+w4i0+j1b4i0+j1).(16)

Submitting h(s,t) into (14), the weights associated with vi2 are equal to [716,1316,1316,1516].

Fig. 3b shows a case where the weights for vi2 are equal to [0,0,0,0]. The non-vanishing basis functions on F (the refined element) are also the non-vanishing basis functions at vi2. Notice that the weights on level 1 are computed such that the basis functions form a partition of unity, so h(s,t)|F≡1, which means the right term in (14) is a zero vector. Since G is invertible, the weights associated with vi2 are thus zeros.

The basis functions correspond to zero weights has no contribution in the entire representation. Notice that zero weights occur if and only if the function g form a partition of unity on the domain occupied by the support mesh of the considered basis vertex. Therefore, to avoid this, the isolated refined elements, like the case in Fig. 3b, are not allowed.

In Fig. 3c, two adjacent elements are refined. The support of vi11 is shaded by yellow color in Figs. 3b and 3c. The support of vi11 in Fig. 3c becomes smaller, thus h(s,t) do not form a partition of unity on F anymore. Therefore, the weights for vi2 are non-zeros.

Based on the above analysis, we made the following assumptions for the hierarchical T-mesh to define positively weighted bases that form a partition of unity.

Definition 2 A hierarchical T-mesh is called an admissible hierarchical T-mesh if there are not any isolated refined elements in the hierarchical T-mesh. If an element on level l is refined, but none of its adjacent elements are refined, it is called an isolated refined element. Here, if two elements of the same level share at least one edge, they are called adjacent elements.

Figs. 3a and 3c are two admissible hierarchical T-meshes, while Fig. 3b is not an admissible hierarchical T-mesh. A hierarchical T-mesh can become an admissible hierarchical T-meshes by performing refinement on adjacent elements of isolated refined elements.

To ensure the refinement is highly localized, we generally require the refinement level between any two neighboring elements cannot be greater than one, which is also required in [18]. Under this requirement, hierarchical T-meshes are admissible.

Proposition 1 There uniquely exists a set of positive weights for the bases constructed in [2] on the admissible hierarchical T-mesh such that the weighted bases form a partition of unity, that is there uniquely exists w1,w2,⋅⋅⋅,wn>0 such that ∑i=1nwibi(u,v)≡1 for (u,v)∈[a,b]×[c,d].

Proof The existence and uniqueness of the weights w4i+jl are implied by the linear system (13). In the following, we only prove w4i+jl>0.

For the basis functions bi0 associated with the basis vertices on level 0, the function h(s,t) defined in (11) is equal to zero, i.e., h(s,t)=0, then the weights solved by (13) are all equal to ones.

For an element F in the initial mesh T0, we denote the B-spline functions associated with the four corners of F by b4i0+j,b4i1+j,b4i2+j,b4i3+j,j=0,1,2,3. Remember that these basis functions are defined over T0, instead of the final hierarchical T-mesh T. Then, the non-vanishing basis functions on F are only these basis functions, which also form a partition of unity on F. As the refinement proceeds, the supports of some of these basis functions are changed, and some new basis functions are added. The PHT-spline space defined over T0 are contained in the PHT-spline defined over T, then b4i0+j,b4i1+j,b4i2+j, b4i3+j can be exactly represented by the basis functions b4i+kl on level l by knot insertion algorithm, that is, for j=0,1,2,3,

b4i0+j=∑l=0N∑i=0nl∑k=03ciklb4i+kl,b4i1+j=∑l=0N∑i=0nl∑k=03diklb4i+kl,b4i2+j=∑l=0N∑i=0nl∑k=03eiklb4i+kl,b4i3+j=∑l=0N∑i=0nl∑k=03fiklb4i+kl.(17)

Since

∑j=03b4i0+j+b4i1+j+b4i2+j+b4i3+j≡1,on F,

then

∑l=0N∑i=0nl∑k=03(cikl+dikl+eikl+fikl)b4i+kl≡1 on F.

This implies that the weights

w4i+kl=cikl+dikl+eikl+fikl≥0,(18)

because cikl,dikl,eikl and fikl are the coefficients in the knot insertion algorithm, thus are nonnegative values. On an admissible hierarchical T-mesh, all the weights are non-zeros, thus w4i+jl>0.

In Fig. 4, we show the PHT surfaces defined over the same hierarchical T-mesh and by the same control points, but using different bases: the proposed weighted bases and non-decay bases [2]. The control points are designed on the plane z=1. We see that the PHT surface defined by weighted basis functions form a partition of unity.

Figure 4: The weighted PHT-splines proposed in this paper form the partition of unity, while the PHT-splines defined by the non-decay bases [2] cannot form a partition of unity

We conclude the method of computing the weights such that the weighted basis functions form a partition of unity in Algorithm 1.

4  Subdivision Scheme for PHT-Splines

Considering a knot vector {s0,s0,s1,s1,s2,s2,s3,s3}, we insert a double knot s′∈[s1,s2] in it. Then, we get the following relationship of B-splines before and after s′ is inserted based on the B-spline knot insertion algorithm,

(N3[s0,s0,s1,s1,s2][s]N3[s0,s1,s1,s2,s2][s]N3[s1,s1,s2,s2,s3][s]N3[s1,s2,s2,s3,s3][s])=W(N3[s0,s0,s1,s1,s′][s]N3[s0,s1,s1,s′,s′][s]N3[s1,s1,s′,s′,s2][s]N3[s1,s′,s′,s2,s2][s]N3[s′,s′,s2,s2,s3][s]N3[s′,s2,s2,s3,s3][s]),

W=(1αsαsβs00001−αs(2−αs−βs)βsβs20000(1−βs)2(1−βs)(βs+γs)γs0000(1−βs)(1−γs)1−γs1),(19)

where αs=s2−s′s2−s0, βs=s2−s′s2−s1 and γs=s3−s′s3−s1.

Based on the above formula, we derive the subdivision scheme for PHT-splines as follows. The support mesh of a basis vertex vik is supposed to be [si0,si,si1]×[tj0,tj,tj1]. The knot intervals of the support mesh are denoted by d0=si−si0,d1=si1−si,e0=tj−tj0,e1=tj1−tj. When elements are refined, new basis vertices are added. The new basis vertices lying in the interior of the refined elements are called face vertices, while the new basis vertices lying on the edge of the mesh are called edge vertices. The old basis vertices are called vertex points. Under the assumption that the level difference in an admissible hierarchical T-mesh is not greater than 1, the subdivision scheme of PHT-splines are deduced as follows.

•   Vertex-points

If the refinement changes the support mesh associated with vik, the control points associated with vik are updated as follows:

(P4ik+1P4i+1k+1P4i+2k+1P4i+3k+1)=(ws1wt1(1−ws1)wt1ws1(1−wt1)(1−ws1)(1−wt1)ws2wt1(1−ws2)wt1ws2(1−wt1)(1−ws2)(1−wt1)ws1wt2(1−ws1)wt2ws1(1−wt2)(1−ws1)(1−wt2)ws2wt2(1−ws2)wt2ws2(1−wt2)(1−ws2)(1−wt2))(P4ikP4i+1kP4i+2kP4i+3k).(20)

where ws1=μ0d0+2μ1d12μ0d0+2μ1d1, ws2=μ1d12μ0d0+2μ1d1, wt1=λ0e0+2λ1e12λ0e0+2λ1e1 and wt2=λ1e12λ0e0+2λ1e1. The parameters μ0,μ1,λ0,λ1 are used to indicate whether the knot intervals d0,d1,e0,e1 are refined, respectively. The parameter is equal to 1 implies the corresponding knot interval is subdivided. Especially when the parameters μ0,μ1,λ0,λ1 are zeros, control points are not updated, that is P4i+jk+1=P4i+jk,j=0,1,2,3. Fig. 5 shows different cases of refinement. In Fig. 5a, four elements in the support mesh of vik are refined, now the parameters are set as [μ0,μ1,λ0,λ1]=[1,1,1,1]. In Figs. 5b and 5c, three elements are refined and the corresponding parameters are set as [μ0,μ1,λ0,λ1]=[0,1,0,1] and [μ0,μ1,λ0,λ1]=[1,0,1,0]. In Figs. 5d and 5e, two elements are refined and the corresponding parameters are set as [μ0,μ1,λ0,λ1]=[0,0,1,0] and [μ0,μ1,λ0,λ1]=[1,0,0,0]. In Fig. 5f, one element is refined. If vik is on the boundary, [μ0,μ1,λ0,λ1]=[0,1,0,1]. Otherwise, [μ0,μ1,λ0,λ1]=[0,0,0,0], which means the old vertex vik is not updated.

Figure 5: Subdivision scheme of vertex points

•   Edge-points

Suppose vik+1 is a edge point lying on the edge with endpoints vi1k and vi2k, referring to Fig. 6. Since the level difference is assumed to be less than 2, vi1k and vi2k are two basis vertices. As shown in Fig. 6a, the knot intervals for vi1k and vi2k are denoted by {d0,d1;e0,e1} and {d1,d2;e0,e1}, respectively. We have the following edge-points updating formula.

P4ik+1=18(d0+d1)(e0+e1)[(3d0+2d1)(e0+2e1)P4i1+1k+d1e0P4i1+2k+d1(2e1+e0)P4i1k+(3d0+2d1)e0P4i1+3k+(d0+d1)(e0+2e1)P4i2k+(d0+d1)e0P4i2+2k],

P4i+1k+1=18(d1+d2)(e0+e1)[(2d1+3d2)(e0+2e1)P4i2k+d1e0P4i2+3k+d1(e0+2e1)P4i2+1k+(2d1+3d2)e0P4i2+2k(d1+d2)(e0+2e1)P4i1+1k+(d1+d2)e0P4i1+3k],

P4i+2k+1=18(d0+d1)(e0+e1)[(3d0+2d1)(2e0+e1)P4i1+3k+d1e1P4i1k+d1(2e0+e1)P4i1+2k+(3d0+2d1)e1P4i1+1k+(d0+d1)e1P4i2k+(d0+d1)(2e0+e1)P4i2+2k],

P4i+3k+1=18(d1+d2)(e0+e1)[(2d1+3d2)(2e0+e1)P4i2+2k+d1e1P4i2+1k+d1(2e0+e1)P4i2+3k+(2d1+3d2)e1P4i2k+(d1+d2)e1P4i1+1k+(d1+d2)(2e0+e1)P4i1+3k].(21)

Figure 6: The subdivision scheme of edge points and face points

When the edge point lies on a vertical edge, as shown in Fig. 6b, the update formula is similar to the above one, except that eis and dis are swapped.

•   Face-points

Since we assume the level difference in the mesh is not greater than 1, so when an element is refined, the four corners are exactly basis vertices. Suppose the four basis vertices of the refined face F are denoted by vi1,vi2,vi3,vi4. The knot intervals of these four basis vertices are denoted by {d0,d1,d2;e0,e1,e2}, see Fig. 6c for a reference. Here, we use the parameter aj=1 and aj=0 to indicate the support mesh of vij is changed and unchanged, respectively. Let P4ij+lk:=ajP4ij+lk. We have the following formulas:

P4ik+1=116(d0+d1)(e0+e1)[d1e1P4i1k+(3d0+2d1)e1P4i1+1k+d1(3e0+2e1)P4i1+2k+(3d0+2d1)(3e0+2e1)P4i1+3k+(d0+d1)e1P4i2k+(d0+d1)(3e0+2e1)P4i2+2k+d1(e0+e1)P4i3k+(3d0+2d1)(e0+e1)P4i3+1k+(d0+d1)(e0+e1)P4i4k],

P4i+1k+1=116(d1+d2)(e0+e1)[(2d1+3d2)e1P4i2k+d1e1P4i2+1k+(2d1+3d2)(3e0+2e1)P4i2+2k+d1(3e0+2e1)P4i2+3k+(d1+d2)e1P4i1+1k+(d1+d2)(3e0+2e1)P4i1+3k+(2d1+3d2)(e0+e1)P4i4k+d1(e0+e1)P4i4+1k+(d1+d2)(e0+e1)P4i3+1k],

P4i+2k+1=116(d0+d1)(e1+e2)[d1(2e1+3e2)P4i3k+d1e1P4i3+2k+(3d0+2d1)(2e1+3e2)P4i3+1k+(3d0+2d1)e1P4i3+3k+d1(e1+e2)P4i1+2k+(3d0+2d1)(e1+e2)P4i1+3k+(d0+d1)(2e1+3e2)P4i4k+(d0+d1)e1P4i4+2k+(d0+d1)(e1+e2)P4i2+2k],(22)