Stable Simplex Spline Bases for C3 Quintics on the Powell–Sabin 12-Split

STABLE SIMPLEX SPLINE BASES

FOR C³ QUINTICS ON THE POWELL-SABIN 12-SPLIT

TOM LYCHE AND GEORG MUNTINGH

Abstract. For the space ofC³ quintics on the Powell-Sabin 12-split of a triangle, we determine explicitly the six symmetric simplex spline bases that reduce to a B-spline basis on each edge, have a positive parti- tion of unity, a Marsden identity that splits into real linear factors, and an intuitive domain mesh. The bases are stable in theL∞norm with a condition number independent of the geometry, have a well-conditioned Lagrange interpolant at the domain points, and a quasi-interpolant with local approximation order 6. We show anh² bound for the distance be- tween the control points and the values of a spline at the corresponding domain points. For one of these bases we deriveC⁰,C¹,C²andC³con- ditions on the control points of two splines on adjacent macrotriangles.

Keywords Stable basis · Powell-Sabin 12-split · Simplex splines · Marsden identity·Quasi-interpolation

Mathematics Subject Classification 41A15·65D07·65D17

1. Introduction

Piecewise polynomials or splines defined over triangulations form an indis- pensable tool in the sciences, with applications ranging from scattered data fitting to finding numerical solutions to partial differential equations. See [3, 14] for comprehensive monographs.

In applications like geometric modelling [4] and solving PDEs by isogeo- metric methods [11], one often desires a low degree spline with C¹,C²orC³ smoothness. For a general triangulation, it was shown in Theorem 1.(ii) of [30] that the minimal degree of a triangularC^relement is 4r+1, e.g., degrees 5,9,13 for the classesC¹,C²,C³. To obtain smooth splines of lower degree one can split each triangle in the triangulation into several subtriangles.

One such split is the Powell-Sabin 6-split of a triangle , which has been used for the construction of B-spline-like bases over arbitrary triangu- lations [7, 26, 28]. These bases form a nonnegative partition of unity, admit a Marsden identity and yield stable quasi-interpolants [15, 29]. Moreover, by imposing local C³ smoothness across the interior edges of , the quin- tic spline space can be reduced while maintaining the approximation order [1, 27].

In this paper we instead focus on the Powell-Sabin 12-split ; see Figure 2.

This split is suitable for the construction of simplex splines, as it is the complete graph of the vertices and midpoints of a triangle. Similar to the 6-split, global C¹ smoothness can be obtained with degree only 2 [20], and C² smoothness with degree only 5 [13, 16, 24] on any (planar) triangulation.

1

Final version available at SpringerLink : http://dx.doi.org/10.1007/s00365-016-9332-8

Once a space is chosen one determines its dimension. The space S₂¹ and S₅³ ofC¹ quadratics andC³ quintics on the 12-split of a single triangle have dimension 12 and 39, respectively. For a general triangulationT of a polyg- onal domain, we can replace each triangle in T by its 12-split, obtaining a refined triangulationT₁₂. The dimensions of the correspondingC¹ quadratic andC² quintic spaces (the latter withC³supersmoothness at the vertices of T and interior edges of the 12-split of each triangle inT) are 3|V|+|E|and 10|V|+ 3|E|, where |V| and |E| are the number of vertices and edges in T. Moreover, in addition to providing C¹ and C² spaces on any triangulation, these spaces are suitable for multiresolution analysis [6, 8, 16, 19].

For a general smoothness, degree, and triangulation, it is a hard problem to find a basis of the corresponding spline space with all the usual properties of the univariate B-spline basis. One reason is that it is difficult to find a sin- gle recipe for the various valences and topologies of the triangulations. In [5]

a basis, called the (quadratic) S-basis, was constructed for theC¹quadratics on the Powell-Sabin 12-split on one triangle. The S-basis consists of sim- plex splines [17, 21] and has all the usual properties of univariate B-splines, including a recurrence relation down to piecewise linear polynomials and a Marsden identity. Moreover, analogous to the Bernstein-B´ezier case, C⁰ and C¹ smoothness conditions were given, tying the S-bases on neighboring triangles together to give C¹ smoothness on the refined triangulation T₁₂.

In the next section we recall some background on splines and the 12-split.

Section 3 introduces dimension formulas for spline spaces on the 12-split. In Section 4 some basic properties of simplex splines are recalled, after which an exhaustive list is derived of the C³ quintic simplex splines on the 12-split that reduce to a B-spline on the boundary. Section 5 introduces a barycentric form of the Marsden identity and describes how the dual basis in [16] can be applied to find the simplex spline bases of S₅³ satisfying a Marsden identity.

Making use of the computer algebra system Sage [22], these techniques are applied in Section 6 to discover six symmetric simplex spline bases that reduce to a B-spline basis on each boundary edge, have a positive partition of unity, a Marsden identity, and domain points with an intuitive control mesh and unisolvent for S₅³. The concise and coordinate independent form of the barycentric Marsden identity makes it possible to state these results in Table 3. Analogous to the Bernstein-B´ezier case, we findC⁰,C¹,C² and even C³ conditions for one of these bases, tying its splines together on a general triangulation. One of the latter smoothness conditions involves just the control points in a single triangle and the geometry of the neighboring triangles, showing that C³ smoothness using S₅³ cannot be achieved on a general refined triangulationT₁₂. Finally a conversion to the Hermite nodal basis from [16] is provided.

2. Background

2.1. Notation. On a triangulation T of a polygonal domain Ω ⊂ R² we define the spline spaces

S_d^r(T) :={f ∈C^r(Ω) :f| ∈Π_d∀ ∈ T }, r, d∈Z, r≥ −1, d≥0,

(a) (b) (c)

Figure 1. Domain mesh for the quintic B´ezier basis (left), the S-basis from [5] (middle) and the basis B_c in this paper (right).

where Πd is the space of bivariate polynomials of total degree at most d.

With the convention Π_d:=∅and dim Π_d:= 0 if d <0, one has dim Π_d= 1

2(d+ 2)(d+ 1)₊, d∈Z, where z₊ := max{0, z}for any real number z.

Any point x in a nondegenerate triangle [v1,v2,v3] can be represented by its barycentric coordinates (β1, β2, β3), which are uniquely defined by x =β₁v₁+β₂v₂+β₃v₃ and β₁+β₂+β₃ = 1. Similarly, each vector u is uniquely described by itsdirectional coordinates, i.e., the triple (β1−β₁⁰, β2− β₂⁰, β₃−β₃⁰) with (β₁, β₂, β₃) and (β₁⁰, β⁰₂, β₃⁰) the barycentric coordinates of two pointsxandx⁰ such thatu=x−x⁰. Sometimes we writexas a linear combination of more than three vertices, in which case these coordinates are no longer unique.

A bivariate polynomial p of total degree ddefined on a triangle ⊂R² is conveniently represented by its B´ezier form

p(x) = X

i+j+k=d i,j,k≥0

c_ijkB_ijk^d (x), B^d_ijk(x) := d!

i!j!k!β₁ⁱβ₂^jβ₃^k,

with (β₁, β₂, β₃) the barycentric coordinates of x with respect to . Here theB^d_ijkare theBernstein basis polynomials of degreedand the coefficients c_ijkare theB´ezier ordinates ofp. We associate each B´ezier ordinatec_ijk∈R to thedomain pointξ_ijk:= _dⁱv₁+_d^jv₂+^k_dv₃ ∈R²and combine them into the control point (ξ_ijk, c_ijk)∈R³. By connecting any two domain pointsξ_i1j1k1 andξ_i2j2k2 by a line segment whenever|i₁−i2|+|j₁−j2|+|k₁−k₂|= 1, one arrives at thedomain mesh ofp; see Figure 1a. Thecontrol mesh is defined similarly by connecting control points.

We consider finite multisets K = {v^m₁¹· · ·v^m_nⁿ} ⊂ R², in which the distinct elements v₁, . . . ,v_n are counted with corresponding multiplicities m1, . . . , mn ≥ 0. Write |K| := m1 +· · · +mn for the total number of elements inK. For any two integers i, j, the Kronecker delta is the symbol

δ_ij :=

1 ifi=j, 0 otherwise.

1 2 3

4

5 6

7 8

9 1 10

2 3

4 5 6 7

8 9 11 10 12

(a)

1

2 3

4 5 6 7 12 8 9 11 10

(b)

Figure 2. The Powell-Sabin 12-split with labelling of ver- tices and faces (left), and a scheme assigning every point in the macrotriangle to a unique face of the 12-split (right).

For any set A⊂R², define the indicator function 1_A:R² −→ {0,1}, 1_A(x) :=

1 ifx∈A, 0 ifx /∈A.

2.2. The Powell-Sabin 12-split. Given a triangle = [v₁,v₂,v₃] with vertices v1,v2,v3 write e1 := [v2,v3], e2 := [v3,v1], and e3 := [v1,v2] for its (nonoriented) edges. Connecting vertices and the edge midpoints v₄ := (v₁+v₂)/2,v₅:= (v₂+v₃)/2 andv₆:= (v₁+v₃)/2, we arrive at the Powell-Sabin 12-split of ; see Figure 2a for the labelling of the vertices v₁, . . . ,v₁₀ and faces ₁, . . . , ₁₂.

To decide to which face of the 12-split points on the interior edges be- long, we follow the convention in [25] shown in Figure 2b, which can be quickly computed by Algorithm 1.1 in [5]. If K is a multiset satisfying K⊂ {v₁, . . . ,v₁₀}as sets, its convex hull [K] is a union of some of the faces of the 12-split, and we define the half-open convex hull [K) as the union of the corresponding half-open faces depicted in Figure 2b.

2.3. A basis for the dual space of a space of C³ quintics. Let be the 12-split of a triangle with vertices V = {v₁,v₂,v₃} and edges E ={[v₁,v2], [v2,v3],[v3,v1]}. For any edge e= [vi,vj]∈ E with opposing vertex v_k, let

q_1,e := 3v_i+v_j

4 , me := v_i+v_j

2 , q_2,e := v_i+ 3v_j 4

be its midpoint and quarterpoints. With εv the point evaluation at v and D_u the directional derivative in the direction u, let

(1) Λ := [

v∈V

[

i+j≤3 i,j≥0

{ε_vDⁱ_xvD^j_y

v} ∪ [

e∈E

{ε_q

1,eD²_ue, ε_meD_ue, ε_q2,eD²_ue};

Figure 3. A basis for the dual space of S₅³ on a single tri- angle. A bullet represents a point evaluation, three circles represent all derivatives up to order three, and a single and double arrow represent a first- and second-order directional derivative. These derivatives are evaluated at the rear end of the arrows, which are located at the midpoints and quarter- points.

see Figure 3. Here, for every vertex v and edge e, the symbols x_v,y_v are any linearly independent vectors and ue is any vector not tangent to e, for example the outside unit normal as shown in Figure 3. It was shown in [16, Theorem 4] that Λ is a basis for the dual space to S₅³( ), which therefore has dimension 39.

3. Dimension formulas

Consider a polygonal domain Ω ⊂ R² with a triangulation T with sets of vertices V, edges E, faces F, and 12-split refinement T₁₂. Let

S₅^2,3(T₁₂) :={s∈ S₅²(T₁₂) :s∈C³( ), ∈ F and s∈C³(v),v∈ V}.

Here s ∈ C³(v) means that all polynomials s| such that is a triangle with vertex at v have common derivatives up to order three at the pointv.

Note that ifT consists of a single triangle, thenS₅^2,3(T₁₂) =S₅³(T₁₂).

Since Λ from (1) specifies the value and partial derivatives up to or- der three at each vertex of T and the value, first- and second-order cross- boundary derivatives at each edge of T, the following theorem is an imme- diate consequence of [16, Theorem 4].

Theorem 1. For any triangulation T with |V| vertices and |E| edges, the set Λ is a basis for the dual space of S₅^2,3(T₁₂). In particular

dimS₅^2,3(T₁₂) = 10|V|+ 3|E|.

Next, let denote the 12-split triangulation of a single triangle . For future reference, we state the following formula for the dimension of the space of C^r splines of degreedon , which is a special case of Theorem 3.1 in [2] (also compare [23]).

dimS_d^r( ) C⁻¹ C⁰ C¹ C² C³ C⁴ C⁵ C⁶ C⁷ C⁸ C⁹

d= 0 12 1

d= 1 36 10 3

d= 2 72 31 12 6

d= 3 120 64 30 16 10

d= 4 180 109 60 34 21 15 d= 5 252 166 102 61 39 27 21 d= 6 336 235 156 100 66 46 34 28 d= 7 432 316 222 151 102 73 54 42 36 d= 8 540 409 300 214 150 109 81 63 51 45 d= 9 660 514 390 289 210 154 117 91 73 61 55

Table 1. Dimensions of S_d^r( ), with (r, d) = (3,5) highlighted.

Theorem 2. For any integers d, r withd≥0 andd≥r≥ −1,

(2)

dimS_d^r( ) = 1

2(r+ 1)(r+ 2) + 9

2(d−r)(d−r+ 1) + 3

2(d−2r−1)(d−2r)₊+

d−r

X

j=1

(r−2j+ 1)₊.

To quickly look up dimS_d^r( ) for small values ofr and d, we have listed these first dimensions in Table 1.

4. Simplex splines

In this section we first recall the definition and some basic properties of the simplex spline, and then proceed to determine the C³ quintic simplex splines on the 12-split that reduce to a B-spline on the boundary. For a comprehensive account of the theory of simplex splines, see [17, 21].

4.1. Definition and properties. The following definition of the simplex spline is convenient for our purposes.

Definition 1. For any finite multiset K= {v^m₁¹· · ·v^m₁₀¹⁰} ⊂ R² composed of vertices of , the (area normalized) simplex spline Q[K] : R² −→ R is recursively defined by

Q[K](x) :=







0 if area([K]) = 0,

1_[K)(x)_area([K])^{area( )} if area([K])6= 0 and |K|= 3, P10

j=1β_jQ[K\v_j](x) if area([K])6= 0 and |K|>3, withx=β₁v₁+· · ·+β₁₀v₁₀, β₁+· · ·+β₁₀= 1, andβ_i = 0 wheneverm_i= 0.

By Theorem 4 in [17] this definition is independent of the choice of the βj. It is well known that Q[K] is a piecewise polynomial with support [K] ⊂ and of total degree at most |K| −3. One shows by induction on |K| that R

R²Q[K](x)dx = area( ) · ^|K|−1₂ ⁻¹

. Although M[K] :=

area( )^{−1 |K|−1}₂

Q[K] has unit integral and is used more frequently, our discussion is simpler in terms of Q[K].

Whenever m₇=m₈ =m₉ =m₁₀= 0, we use the graphical notation i j

k lm

n :=Q[vⁱ₁v^j₂v^k₃v^l₄v^m₅ vⁿ₆].

Example 1. The linear S-spline basis in [5] only uses vertices v1, . . . ,v6, v₁₀, while the quadratic S-spline basis only uses v₁, . . . ,v₆. It is given by (3) S_j,2= area([K_j,2])

6 M[K_j,2] = area([K_j,2])

area( ) Q[K_j,2], j = 1, . . . ,12, where by (2.5) in [5], as (unordered) sets,

(4)

{Q[K_1,2], . . . , Q[K12,2]}= n

3 11 , ₁¹₃, ¹³¹ , ₂¹_{1 1}, ¹₁¹₂, ₁¹²¹ , 11111, ¹¹₁¹₁, ₁¹¹₁¹ , _{1 1}¹₂, ₂¹¹₁ , ¹²¹₁o

. Example 2. If K = {v^µ_iⁱ⁺¹v^µ_j^j+1v^µ_k^k+1} has three distinct elements with area([K])>0, then, with (βi, βj, β_k) the barycentric coordinates of x with respect to [v_i,v_j,v_k], it follows by induction that

Q[K](x) =1_[K)(x) area( ) area([K])

(µi+µj+µk)!

µ_i!µ_j!µ_k! β_i^µiβ_j^µjβ_k^µk, which is, up to a scalar, a Bernstein polynomial on [K).

Continuity. For any edgeeof , ifKhas at mostmknots (counting multi- plicities) along the affine hull ofe, thenQ[K] is|K|−m−2 times continuously differentiable overe. For instance, everyC³ quintic simplex spline on has at most three knots along the affine hull of any interior edge e.

Differentiation. LetK={v^m₁¹· · ·v^m₁₀¹⁰} be a finite multiset. Ifu=α1v1+

· · ·+α₁₀v₁₀ is such that α₁+· · ·+α₁₀= 0 and α_j = 0 whenever m_j = 0, then one has a differentiation formula

(5) DuQ[K] = (|K| −3)

10

X

j=1

αjQ[K\v_j].

Knot insertion. If y =β₁v₁+· · ·+β₁₀v₁₀ is such that β₁+· · ·+β₁₀ = 1 and β_j = 0 wheneverm_j = 0, then one has a knot insertion formula

(6) Q[K] =

10

X

j=1

β_jQ[(Kty)\v_j].

For instance, if v1,v2 ∈K, then, since v4 = ¹₂v1+¹₂v2, Q[K] = 1

2Q[(Ktv₄)\v₁] +1

2Q[(Ktv₄)\v₂], and for example

1 4

1 11 = 1

2 ⁴

1 21 + 1

2 ^{1 3} 1 21 .

Restriction to an edge. Let e= [v_i,v_k] be an edge of with midpoint v_j and let ϕ_ik(t) := (1−t)vi+tv_k. By induction on|K|,

(7) Q[K]◦ϕ_ik(t) =

( 0 ifm_i+m_j+m_k<|K| −1,

area( )

area([K])B(t) ifm_i+m_j+m_k=|K| −1, where B is the univariate B-spline with knot multiset {0^mi0.5^mj1^mk}.

We say that Q[K] reduces to a B-spline on the boundary when B is one of the consecutive univariate quintic B-splines B₁⁵, . . . , B₈⁵ on the open knot multiset {0⁶0.5²1⁶}; see Table 5. Similarly a basis B = {S₁, . . . , S₃₉} of S₅³( ) reduces to a B-spline basis on the boundary when

S₁◦ϕ_ik, . . . , S₃₉◦ϕ_ik ={ B₁⁵1

· · · B₈⁵1

0³¹}, 1≤i < k≤3, as multisets. This scaling of B ensures simple C⁰ conditions for connecting two adjacent patches expressed in terms of B.

Symmetries. The dihedral group S3 of the equilateral triangle consists of the identity, two rotations and three reflections, i.e.,

1 2

3

2 3

1

3 1

2

1 3

2

3 2

1

2 1

3

The affine bijection sending v_k to cos(2πk/3),sin(2πk/3)

, fork= 1,2,3, maps to the 12-split of an equilateral triangle. Through this correspon- dence, the dihedral group permutes the vertices v₁, . . . ,v₁₀ of . Every such permutation σ induces a bijection Q[v^m₁¹· · ·v^m₁₀¹⁰] 7−→ Q[σ(v₁)^m1· · · σ(v10)^m10] on the set of all simplex splines on . For any set Bof simplex splines, we write

[B]_S3 :={Q[σ(K)] : Q[K]∈ B, σ∈S3}

for the S3 equivalence class ofB, i.e., theset of simplex splines related toB by a symmetry in S₃. For example,

6 11

S3

=

6 11 , ₁¹₆, ¹⁶¹

,

5 21

S3

=

5 21 , _{5 1}² , ₂¹₅, ₁²₅, ¹⁵² , ²⁵¹

, and (4) takes the compact form

3 11 , ₂¹_{1 1}, ₁¹₁¹₁

S3

.

We say that B is S₃-invariant whenever [B]_S3 =B.

4.2. C³ quintic simplex splines on the 12-split. Any simplex spline Q[K] of degree d = 5 on is specified by a multiset K = {v^m₁¹· · ·v^m₁₀¹⁰} satisfying

(8) m1+m2+· · ·+m10=d+ 3 = 8.

Lemma 1. Suppose a quintic simplex spline Q[K] on is of class C³. Then

(9) m7 =m8=m9 =m10= 0, and

(10) m1+m5, m3+m4, m2+m6, m4+m6, m4+m5, m5+m6 ≤3, whenever both multiplicities are nonzero.

Proof. In the 12-split certain knot lines do not appear, leading to the con- ditions

(11)

m₁m₈=m₁m₉ =m₈m₉ = 0, m2m7=m2m9 =m7m9 = 0, m3m7=m3m8 =m7m8 = 0.

To achieve C³ smoothness over the remaining knot lines, the sum of the multiplicities along each line in the 12-split has to be at most three,

(12)

m₁+m₅+m₇+m₁₀≤3, m₄+m₆+m₇ ≤3, m₂+m₆+m₈+m₁₀≤3, m₄+m₅+m₈ ≤3, m3+m4+m9+m10≤3, m5+m6+m9 ≤3,

whenever the line contains at least two knots with positive multiplicities.

Suppose m₇ ≥1. Then (11) implies m₂ =m₃ = m₈ = m₉ = 0, and by (8) and (12), 8 = (m1+m5+m7+m10) + (m4+m6)≤3 + 2, which is a contradiction. It follows that m7 (and similarly m8 and m9) must be equal to zero. Moreover, by (8) and (12), 8 = (m₁+m₅+m₇+m₁₀) + (m₂+m₆+ m8+m10) + (m3+m4+m9+m10)−2m10≤9−2m10. Therefore m10= 0 and (9) holds, and (10) follows immediately from (12).

In addition we demand thatQ[K] reduces to a B-spline on the boundary.

By (7), if m_i+m_j+m_k <7, thenQ[K]|_e = 0 and this condition is satisfied.

The remaining case mi+mj+mk= 7 yields the conditions

(13)

not(m1+m4+m2 = 7 and m4 ≥3),

not(m₁+m₄+m₂ = 7 and m₁ ≥1 and m₂ ≥1 andm₄6= 2), not(m₂+m₅+m₃ = 7 and m₅ ≥3),

not(m₂+m₅+m₃ = 7 and m₂ ≥1 and m₃ ≥1 andm₅6= 2), not(m1+m6+m3 = 7 and m6 ≥3),

not(m1+m6+m3 = 7 and m1 ≥1 and m3 ≥1 andm66= 2).

Adding the entries in the second column of (12) and using (9) gives 2m4+ 2m₅+ 2m₆≤9, implying, by (8),

(14) m₁+m₂+m₃ ≥4, m₄+m₅+m₆≤4.

Theorem 3. With one representative for each S3 equivalence class, Table 2 is an exhaustive list of the C³ quintic simplex splines on that reduce to a B-spline on the boundary.

Q^a Q^b Q^c Q^d Q^e Q^f Q^g Q^h Qⁱ Q^j 6 11

5 21 5

1

2 421 1

412 1

312 2

21212

4 2

2

3 3

2

4 1

2 1

Q^k Q^l Q^m Qⁿ Q^o Q^p Q^q Q^r Q^s Q^t

3 2

2

1 1 4

1 11

1 3

2 11

2 2

2 11

2 2

1 11

1 4 1

1

2 3 2

1

2 1 3

1 21

2 2

1 21

1 2

1 21 1 Table 2. The C³ quintic simplex splines on , one rep- resentative for each S3 equivalence class, that reduce to a B-spline on the boundary.

Proof. Recall from (9) that m7 = m8 = m9 = m10 = 0. We distinguish cases according to the support [K] of Q[K], up to a symmetry in S3.

Case 0, no corner included,[K] = [v₄,v₅,v₆]: By (8),m₄+m₅+m₆= 8, so that the sum of two of these multiplicities will be at least 5, contradict- ing (10). Therefore this case does not occur.

Case 1a, 1 corner included, [K] = [v1,v4,v6]: For a positive support m₁, m₄, m₆ ≥1, and sincem₄+m₆ ≤3 by (10), we obtain

[Q^a]S3 =

6 11

S3

, [Q^b]S3 =

5 21

S3

.

Case 1b, 1 corner included, [K] = [v₁,v₄,v₅,v₆]: By (8) and (14) one has m1 = 8−m4−m5−m6 ≥ 4, contradicting m1+m5 ≤ 3 from (10).

Therefore this case does not occur.

Case 2a, 2 corners included, [K] = [v1,v2,v6]: If m4 ∈ {0,1}, then m₆≥2 by the second line in (13), and sincem₆≤3−m₂ by (10), it follows that m2 = 1 and m6 = 2. We obtain

[Q^c]_S

3 =

52 1

S3

, [Q^d]_S3 =

421 1

S3

.

If m₄ ≥2, then, since m₄+m₆ ≤ 3 by (10), one hasm₄ = 2 and m₆ = 1.

Since m2≤3−m6 = 2 by (10), we obtain [Q^e]_S3 =

412 1

S3

, [Q^f]_S3 =

312 2

S3

.

Case 2b, 2 corners included, [K] = [v₁,v₂,v₅,v₆]: Since m₅+m₆ ≤3 by (10), one has 1 ≤ m₅, m₆ ≤ 2. Suppose m₆ = 2. Then m₂ =m₅ = 1 and m4≤1 by (10), andm1+m5= 8−m2−m4−m6 ≥4, contradicting (10).

We conclude thatm₆= 1 and similarly thatm₅= 1. Then m₁+m₂+m₄ = 8−m5−m6 = 6, and since (10) implies m1, m2, m4 ≤2, one obtains

[Q^g]_S3 =

21212

S3

.

Case 3, 3 corners included, [K] = [v₁,v₂,v₃]: We distinguish cases for (m4, m5, m6), withm4 ≥m5 ≥m6, first bym4, and then bym4+m5+m6, which is at most 4 by (14).

(0,0,0) One has m₁, m₂, m₃≥2 by (13), and we obtain [Q^h]_S3 =

4 2

2

S3

, [Qⁱ]_S3 =

3 3

2

S3

.

(1,0,0) One has 2≤m₃ ≤3−m₄ by (13) and (10), yielding [Q^j]S3 =

4 1

2 1

S3

, [Q^k]S3 =

3 2

2 1

S3

. (1,1,0) One has m1, m3≤2 by (10), yielding

[Q^l]S3 =

1 4

1 11

S3

, [Q^m]S3 =

1 3

2 11

S3

, [Qⁿ]S3 =

2 2

2 11

S3

. (1,1,1) One has m1, m2, m3≤2 by (10), and we obtain

[Q^o]_S3 =

2 2

1 11

1

S3

.

(2,0,0) One has m₃ = 1 by (10), yielding [Q^p]_S3 =

4 1

1 2

S3

, [Q^q]_S3 =

3 2

1 2

S3

. (2,1,0) One has m₃ = 1 andm₁≤2 by (10), yielding

[Q^r]S3 =

1 3

1 21

S3

, [Q^s]S3 =

2 2

1 21

S3

.

(2,1,1) One has m3 = 1 andm1, m2≤2 by (10), yielding [Q^t]_S3 =

1 2

1 21 1

S3

.

5. Simplex spline bases for S₅³

Let Λ as in (1) be a basis of the dual space of S₅³( ). In this section we describe a recipe for determining the S3-invariant simplex spline bases that reduce to a B-spline basis on the boundary, having a positive partition of unity and a Marsden identity with only real linear factors.

5.1. Potential bases. Of the C³ quintic simplex splines in Table 2, only

Q^a Q^b Q^c Q^e Q^f Q^p Q^q

6 11 , _{5 2}¹ , ₅ ¹₂ , ₄¹_{2 1}, ₃¹_{2 2}, ₄ ¹_{2 1}, ₃ ¹_{2 2}

are nonzero on [v1,v2]. Any S3-invariant basis reducing to a B-spline basis on the boundary should therefore contain theS₃ equivalence class of one of the eight possible combinations in the Cartesian product

6 11

×

5 21 , ₅ ¹₂

×

412 1, ₄ ¹_{2 1}

×

312 2, ₃ ¹_{2 2}

. This determines 3 + 6 + 6 + 6 = 21 of the 39 basis elements.

Of the 13 remaining S3 equivalence classes in Table 2 of simplex splines that are zero on [v1,v2], there are 6 of size 3,

Q^g Q^h Qⁱ Q^l Qⁿ Q^o h

21212, ₄ ² ₂, ₃ ² ₃, ₁ ¹₁1₄, ₂ ²₁1₂, ₂1¹₁1₂i

S3

and 7 of size 6,

Q^d Q^j Q^k Q^m Q^r Q^s Q^t

h

421 1, ₄ ²_{1 1}, ₃ ²_{1 2}, ₁ ²₁¹₃, ₁ ¹₂¹₃, ₂ ¹₂¹₂, ₁¹¹₂¹₂i

S3

.

For each of the above 8 choices of 21 simplex splines that are nontrivial on the boundary of , we complete the basis by addingS3 equivalence classes with together 18 elements, resulting in

8 7

3

+ 7

2 6

2

+ 7

1 6

4

+ 6

6

= 3648

potential S3-invariant simplex spline bases B of S₅³ that reduce to a B- spline basis on the boundary. One selects the linearly independent sets B by checking that the collocation-like matrix {λ(Q)}_{Q∈B,λ∈Λ} has full rank.

5.2. Positive partition of unity. Let B = {Q₁, . . . , Q39} be an ordered simplex spline basis ofS₅³( ) that reduces to a B-spline basis on the bound- ary. We desire that B has a positive partition of unity, i.e., that there exist weights w1, . . . , w39>0 for which

(15)

39

X

i=1

wiQi(x) = 1

holds identically. Applying the functionals in Λ ={λ₁, . . . , λ39} yields

39

X

i=1

wiλj(Qi) =λj(1), j= 1, . . . ,39,

which has a unique solution (w₁, . . . , w₃₉)∈R³⁹ by the linear independence of B and Λ. One then checks whether the weightswi are all positive.

5.3. Marsden identity. More generally, we would likeBto satisfy aMars- den identity

(16) (1 +x^Ty)⁵=

39

X

i=1

wiQi(x)ψi(y), x∈ , y∈R², for certain dual polynomials ψ_i and dual points p^∗_i,r of the form

ψi(y) :=

5

Y

r=1

(1+p^∗T_i,ry), p^∗_i,r ∈R², r= 1, . . . ,5, i= 1, . . . ,39.

In particular one recovers the partition of unity by setting y= 0. Similarly one can generate all other polynomials of degree at most 5. For instance, differentiating (16) with respect to y and evaluating aty= 0 yields

(17) x=

39

X

i=1

ξ_iwiQi(x), ξ_i := 1 5∇ψ_i

y=0

= p^∗_i,1+· · ·+p^∗_i,5

5 ,

where ξ_i is the domain point associated to Qi.

While the Marsden identity (16) is the form commonly encountered in the literature [9, 10], we instead present a barycentric form that is independent of the vertices of the macrotriangle.

Theorem 4 (Barycentric Marsden identity). Let βj =βj(x), j = 1,2,3, be the barycentric coordinates of x∈R² with respect to = [v1,v2,v3]. Then (16) is equivalent to

(18) (β₁c₁+β₂c₂+β₃c₃)⁵ =

39

X

i=1

w_iQ_i(β₁v₁+β₂v₂+β₃v₃)Ψ_i(c₁, c₂, c₃), where x∈ , c1, c2, c3∈R, and, for i= 1, . . . ,39,

Ψi(c1, c2, c3) :=

5

Y

r=1

β1(p^∗_i,r)c1+β2(p^∗_i,r)c2+β3(p^∗_i,r)c3

.

Proof. Let Xj := (1,vj)^T ∈ R³, j = 1,2,3. Since {X₁,X2,X3} is linearly independent, there are Yi ∈ R³, i= 1,2,3, such that X^T_jYi = δij. Given x ∈ and c₁, c₂, c₃ ∈ R we define X,Y ∈ R³ by X := (1,x)^T and Y = (1,y)^T:=c₁Y₁+c₂Y₂+c₃Y₃. SinceX=β₁(x)X₁+β₂(x)X₂+β₃(x)X₃,

1 +x^Ty=X^TY=β₁(x)c₁+β₂(x)c₂+β₃(x)c₃, and in particular with x=p^∗_i,r and P^∗_i,r := (1,p^∗_i,r)^T,

1 +p^∗T_i,ry=P^∗T_i,rY =β₁(p^∗_i,r)c₁+β₂(p^∗_i,r)c₂+β₃(p^∗_i,r)c₃.

It follows that (18) is equivalent to (16).

For a compact representation of the dual polynomials, we introduce the shorthands (compare Figure 2a)

(19)

c4:= c1+c2

2 , c5:= c2+c3

2 , c6:= c1+c3

2 , c7:= c4+c6

2 , c8:= c4+c5

2 , c9:= c5+c6

2 , c10:= c1+c2+c3

3 .

Example 3. The barycentric Marsden identity for the quadratic S-basis in [5, Theorem 3.1] (cf. Example 1) is

(c1β1+c2β2+c3β3)² =

c12 1

1 +c21 2

1 +c31 1 2 2

= +1

4c²_{13 1}¹ +3

4c₄c₁₀₁¹₁¹₁ +1

2c₁c₄₂¹_{1 1} +1

2c₂c_{41 1}¹₂ +1

4c²_{2 1}¹₃ +3

4c5c10 1 1 11

1 +1

2c2c5 2 1 11 +1

2c3c5 1 21 1 +1

4c²₃ ¹³¹ +3

4c₆c₁₀₁¹¹₁¹ +1

2c₃c_{6 1}¹²¹ +1

2c₁c₆₂¹¹₁ , where the weights follow from (3).

For every basis B={Q₁, . . . , Q₃₉}with a positive partition of unity (15), we apply the functionals in Λ to (16) with respect to x, which gives







λ1(Q1) · · · λ1(Q39) ... . .. ... λ39(Q1) · · · λ39(Q39)













w1ψ1(y) ... w39ψ39(y)





=







λ1 (1 +x^Ty)⁵ ...

λ₃₉ (1 +x^Ty)⁵





. This system has a unique solution in the module R[y]³⁹, with each compo- nent wiψi(y) a polynomial of degree at most 5 by Cramer’s rule. The basis Bhas a Marsden identity (16) if and only if these polynomials split into real linear factors.

6. Main results

Some of the remaining computations are too large to carry out by hand. We have therefore implemented the above computations in Sage[22]. From the website [18] of the second author, the resulting worksheet can be downloaded and tried out online in SageMathCloud.

6.1. Six bases for S₅³. Checking linear independence for each of the 3648 potential bases, we discover that there are 1024 S3-invariant simplex spline bases that reduce to a B-spline basis on the boundary. There are 243 such bases with a nonnegative partition of unity, of which there are 47 bases with a positive partition of unity. Of these there are 9 with all domain points inside the macrotriangle, of which there are 7 with precisely 8 domain points on each boundary edge. Of these there are 6 bases B_a,B_b,B_c,B_d,B_e,B_f for which each dual polynomial only has real linear factors. These bases, together with their weightsw?, dual polynomials Ψ?, and domain pointsξ_?, are listed in Table 3. For instance, the highlighted rows in the table yield the set

(B_c) 1

4 ^{6 11} ,1

4 ^{5 21} ,1

2 ^{412 1},1

2 ^{312 2},3

4 ²¹²¹², ₁ ¹₁¹₄,1 2^{1 3}

1 21 ,3

4 ^{1 2} 1 21

1

S3

. We summarize these results in the following theorem.

Theorem 5. The setsB=B_a,B_b,B_c,B_d,B_e,B_f are the only sets satisfying:

(1) B is a basis of S₅³( ) consisting of simplex splines.

(2) B is S3-invariant.

(3) B reduces to a B-spline basis on the boundary.

(4) B has a positive partition of unity and a Marsden identity (18), for which the dual polynomials have only real linear factors; see Table 3.

(5) Bhas all its domain points inside the macrotriangle , with precisely 8 domain points on each edge of .

Remark 6. Table 3 shows that, for the simplex splines Q^g, Q^s, Q^t

S3 =

21212, ₂ ¹₂1₂, ₁1¹₂1₂

S3

,

the weights and dual polynomials depend on the entire basis, and they cannot be determined directly from the corresponding knot multiset.