Tom Lyche

Jean-Louis Merrien
email: jmerrien@insa-rennes.fr

Tomas Sauer
email: tomas.sauer@uni-passau.de

Simplex-Splines on the Clough-Tocher Split with Arbitrary Smoothness

Keywords: Triangle Mesh, Piecewise polynomials, Interpolation, Simplex Splines, Marsden-like Identity

The space of piecewise polynomials of smoothness r and degree 3r is considered on the Clough-Tocher split of a triangle. For any r ≥ 1 we give a basis of simplex splines for this space, then a Marsden-like identity, which is proved explicitly for r ≤ 3 and symbolically for 4 ≤ r ≤ 6. In addition, generalizing results for r = 1, we prove for r = 2, 3 a geometry independent bound for the condition number in the infinity norm of this basis, and conditions to connect two triangles with smoothness r. For odd values of r, and especially r = 3, using the previous constructions on each triangle, and the connecting conditions, we obtain global Cr continuity on any triangulation.

Introduction

Splines over triangulations have applications in several branches of the sciences ranging from finite element analysis, surfaces in computer aided design and other engineering problems, see for example [START_REF] Cohen | Geometric Modeling with Splines[END_REF][START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF][START_REF]Spline Functions on Triangulations[END_REF]. For many of these applications, piecewise linear C 0 surfaces do not suffice. In some cases, we need smoother elements for modeling, or higher polynomial degrees to increase the approximation order.

In this paper, we are interested in the space of polynomial splines over a triangulation ∆ of a polygonal domain Ω of R 2 , S r d (∆) := { f ∈ C r (Ω) : f | T ∈ P d , for all T ∈ ∆} , where d > r > 0 are given integers, and P d is the space of bivariate polynomials of total degree ≤ d. The dimension of this finite dimensional vector space is difficult to determine in general [START_REF]Spline Functions on Triangulations[END_REF], but with the restriction d ≥ 3r + 2 its dimension can be expressed solely in terms of d and r, see [START_REF] Ibrahim | Super spline spaces of smoothness r and degree d ≥ 3r + 2[END_REF].

We can use lower degrees if we are willing to split each triangle into a number of subtriangles. The most well known examples are the Clough-Tocher split [START_REF] Clough | Finite element stiffness matrices for analysis of plates in bending[END_REF], and the Powell-Sabin 6 and 12 splits [START_REF] Powell | Piecewise quadratic approximations on triangles[END_REF]. For these splits each triangle is divided into 3,6 and 12 subtriangles, respectively. For material on these splits and B-spline like bases for splines on triangulations see [START_REF] Alfeld | Smooth macro-elements based on Powell-Sabin triangle splits[END_REF][START_REF] Cao | A finite element framework based on bivariate simplex splines on triangle configurations[END_REF][START_REF] Cohen | A B-spline-like basis for the Powell-Sabin 12-split based on simplex splines[END_REF][START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF][START_REF] Davydov | Refinable C 2 piecewise quintic polynomials on Powell-Sabin-12 triangulations[END_REF][START_REF] Grošelj | A normalized representation of super splines of arbitrary degree on Powell-Sabin triangulations[END_REF][START_REF] Grošelj | Construction and analysis of cubic Powell-Sabin B-splines[END_REF][START_REF] Lai | Macro-elements and stable local bases for splines on Clough-Tocher triangulations[END_REF][START_REF]Macro-elements and stable local bases for splines on Powell-Sabin triangulations[END_REF][START_REF]Spline Functions on Triangulations[END_REF][START_REF] Laghchim-Lahlou | Triangular finite elements of HCT type and class C p[END_REF][START_REF] Lyche | Simplex-splines on the Clough-Tocher element[END_REF][START_REF] Lyche | Stable simplex spline bases for C 3 quintics on the Powell-Sabin 12-split[END_REF][START_REF] Sablonnière | Composite finite elements of class C k[END_REF][START_REF] Schumaker | Smooth macro-elements on Powell-Sabin-12 splits[END_REF][START_REF] Speleers | A normalized basis for quintic Powell-Sabin splines[END_REF][START_REF]A normalized basis for reduced Clough-Tocher splines[END_REF][START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusionreaction problems[END_REF][START_REF] Speleers | Construction of normalized B-splines for a family of smooth spline spaces over Powell-Sabin triangulations[END_REF][START_REF]A family of smooth quasi-interpolants defined over Powell-Sabin triangulations[END_REF][START_REF]A new B-spline representation for cubic splines over Powell-Sabin triangulations[END_REF][START_REF] Speleers | Optimizing domain parameterization in isogeometric analysis based on Powell-Sabin splines[END_REF][START_REF] Hughes | Polynomial splines of non-uniform degree on triangulations: Combinatorial bounds on the dimension[END_REF][START_REF] Wang | S µ µ+1 surface interpolations over triangulations[END_REF][START_REF] Wang | Isogeometric shape optimization on triangulations[END_REF][START_REF] Ženíšek | A general theorem on triangular finite C (m) -elements[END_REF].

Here we consider the Clough-Tocher triangulation ∆ CT , see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Clough | Finite element stiffness matrices for analysis of plates in bending[END_REF], where each triangle in the original triangulation ∆ is split into 3 subtriangles by connecting the vertices of each triangle to its barycenter , see Figure 1.1. In [START_REF] Laghchim-Lahlou | Triangular finite elements of HCT type and class C p[END_REF], Hermite interpolation problems were considered for super-spline subspaces of S r 3r+1 (∆ CT ) and S r 3r (∆ CT ) for r even and odd, respectively. It was also stated that these degrees are minimal for global C r . See also [START_REF] Sablonnière | Composite finite elements of class C k[END_REF]. In [START_REF] Lai | Macro-elements and stable local bases for splines on Clough-Tocher triangulations[END_REF] stable local bases were constructed for even smaller super-spline subspaces of S r 3r+1 (∆ CT ) and S r 3r (∆ CT ). In this paper, we consider for any r ∈ N, the spaces S r 3r ( ) on a single triangle T = in ∆ CT for which we construct a B-spline like basis made out of simplex splines. This extends and generalizes the case r = 1 that was considered in [START_REF] Lyche | Simplex-splines on the Clough-Tocher element[END_REF]. For more on the Clough-Tocher split see [START_REF] Alfeld | A bivariate C 2 Clough-Tocher scheme[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Farin | A modified Clough-Tocher interpolant[END_REF][START_REF] Grošelj | Construction and analysis of cubic Powell-Sabin B-splines[END_REF][START_REF] Lai | Macro-elements and stable local bases for splines on Clough-Tocher triangulations[END_REF][START_REF]Spline Functions on Triangulations[END_REF][START_REF] Laghchim-Lahlou | Triangular finite elements of HCT type and class C p[END_REF][START_REF] Mann | Cubic precision Clough-Tocher interpolation[END_REF][START_REF]A normalized basis for reduced Clough-Tocher splines[END_REF][START_REF]A new B-spline representation for cubic splines over Powell-Sabin triangulations[END_REF]. Looking in more detail at the cases r = 2, 3, corresponding to degrees d = 6, 9, we moreover give explicit formulas for connecting two neighboring triangles in a C r fashion across an edge using Bernstein-Bézier techniques, and give an upper bound for the L ∞ condition number of the basis. This upper bound is independent of the shape of the triangle T . We also give a Marsden-like identity for the reproduction of polynomials which is proved for r ≤ 3 and shown symbolically for r ≤ 6. We conjecture that it holds for any r.

It was shown in [START_REF] Lai | Macro-elements and stable local bases for splines on Clough-Tocher triangulations[END_REF][START_REF] Laghchim-Lahlou | Triangular finite elements of HCT type and class C p[END_REF] that global C 2 continuity cannot be achieved for d = 6 for a general triangulation refined by Clough-Tocher splits into ∆ CT . However, for r odd it follows, again from [START_REF] Lai | Macro-elements and stable local bases for splines on Clough-Tocher triangulations[END_REF][START_REF] Laghchim-Lahlou | Triangular finite elements of HCT type and class C p[END_REF], that global C r continuity holds for d = 3r. This means in particular that the formulas for C 3 continuity across an edge in Section 1.5.2 can be used to compute with elements in S 3 9 (∆ CT ), using the simplex spline basis on each triangle in S 3 9 (∆) in the usual Bernstein-Bézier fashion.

The paper is organized as follows. Since it depends heavily on properties of Bernstein polynomials and simplex splines, we recall some well known facts about these functions in the next section. Section 1.3 introduces the Clough-Tocher split on a triangle and gives some basic facts about the dimension of the associated spline spaces. In Section 1.4, we derive a local basis for the spline space and prove its basis property, mainly by a combinatorial argument. In Section 1.5, we consider the cases of global C 2 and C 3 regularity in more detail and state and prove a Marsden-like identity.

Preliminaries

In this section we recall some properties of Bernstein polynomials on a triangle and bivariate simplex splines. Here we use the notation d ∈ N 0 := N ∪ {0}, and let S denote the convex hull of the set S ⊂ R 2 .

Bernstein Polynomials

For a given nondegenerate triangle T := { p 1 , p 2 , p 3 } ∈ R 2 , and i, j, k ∈ N 0 , the Bernstein polynomial B d i jk : R 2 → R of degree d := i + j + k ∈ N 0 , is defined by

B d i jk (x, y) = B d i jk (β 1 , β 2 , β 3 ) := d! i! j!k! β i 1 β j 2 β k 3 , (1.1)
where β = (β 1 , β 2 , β 3 ) are the barycentric coordinates of x = (x, y) ∈ R 2 with respect to T , i. e.,

x = β 1 p 1 + β 2 p 2 + β 3 p 3 , β 1 + β 2 + β 3 = 1. (1.2)
The barycentric form of Marsden's identity for Bernstein polynomials then is simply the multinomial expansion

(u 1 β 1 + u 2 β 2 + u 3 β 3 ) d = i+j+k=d u i 1 u j 2 u k 3 B d i jk (β 1 , β 2 , β 3 ), (u 1 , u 2 , u 3 ) ∈ R 3 , β 1 + β 2 + β 3 = 1, (1.3)
where in the expression i+j+k=d it is understood that i, j, k ∈ N 0 , which is consistent with the convention that B d i jk = 0 if one of the indices becomes negative.

Taking partial derivatives of order ν, µ, κ ∈ N 0 with respect to u 1 , u 2 , u 3 , respectively, in (1.3) and setting u = (1, 1, 1) we obtain

β ν 1 β µ 2 β κ 3 = i+j+k=d γ i jk (β ν 1 β µ 2 β κ 3 )B d i jk (β 1 , β 2 , β 3 ), ν + µ + κ ≤ d, γ i jk (β ν 1 β µ 2 β κ 3 ) ∈ R. (1.4) We find 1 = i+j+k=d B d i jk (β 1 , β 2 , β 3 ) (β 1 , β 2 , β 3 ) = i+j+k=d b * i jk B d i jk (β 1 , β 2 , β 3 ), b * i jk = i d , j d , k d .
(1.5)

The vector b * i jk is called the barycentric form of the domain point of B d i jk . From (1.4) and (1.5) it follows that the elements in the set

B d := {B d i jk : i, j, k ≥ 0, i + j + k = d} (1.6)
form a partition of unity basis for P d . Indeed, the number of elements #B d of B d equals d+2 2 , the dimension of P d . We refer to [START_REF]Spline Functions on Triangulations[END_REF] for further properties of B d i jk .

Bivariate Simplex Splines

For our purpose it is convenient to work with area normalized bivariate simplex splines [START_REF] Lyche | Stable simplex spline bases for C 3 quintics on the Powell-Sabin 12-split[END_REF] of degree d ≥ 0 with knots

K := {k 1 , . . . , k d+3 }, k j ∈ R 2 \ {0}, j = 1, . . . , d + 3.
We can consider K either as a multiset or as a matrix

K ∈ R 2×(d+3) . The simplex spline Q[K] : R 2 → R, is now defined as Q[K](x) = 0 for all x ∈ R 2 if rank K < 2
, and otherwise

Q[K] := area (T ) d+2 2 M[K], (1.7)
where area (T ) is the positive area of a fixed reference triangle in the original triangulation ∆, and M[K] is a unit integral bivariate normalized simplex spline, defined as a linear functional

M[K] : C(R 2 ) → R given by M[K](ϕ) := (d + 2)! ∫ S d+2 ϕ d+3 j=1 k j t j dt 1 • • • dt d+2 , ϕ ∈ C(R 2 ), (1.8)
with

S n := {(t 1 , . . . , t n+1 ) ∈ R n+1 : t i ≥ 0, n+1 i=1 t i = 1}, the unit simplex in R n , n ∈ N. If K has full rank, rank(K) = 2, then M[K] can be identified with a function M[K] : R 2 → R, and we write (1.8) in the form ∫ R s M[K](x)ϕ(x)dx = (d + 2)! ∫ S d+2 ϕ d+3 j=1 k j t j dt 1 • • • dt d+2 ,, ϕ ∈ C(R 2 ). (1.9)
We mention the following well-known properties of M[K] [START_REF] Micchelli | On a numerically efficient method for computing multivariate B-splines[END_REF][START_REF] Prautzsch | Bézier and B-Spline Techniques, Mathematics and Visualization[END_REF] and

Q[K]. 1. Q[K] and M[K] are piecewise polynomials of degree d = #K -3 witht support K . 2.
Local smoothness: Across a knot line, which is a line in the complete graph associated with K, we have that

M[K], Q[K] ∈ C d+1-µ
, where µ is the number of knots on that knot line, including multiplicities. 3. Differentiation formula: For u = (u 1 , u 2 ) ∈ R 2 and any choice of a 1 , . . . , a d+3 such that j a j k j = u, j a j = 0, one has

D u M[K] = (d + 2) d+3 j=1 a j M[K \ k j ], D u Q[K] = d d+3 j=1 a j Q[K \ k j ], (1.10)
where

D u := u 1 D 1 + u 2 D 2 and D 1 , D 2 denotes partial derivatives. (A-recurrence) 4. Recurrence relation: For any x ∈ R 2 and any b 1 , . . . , b d+3 such that j b j k j = x, j b j = 1, one has M[K](x) = d + 2 d d+3 j=1 b j M[K \ k j ](x), Q[K](x) = d+3 j=1 b j Q[K \ k j ](x). (1.11)
(B-recurrence) 5. Knot insertion formula: For any y ∈ R 2 and any c 1 , . . . , c d+3 such that j c j k j = y, j c j = 1, one has

M[K] = d+3 j=1 c j M[K ∪ y \ k j ], Q[K] = d+3 j=1 c j Q[K ∪ y \ k j ].
(1.12) (C-recurrence) 6. Degree zero:

For K = {k 1 , k 2 , k 3 }
T 3 := p T , p 1 , p 2 M[K](x) := 1 area( K ) , Q[K](x) := area(T ) area( K ) , x ∈ K o , M[K](x) := Q[K](x) = 0, x K , (1.13)
where S o is the interior of the set S. The values of M[K] and Q[K] on the boundary of K has to be dealt with separately, see below.

We refer to [START_REF] Micchelli | On a numerically efficient method for computing multivariate B-splines[END_REF][START_REF] Prautzsch | Bézier and B-Spline Techniques, Mathematics and Visualization[END_REF] for further properties of M[K].

Clough-Tocher split, dimension and smoothness

Given a nondegenerate triangle T in R 2 , we connect the vertices p 1 , p 2 , p 3 to the barycenter p T := ( p 1 + p 2 + p 3 )/3. With this construction, known as the Clough-Tocher split , we obtain three subtriangles T 1 := p T , p 2 , p 3 , T 2 := p T , p 3 , p 1 and T 3 := p T , p 1 , p 2 , see Figure 1.1. We consider the spline space S r d ( ) with respect to these three subtriangles on T . To obtain a unique function value at each point in T we associate the half open edges

p i , p T ) := {(1 -t) p i + t p T : 0 ≤ t < 1}, i = 1, 2, 3,
to the three subtriangles of T as follows

p 1 , p T ∈ T 2 , p 2 , p T ∈ T 3 , p 3 , p T ∈ T 1 , (1.14)
and we somewhat arbitrarily associate the point p T to T 2 . We first give the dimension of S r d ( ), and a set of conditions cs r d,i induced by requesting smoothness C r across an edge p i , p i+1 of two neighbouring triangles in ∆. For 1 ≤ m ≤ d and s ∈ S m d ( ) let D n i s( p i j,ℓ ) denote a derivative across the edge p i , p i+1 at the point p i j,ℓ , where

p i j,ℓ := ( j -ℓ + 1) p i + ℓ p i+1 j + 1 , for ℓ = 1, . . . , j, j ∈ N, and i = 1, 2, 3, p 4 := p 1 . For i = 1, 2, 3 we let, for m ∈ N, V m,i (s) := {D α 1 D β 2 s( p i ) : 0 ≤ α + β ≤ m}, E even m,i (s) := {D j-1 n i s( p i j,ℓ ) : ( j, ℓ) ∈ I m }, E odd m,i (s) := {D j n i s( p i j,ℓ ) : ( j, ℓ) ∈ I m }, I m := {( j, ℓ) : ℓ = 1, .
. . , j and j = 1, . . . , m}.

(1.15) We define for s ∈ S r d ( ) and i = 1, 2, 3 the following independent degrees of freedom (interpolation conditions),

cs r d,i := V q-1,i (s) ∪ E even r+1,i (s), if d = 2q is even, V q,i (s) ∪ E odd r ,i (s), if d = 2q + 1 is odd. (1.16)
These conditions guarantee smoothness C r across the edge p i , p i+1 .

Proposition 1

We have

dim S r d ( ) = r + 2 2 + 3 d -r + 1 2 + d-r j=1 (r + 1 -2 j) + (1.17)
and

#cs r 2q,i = q + 1 2 + r + 2 2 , #cs r 2q+1,i = q + 2 2 + r + 1 2 . (1.18)
Moreover, dim S r 3r ( ) ≥ #cs r 3r := 3 i=1 #cs r 3r ,i for all r ∈ N. Proof The formula (1.17) for dim S r d ( ) follows from Theorem 9.3 in [START_REF]Spline Functions on Triangulations[END_REF] with n = m v = 3. Since V m,i (s) and I m given by (1.15) contain m+2 2 and m+1 2 linearly independent, hence nonredundant, elements, respectively, (1.18) follows.

In the even and odd case (1.17) and (1.18) .

In Figure 1.2 we display the degrees freedom cs r 3r for r = 1, 2, 3.

A basis for S r 3r

We now focus on a special collection of simplex splines on the Clough-Tocher split of one single triangle in the following way: for integers i, j, k, ℓ, we consider the simplex spline

Q p {i } 1 , p { j } 2 , p {k } 3 , p {ℓ } T
, where p {i } j , means that the vertex p j has multiplicity i, i. e., is repeated i times. For brevity, we will use the alternative notations, a more compact and a more illustrative one, for simplex splines of degree d = i + j + k + ℓ -3 from now on, namely,

Q p {i } 1 , p { j } 2 , p {k } 3 , p {ℓ } T =: [i, j, k; ℓ] =: i jkℓ =: k ℓ j i , i, j, k ∈ N, ℓ ∈ N 0 := N ∪ {0}.
(1.20)

By the local smoothness property for simplex splines we note that i jkℓ has smoothness d + 1iℓ across the knotline p 1 , p T for ℓ > 0. In this notation, the Bernstein polynomials B i jk (u, v, w) := (i+j+k)! i!j!k! u i v j w k of degree

d = i + j + k have the form B i jk = i+1, j+1,k+1,0 = k + 1 0 j + 1 i + 1 , i + j + k = d. (1.21)
We will use the more graphic form on the right hand side of (1.20) whenever possible to make the basic ideas more accessible, but it is convenient to use the more compact notations in computations. The well-known knot insertion formula for the insertion of a knot at the barycenter, as a special case of the C-recurrence (1.12), for example, akes the form

k ℓ j i = 1 3 k ℓ + 1 j i -1 + k ℓ + 1 j -1 i + k -1 ℓ + 1 j i (1.22)
in the graphical notation.

Definition 1

The number µ := d + 1r is called the maximum multiplicity of an interior knot line of the simplex spline i jkℓ . Some obvious properties for simplex splines to be in S r d ( ) are listed in the following lemma.

Lemma 1 If i jkℓ ∈ S r d ( ) then 1. i + j + k + ℓ = d + 3 = µ + r + 2, 2. if, in addition, ℓ 0 then regularity C r implies max{i, j, k} + ℓ ≤ µ.
Next, let us focus on specific elements of S r 3r ( ) that will turn out to be useful. The three types are categorized by whether max{i, j, k} + ℓ equals µ, exceeds it or is strictly less than µ, where the second case occurs only for ℓ = 0. Definition 2 Let Σ r 3r denote the set of elements i jkℓ of S r 3r ( ) that consists of the elements of the following three types: Type (1): max{i, j, k} + ℓ = µ and min{i, j, k} ≥ 1, Type (2): ℓ = 0 and max{i, j, k} > µ, Type (3): max{i, j, k} + ℓ < µ and min{i, j, k} = 1.

Remark 1

The types are symmetric with respect to i, j, k, i. e., all subclasses are closed under permutation of the indices.

For r = 1, 2, 3, 4, Definition 2 results in the following set of splines. One central result of this paper, stated in Theorem 1, is that Σ r 3r is in fact a basis for S r 3r ( ). The proof of this fact will consist of showing that Σ r 3r is a subset of dim(S r 3r ( )) linearly indenpendent elements of the space S r 3r ( ).

To count the number of elements of Σ r 3r , we start with some bounds of ℓ with respect to the different types of functions in Σ r 3r .

Removed Type (1)

Type ( 2)

Type ( 1)

Type ( 2 Lemma 2 For Type [START_REF] Alfeld | A bivariate C 2 Clough-Tocher scheme[END_REF] we have ℓ ≤ d -3r/2 while for Type [START_REF] Cao | A finite element framework based on bivariate simplex splines on triangle configurations[END_REF]

ℓ ≤ d -2r -2 holds.
Proof For a Type (1) element, assume that i

+ ℓ = µ, then j + ℓ ≤ µ and k + ℓ ≤ µ, hence i + j + k + 3ℓ ≤ 3µ. Since i + j + k + ℓ = d + 3 and µ = d -r + 1, we thus deduce that d + 3 + 2ℓ ≤ 3d + 3 -3r or ℓ ≤ d -3r/2. For Type (3), we assume that k = 1, hence i + ℓ ≤ µ -1 and j + ℓ ≤ µ -1, then i + j + k + 2ℓ ≤ 2µ -1. Since i + j + k + ℓ = d + 3 and µ = d -r + 1, it follows that d + 3 + ℓ ≤ 2d -2r + 2 -1 or ℓ ≤ d -2r -2.
Next we prove that the sum of the numbers of elements of the three types is exactly the dimension of dim S r 3r ( ).

Proposition 2 For r ≥ 1 and d = 3r, we have the following table according to the parity of r.

r 2s, s > 0 2s + 1, s ≥ 0 d 6s 6s + 3 dim S r d ( ) 27s 2 + 9s + 1 3(9s 2 + 12s + 4) #T ype(1) 3s(5s + 3) + 1 3(s + 1)(5s + 3) #T ype(2) 3s(2s + 1) 3(s + 1)(2s + 1) #T ype(3) 3s(2s -1) 3s(2s + 1) (1.23) Hence, #T ype(1) + #T ype(2) + #T ype(3) = dim S r 3r ( ), r ≥ 1. (1.24)
Proof We begin with r of even parity, i. e., r = 2s, d = 3r = 6s, µ = 4s + 1, i + j + k + ℓ = 6s + 3 and count the basis elements of different types.

1. Type (1): According to Lemma 2 we have ℓ ≤ d -3r/2 = 3s and a generic element of Type ( 1) is of the form

[i = s + p + 1, j, k = 2s + 2 -j; 3s -p], 1 ≤ j ≤ i, 1 ≤ k ≤ i, (1.25)
from which it follows that sp + 1 ≤ j ≤ 2s + 1 and p = 0, . . . , 3s. We count the number of elements with respect to p:

a. For p = 0: 1 element. The only possible choice is [s + 1, s + 1, s + 1; 3s]. b. For 1 ≤ p ≤ s: 6p elements.
The elements [i, j, k; 3sp] are of the form

[s + p + 1, s -p + 1 + q, s + p + 1 -q; 3s -p], q = 0, . . . , 2p -1, (1.26)
i. e., i = s + p + 1, j = sp + 1 + q, s + p + 1q, and with the permutations

[i, j, k], [ j, k, i] and [k, j, i],
which gives 3 elements for any q in (1.26). We notice that in two cases we obtain the same three elements up to the permutation, namely, for q = 0 the values i = k = s + p + 1, j = sp + 1, and for q = 2p the values i

= j = s + p + 1, k = s -p + 1.
In the other cases, 0 < q < 2p, we have

j = s -p + 1 + q < i = s + p + 1 and k = s + p + 1 -q < i = s + p + 1, so that the elements are different. c. For s + 1 ≤ p ≤ 3s: 6s + 3 elements consisting of [s + p + 1, q + 1, 2s -q + 1; 3s -p], q = 0, . . . , 2s, (1.27)
and again the permutations

[i, j, k], [ j, k, i], [k, j, i].
Hence the total number of elements of Type ( 1) is

1 + s p=1 6p + 3s p=s+1 (6s + 3) = 1 + 3s(s + 1) + 2s(6s + 3) = 1 + 3s(5s + 3)
as listed in (1.23). 2. Type (2) requires ℓ = 0 by definition and the generic elements of the form

[i = 4s + 2 + p, j, k = s + 1 -j -p; 0], 1 ≤ j ≤ i, i ≤ k ≤ i, (1.28)
with the additional constraint i ≥ µ + 1 = 4s + 2 that leads to 1 ≤ j ≤ 2sp, p = 0, . . . , 2s -1. This gives 3(2sp) elements

[4s + 2 + p, q + 1, 2s -p -q; 0], q = 0, . . . , 2s -p -1, (1.29)
and the respective permutations so that the total number of elements of Type ( 2) is

2s-1 p=1 3(2s -p) = 2s q=1 q = 3s(2s + 1)
in this case. 3. Type (3): up to symmetry we can assume that k = 1 and, by Lemma 2 that ℓ ≤ d -2r -2 = 2s -2. The generic element is of the form

[i, j = 4s + 4 + p -i, k = 1; ℓ = 2s -2 -p], i + ℓ ≤ µ -1 = 4s, j + ℓ ≤ 4sp = 0, . . . , 2s -2, (1.30) since i + j + k + ℓ = d + 3 = 6s + 3. Hence, for 0 ≤ p ≤ 2s -2 we get 3(p + 1) elements [2s + 2 + p -q, 2s + 2 + q, 1; 2s -2 -p], q = 0, . . . , p, (1.31)
and the permutations

[i, j, k], [ j, k, i], [k, j, i], leading to a total of 2s-2 p=0 3(p + 1) = 3s(2s -1)
elements of Type (3).

In the case of odd parity, i. e., r = 2s + 1, d = 3r = 6s + 3, µ = 4s + 3, i + j + k + ℓ = 6s + 6, we proceed in the same way and distinguish by types.

1. For Type [START_REF] Alfeld | A bivariate C 2 Clough-Tocher scheme[END_REF] we have the bound ℓ ≤ 3s + 1 and the generic element

[i = s + p + 2, j, k = 2s + 3 -j; 3s + 1 -p], 1 ≤ j ≤ i 1 ≤ k ≤ i, p = 0, . . . , 3s + 1. (1.32)
Again, Type (1) request the distinction of several cases according to p.

a. For 0 ≤ p ≤ s: 6p + 3 elements The generic elements are

[s + p + 2, s + 1 -p + q, s + 2 + p -q; 3s + 1 -p], q = 0, . . . , 2p, (1.33)
and the permutations

[i, j, k], [ j, k, i], [k, j, i].
Again, we notice that we obtain the same three elements up to the permutations for q = 0, namely i = k = s + p+2, j = s -p+1, and for q = 2p+1, namely i = j = s + p+2, k = sp + 1, respectively. For 0 < q < 2p, on the other hand, we have j = s + 1p + q < i = s + p + 2 and k = s + 2 + pq < i = s + p + 2 so that all the elements are different again, just like in the case of even r. b. For s + 1 ≤ p ≤ 3s + 1: 6(s + 1) elements based on the generic element

[s + p + 2, q + 1, 2s -q + 2; 3s + 1 -p], q = 0, . . . , 2s + 1, (1.34)
and its permutations.

Therefore, the total number of elements of Type ( 1) is s p=0

(6p + 3) + 3s+1 p=s+1 6(s + 1) = 3(s + 1)(5s + 3).
2. Type (2) again requests ℓ = 0 and leads to 3(2s + 1p) elements based on the generic element [4s + 4 + p, q + 1, 2s + 1pq; 0], q = 0, . . . , 2sp, (1.35) and its permutations, so that the total number of elements of Type ( 2) is 2s p=0

3(2s + 1p) = 3(s + 1)(2s + 1).

3. For Type (3) we again assume that k = 1, note ℓ ≤ d -2r -2 = 2s -1 and obtain 3(p + 1) elements from the generic element

[2s + 3 + p -q, 2s + 3 + q, 1; 2s -1 -p], q = 0, . . . , p (1.36)
and its permutations totalling up to 2s-1

p=0 3(p + 1) = 3s(2s + 1)
elements of Type (3).

Having completed the table in (1.23) the claim (1.24) follows from summing up the columns of the table.

Theorem 1 Σ r 3r is a basis of S r 3r .

We prove Theorem 1 by verifying in Proposition 3 that the elements of Σ r 3r are linearly independent. Since we already know from Proposition 2 that #Σ r 3r = dim S r 3r , this indeed shows that they are a basis of the spline space. Consequently, Σ r 3r spans the space of all simplex splines contained S r 3r . We give an independent proof of this fact in Proposition 8 in the appendix as it may be of independent interest and motivates the classification of the simplex splines according to the three types.

To prove linear independence, we need the following technical tool concerning particular derivatives of simplex splines.

Lemma 3 For i > 0 We have that

D p 1 -x [i, j, k; ℓ](x) = d [i -1, j, k; ℓ](x) -[i, j, k; ℓ](x) , (1.37)
while

D p 1 -x [0, j, k; ℓ](x) = d 3 [0, j, k; ℓ -1](x) -[0, j -1, k; ℓ](x) -[0, j, k -1; ℓ](x) -[0, j, k; ℓ](x) . (1.38)
Since for i, j, k, l ∈ N 0 , and x ∈ T with barycentric coordinates β 1 , β 2 , β 3 with respect to p 1 , p 2 , p 3 , we have

[0, j + 1, k + 1; ℓ + 1](x) = ( j + k + ℓ)! j!k!ℓ! (β 2 -β 1 ) j (β 3 -β 1 ) k (3β 1 ) ℓ , x ∈ T 1 , [i + 1, 0, k + 1; ℓ + 1](x) = (i + k + ℓ)! i!k!ℓ! (β 1 -β 2 ) i (β 3 -β 2 ) k (3β 2 ) ℓ , x ∈ T 2 , [i + 1, j + 1, 0; ℓ + 1]x) = (i + j + ℓ)! i! j!ℓ! (β 1 -β 3 ) i (β 2 -β 3 ) j (3β 3 ) ℓ , x ∈ T 3 , [i + 1, j + 1, k + 1; 0](x) = (i + j + k)! i! j!k! β i 1 β j 2 β k 3 = B d i jk (x), x ∈ T , (1.39)
we observe that the formula (1.38) in fact corresponds to taking a partial derivative with respect to β 1 in (1.39).

Proof Write x = α j p j , with α j = 1. Then the derivative formula and the recurrence yield that

D p 1 -x [i, j, k; ℓ](x) = d (1 -α 1 ) [i -1, j, k; ℓ](x) -α 2 [i, j -1, k; ℓ](x) -α 3 [i, j, k -1; ℓ](x) = d [i -1, j, k; ℓ](x) -d α 1 [i -1, j, k; ℓ](x) + α 2 [i, j -1, k; ℓ](x) + α 3 [i, j, k -1; ℓ](x) = d [i -1, j, k; ℓ](x) -[i, j, k; ℓ](x) ,
which is (1.37). For the second identity we note that p T = 1 3 ( p 1 + p 2 + p 3 ) implies p 1 = 3p Tp 2p 3 and hence, writing x = α T p T + α 2 p 2 + α 3 p 3 with α j = 1, we obtain that

D p 1 -x [0, j, k; ℓ](x) = d (3 -α T ) [0, j, k; ℓ -1](x) -(1 + α 2 ) [0, j -1, k; ℓ](x) -(1 + α 3 ) [0, j, k -1; ℓ](x)
which can be recombined as above to yield (1.38).

Proposition 3

The elements of Σ r 3r are linearly independent. Proof Denote by

I r := I(Σ r 3r ) = (i, j, k, ℓ) : [i, j, k; ℓ] ∈ Σ r 3r (1.40)
the set of all knot multiplicities of splines in Σ r 3r . Now assume that there exist coefficients a i jkℓ such that

s := (i, j,k,ℓ)∈I r a i jkℓ [i, j, k; ℓ] = 0.
On the boundary p 2 , p 3 we have [i, j, k; ℓ] 0 if and only if i = 1 and ℓ = 0 and the splines [1, j, k; 0] reduce to univariate Bernstein polynomials that can be classified as follows:

1. Type (1): [1, µ -ℓ, ℓ + r + 1; 0], 2. Type (2): [1, µ + 1, r + 1; 0], . . . , [1, µ + r + 1, 1; 0], 3. Type (3): [1, µ -ℓ -m, ℓ + r + 1 + m; 0], 1 ≤ m ≤ µ-r-1 2 -ℓ,
together with their symmetric elements where j and k are interchanged. Recall that if these symmetries coincide they are considered as only one element in Σ r 3r . These Bernstein polynomials are linearly independent within the same type by construction and between types since the maximal multiplicity is = µℓ in Type (1), > µℓ in Type (2) and < µℓ in Type [START_REF] Cao | A finite element framework based on bivariate simplex splines on triangle configurations[END_REF]. Therefore, restricting s to the boundary p 2 , p 3 , it follows that a i jkℓ = 0 for (i, j, k, ℓ) = (1, j, k, 0). Considering the other boundaries of T , we can thus conclude that a i jkℓ = 0 whenever min{i, j, k} = 1 and ℓ = 0.

Starting from this observation, we prove by induction on m = 1, 2, . . . that

a i jkℓ = 0, min{i, j, k} + ℓ = m, (1.41)
where the case m = 1 has been treated in the first part of this proof. We will treat the case m = 2 explicitly as the general procedure will become clear by then. We assume that i is the minimal value, consider the identity

0 = D p 1 -x s(x), x ∈ p 2 , p 3
and find by Lemma 3 that there are only two types of splines which are nonzero on the boundary. The first is

D p 1 -x [2, j, k; 0](x) = [1, j, k; 0](x) -[2, j, k; 0](x)
which coincides with [1, j, k; 0](x) on p 2 , p 3 . The second is

D p 1 -x [1, j, k; 1](x) = [0, j, k; 1](x) -[1, j, k; 1](x),
which by (1.39) coincides with [0, j, k; 1] = [1, j, k; 0] on p 2 , p 3 . Again, these splines are linearly independent within the types by construction and across the types by the different values of the maximal multiplicity, and (1.41) for m = 2 follows by considering all three faces of T by a symmetry argument. The general induction step proceeds in exactly this way by assuming that i = min{i, j, k} and i + ℓ = m + 1 and applying (1.37) i times we find from (1.39) with x = (β 1 , β 2 , β 3 ) and ℓ = m + 1i

D i p 1 -x [i, j, k; ℓ](x) = K [0, j, k; ℓ](x) = ( j + k + ℓ)! j!k!ℓ! (β 2 -β 1 ) j (β 3 -β 1 ) k (3β 1 ) ℓ
for some connstant K. Differentiating ℓ times with respect to β 1 and setting β 1 = 0 we find

∂ ℓ ∂ β ℓ 1 D i p 1 -x [i, j, k; ℓ](x) = K β j 2 β k 3 ,
where again K is some nonzero constant. where the admissible values for j, k lead to linearly independent polynomials on the boundary. This completes the proof of Proposition 3. on the "lower" edge of the triangle where alle other simplex splines from our list vanish. Note that the three types are then distinguished by whether the maximum equals µ = 7 (the element on the left), exceeds this values (the three elements in the middle) or is strictly smaller (the element on the right). Therefore, the coefficients of these splines and all their symmetries have to vanish which deals with the first "ring" of coefficients on the boundary.

Applying the differential operator D p 1 -x once, gives us the Bernstein polynomials , respectively. Observe that in each step of the process always Bernstein polynomials of the same fixed degree are considered on the boundary.

Marsden-like Identity, Domain points, Stability and C r -connection

A Marsden-like identity allows us to derive explicit formulas for the representation of polynomials of degree up to 3r in the simplex spline basis of S r 3r ( ). Even if we have identified a basis of S r 3r ( ), it turns out that for the partition of unity and the Marsden-like identity alone, i.e., for the generation of polynomials and espcially the constant function, we do not need all the elements of that basis. Hence, instead of having redundancy or a null component in the sums, we remove one further element according to the following definition.

Definition 3 Depending on the parity of r, we define the following sets and spaces:

1. If r = 2s + 1 with s ≥ 0, then Σr 3r := Σ r 3r , Sr 3r ( ) = S r 3r ( ).

If r = 2s then

Σr 3r := Σ r 3r \ (s + 1, s + 1, s + 1; 3s) , Sr 3r ( ) := span Σr 3r . In analogy with (1.40), the set of indices (i, j, k, ℓ) of the basis Σr 3r will be written as Īr . Moreover, for the Bernstein polynomials in Σr 3r (cf. Figure 1.3), we define the index sets

I r removed := {(i, j, k) : i + j + k = 3r, max{i, j, k} ≤ 2r -1, min{i, j, k} ≥ 1}, I r Bernstein := {(i, j, k) : i, j, k ∈ N 0 , i + j + k = 3r } \ I r removed .
(1.42)

It immediately follows that

#I r Bernstein = 3 2 r(r + 5), #I r removed = 3r(r -1) + 1 #I r Bernstein + #I r removed = 3r + 2 2 = dim P 3r . (1.43)
We now consider the cases r = 2, 3 in more detail. We prove, using knot insertion, a Marsden-like identity for r = 2, and state it in a form valid for any r ≥ 1. It is verified symbolically for r ≥ 6. Shouldn't it be r ≤ 6?

1.5.1 The C 2 elements, Σ2

+ 6β 3 β 5 1 -15β 2 3 β 4 1 + 20β 3 3 β 3 1 {4, 1, 3, 1} β 6 1 -6β 3 β 5 1 + 15β 2 3 β 4 1 {4, 2, 2, 1} 4β 6 1 -12β 2 β 5 1 -12β 3 β 5 1 + 60β 2 β 3 β 4 1 {3, 3, 1, 2} 9β 6 1 -36β 2 β 5 1 + 45β 2 2 β 4 1 {2, 4, 2, 1} 8β 6 1 -36β 2 β 5 1 -12β 3 β 5 1 + 60β 2 2 β 4 1 + 60β 2 β 3 β 4 1 -40β 3 2 β 3 1 -120β 2 2 β 3 β 3 1 + 120β 3 2 β 3 β 2 1 {1, 3, 3, 2} S 6 31 (1) {2, 2, 4, 1} 8β 6 1 -12β 2 β 5 1 -36β 3 β 5 1 + 60β 2 3 β 4 1 + 60β 2 β 3 β 4 1 -40β 3 3 β 3 1 -120β 2 β 2 3 β 3 1 + 120β 2 β 3 3 β 2 1 {3, 1, 3, 2} 9β 6 1 -36β 3 β 5 1 + 45β 2 3 β 4 1 {3, 2, 2, 2} 36β 6 1 -72β 2 β 5 1 -72β 3 β 5 1 + 180β 2 β 3 β 4 1 {2, 3, 2, 2} -72β 6 1 + 216β 2 β 5 1 + 108β 3 β 5 1 -180β 2 2 β 4 1 -360β 2 β 3 β 4 1 + 360β 2 2 β 3 β 3 1 {2, 2, 3, 2} -72β 6 1 + 108β 2 β 5 1 + 216β 3 β 5 1 -180β 2 3 β 4 1 -360β 2 β 3 β 4 1 + 360β 2 β 2 3 β 3 1 (1.46)
where the first entry denotes the position of the spline with respect to aforementioned ordering, the second the multiplicity of the knots and the third one the explicit expression on T 1 , using

S 6 24 (1) = -10β 6 1 + 36β 2 β 5 1 + 24β 3 β 5 1 -45β 2 2 β 4 1 -15β 2 3 β 4 1 -90β 2 β 3 β 4 1 + 20β 3 2 β 3 1 + 60β 2 β 2 3 β 3 1 + 120β 2 2 β 3 β 3 1 -90β 2 2 β 2 3 β 2 1 -60β 3 2 β 3 β 2 1 + 60β 3 2 β 2 3 β 1 S 6 25 (1) = -10β 6 1 + 24β 2 β 5 1 + 36β 3 β 5 1 -15β 2 2 β 4 1 -45β 2 3 β 4 1 -90β 2 β 3 β 4 1 + 20β 3 3 β 3 1 + 120β 2 β 2 3 β 3 1 + 60β 2 2 β 3 β 3 1 -60β 2 β 3 3 β 2 1 -90β 2 2 β 2 3 β 2 1 + 60β 2 2 β 3 3 β 1 S 6 31 (1) = 90β 6 1 -216β 2 β 5 1 -216β 3 β 5 1 + 135β 2 2 β 4 1 + 135β 2 3 β 4 1 + 540β 2 β 3 β 4 1 -360β 2 β 2 3 β 3 1 -360β 2 2 β 3 β 3 1 + 270β 2 2 β 2 3 β 2 1 (1.47)
Here, S 6 k (1) is used to indicate the restriction of S 6 k to the triangle T 1 . On T 2 we have, in the same fashion, Proof By the Marsden-like identity we have Next, we define, as usually, the ∞-norm condition number of the Σr 3r bases for Sr 3r ( ) by

{4, 3, 1, 1} -β 6 2 + 6β 1 β 5 2 -15β 2 1 β 4 2 + 20β 3 1 β 3 2 {3, 4, 1, 1} β 6 2 -6β 1 β 5 2 + 15β 2 1 β 4 2 {1, 4, 3, 1} β 6 2 -6β 3 β 5 2 + 15β 2 3 β 4 2 {1, 3, 4, 1} -β 6 2 + 6β 3 β 5 2 -15β 2 3 β 4 2 + 20β 3 3 β 3 2 {3, 1, 4, 1} S 6 26 (2) {4, 1, 3, 1} S 6 27 (2) {4, 2, 2, 1} 8β 6 2 -36β 1 β 5 2 -12β 3 β 5 2 + 60β 2 1 β 4 2 + 60β 1 β 3 β 4 2 -40β 3 1 β 3 2 -120β 2 1 β 3 β 3 2 + 120β 3 1 β 3 β 2 2 {3, 3, 1, 2} 9β 6 2 -36β 1 β 5 2 + 45β 2 1 β 4 2 {2, 4, 2, 1} 4β 6 2 -12β 1 β 5 2 -12β 3 β 5 2 + 60β 1 β 3 β 4 2 {1, 3, 3, 2} 9β 6 2 -36β 3 β 5 2 + 45β 2 3 β 4 2 {2, 2, 4, 1} 8β 6 2 -12β 1 β 5 2 -36β 3 β 5 2 + 60β 2 3 β 4 2 + 60β 1 β 3 β 4 2 -40β 3 3 β 3 2 -120β 1 β 2 3 β 3 2 + 120β 1 β 3 3 β 2 2 {3, 1, 3, 2} S 6 33 (2) {3, 2, 2, 2} -72β 6 2 + 216β 1 β 5 2 + 108β 3 β 5 2 -180β 2 1 β 4 2 -360β 1 β 3 β 4 2 + 360β 2 1 β 3 β 3 2 {2, 3, 2, 2} 36β 6 2 -72β 1 β 5 2 -72β 3 β 5 2 + 180β 1 β 3 β 4 2 {2, 2, 3, 2} -72β 6 2 + 108β 1 β 5 2 + 216β 3 β 5 2 -180β 2 3 β 4 2 -360β 1 β 3 β 4 2 + 360β 1 β 2 3 β 3 2 (1.48) where S 6 26 (2) = -10β 6 2 + 24β 1 β 5 2 + 36β 3 β 5 2 -15β 2 1 β 4 2 -45β 2 3 β 4 2 -90β 1 β 3 β 4 2 + 20β 3 3 β 3 2 + 120β 1 β 2 3 β 3 2 + 60β 2 1 β 3 β 3 2 -60β 1 β 3 3 β 2 2 -90β 2 1 β 2 3 β 2 2 + 60β 2 1 β 3 3 β 2 S 6 27 (2) = -10β 6 2 + 36β 1 β 5 2 + 24β 3 β 5 2 -45β 2 1 β 4 2 -15β 2 3 β 4 2 -90β 1 β 3 β 4 2 + 20β 3 1 β 3 2 + 60β 1 β 2 3 β 3 2 + 120β 2 1 β 3 β 3 2 -90β 2 1 β 2 3 β 2 2 -60β 3 1 β 3 β 2 2 + 60β 3 1 β 2 3 β 2 S 6 33 (2) = 90β 6 2 -216β 1 β 5 2 -216β 3 β 5 2 + 135β 2 1 β 4 2 + 135β 2 3 β 4 2 + 540β 1 β 3 β 4 2 -360β 1 β 2 3 β 3 2 -360β 2 1 β 3 β 3 2 + 270β 2 1 β 2 3 β 2 2 , (1
β m = (i, j,k,ℓ)∈ Ī2 ∂ ∂u m ρ i jkℓ (1, 1, 1)S 6 i jkℓ (β), m = 1,
κ d,∞ (T ) := max c 0 b T c L ∞ (T) c ∞ max c 0 c ∞ b T c L ∞ (T) ,
where d = 3r and b T c := (i, j,k,ℓ)∈ Īr c i jkℓ S d i jkℓ ∈ Sr d ( ). This number turns out to be bounded by a moderate number independently of the shape of the basis triangle T .

Proposition 5 (Stability)

For any triangle T we have κ 6,∞ (T ) < 1350.

Proof Since the S 6

i jkℓ , (i, j, k, ℓ) ∈ Ī2 form a nonnegative partition of unity it follows that for any

x ∈ |b(x) T c| ≤ (i, j,k,ℓ)∈ Ī2 |c i jkℓ | |S 6 i jkℓ (x)| ≤ c ∞ (i, j,k,ℓ)∈ Ī2 S 6 i jkℓ (x) = c ∞ with an equality if all the c i jkℓ = 1, so that max c 0 b T c L ∞ (T) / c ∞ = 1.
To bound the second part of κ 6,∞ (T ), we consider a spline s ∈ S2 6 interpolating given data at the 36 domain points. Using the ordering of domain points p * m and basis functions S 6 n , m, n = 1, . . . , 36, shown in Figure 1.5 we obtain a linear system Finally, we consider next the problem of obtaining C 2 -continuity across an edge between two adjacent triangles in a global triangulation, making use of the local properties of our simplex spline basis.

Ac =        S 6 1 ( p * 1 ) • • • S 6 36 ( p * 1 ) . . . . . . S 6 1 ( p * 36 ) • • • S 6 36 ( p * 36 )               c 1 . . . c 36        =        s( p * 1 ) . . . s( p * 36 )        = s *   1 2 3 4 5 6
A -1 = B 1 0 B 2 B 3 , and
B 1 = A -1 1 , B 3 = A -1 3 , B 2 = -B 3 A 2 B 1 . We compute A -1 ∞ = 12
We begin with a technical observation.

Proposition 6 Let s = (i, j,k,ℓ)∈ Ī2 c i jkℓ S 6 i jkℓ or s = (i, j,k,ℓ)∈ Ī2 cijkℓ S6 i jkℓ ,
respectively, be defined on the triangle T := p 1 , p 2 , p 3 (resp. T := p 1 , p 2 , p3 ). We suppose that p3 =

λ 1 p 1 + λ 2 p 2 + λ 3 p 3 with λ 1 + λ 2 + λ 3 = 1. The function s + = s on T s on T is C r with r ≤ 2 if and only if        c0 . . . cr        =        C 00 . . . C 0r . . . . . . C r0 . . . C rr               c 0 . . . c r        , (1.58)
where the c0 , c1 , c2 are the coefficients of s corresponding to (1.57), and the matrices C mn are defined by

C 00 = I ∈ R 7×7 , C 01 = 0 ∈ R 7×6 , C 02 = 0 ∈ R 7×5 , (1.59)
C 10 =           λ 1 λ 2 0 . . . 0 0 λ 1 λ 2 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 λ 1 λ 2           ∈ R 6×7 , C 11 = λ 3 I ∈ R 6×6 , C 12 = 0 ∈ R 6×5 , (1.60)
C 20 =           λ 2 1 2λ 1 λ 2 λ 2 2 0 0 0 0 0 λ 2 1 /2 λ 1 (1/2 + λ 2 ) λ 2 (1 + λ 2 )/2 0 0 0 0 0 λ 1 (1 + λ 1 )/3 (λ 1 + λ 2 + 2λ 1 λ 2 )/3 λ 2 (1 + λ 2 )/3 0 0 0 0 0 λ 1 (1 + λ 1 )/2 λ 2 (1/2 + λ 1 ) λ 2 2 /2 0 0 0 0 0 λ 2 1 2λ 1 λ 2 λ 2 2           ∈ R 5×7 , C 21 =          2λ 1 λ 3 2λ 2 λ 3 0 0 0 0 0 λ 1 λ 3 λ 3 (1/2 -λ 3 /2 + λ 2 ) 0 0 0 0 0 λ 3 (1 -λ 3 + 2λ 1 )/3 λ 3 (1 -λ 3 + 2λ 2 )/3 0 0 0 0 0 λ 3 (1/2 -λ 3 /2 + λ 1 ) λ 2 λ 3 0 0 0 0 0 2λ 1 λ 3 2λ 2 λ 3          ∈ R 5×6 , C 22 = λ 2 3 I ∈ R 5×5 . (1.61)
Proof We begin by the C r -continuity for the Bézier surfaces using the Bernstein basis. Let σ = ν+µ+κ=6 γ νµκ B 6 νµκ ∈ P 6 be defined on the triangle T := p 1 , p 2 , p 3 (respectively σ = ν+µ+κ=6 γνµκ B6 νµκ ∈ P 6 on T := p 1 , p 2 , p3 ) where B6

νµκ are the Bernstein polynomials with barycentric variables β1 , β2 , β3 .

We recall, see [START_REF]Spline Functions on Triangulations[END_REF]Theorem 2.28], that the function σ + = σ on T σ on T is C r if and only if

[ γm ] m=0,...,r = [Γ mn ] m,n=0,...,r [γ n ] n=0,...,r , (1.62)
where

γ 0 = [γ 600 , γ 510 , γ 420 , γ 330 , γ 240 , γ 150 , γ 160 ] T , γ 1 = [γ 501 , γ 411 , γ 321 , γ 231 , γ 141 , γ 051 ] T , γ 2 = [γ 402 , γ 312 , γ 222 , γ 132 , γ 042 ] T ,
similarly for the γm and the matrices are defined by

Γ 00 = I ∈ R 7×7 , Γ 01 = 0 ∈ R 7×6 , Γ 02 = 0 ∈ R 7×5 , (1.63)
Γ 10 =           λ 1 λ 2 0 . . . 0 0 λ 1 λ 2 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 λ 1 λ 2           ∈ R 6×7 , Γ 11 = λ 3 I ∈ R 6×6 , Γ 12 = 0 ∈ R 6×5 , (1.64)
Γ 20 =            λ 2 1 2λ 1 λ 2 λ 2 2 0 . . . 0 0 λ 2 1 2λ 1 λ 2 λ 2 2 . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 λ 2 1 2λ 1 λ 2 λ 2 2            ∈ R 5×7 , Γ 21 =           2λ 1 λ 3 2λ 2 λ 3 0 . . . 0 0 2λ 1 λ 3 2λ 2 λ 3 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 2λ 1 λ 3 2λ 2 λ 3           ∈ R 5×6 , Γ 22 = λ 2 3 I ∈ R 5×5 .
(1.65)

See also Figure 1.4.

Consider now the function s + = s on T s on T on T ∪ T . To study the C r -continuity through the edge p 1 p 2 , it is sufficient to consider s on T 3 and s on T3 . On T 3 (resp. T3 ), any S 6 i jkℓ of the basis (resp. S6 i jkℓ ) is a polynomial of degree at most 6,

S 6 i jkℓ = ν+µ+κ=6 w i jkℓ νµκ B 6 νµκ (resp. S6 i jkℓ = ν+µ+κ=6 w i jkℓ νµκ B6 νµκ )
. So that s| T 3 and s| T3 can also be written in the Bernstein bases

s| T 3 = (i, j,k,ℓ)∈ Ī2 c i jkℓ S 6 i jkℓ | T 3 = ν+µ+κ=6 γ νµκ B 6 νµκ with γ νµκ = (i, j,k,ℓ)∈ Ī2 c i jkℓ w i jkℓ νµκ s| T3 = (i, j,k,ℓ)∈ Ī2 cijkℓ S6 i jkℓ | T 3 = ν+µ+κ=6 γνµκ B6 νµκ with γνµκ = (i, j,k,ℓ)∈ Ī2 cijkℓ wijkℓ νµκ .
From (1.56), we deduce the components wijkℓ νµκ for κ = 0, 1, 2 and we put forward the corresponding components

s| T 3 =
+O(β 3 3 )
and a similar expression for s| T3

The conditions for the regularity C 0 of s + , (1.59), is a consequence of (1.63) and similarly for C 1 with also (1.60) coming from (1.64). To obtain C 2 , we add the conditions (1.65). They can be rewritten

c5130 = c 7110 λ 2 1 + 2c 6210 λ 1 λ 2 + c 5310 λ 2 2 + 2c 5220 λ 2 λ 3 + c 5130 λ 2 3 + 2c 6120 λ 1 λ 3 , 2 c4221 -c4311 = c 6210 λ 2 1 + 2c 5310 λ 1 λ 2 + c 4410 λ 2 2 + 2c 4311 λ 2 λ 3 +(2c 4221 -c 4311 )λ 2 3 + 2c 5220 λ 1 λ 3 , 3 c3312 -c4311 -c3411 = c 5310 λ 2 1 + 2c 4410 λ 1 λ 2 + c 3510 λ 2 2 + 2c 3411 λ 2 λ 3 +(3c 3312 -c 4311 -c 3411 )λ 2 3 + 2c 4311 λ 1 λ 3 2 c2421 -c3411 = c 4410 λ 2 1 + 2c 3510 λ 1 λ 2 + c 2610 λ 2 2 + 2c 2520 λ 2 λ 3 +(2c 2421 -c 3411 )λ 2 3 + 2c 3411 λ 1 λ 3 , c1530 = c 3510 λ 2 1 + 2c 2610 λ 1 λ 2 + c 1710 λ 2 2 + c 1620 λ 2 λ 3 + c 1530 λ 2 3 + 2c 2520 λ 1 λ 3 .
We already have c4311 and c3411 in (1.60) so that we deduce that the last components have to satisfy (1.61).

Several examples have been considered for scattered data on the CT-split, see for example [START_REF] Farin | A modified Clough-Tocher interpolant[END_REF][START_REF] Mann | Cubic precision Clough-Tocher interpolation[END_REF]. Here, we consider a surface on two triangles, see Figure 1.6. With the 18 conditions from Proposition 6, we obtain a C 2 surface on the two triangles.

For the stability, the computation is similar to the proof of Proposition 5, except that since S 9 3324 , S 9 3234 , S 9 2334 have the same domain points (1/3, 1/3, 1/3). We replace this triple point by the three points (3, 3, 1)/7, (3, 1, 3)/7, (1, 3, 3)/7 and find Let s = (i, j,k,ℓ)∈ Ī3 c i jkℓ S 9 i jkℓ (respectively s = (i, j,k,ℓ)∈ Ī3 cijkℓ S9 i jkℓ ) be defined on the triangle T := p 1 , p 2 , p 3 (resp. T := p 1 , p 2 , p3 ). We suppose that p3 = λ 1 p 1 + λ 2 p 2 + λ 3 p 3 with λ 1 + λ 2 + λ 3 = 1.

κ 9,∞ (T ) := max c 0 b T c L ∞ (T) c ∞ max c 0 c ∞ b T c L ∞ (T) ≃ 159
With an extension of the notations of Proposition 6, we define For r ≤ 3, we connect with smoothness C r two adjacent triangles in the following proposition. Proposition 7 Let s = (i, j,k,ℓ)∈ Ī2 c i jkℓ S 9 i jkℓ (respectively s = (i, j,k,ℓ)∈ Ī2 cijkℓ S9 i jkℓ ) be defined on the triangle T := p 1 , p 2 , p 3 (resp. T := p 1 , p 2 , p3 ). We suppose that p3 = λ 1 p 1 + λ 2 p 2 + λ 3 p 3 with λ 1 + λ 2 + λ 3 = 1.

The function s + = s on T s on T is C r with r ≤ 3 if and only if

[ cm ] m=0,...,r = [C mn ] m,n=0,...,r [c n ] n=0,...,n , C mn =∈ R (10-m)×(10-n) (1.72)
where the nonzero submatrices or components are written below in (1.73) (1.74), (1.75) and (1.76).

The proof of the proposition is a reproduction of the one of Proposition 6, firstly by connecting two polynomials written in the two corresponding Bernstein basis, then computing the Bernstein polynomials in the Simplex-Splines basis.

C 00 = I ∈ R 10×10 , (1.73)
C 10 =           λ 1 λ 2 0 . . . 0 0 λ 1 λ 2 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 λ 1 λ 2           ∈ R 9×10 , C 11 = λ 3 I ∈ R 9×9 , (1.74)
C 20 (1, 1 : 3) = λ 2 1 2λ 1 λ 2 λ 2 2 C 20 (2, 2 : 4) = λ 2 1 2λ 1 λ 2 λ 2 2 C 20 (3, 3 : 5) = λ 2 1 2 λ 1 λ 2 + λ1
+ λ3 2 λ 2 λ 3 C 21 (7, 7 : 8) = 2λ 1 λ 3 2λ 2 λ 3 C 21 (8, 8 : 9) = 2λ 1 λ 3 2λ 2 λ 3 , C 22 = λ 2 3 I ∈ R 8×8 , (1.75)
C 30 (1, 1 : 4) = λ 3 1 3λ 2 1 λ 2 3λ 1 λ 2 2 λ 3 2 C 30 (2, 2 : 5) = λ 3 1 2 3λ 2 1 λ2 2 + λ 2 1 2 3λ1 λ 2 2 2 + λ 1 λ 2 λ 3 2 2 + λ 2 C 30 (3, 3 : 6) = λ 3 1 6 + λ 2 1 6 λ 2 1 λ2 2 + λ1 λ2 3 + λ 2 1 3 + λ1
+ λ 1 λ 2 3λ1 λ 2 2 2 + λ 2 2 2 λ 3 2 2 C 30 (7, 7 : 10) = λ 3 1 3λ 2 1 λ 2 3λ 1 λ 2 2 λ 3 2 (1.76) C 31 (1, 1 : 2) = 3λ 2 1 λ 3 6λ 1 λ 2 λ 3 C 31 (1, 3) = 3λ 2 2 λ 3 C 31 (2, 2 : 3) = 3λ 2 1 λ3 2 3λ 1 λ 2 λ 3 + λ 1 λ 3 C 31 (2, 4) = λ 3 3 -3λ2 λ 2 3 +3λ 2 2 λ3 2 + 2λ2 λ3-λ 2 3 2 C 31 (3, 3 : 4) = λ 2 1 λ3 2 + λ1 λ3 3 2λ 3 3 -3λ2 λ 2 3 -3λ1 λ 2 3 +6λ1 λ2 λ3 6 + 2λ2 λ3-λ 2 3 6 + 2λ1 λ3-λ 3 + λ3 6 C 31 (3, 5) = λ 3 3 -3λ2 λ 2 3 +3λ 2 2 λ3 6 + 2λ2 λ3-λ 2 3 3 + λ3 6 C 31 (4, 4 : 5) = λ 3 3 -3λ1 λ 2 3 +3λ 2 1 λ3 9 + 2(2λ1 λ3-λ 2 3 ) 9 + λ3 9 2λ 3 3 -3λ2 λ 2 3 -3λ1 λ 2 3 +6λ1 λ2 λ3 9 + 2(2λ2 λ3-λ 2 3 ) 9 + 2(2λ1 λ3-λ 2 3 ) 9 + 2λ3 9 C 31 (4, 6) = λ 3 3 -3λ2 λ 2 3 +3λ 2 2 λ3 9 + 2(2λ2 λ3-λ 2 3 ) 9 + λ3 9 C 31 (5, 5 : 6) = λ 3 3 -3λ1 λ 2 3 +3λ 2 1 λ3 6 + 2λ1 λ3-λ 2 3 3 + λ3 6 2λ 3 3 -3λ2 λ 2 3 -3λ1 λ 2 3 +6λ1 λ2 λ3 6 + 2λ2 λ3-λ 2 3 3 + 2λ1 λ3-λ 6 + λ3 6 C 31 (5, 7) = λ 2 2 λ3 2 + λ2 λ3 3 C 31 (6, 6 : 7) = λ 3 3 -3λ1 λ 2 3 +3λ 2 1 λ3 2 + 2λ1 λ3-λ 2 3 2 C 31 (6, 8) = 3λ 1 λ 2 λ 3 + λ 2 λ 3 3λ 2 2 λ3 2 C 31 (7, 7 : 8) = 3λ 2 1 λ 3 6λ 1 λ 2 λ 3 C 31 (7, 9) = 3λ 2 2 λ 3 C 32 (1, 1 : 2) = 3λ 1 λ 2 3 3λ 2 λ 2 3 C 32 (2, 2 : 3) = 3λ1 λ 2 3 2 6λ2 λ 2 3 -2λ 3 3 2 + λ 2 3 C 32 (3, 3 : 4) = 6λ1 λ 2 3 -2λ 3

Higher smoothness

For any r ≥ 1, we have built a basis for S r 3r ( ). For even values of r, we remove one of the elements to obtain the partition of unity and the Marsden-like identity. The following theorem is proved for r = 1, 2, 3 and symbolically for r = 4, 5, 6. It is a conjecture for r > 6.

Theorem 4 (The barycentric Marsden-like identity for degree 3r)

For r ∈ N, d = 3r, u i , β i ∈ R, with β i ≥ 0, i = 1, 2, 3, and β 1 + β 2 + β 3 = 1 we have

(u 1 β 1 + u 2 β 2 + u 3 β 3 ) d = (i, j,k,ℓ)∈ Īr ρ i jkℓ (u 1 , u 2 , u 3 )S d i jkℓ (β 1 , β 2 , β 3 ),
where the index set Īr is given in Definition 3 and

S 3r i jkℓ := c i jkℓ i jkℓ , where c i jkℓ :=          1, if ℓ = 0, 2/3 if max(ǫ 1 , ǫ 2 , ǫ 3 ) = 1/2, 1/3 otherwise, ρ i jkℓ (u 1 , u 2 , u 3 ) := u i-1 1 u j-1 2 u k-1 3 ūλ-1 123 (ǫ 1 u 1 + ǫ 2 u 2 + ǫ 3 u 3 + δ 0ℓ ) 1-δ 0ℓ .
(1.77)

Here, λ := max(ℓ, 1), ū123 := (u 1 + u 2 + u 3 )/3, ν := max(i, j, k), and

ǫ 1 := 1/γ 1 , i < ν & j 1 & k 1, 0, otherwise , γ 1 := 2, i 1 & j k, 1, otherwise , ǫ 2 := 1/γ 2 , j < ν & i 1 & k 1, 0, otherwise , γ 2 := 2, j 1 & i k, 1, otherwise , ǫ 3 := 1/γ 3 , k < ν & i 1 & j 1, 0, otherwise , γ 3 := 2, k 1 & i j, 1 , otherwise .
(1.78)

Proof For r = 1, 2, 3 this is an alternative way of formulating Theorem 5 in [START_REF] Lyche | Simplex-splines on the Clough-Tocher element[END_REF] for r = 1, and Theorems 2,3. To see this consider first r = 1. In [START_REF] Lyche | Simplex-splines on the Clough-Tocher element[END_REF] it was shown that hence belong to Σ r 3r , thus to M ℓ+1 . The same holds true for [i -1, j, k; ℓ + 1] if j = i or k = i. In the remaining case, i > max{ j, k}, we decompose

[i -1, j, k; ℓ + 1] = 1 3 [i -2, j, k; ℓ + 2] + [i -1, j -1, k; ℓ + 2] + [i -1, j, k -1; ℓ + 2]
and note that the last two functions on the right hand side again belong to Σ r 3r , thus to M ℓ+2 ⊂ M ℓ+1 , so that again one only has to look at the first term. This procedure is repeated n := i -max{ j, k} times when in = max{ j, k} and thus [in, j, k; ℓ + n] ∈ M ℓ+n ⊂ M ℓ+1 .

To advance the induction, suppose that the result has been verified for some defect ν ≥ 1 and again apply the decomposition (1.80) to a spline in Σ r 3r of defect ν +1. Suppose again without loss of generality that i = max{i, j, k}. We begin by looking at the first element of the decomposition where three things can happen:

1. i = 2, then immediately [i -1, j, k; ℓ + 1] ∈ W ℓ+1 , 2. the defect is ν if i = j or i = k,
and the hypothesis yields that [i -1, j, k; ℓ + 1] ∈ M ℓ+1 + W ℓ+1 , 3. still we have ν+1 which happens if i > max{ j, k}. In this case, we repeat the above argument of n := i-max{ j, k} iterated decompositions until, eventually,

[i -n, j, k; ℓ + n] ∈ M ℓ+n + W ℓ+n ⊂ M ℓ+1 + W ℓ+1 .
For the two other elements of the decomposition, the defect of

[i, j -1, k; ℓ + 1] and [i, j, k -1; ℓ + 1] is µ -(ℓ + 1) -max{i, j -1, j, k -1, k} = µ -(ℓ + 1) -i = (µ -ℓ -max{i, j, k}) -1 = ν + 1 -1 = ν
and the induction hypothesis yields that

[i, j -1, k; ℓ + 1] [i, j, k -1; ℓ + 1] ∈ M ℓ+2 + W ℓ+2 ⊂ M ℓ+1 + W ℓ+1 .
This advances the induction hypothesis and completes the proof of the lemma.

Proposition 8 Σ r 3r generates the space of all simplex splines in S r 3r .

Proof First recall that, by assumption, µ = 3r + 1r = 2r + 1 and let [i, j, k; ℓ] be one of the simplex spline generating S r 3r ( ), which implies that 1. min{i, j, k, ℓ} ≥ 0 and i + j + k + ℓ = 3r + 3, 2. ℓ > 0 implies max{i, j, k} ≤ µℓ = 2r + 1ℓ for C r -smoothness, 3. if min{i, j, k} = 0, then ℓ > 0 since the Bernstein polynomials with zero multiplicity correspond to distributions defined only on the boundary of the simplex.

We distinguish three cases:

1. max{i, j, k} > µℓ, then ℓ = 0 and [i, j, k; ℓ] is already an element of Type (2), 2. max{i, j, k} = µℓ then a. if i := min{i, j, k} = 0 and j := max{i, j, k} = µℓ = 2r + 1ℓ. Since i + j + k + ℓ = 3r + 3, we deduce that k = 2 + r ≥ 2. By knot insertion,

[0, j, k; ℓ] = 1 3 [1, j -1, k; ℓ] + [1, j, k -1; ℓ] + [1, j, k; ℓ -1] .
The first element in this expression, [1, j -1, k; ℓ] is of Type (3) if k < j and of Type (1) if k = j, the second one, [1, j, k -1; ℓ] is always of Type ( 1) and the third one, [1, j, k; ℓ -1], of Type (3). Consequently, [0, j, k; ℓ] ∈ span(Σ r 3r ). b. if min{i, j, k} ≥ 1, then [i, j, k; ℓ] is already an element of Type (1), hence in Σ r 3r . 3. max{i, j, k} < µℓ a. if i := min{i, j, k} = 0, assume that j := max{i, j, k}. Since i + j + k + ℓ = 3r + 3, we deduce that j ≥ k > 2 + r ≥ 2 and knot insertion yields where the three elements on the right hand side are all of Type (3), hence [0, j, k; ℓ] ∈ span(Σ r 3r ). b. if min{i, j, k} = 1, then [i, j, k; ℓ] is immediately an element of Type (3). c. if min{i, j, k} > 1 we refer to Lemma 5 and find that [i, j, k; ℓ] ∈ M ℓ+1 + W ℓ+1 ⊂ M 0 + W 0 .

To finish the proof and to complete the last of these cases, we have to show that M 0 + W 0 ⊆ span(Σ r 3r ). In fact, M 0 ⊂ span(Σ r 3r ) was exactly proved in case 2. To complete the proof, let [i, j, k; ℓ] be a generator of W 0 . If max{i, j, k} > µℓ then once more ℓ = 0 and [i, j, k; ℓ] = [i, j, k; 0] is an element of Type [START_REF] Alfeld | Smooth macro-elements based on Powell-Sabin triangle splits[END_REF]. If, on the other hand, max{i, j, k} = µℓ, then [i, j, k; ℓ] is an element of Type (1), and if max{i, j, k} < µℓ, then [i, j, k; ℓ] is an element of Type (3). Consequently, [i, j, k; ℓ] ∈ Σ r 3r and thus W 0 ⊂ span(Σ r 3r ) which completes the proof that M 0 + W 0 ⊆ span(Σ r 3r ) (1.81) and the proof of the theorem.

A special case of linear independence

Here we give a direct proof of the fact stated in Proposition 4 that the functions in Σ2 6 are linearly independent on . Proof Suppose for some real numbers {c j } that 36 j=1 c j S 6 j (x) = 0 for all x ∈ . We first show that c j = 0 for j = 1, . . . , 18. These corresponds to Bernstein polynomials with domainpoints on the boundary of . Consider the edge p 1 , p 2 corresponding to β 3 = 0. Looking at Figures 1.4, 1.5 we see that i + j ≤ 7 for S 6 8 , . . . , S 6 36 . By the local smoothness property only S 6 1 , . . . , S 6 7 can be nonzero on this edge. Moreover, on the same edge, these functions reduce to linearly independent univariate Bernstein polynomials B 6 i j for i + j = 6. It follows that c 1 = • • • = c 7 = 0. With similar arguments on the edges p 2 , p 3 and p 3 , p 1 we conclude that c j = 0 for j ≤ 18.

The remaining simplex splines S 6 j j = 19, . . . , 36 are located on 3 rings. On ring k we find S 6 j for j = 19, . . . , 27 for k = 1, j = 28, . . . , 33 for k = 2, and j = 34, 35, 36 for k = 3, see Figure 1.5. On the horizontal part of these rings we take partial derivatives of order k with respect to β 3 and evaluate at β 3 = 0. On the parts parallel to p 2 , p 3 we take partial derivatives of order k with respect to β 1 and evaluate at β 1 = 0. Similarly we use β 2 on the last parts. The details are as follows

The horizontal part of the first inner ring contains the functions S

Fig. 1 . 1

 11 Fig. 1.1 The Clough-Tocher split, T = p 1 , p 2 , p 3 , T 1 := p T , p 2 , p 3 , T 2 := p T , p 3 , p 1 and T 3 := p T , p 1 , p 2

Fig. 1 . 2

 12 Fig. 1.2 The degrees of freedom cs r 3r for C 1 , C 2 and C 3 . For the order of a cross boundary derivative D α ni we use a dot for α = 0, and colors red, blue and green for α = 1, 2, 3.

Example 1 Example 2 Example 3 9 (Example 4 Remark 2

 123942 For r = 1, d = 3, µ = 3, dim S 1 3 ( ) = 12 we have 1. Type (1)For r = 2, d = 6, µ = 5, dim S 2 6 ( ) = 37 the different types look as follows: 1. Type (1): 25 elements In the case r = 3, d = 9, µ = 7, dim S 3 For r = 4, d = 12, dim S 4 12 ( ) = 127, the elements look as follows: 1Some Bernstein polynomials are not of any of the three types and thus are not in the basis. A Bernstein polynomial in the basis can be of Type (1), (2) or (3), see Figure 1.3 and the previous examples.

Fig. 1 . 3

 13 Fig. 1.3 Domain points of the Bernstein polynomials used in the CTS basis, left: r = 1, d = 3, middle: r = 2, d = 6, right: r = 3, d = 9. "Removed" means that the respective Bernstein polynomials are of none of the Types (1), (2) or (3).

Example 5

 5 We illustrate the elimination procedure of Proposition 3 for the case r = 3. The elimination proceudre starts by considering the Bernstein polynomials

2 p 1 1 3 6 .

 2116 [START_REF] Alfeld | A bivariate C 2 Clough-Tocher scheme[END_REF] and[START_REF] Alfeld | Smooth macro-elements based on Powell-Sabin triangle splits[END_REF]. As the list for r = 4 shows, Type (3) elements are not excluded in principle. In addition, we get to nonvanishing univariate Bernstein polynomials on the boundary. This finishes off all elements with min{i, j, k} + ℓ = 2 and especially all elements of the Types (2) and (3) for r = 2.The application of D 2 p 1 -x gives us D Berstein polynomias which are subtracted vanish on the boundary. The other nonvanishing element is 0 Now it should be clear how the process is completed: D 3 p 1 -x extracts the nonzero polynomials with their symmetries and applying D 4 p 1 -x and D 5 p 1 -x we get

2 , 3 ,Fig. 1 . 4

 2314 Fig. 1.4 Bernstein domain points on the left, Simplex splines domain points on the right.

36 Fig. 1 . 5

 3615 Fig. 1.5 The position, left, and ordering of the domainpoints and corresponding basis functions in Figure 1.4, right.

| T 3 = O(β 3 3 )

 33 for the remaining splines. Proof The set S 6 i+1, j+1,k+1,0 with {(i, j, k) : i + j + k = 6, i ≥ 4 or j ≥ 4} are 2 × 6 = 12 classical Bernstein polynomials on the triangle T , located on the lower left and right corners in Figure 1.4. Using the explicit forms (1.50) and (1.51), (1.56) follows by inspection. For s = (i, j,k,ℓ)∈ Ī2 c i jkℓ S 6 i jkℓ let c 0 = [c 7110 , c 6210 , c 5310 , c 4410 , c 3510 , c 2610 , c 1710 ] T ∈ R 7 , c 1 = [c 6120 , c 5220 , c 4311 , c 3411 , c 2520 , c 1620 ] T ∈ R 6 , c 2 = [c 5130 , c 4221 , c 3312 , c 2421 , c 1530 ] T ∈ R 5 , (1.57) be the coefficients involved in obtaining C 2 continuity across an edge between two triangles, see Figure 1.4.

Fig. 1 . 7 Fig. 1 . 8

 1718 Fig. 1.7 Simplex splines domain points on the right, with their positions on the left. S 9 3324 , S 9 3234 , S 9 2334 have the same domain point (1/3, 1/3, 1/3) as indicated by a *

c 0 =

 0 [c 10-i,i+1,1,0 ] T i=0,...,9 , c 1 = [c 9120 , c 8220 , c 7320 , c 6411 , c 5511 , c 4611 , c 3720 , c 2820 , c 1920 ] T , c 2 = [c 8130 , c 7230 , c 6321 , c 7412 , c 4512 , c 3621 , c 2730 , c 1830 ] T , c 3 = [c 7140 , c 6231 , c 5322 , c 4413 , c 43522 , c 2631 , c 1740 ] T .

2 C 20 ( 8 , 8 : 10 ) = λ 2 1 2λ 1 λ 2 λ 2 2 Cλ 3 2λ 2 λ 3 C 21 ( 3 , 3 : 4 ) = λ 1 λ 3 2λ2 λ3-λ 2 3 2 + λ3 2 C 21 ( 4 , 4

 2208810123213343222144 [START_REF] Mann | Cubic precision Clough-Tocher interpolation[END_REF] (1,1 : 2) = 2λ 1 λ 3 2λ 2 λ 3 C 21 (2, 2 : 3) = 2λ 1

3 C

 3 33 = λ 3 3 I ∈ R 7×7

[ 0 ,

 0 j, k; ℓ] = 1 3 [1, j -1, k; ℓ] + [1, j, k -1; ℓ] + [1, j, k; ℓ -1]

β 3 + O(β 2 3 ), S 6 20 | T 3 = 30β 1 β 4 2 β 3 . 19 += 30c 19 β 4 1 β 2 + 60c 22 β 3 1 β 2 2 + 60c 23 β 2 1 β 3 2 + 30c 20 β 1 β 4 2 ,= 120 * 4c 28 β 3 1 β 2 + 270 * 4c 29 β 2 1 β 2 2 + 120 * 4c 30 β 1 β 3 2 , 3 c 34 S 6 34 +β 2 1 β 2 c 34 + β 1 β 2 2 c 35 .

 332033192222222334343435 With x = (β 1 , β 2 , 0) we then find c 23 S 6 22 + c 23 S 6 23 + c 20 S 6 20 (x) a linear combination of linearly independent univariate Bernstein polynomials of degree 5. It follows that c 19 = c 22 = c 23 = c 20 = 0. With a similar argument with S 6 24 , S 6 25 , S 6 21 on T 1 and S 6 26 , S 6 27 on T 2 we conclude that c j = 0 for j ≤ 27. Moving to the next ring we consider S 6 28 , S 6 29 , S 6 30 on the horizontal part and obtain from (1.50) and (1.51) with x = (β 1 , β 2 , a linear combination of linearly independent univariate Bernstein polynomials of degree 4. It follows that c 28 = c 29 = c 30 = 0. Moving around this ring we conclude that c j = 0, j ≤ 33. Finally, by taking third derivatives with respect to β 3 we obtain by (1.50) c 35 S 6 35 (x) = 360 * 6 This implies that c 34 = c 35 = 0 and then c 36 = 0 since by (1.50) S 6

6

  According to Definition 3, we use only 36 elements of Σ2 6 with indices from the set Ī2 , which corresponds to removing 2223 from Σ 2 6 . In Definition 4 we list these elements, normalized to ensure partition of unity. The set Σ2 6 consists of Bernstein polynomials as depicted in Figure1.3, and 15 other simplex splines with at least one knot at the barycenter of the triangle.For the numbering of the functions S6 1 , . . . , S 6 21 see Figures 1.4 and the new figure 1.5 below. In the following we give explicit formulas for the simplex splines in (1.45). On T 1 we have

	{4, 3, 1, 1} {3, 4, 1, 1} {1, 4, 3, 1} {1, 3, 4, 1} {3, 1, 4, 1}		β 6 1 -6β 2 β 5 1 + 15β 2 2 β 4 1 1 + 6β 2 β 5 -β 6 1 -15β 2 2 β 4 1 + 20β 3 2 β 3 1 S 6 24 (1) S 6 25 (1) -β 6 1
	Definition 4 (The functions Σ2 6 )
	The functions Σ2 6 consists of the 21 Bernstein polynomials
	S 6 i+1, j+1,k+1,0 := i+1, j+1,k+1;0 = B 6 i jk ,	(i, j, k) ∈ I 2 Bernstein ,		(1.44)
	and the additional simplex splines
	S 6 22 := S 6 4311 :=	1 3	4311 , S 6 23 := S 6 3411 :=	1 3	3411 , S 6 24 := S 6 1431 :=	1 3	1431 ,
	S 6 25 := S 6 1341 :=	1 3	1341 , S 6 26 := S 6 3141 :=	1 3	3141 , S 6 27 := S 6 4131 :=	1 3	4131 ,
	S 6 28 := S 6 4221 :=	2 3	4221 , S 6 30 := S 6 2421 :=	2 3	2421 , S 6 32 := S 6 2241 :=	2 3	2241 ,	(1.45)
	S 6 29 := S 6 3312 :=	1 3	3312 , S 6 31 := S 6 1332 :=	1 3	1332 , S 6 33 := S 6 3132 :=	1 3	3132 ,
	S 6 34 := S 6 3222 :=	2 3	3222 , S 6 35 := S 6 2322 :=	2 3	2322 , S 6 36 := S 6 2232 :=	2 3	2232 .

  c 7110 B 6 600 + c 6210 B 6 510 + c 5310 B 6 420 + c 4410 B 6 330 + c 3510 B 6 240 + c 2610 B 6 150 + c 1710 B 6

	060
	+c 6120 B 6 501 + c 5220 B 6 411 + c 4311 B 6 321 + c 3411 B 6 231 + c 2520 B 6 141 + c 1620 B 6 051
	+c 5130 B 6 402 + (2c 4221 -c 4311 )B 6 312 + (3c 3312 -c 4311 -c 3411 )B 6 222
	+(2c 2421 -c 3411 )B 6 132 + c 1530 B 6 042

  6 19 , S6 22 , S 6 23 , S 6 20 corresponding to S 6 i jkℓ with i + j = 7. By (1.51) S 6 j | T 3 = O(β 2 3 ) for j = 19, . . . , 36, j 19, 20, 22, 23 and

Simplex-Splines on the Clough-Tocher Split with Arbitrary Smoothness.

Acknowledgements

We thank the referees for their extremely thorough and careful review and for pointing out a flaw in the first version of the paper.

The restrictions of the simplex splines can even be written in terms of Bernstein polynomials on the three subtriangles. For

where, with ūm,n := (u m + u n )/2 for m, n = 1, 2, 3, and ū123 := (u 1 + u 2 + u 3 )/3,

and

(1.53)

Proof

The barycentric Marsden-like identity follows from the barycentric form (1.3) of the Marsden identity for Bernstein polynomials by expressing the removed Bernstein polynomials in terms of the elements in Σ2 6 . Here are some details. For (i, j, k) ∈ I 2 removed , where

we insert knots at the barycenter using (1.22), (1.45), and find (1.54)

For (i, j, k) ∈ I 2 Bernstein we have

In the second sum we insert the expressions in (1.54) for B 6 i jk , and collect terms for each i jkℓ to obtain (1.53). We show this for three typical cases.

and (1.53) follows.

Remark 3 Since ρ i jkℓ (1, 1, 1) = 1 for any i j kℓ it follows that In this section we repeat the process of the preceding one for r = 3. The methods are the same, but the expression become more lengthy. For that reason, we essentially list the results. The partition of unity basis Σ3 9 is constructed from the 75 elements i jkℓℓ defined in Example 3 as follows,

(1.66)

Theorem 3 (Barycentric Marsden-like identity for d = 9)

where

ūℓ-1 123 , otherwise .

(1.68)

Here ū123 := (u 1 + u 2 + u 3 )/3, and

Proof As for r = 2 the barycentric Marsden-like identity follows from the barycentric form (1.3) of the Marsden identity for Bernstein polynomials by expressing the 19 removed Bernstein polynomials B 9 i jk in terms of the elements in Σ3 9 . For (i, j, k) ∈ I

For (i, j, k) ∈ I 3 Bernstein we have B 9 i, j,k (β) = S 9 i+1, j+1,k+1,0 (β) and hence ρ i+1, j+1, k+1,0 (u

In the second sum we insert the expressions in (1.69)-(1.71) for B 9 i jk , and collect terms for each i jkℓ to obtain (1.67). We show this for 7 typical cases. Let ūr,s := (u r + u s )/2 for r, s = 1, 2, 3. Then

where β := (β 1 , β 2 , β 3 ) and u := (u 1 , u 2 , u 3 ). For the first 9 (Bernstein) dual functions we find ǫ 1 = ǫ 2 = ǫ 3 = 0, δ 0ℓ = 1, and (1.77) holds. For the dual function ρ[1221] corresponding to S 10 = 1221 /3 we find For r > 3 the domain points can be computed as for r ≤ 3 using the Marsden like-identity, which gives interpolation points to study the stability . The coefficients to obtain the C r -connection between two triangles T := p 1 , p 2 , p 3 and T := p 1 , p 2 , p3 can be found by a computation in the Bernstein basis of the polynomials in p 1 , p 2 , p T and p 1 , p 2 , pT .

Appendix

We provide the reader with free extra entertainment by giving alternative proofs, not facts, for some properties of our spline basis.

Generating system for the spline space

To show that Σ r 3r generates the space of all simplex splines in, S r 3r ( ), we begin with a definition and a lemma. Proof We prove the lemma by induction on the defect ν = µℓ -max{i, j, k} which starts at ν = 1 since max{i, j, k} < µℓ.

If ν = 1, suppose without loss of generality that max{i, j, k} = i, otherwise one can permute the knots appropriately. The knot insertion formula (1.22) yields that