Construction of C 2 cubic splines on arbitrary triangulations

In this paper, we address the problem of constructing C 2 cubic spline functions on a given arbitrary triangulation T . To this end, we endow every triangle of T with a Wang–Shi macro-structure. The C 2 cubic space on such a reﬁned triangulation has a stable dimension and optimal approximation power. Moreover, any spline function in such space can be locally built on each of the macro-triangles independently via Hermite interpolation. We provide a simplex spline basis for the space of C 2 cubics deﬁned on a single macro-triangle which behaves like a Bernstein/B-spline basis over the triangle. The basis functions inherit recurrence relations and diﬀerentiation formulas from the simplex spline construction, they form a nonnegative partition of unity, they admit simple conditions for C 2 joins across the edges of neighboring triangles, and they enjoy a Marsden-like identity. Also, there is a single control net to facilitate control and early visualization of a spline function over the macro-triangle. Thanks to these properties, the complex geometry of the Wang–Shi macro-structure is transparent to the user. Stable global bases for the full space of C 2 cubics on the Wang–Shi reﬁned triangulation T are deduced from the local simplex spline basis by extending the concept of minimal determining sets.


Introduction
Piecewise polynomial spaces defined over polygonal partitions, usually triangulations, have applications in several branches of the sciences including geometric modeling, signal processing, data analysis, visualization, and numerical simulation; we refer the reader to [10,13,26] and Section 5.2 for some examples.For many of these applications, a smooth join between the different pieces is beneficial or even required; C 2 smoothness is often preferred.Such spaces are commonly referred to as (bivariate) spline spaces.According to [26, page 197], in general we would like to work with low degree splines: they involve fewer coefficients, and have less tendency to oscillate.
An indispensable feature for a spline space to be useful in practice is having a stable dimension that only depends on the degree (d), the order of smoothness (r), and combinatorial -or other easy to check -properties of the partition (T ).When T is a triangulation, the dimension can be expressed in terms of the above quantities for spline spaces with d ≥ 3r + 2; see [20] and [26,Chapter 9].On the other hand, instability in the dimension has been illustrated for d = 2r in [14].We refer the reader to [52] for recent results on the dimension of spline spaces on triangulations with nonuniform degrees.Similar results are known for spline spaces over general rectilinear partitions; see [7,31] and references therein.
Spline spaces with too low degree compared to the smoothness are also exposed to several other shortcomings.In particular, they might lack optimal approximation power, a property strongly related to the possibility of constructing stable bases with local support for the considered spaces [26].In this perspective, the bound d ≥ 3r+2 plays again an important role in identifying the spline spaces with optimal approximation power on a given triangulation [26,Chapter 10].Furthermore, the possibility of constructing any function of the spline space locally on each of the elements of T is often seen as a desirable, if not imperative, property for practical purposes.On a triangulation, a degree d ≥ 4r + 1 is necessary to admit such a local construction [8,26,56].
The above lower bounds on the degree can be alleviated by considering so-called macro-elements, where the partition T is further refined in a specific manner (often referred to as splits).In case T is a triangulation, the most famous examples are the Clough-Tocher (CT) split [8,9,26,40] and the Powell-Sabin (PS) 6 and 12 splits [1,26,36,40,43].They subdivide each triangle of T into 3, 6, and 12 subtriangles, respectively.To achieve global C 2 smoothness, polynomial pieces of at least degree d = 7 are necessary for the CT split, while at least degree d = 5 is required for both PS splits of T .All these spline spaces have a stable dimension and possess optimal approximation power [24][25][26].Other common macro-elements also require at least degree d = 5 to realize C 2 splines with the above properties on a refined partition only consisting of triangles [26,Section 7.7].
The Bernstein polynomial basis is the most common tool for the construction and analysis of splines on a given triangulation T [26], as it helps in localizing imposition of smoothness conditions across edges of (the refinement of) T .Interesting alternatives have been developed for CT and PS splits in [11,[28][29][30], where a simplex spline basis for the local spline space over a triangle of T has been considered.Such a basis behaves like a Bernstein polynomial basis for imposing smoothness across edges of T and like a B-spline basis internal to each triangle of T .Neither the Bernstein polynomial basis nor the simplex spline basis provide a global basis for the full spline space on (the refinement of) T .To achieve a global basis, one may apply the general framework of minimal determining sets via the local Bernstein basis; see [26].Global B-spline bases have been constructed for C 1 PS spline spaces on triangulations [12,18,19,49], for PS spline spaces with higher smoothness [17,45,47], and for CT spline spaces [46].
While in the univariate case C 2 cubics are probably the best known and most used splines, the above discussion shows that dealing with C 2 cubics in the bivariate setting is an arduous task.In this paper, we address the problem of building and handling C 2 cubic splines on a suitable refinement of any given triangulation T .Our wish list for the spline space consists of stable dimension, optimal approximation power, and local construction on any (refined) triangle of T .Moreover, we want a practical construction of a stable global basis for the space at our disposal.
Despite the high smoothness and the minimum gap between degree and smoothness, a C 2 cubic space can be obtained by splitting any triangle ∆ in T according to a degree-dependent scheme introduced by Wang and Shi [54].Contrarily to the well-known splits mentioned above, the family of Wang-Shi (WS) splits generates a very large number of polygonal pieces in each ∆; for cubics we get a set of 75 polygons which includes triangles, quadrilaterals, and pentagons.In practice, this complex geometry hampers a piecewise treatment -in terms of a local polynomial basisof spline functions on WS splits and discourages the use of such an interesting space.
To overcome this issue, we propose a simplex spline basis for the local space of C 2 cubics on the (cubic) WS split of any ∆ in T .The basis functions enjoy the following properties: • they form a nonnegative partition of unity; • they inherit recurrence relations and differentiation formulas from the simplex spline structure; • for each of them, the restriction to a boundary edge of ∆ reduces to a classical C 2 cubic univariate B-spline; • they admit simple conditions for C 2 joins to neighboring triangles in T ; • cubic polynomials can be represented through a Marsden-like identity; • they lead to well-conditioned collocation matrices for Lagrange and Hermite interpolation using certain sites; • a control net can be formed that mimics the shape of the spline function and exhibits distance O(h 2 ) to any one of its control points from its surface, where h is the length of the longest edge.
Thanks to the characteristics of the simplex spline basis, one can avoid to consider separate polynomial representations on each of the polygonal subelements of ∆.Instead, there is a single control net to facilitate control and early visualization of a spline function over each element ∆ in T .This makes that the complex geometry of the WS split is transparent to the user.However, the simplex spline basis is a local basis and does not provide a global basis for the full space of C 2 cubics on the (cubic) WS refinement of T .To this end, we extend the concept of minimal determining sets and use the simplex spline basis as a stepping stone to the construction of a stable global basis for the full space.The remainder of this paper is divided into four sections.In Section 2, we summarize the definition and some properties of simplex splines and describe the family of WS splits.In Section 3, we present a local simplex spline basis for the refinement of a single triangle and discuss some of its properties.To simplify some computations, an alternative basis is also provided.Smoothness conditions across the edges of the given triangulation and stable global bases for the C 2 cubic space on the (cubic) WS refinement of any triangulation are considered in Section 4. Section 5 collects some concluding remarks about implementation aspects, possible application areas, and a higher-order extension of the basis.Finally, the appendix aggregates data related to the presented simplex spline bases that might be useful for practical computations.
Throughout the paper, we use small boldface letters for vectors and capital boldface letters for matrices.Calligraphic letters like B indicate sets, and we write #B for the cardinality of B. Function spaces are denoted by symbols like S. In particular, P d stands for the space of bivariate polynomials with real coefficients of total degree ≤ d.The partial derivatives in x and y are denoted by D x and D y , respectively.Given a vector u, the associated directional derivative is denoted by D u .The directional derivative in the direction of the vector from point p 1 to p 2 is denoted by D p 1 p 2 .

Preliminaries
This section contains some preliminary material about simplex splines and the splits of interest in the rest of the paper.

A short summary of simplex splines
For e ∈ N, d ∈ N 0 , let n := d + e and Ξ := {ξ 1 , . . ., ξ n+1 } be a sequence of possibly repeated points in R e called knots.The multiplicity of a knot is the number of times it occurs in the sequence.Let • denote the convex hull of a sequence of points.For the sake of simplicity, we assume Ξ is nondegenerate, i.e., vol e ( Ξ ) > 0. Let σ = ξ 1 , . . ., ξ n+1 be any simplex in R n with vol n (σ) > 0, whose projection π : R n → R e onto the first e coordinates satisfies π(ξ i ) = ξ i for i = 1, . . ., n + 1.
The simplex spline M Ξ can be defined geometrically by For d = 0 we have and the value of M Ξ on the boundary of Ξ has to be dealt with separately.For properties of M Ξ and proofs, we refer the reader to, e.g., [32,37].Here, we mention: • Knot dependence: M Ξ only depends on Ξ; in particular, it is independent of the choice of σ and the ordering of the knots.
• Nonnegativity: M Ξ is a nonnegative piecewise polynomial of total degree d.
• Differentiation formula (A-recurrence): For any u ∈ R e and any a 1 , . . ., a d+e+1 such that i a i ξ i = u, i a i = 0, we have • Recurrence relation (B-recurrence): For any x ∈ R e and any b 1 , . . ., b d+e+1 such that • Knot insertion formula (C-recurrence): For any y ∈ R e and any c 1 , . . ., c d+e+1 such that i c i ξ i = y, i c i = 1, we have .
If e = 1 then M Ξ is the univariate B-spline of degree d with knots Ξ, normalized to have its integral equal to one.In the bivariate case, e = 2, the lines in the complete graph of Ξ are called knot lines.They provide a partition of Ξ into polygonal elements.The simplex spline M Ξ is a polynomial of degree d = #Ξ − 3 in each region of this partition, and across a knot line where µ is the number of knots on that knot line, including multiplicities.

The Wang-Shi splits
Given three noncollinear points p 1 , p 2 , p 3 in R 2 , the triangle ∆ := p 1 , p 2 , p 3 with vertices p 1 , p 2 , p 3 will serve as our macro-triangle.Given a degree d ∈ N, we divide each edge of ∆ into d equal segments, respectively, resulting into 3d boundary points.Then, we refine ∆ into a number of subelements delineated by the complete graph connecting those boundary points.This is called the WS d split of ∆ as it was originally proposed by Wang and Shi [54].We denote by ∆ WS d the obtained mesh structure, and by P d the set of polygons in ∆ WS d .All the possible intersections of the various lines connecting the boundary points are called vertices of ∆ WS d .In particular, the boundary points are vertices of ∆ WS d .The cases d = 2, 3, 4 are shown in Figure 1.For d = 2 we obtain the well-known PS-12 split [36], while for d = 1 we have P 1 = {∆}.Note that for d > 2 not all elements of P d are triangles.We consider the space When the degree increases, the complexity of the mesh grows quickly.There are • 3d boundary points and 3d(d − 1) interior lines in the complete graph; • the maximum number of lines intersecting at an interior vertex is 3, 3, 4, 5, 6, 7, 6 for d = 2, 3, . . ., 8; • the number of vertices of ∆ WS d is 10, 58, 178, 558, 1255, 2532, 4786 for d = 2, 3, . . ., 8.
The dimension of S d−1 d (∆ WS d ) can be computed using the general dimension formula for spline spaces over cross-cut partitions from [7,Theorem 3.1].A partition T c of a domain Ω is called a cross-cut partition if it is obtained by drawing lines across Ω.Let S r d (T c ) be the space of functions in C r (Ω) which belong to P d when restricted to any polygon of T c .

Theorem 1.
Let Ω be a simply connected domain in R 2 .Let T c be a cross-cut partition of Ω, with m cross-cuts, n interior vertices v 1 , . . ., v n , and m k cross-cuts intersecting at v k , k = 1, . . ., n.Then, the dimension of the spline space where As usual, x denotes the largest integer smaller than or equal to x, and (x) + := max{x, 0}.Proof.We make use of the dimension formula (2) in Theorem 1.In our case we have r = d − 1, and it is easy to check that ς(l) = 0, l = 1, . . ., d + 1.
Since at most d + 1 lines cross at each interior vertex, it immediately follows from (2) that which completes the proof.
3 Simplex spline bases for S 2 3 (∆ WS 3 ) In this section, we focus on the case d = 3, provide two (scaled) simplex spline bases for the space S 2 3 (∆ WS 3 ) in (1), and prove some properties of these bases.With a slight abuse of notation we also refer to the corresponding basis functions as simplex splines.

A simplex spline basis
For a given triangle ∆ = p 1 , p 2 , p 3 , the WS 3 split is shown in the middle plot of Figure 1.From Theorem 2 we know that the dimension of S 2 3 (∆ WS 3 ) is 28.In order to construct a basis for this space, we first specify nine points along the boundary of the triangle (see Figure 2): the three vertices p 1 , p 2 , p 3 and Figure 3: Sequences of knots for a set of simplex spline basis functions for S 2 3 (∆ WS 3 ).Each black disc shows the position of a knot and the number inside indicates its multiplicity.
Note that these points are part of the WS 3 split.We then consider the cubic simplex splines M 1 , . . ., M 28 as schematically illustrated in Figure 3, where each simplex spline has six (including multiplicity) knots chosen among the nine points above.For instance, M 4 is defined by the sequence {ξ 1 , . . ., ξ 6 } = {p 1 , p 1 , p 1 , p 3,1 , p 3,2 , p 2,1 }.Each of them can be computed using the B-recurrence relation.We define the following set of 28 (scaled) simplex splines: where, denoting by |∆| the area of ∆, the scaling factors are given by Theorem 3. The simplex splines {B 1 , . . ., B 28 } in (4) form a nonnegative partition of unity basis for the space S 2 3 (∆ WS 3 ).Proof.Let B be one of the functions B i .We first prove that B ∈ S 2 3 (∆ WS 3 ).Since B has six knots, it is a piecewise cubic polynomial.Moreover, the knots of B are a subset of the knots shown in Figure 2. Thus, the knot lines of B are a subset of the knot lines in the complete graph; see Figure 1.Since each interior knot line contains exactly two knots, B has C 2 smoothness according to the smoothness property of simplex splines.It follows that B ∈ S 2 3 (∆ WS 3 ).We now consider linear independence.Using the recurrence relation and differentiation formula for simplex splines and the scaling factors w, we compute values and derivatives of B corresponding to the following 28 operators: ρ 1 , . . ., ρ 18 are related to the vertices, (5) ρ 19 , . . ., ρ 27 are related to the edges, where and the final ρ 28 is related to the triangle, The computed values are shown in Table 2 in the appendix.Since the matrix [ρ j (B i )] ∈ R 28×28 is block upper triangular with nonsingular 2 × 2 blocks, linear independence of the set of functions {B 1 , . . ., B 28 } follows.From Theorem 2 we know that the dimension of S 2 3 (∆ WS 3 ) is 28, and thus the 28 linearly independent functions in (4) form a basis of this space.At the same time, we may conclude linear independence of the set of operators {ρ 1 , . . ., ρ 28 } defined on S 2 3 (∆ WS 3 ).Simplex splines are nonnegative, so it only remains to prove that the functions in ( 4) sum up to one on ∆.An inspection of Table 2 shows that By linear independence of the operators ρ j , 28 i=1 B i must be equal to the unity function which belongs to S 2 3 (∆ WS 3 ).
From the proof of Theorem 3 it follows that we can formulate a Hermite interpolation problem to characterize any spline in S 2 3 (∆ WS 3 ).Corollary 4. For given data f k,α,β , g k , g k,l , and h 0 , there exists a unique spline s ∈ S 2 3 (∆ WS 3 ) such that where n k is the normal direction of the edge opposite to vertex p k , and the points p k,l , q k , and q are defined in (3), (7), and (8), respectively.
The Hermite degrees of freedom specified in Corollary 4 are schematically visualized in Figure 4 using graphical symbols that are common in finite element literature; see, e.g., [8].

Domain points and condition number
We now associate a special point in ∆ with each basis function B i in (4), that plays an important role in geometric modeling.We solve the system ρ j ( 28 for the two functions f 1 (x, y) := x and f 2 (x, y) := y.The points are called the domain points of the basis (4).Together with the partition of unity the domain points provide an explicit representation of any affine function with respect to the basis (4).The barycentric coordinates with respect to the triangle ∆ of the domain points ( 9) are given by b (10) ).Left: the domain points (10).Right: the domain points (16).Figure 6: A simplex spline surface and its control net for the domain points (10).
These points are visualized in Figure 5 (left).When representing a spline s ∈ S 2 3 (∆ WS 3 ) in the basis (4), it is common to organize the coefficients in terms of control points (b * i , b i ), i = 1, . . ., 28.There are several possibilities to connect these points into a control net.Such a control net forms a caricatural approximation for the graph of the function that is useful for geometric modeling.A viable option for connecting these points is shown in Figure 5 (left); the configuration consists of a small number of regions but both triangles and quadrilaterals are involved.As shown in Section 4.1, this choice allows for a geometric interpretation of C 1 smoothness conditions analogous to the classical Bernstein representation for polynomial triangular patches.An example spline and its corresponding control net is illustrated in Figure 6.
In many applications it is of interest to have a bound on the condition number of the basis we are dealing with.We consider the infinity norm, and we look for constants The condition number of the basis is then defined by The condition number of the simplex spline basis (4) is bounded by Proof.Since the simplex spline basis (4) forms a nonnegative partition of unity, it is clear that K + ∞ = 1 satisfies (12).Let A be the matrix for any unisolvent Lagrange interpolation problem with respect to the basis (4).Then, Considering interpolation at the domain points (10), a direct computation gives A −1 ∞ < 37.This implies that K − ∞ = 1/37 satisfies (12).
Note that the bound on the condition number in Proposition 5 is independent of the shape of the triangle ∆.From the proof we also deduce that the condition number of B in (4) can be computed as κ ∞ (B) = inf{1/K − ∞ : K − ∞ satisfies (12)}.By means of this number we can easily obtain the following distance result.Proposition 6.Let s ∈ S 2 3 (∆ WS 3 ) be given as in (11).Then, where h is the length of the largest edge of ∆.
Proof.Let Π i ∈ P 1 be the linear Taylor approximation to s at the domain point b Then, using the definition of domain points, we have so that from (12) we obtain Furthermore, since s ∈ C 2 (∆), Taylor approximation error analysis tells us that which completes the proof.
Proposition 6 implies that the control points of the spline s in (11) converge like O(h 2 ) to the graph of s.

A Marsden-like identity
In Section 3.2, we have provided the representation of any affine function with respect to the basis (4).We now extend this result by providing a Marsden-like identity which allows us to represent any cubic polynomial.In the univariate B-spline case, the Marsden identity is given by where B j,d is a normalized B-spline of degree d defined by the knots ξ j , ξ j+1 , . . ., ξ j+d+1 ; see, e.g., [27,Theorem 2].Dividing both sides by z d and setting y := z −1 we obtain a form more amenable to multivariate generalization The functions ψ j,d are polynomials of degree d and are called dual polynomials.The following result is obtained by a direct computation.
Theorem 7. We have where the dual polynomials ψ i , i = 1, . . ., 28 are defined as follows.Recalling the points in  It is remarkable that the dual polynomials ψ i , i = 1, . . ., 27, can be written as products of three linear polynomials and mimic the classical univariate Marsden identity (13).It is worth noting that these functions are the polar forms (or blossoms) of (1+y T x) 3 evaluated at appropriate points [38].Similarly, polar forms were used in [33] for the representation of polynomials in terms of simplex splines whose knots are in generic position.This is not the case for the knots in Figures 3 and 7.

An alternative simplex spline basis
An alternative basis for the space S 2 3 (∆ WS 3 ) is provided by the simplex splines identified by the knot sequences in Figure 7 and scaled to form a partition of unity.We denote this basis by It can be checked that ( These points are depicted in Figure 5 (right).Note that Theorem 3 in combination with (15) confirms that the functions B 1 , . . ., B 28 are linearly independent and form a partition of unity.On the other hand, they are not all nonnegative.For subsequent use, some Hermite data of these basis functions are collected in the appendix (Tables 2 and 3).Any spline s ∈ S 2 3 (∆ WS 3 ) can be represented in terms of this alternative basis, so Note that b = Cb, where C is the conversion matrix already used to obtain (16).It can be easily checked that C ∞ = C −1 ∞ = 3.Therefore, from ( 12) and ( 17) we immediately deduce and the condition number can be bounded as κ ∞ ( B) < 333.We can also formulate the analogue of Proposition 6 as well as a Marsden-like identity for the basis (14).We omit the details for the sake of brevity.The simplex spline basis (4) forms a convex partition of unity and so it is particularly useful for geometric modeling.On the other hand, as we will show in Section 4.1, the simplex spline basis ( 14) is more suited to handle C 2 smoothness conditions between spline functions on adjacent macrotriangles.Of course, there are many more alternative sets of simplex spline basis functions.One could, for instance, take the 10 cubic Bernstein polynomials defined on ∆ (they are special simplex splines; see [37]) and enrich them with 18 more simplex splines that are linearly independent.

C 2 cubic splines on the WS 3 refinement of a triangulation
In the previous section, we have provided simplex spline bases for the spline space S 2 3 (∆ WS 3 ) of cubic C 2 splines on the WS 3 split of a given triangle ∆.Let T be a triangulation of a polygonal domain Ω and let T WS 3 denote its refinement obtained by taking the WS 3 split of each of its triangles.In this section, we consider the spline space of C 2 cubic splines on T WS 3 , i.e., The unisolvency of the Hermite interpolation problem stated in Corollary 4 implies that the dimension of the space only depends on combinatorial properties of the triangulation, and so it is stable.From the corollary we directly deduce that (see also [54]) where n V , n E , n T are the number of vertices, edges, and triangles of T , respectively.Moreover, any spline function of S 2 3 (T WS 3 ) can be locally constructed on each (macro-)triangle ∆ of T via the Hermite data, and the corresponding spline piece on ∆ can be represented in the form (17). Conversely, any function, which is represented locally in the form (17) on each ∆ of T , is C 2 smooth over each ∆ of T but not necessary C 2 smooth across the edges of T .First we derive conditions on the local spline coefficients to ensure global C 2 smoothness, and then we describe a stable global basis with local support for S 2 3 (T WS 3 ).

Smoothness conditions
Let T be a triangulation of a polygonal domain Ω ⊂ R 2 .We seek conditions on the local spline coefficients in (17) to guarantee C r smoothness across a common edge of two adjacent triangles of T for r = 0, 1, 2.
Theorem 8. Suppose the triangles ∆ L := p 1 , p 2 , p 3 and ∆ R := p 1 , p 2 , p 4 share the common edge with vertices p 1 , p 2 , and let Let {B L i , i = 1, . . ., 28} and {B R i , i = 1, . . ., 28} be the scaled simplex spline basis defined by the knot sequences in Figure 3 on ∆ L and ∆ R , respectively.We assume the numbering of the basis functions in agreement with Figure 3. Let us consider the spline functions We have • s L , s R join C 1 across the common edge if and only if they join C 0 and in addition Proof.Let us first discuss C 0 smoothness.Along the common edge p 1 p 2 the two functions s L and s R are univariate cubic C 2 splines with (interior) knots at the points p 3,1 and p 3,2 .Considering the restriction onto the edge of the basis functions {B L i , i = 1, . . ., 28} and {B R i , i = 1, . . ., 28}, we obtain that the only nonzero elements are Since they are linearly independent, C 0 smoothness is equivalent to agreement of the corresponding coefficients.This proves (21).We now discuss C 1 smoothness across the common edge.It suffices to prove that along the edge p 1 p 2 the functions D q 3 p 3 s L and D q 3 p 3 s R agree.These functions are univariate C 1 quadratic splines with (interior) knots at the points p 3,1 and p 3,2 .Therefore, each of them is uniquely determined by its value and first derivative at the two endpoints of the edge and by the value at the midpoint q 3 .From (20) we obtain and so Then, by employing the C 0 smoothness conditions and the values in Table 2 we get Equating the above expressions results in the first condition of (22).With the same line of arguments we deduce the remaining four conditions.
From the relations in (15) it is clear that the conditions in ( 21) and ( 22) also ensure C 0 and C 1 smoothness, respectively, for local spline representations in the alternative basis (14).
Corollary 9. Consider the same assumptions as in Theorem 8.The C 0 smoothness conditions for the control points can be written as and the C 1 smoothness conditions for the control points can be written as Proof.The statements follow by direct computation from (22) and from the expressions of the domain points in (10).
The C 1 smoothness conditions in Corollary 9 have a nice geometric interpretation.There are five sets of four control points that need to be coplanar.In terms of our control net configuration in Figure 5 (left), that means that the five triangles in both control nets along the common edge must be all pairwise coplanar.This is illustrated in Figure 8.
Proof.Assume s L and s R join C 1 across the common edge p 1 p 2 .To prove C 2 smoothness across the same edge, it suffices to prove that along p 1 p 2 the functions D 2 p 1 p 3 s L and D 2 p 1 p 3 s R agree.Along the edge p 1 p 2 , these functions are univariate C 0 linear splines with (interior) knots at the points p 3,1 and p 3,2 .Therefore, each of them is uniquely determined by its value at the two endpoints of the edge and by the value at the points p 3,1 and p 3,2 .From (20) we know that so that Then, by employing the values in Table 2 we get By equating the above expressions and taking into account the C 1 smoothness conditions, we obtain the first condition of (23).With the same line of arguments, taking into account (15) and the additional values in Table 3, we deduce the remaining three conditions.
Corollary 11.Consider the same assumptions as in Theorem 10.The C 2 smoothness conditions for the control points can be written as Proof.The statements follow by direct computation from (23) and from the expressions of the domain points in (10) and (16).
Using the relations between the domain points b * i and b * i in ( 16), we can immediately rewrite the C 2 smoothness conditions solely in terms of the control points of the basis (4).For instance, the third condition in Corollary 11 reads as and the fourth condition as The smoothness conditions in Corollary 11 show a structural similarity with the C 2 join of two adjacent triangular Bernstein-Bézier patches.We refer to [23,Example 2] for a geometric interpretation.

Stable bases for S 2 3 (T WS 3 )
In Sections 3.1 and 3.4, we have constructed two simplex spline bases for the space of C 2 cubic splines on the WS 3 split of a given triangle ∆.We have also shown that these bases enjoy similar properties to the Bernstein polynomial basis defined on a triangle.Here, we consider the space S 2 3 (T WS 3 ) and we take the similarity between these bases one step further: we extend the concept of minimal determining sets developed for the Bernstein polynomial basis [26] to our simplex spline bases, with the aim of constructing stable bases with local support for the space S 2 3 (T WS 3 ).For the sake of simplicity we focus on the simplex spline basis (14).
To stress their dependence on a specific triangle ∆, from now on we denote the simplex spline basis ( 14) by { B i,∆ , i = 1, . . ., 28}.Any spline s ∈ S 2 3 (T WS 3 ) can be identified by its coefficients { bi,∆ , i = 1, . . ., 28} with respect to the above basis over any triangle ∆ of T , i.e., where the coefficients in (24) must satisfy the smoothness conditions derived in Section 4.1 to ensure the C 2 joins across the edges of T .Following [26,Chapter 5], for any triangle ∆ of T we denote by D ∆ the set of domain points specified in (16).Then, we define a (minimal) determining set as follows.
Definition 12. Assume a set D ⊆ (∪ ∆∈T D ∆ ) is such that if s ∈ S 2 3 (T WS 3 ) has all the coefficients corresponding to elements in D equal to zero then s ≡ 0.Then, D is a determining set for S 2 3 (T WS 3 ).A determining set is a minimal determining set in case it has the smallest possible cardinality.3 (T WS 3 ).Filled circles: the domain points associated with a vertex.Filled squares: the domain points associated with an edge.Filled triangle: the domain points associated with a triangle.Empty circles: the domain points associated with coefficients obtained from those in the minimal determining set by imposing the smoothness conditions.By using the same line of arguments as the proof of [26,Theorem 5.13] we infer that the cardinality of a minimal determining set for S 2 3 (T WS 3 ) equals the dimension of the space.With the aim of specifying such a minimal determining set, let us first introduce some terminology regarding the domain points in a triangle ∆ = p 1 , p 2 , p 3 of T : • the domain points associated with the vertex p i are the six points in (16)  • the domain point associated with the triangle ∆ is the point b * 28 in (16).
We can construct a minimal determining set for S 2 3 (T WS 3 ) as follows; see also Figure 9.
Theorem 13.For a given triangulation T , let M be the set consisting of the following domain points: • for each vertex p of T , choose a triangle ∆ of T such that p is a vertex of ∆ and select the six domain points in ∆ associated with p; • for each edge of T , choose a triangle ∆ of T sharing this edge and select the three domain points in ∆ associated with the edge; • for each triangle ∆ of T , select the domain point associated with it.
Then, M is a minimal determining set for S 2 3 (T WS 3 ).
Proof.Let s ∈ S 2 3 (T WS 3 ).Assume that all its coefficients associated with the domain points in M are set equal to 0. Let p be any vertex of T .All the coefficients of s corresponding to the domain points associated with p (in any triangle of T surrounding p) are zero because either they belong to M or they are uniquely determined by the conditions in (21), the first two conditions in (22) and the first condition in (23) for C 2 smoothness across the edges emanating from p. Let pq be any edge of T .All the coefficients of s corresponding to the domain points associated with the edge (in any of the two triangles of T sharing the edge pq) are zero because either they belong to M or they are uniquely determined by the third condition in (22) and the third and fourth conditions in (23) for C 2 smoothness across the edge.Finally, the domain point associated with any triangle in T belongs to M and so the corresponding coefficient is 0. Hence, for any triangle ∆ of T all the coefficients in ( 24) are 0 and so s ≡ 0, i.e., M is a determining set.Moreover, the cardinality of M clearly equals the dimension of the space, see (19), and so M is a minimal determining set.Given a minimal determining set M for S 2 3 (T WS 3 ), suppose we assign values to all the coefficients corresponding to the domain points in it.The proof of Theorem 13 shows that these coefficients uniquely identify a spline function of S 2 3 (T WS 3 ) because all the remaining coefficients in the representation ( 24) can be deduced from the smoothness conditions.Therefore, any minimal determining set enables us to built a basis for S 2 3 (T WS 3 ).Let us consider the set of functions where B b * ,T is the unique function of S 2 3 (T WS 3 ) obtained by zeroing all the coefficients corresponding to the domain points in M except the one related to b * which is set equal to 1.By construction, the functions in ( 25) are clearly linearly independent (see also [26,Theorem 5.20]) and their number agrees with the dimension of the space because M is a minimal determining set.Thus, ( 25) is a basis for S 2 3 (T WS 3 ).The following proposition ensures that the elements of the above basis are uniformly bounded and have local support.Moreover, there exists a constant K only depending on the minimal angle of T such that Proof.Let b * be fixed.A direct inspection of the smoothness conditions in Theorems 8 and 10 immediately gives that the coefficients in (24) for B b * ,T are 0 whenever ∆ is not a triangle listed in the items a)-b)-c) and so B b * ,T vanishes outside the union of those triangles.In order to prove (26), we first note that the number of triangles surrounding a vertex of T and the (absolute value of the) barycentric coordinates of a point of a triangle with respect to an adjacent triangle are bounded in terms of the minimum angle of T (see [26,proof of Lemma 2.29]).Let ∆ be a triangle of T belonging to the support of B b * ,T .We have that (B b * ,T ) |∆ can be represented in the form (24) where the coefficients bi,∆ are obtained from the value 1 corresponding to the domain point b * ∈ M by applying the smoothness conditions in Theorems 8 and 10; these conditions consist of linear or quadratic relations involving barycentric coordinates of points in adjacent triangles.Therefore, denoting by b∆ the vector of these 28 coefficients, we get b∆ ∞ ≤ K , where K is a constant only depending on the minimum angle of T .Hence, from (18) we get Taking the maximum over all the triangles of T we arrive at (26) with K = 3K .
A basis with properties as in Proposition 14 is called a stable basis with local support.For such a basis, using the same line of arguments as the proof of [26,Theorem 5.22], we can show that for all c : where K − and K + are positive constants depending only on the smallest angle of T .The inequalities in (27) extend the local stability result in (18) to the full spline space S 2 3 (T WS 3 ).Similar stability results in any L q -norm for (a properly scaled version of) the basis can also be achieved; see again the proof of [26,Theorem 5.22].Furthermore, from the partition of unity property of the local basis in (14), we directly deduce that the global basis in (25) forms a partition of unity as well.
Note that the determining set in Theorem 13 is a stable local determining set in the sense of [26,Definition 5.16].This feature is, roughly speaking, the key ingredient in the proof of Proposition 14. Similarly to [26,Section 5.7], it also ensures that the full spline space S 2 3 (T WS 3 ) has optimal approximation power.More precisely, Theorems 5.18 and 5.19 in [26] hold true for S 2 3 (T WS 3 ).The global basis in (25) can be easily expressed over any triangle ∆ with respect to the local basis (4) through the conversion (15).Contrarily to this local basis, the functions in (25) are in general not nonnegative.However, paraphrasing [26,Section 5.8], we observe that the explicit basis in (25) has mainly a theoretical interest.For computation with splines belonging to S 2 3 (T WS 3 ) it is more convenient to work directly with the local representations provided by the bases (4) or (14), rather than with the basis for the full spline space.

Concluding remarks
In [11], a simplex spline basis was described for the C 1 quadratic spline space on the Powell-Sabin 12 split, which is the quadratic member of the Wang-Shi split family.In this paper, we have addressed the C 2 cubic case and constructed two simplex spline bases for the WS 3 split.The characteristics of the C 2 cubic simplex spline bases make it unnecessary to consider separate polynomial representations on each of the numerous polygonal regions of the partitioned macrotriangle.This paves the path for a practical construction of globally C 2 cubic splines on any triangulation by extending the concept of minimal determining sets.
In the following, we outline some implementation aspects and we identify few problems where the provided simplex spline bases may be prosperous, in order to complement the theoretical interest of our investigation with an application-oriented perspective.We end with a discussion on a higher-order extension of the construction.

Implementation aspects
For computation with splines belonging to S 2 3 (T WS 3 ), it is convenient to work directly with the local representations provided by the bases (4) or (14).On the one hand, evaluation of the simplex spline basis functions can be achieved by applying the recurrence relation (B-recurrence) of simplex splines; see Section 2.1.On the other hand, it might be more convenient to use the explicit expressions of the simplex spline basis functions (4) given in Table 1 in the appendix.The alternative basis functions ( 14) can be immediately deduced from the previous ones by means of the linear relation (15).
Having at our disposal such a table, evaluation of any spline in S 2 3 (∆ WS 3 ) can be efficiently performed by combining a lookup-table process with a search algorithm based on boolean vectors.Given (the barycentric coordinates of) any point p in ∆, the values of the simplex spline basis functions (4) at p can be directly obtained from Table 1 once we have figured out which polygonal region of the macro-triangle the evaluation point belongs to.Since the WS 3 split is a cross-cut partition, any polygonal region in ∆ WS 3 is uniquely identified by the sign of the linear expressions of the 18 interior lines in the split.These signs can be interpreted as binary digits of an integer belonging to [0, 2 18 − 1].Therefore, in order to detect which polygonal region of the macro-triangle a given point p belongs to, it suffices to evaluate all the 18 interior lines at p, to collect the resulting signs in a boolean vector, and to interpret such a vector as binary digits of an integer.A similar search algorithm has been described in [11,Algorithm 1.1].
It is important to remark that the selection of the different polynomial pieces is just an implementation aspect.Thanks to the characteristics of the simplex spline representation, there is a single control net to facilitate control and early visualization of a spline function over each element ∆ in T .This single control net makes that the complex geometry of the WS 3 split (consisting of 75 polygons including triangles, quadrilaterals, and pentagons) is transparent to the user.In this perspective, an interesting topic of possible future research is to investigate whether the control net introduced in the paper can give rise to a de Casteljau/de Boor-type algorithm for evaluation of splines in S 2 3 (∆ WS 3 ).

Application areas
Splines on (refined) triangulations are valuable in several application areas.When dealing with bivariate/multivariate problems, the straightforward approach is to rely on tensor-product structures, and in particular tensor-product splines.Tensor-product structures offer several advantages, mainly the simplicity of their use and the inheritance of univariate properties.Major drawbacks, however, are the lack of adequate local refinement and the struggle to represent geometries with complicated shapes.Although there are several appealing extensions of tensor-product splines towards local refinement (see, e.g., [15,16,44]) and complex geometries (see, e.g., [2,34,39]), splines on triangulations emerge as the natural tool to efficiently deal with problems where local features has to be detected, modeled, or simulated.As mentioned in the introduction, low degree splines are preferable due to their stable behavior and their low computational complexity.In particular, univariate C 2 cubic splines are one of the most used tools in modeling, approximation, and simulation.Constructing C 2 splines of low degree on triangulations is a difficult task, but their interest remains unquestionable in the bivariate setting.We limit ourselves to mention two important application areas: computer aided surface modeling and numerical simulation.
In computer aided design/manufacturing (CAD/CAM) high quality free-form surfaces are of utmost importance.The quality of the surfaces can be checked by different techniques, such as the well-established isophotes [35], to detect irregularities of intrinsic measures of surface smoothness like the Gaussian curvature or the distribution of the surface normals.For milling surfaces by five axis machines, second derivatives should not jump too much across edges and C 2 smoothness is desirable.In this context, our C 2 cubic simplex spline representations on triangulations could be beneficial.The spline surfaces could be constructed by direct (interactive) modeling via the control net or by data fitting using quasi-interpolation schemes based on the Marsden-like identity, similar to [48].In CAD/CAM systems it is common to rely on general parametric surfaces; in our case such surfaces are specified on each macro-triangle by a control net consisting of triangles and quadrilaterals.See also [57] for the use of simplex splines in the context of parametric surface reconstruction.As a possible future work, it is of interest to investigate the interplay with tensorproduct (piecewise) bicubic parametric surfaces in Bernstein-Bézier (or B-spline) form, which are ubiquitous in industrial applications.In particular, an important question is whether one can blend standard bicubic Bernstein-Bézier patches with parametric triangular patches whose components are C 2 cubic splines represented in terms of the simplex spline bases introduced in the paper.
Isogeometric analyis (IgA) is a numerical simulation paradigm that extends finite element analysis (FEA) by providing a true design-through-analysis methodology [13].The isogeometric paradigm has some important advantages over traditional FEA.The geometry of the physical domain is exactly described, so the interaction with the CAD system during any further refinement process in the analysis phase is eliminated.Moreover, the discretization spaces possess an inherent higher smoothness (with respect to the polynomial degree) than classical FEA spaces, leading to a higher accuracy per degree of freedom [4,41].The success of IgA roots in the above two properties, the latter being even more relevant.Besides the use of spline spaces based on (local) tensor-product structures and rather involved multipatch constructions (see, e.g., [3,6,22,42,53]), a powerful IgA formulation has been obtained by considering spline spaces on triangulations (see, e.g., [5,21,50,51,55]).In particular, spline representations obtained from local Bernstein representations by means of minimal determining sets have been profitably applied and efficiently implemented via Bézier-extraction [21].In this context, the space of C 2 cubic splines defined on the WS 3 refinement of a given triangulation is appealing because it combines low degree and high smoothness.Our simplex spline bases are the natural counterpart of Bernstein polynomials to define stable global bases by means of minimal determining sets (see Section 4), and allow for a straightforward extension of the Bézier-extraction procedure for practical implementation.Of course, in order to efficiently exploit the potential of the space S 2 3 (T WS 3 ) and its local representations in terms of simplex spline bases in the context of IgA, several steps are still missing, for instance, the need for tailored quadrature rules.and n k is the normal direction of the edge opposite to vertex p k .Given a general triangulation T , this scheme can be used to construct a globally C d−1 spline of degree d on T where every triangle is refined with the WS d split.Such a construction is local, in the sense that the spline can be built on each macro-triangle ∆ of T separately, and the simplex spline basis would then be useful to represent the corresponding spline piece on ∆, without considering explicitly the complicated geometry in the WS d split.

2 p 1 , 3 Figure 2 :Theorem 2 .
Figure 2: Labeling of the knots on the boundary of the triangle ∆.

Figure 4 :
Figure 4: Hermite degrees of freedom on the WS 3 split.

Figure 7 :
Figure 7: Sequences of knots for an alternative set of simplex spline basis functions for S 2 3 (∆ WS 3 ).Each black disc shows the position of a knot and the number inside indicates its multiplicity.

Figure 8 :
Figure 8: A C 1 spline surface on two adjacent domain triangles.The five pairs of triangles in the control nets that must be coplanar according to the smoothness conditions are colored.

28 join C 2
across the common edge if and only if they join C 1 and in addition b

Figure 9 :
Figure 9: A minimal determining set for S 23 (T WS 3 ).Filled circles: the domain points associated with a vertex.Filled squares: the domain points associated with an edge.Filled triangle: the domain points associated with a triangle.Empty circles: the domain points associated with coefficients obtained from those in the minimal determining set by imposing the smoothness conditions.

•••
with the i-th barycentric coordinate ≥ 2/3, i = 1, 2, 3 (for instance, the domain points associated with p 1 the domain points associated with the edge p 1 p 2 are the three points b the domain points associated with the edge p 1 p 3 are the three points b the domain points associated with the edge p 2 p 3 are the three points b

Proposition 14 .
For a given triangulation T , let M be a minimal determining set for S 2 3 (T WS 3 ) as specified in Theorem 13.For all b * ∈ M, the support of B b * ,T is contained in a) the union of the triangles of T sharing the vertex p, if b * is a domain point in M associated with the vertex p; b) the union of the two triangles of T sharing the edge pq, if b * is a domain point in M associated with the edge pq; c) the triangle ∆ of T , if b * is the domain point in M associated with the triangle ∆.

Figure 11 :
Figure 11: The simplex spline basis function B 4 and its support.

Figure 12 :
Figure 12: The simplex spline basis function B 10 and its support.

Figure 13 :
Figure 13: The simplex spline basis function B 16 and its support.

Figure 14 :
Figure 14: The simplex spline basis function B 19 and its support.

Figure 15 :
Figure 15: The simplex spline basis function B 22 and its support.

Figure 16 :
Figure 16: The simplex spline basis function B 28 and its support.

Figure 17 :
Figure 17: Numbering of the 58 intersection points v k in the WS 3 split.
Note that the scaling factors sum up to |∆|.There are seven different types of simplex splines in B. For each type, a representative B i is depicted in Figures10-16in the appendix.Explicit expressions of their polynomial pieces are given in Table1in the appendix; the remaining ones can be obtained by symmetry.On any edge of ∆, there are six basis functions nonzero.Their restrictions to that edge are nothing but the set of univariate C 2 cubic B-splines defined on a uniform open-knot sequence with two interior knots.For instance, for the edge p 1 p 2 , they correspond to the univariate cubic B-splines on the knot sequence specified by {p 1 , p 1 , p 1 , p 1 , p 3,1 , p 3,2 , p 2 , p 2 , p 2 , p 2 }.