10 December 2013

## Multivariate Bernstein operators and redundant systems

### Tan Duc Do and Shayne Waldron

Department of Mathematics, University of Auckland, Private Bag 92019, Auckland, New Zealand e–mail: [email protected] (http://www.math.auckland.ac.nz/˜waldron)

ABSTRACT

The Bernstein operator B_{n} for a simplex in IR^{d} is naturally defined via the Bernstein
basis obtained from the barycentric coordinates given by its vertices. Here we consider a
generalisation of this basis and the Bernstein operator, which is obtained from generalised
barycentric coordinates that are given for more general configurations of points in IR^{d}.
We call the associated polynomials a Bernstein frame, as they span the polynomials of
degree ≤ n, but may not be a basis. By using this redundant system we are able to give
geometrically motivated proofs of some basic properties of the corresponding generalised
Bernstein operator, such as the fact it is degree reducing and converges for all polynomials.

We also consider the conditions for polynomials in this Bernstein form to join smoothly.

Key Words: multivariate Bernstein operator, barycentric coordinates, frames, redundant system Stirling numbers, blossoming, de Casteljau algorithm

AMS (MOS) Subject Classifications: primary 41A10, 41A36 65D17, secondary 15A18, 42C15,

### 1. Introduction

The Bernstein operator [L53] and its variants [AC94] have important applications, e.g., B´ezier curves and surfaces [F02]. Here we consider a generalisation of this operator, which is based on a redundant “Bernstein basis”.

The Bernstein operator Bn for a simplex in IR^{d} is defined via the Bernstein basis for
Πn(IR^{d}) (thed–variate polynomials of degree≤n). This basis is obtained by taking powers
of the barycentric coordinates given by the vertices of the simplex. In the next section, we
outline the basic properties of generalised barycentric coordinates which are given by more
general configurations of points in IR^{d}, e.g., the vertices of a convex polygon. These lead
naturally to an analogue of the Bernstein basis, a set of polynomials of degree n which
span Πn(IR^{d}). These are not a basis if they are given by more than d+ 1 points, and so
we refer to this possibly redundant system as a Bernstein frame (cf. [C03]).

In Section 3, we define the generalised Bernstein operator given by a Bernstein frame.

We give geometrically motivated proofs of some basic properties of it. These include show- ing that it is degree reducing and converges for all polynomials, that it reproduces the linear polynomials, and more generally has the same spectral structure as the classical Bernstein operator. Similar arguments in terms of a basis would be far more cumbersome. Finally, we explore some applications of our generalised Bernstein operator. These include a de Casteljau algorithm, shape preservation properties (Section 4), and smoothness conditions in terms of the control points of the associated B´ezier surfaces (Section 5).

### 2. The Bernstein frame

Let V consist of d+ 1 affinely independent points in IR^{d}, i.e., be the vertices of a d–

simplex. Thebarycentric coordinates(cf. [B87],[LS07]) of a pointx∈IR^{d} with respect
to V are the unique coefficients (ξv(x))v∈V ∈ IR^{V} for which x can be written as an affine
combination of the points in V, i.e.,

x= X

v∈V

ξ_{v}(x)v, X

v∈V

ξ_{v}(x) = 1. (2.1)

We follow [B87] and index the barycentric coordinates by the points v ∈ V that they correspond to, and use standard multiindex notation. It follows, from (2.1), that the ξv

are linear polynomials which are a basis for Π1(IR^{d}). More generally, for any n ≥ 1, the
polynomials

B_{α} :=

|α|

α

ξ^{α}, |α|=n (α ∈ZZ^{V}_{+})
are a basis for Πn(IR^{d}). Here |α|=P

vαv, ^{n}_{α}

= ^{n!}_{α!}, and ξ^{α} =Q

vξ_{v}^{α}^{v}.

From now on, letV be a sequence (or multiset) ofm=|V|points with affine hull IR^{d},
so that each point x∈IR^{d} can be written as an affine combination

x= X

v∈V

avv, X

v∈V

av = 1, (2.2)

where the coefficients a = (av)v∈V are unique if and only if V consists of d+ 1 points.

Following [W11_{2}], we call the unique minimal ℓ_{2}–norm coefficientsa ∈IR^{V} satisfying (2.2)
the (generalised barycentric)coordinatesgiven byV, and denote them by (ξv(x))v∈V.
By construction, they satisfy (2.1). Each ξv is a linear polynomial, and they span Π_{1}(IR^{d}),
since (2.1) gives

1 = X

v∈V

ξ_{v}(x), x_{j} = X

v∈V

ξ_{v}(x)vj, j = 1, . . . , d, (2.3)

which is equivalent to the following reproduction formula for affine functions f = X

v∈V

f(v)ξv, ∀f ∈Π1(IR^{d}). (2.4)
From the formula for ξv given in [W11_{2}] it is easy to see:

• The coordinates of the barycentre c:= _{m}^{1} P

v∈V v of V are ξ_{v}(c) = _{m}^{1} , ∀v.

• ξ_{v} is constant (equal to _{m}^{1}) if and only if v is the barycentre c.

• ξ_{v} =ξ_{w} if and only if v=w.

• The ξ_{v} are continuous functions of the points v ∈V (with affine hull IR^{d}).

These imply that the set of points where the coordinates are nonnegative

NV :={x∈IR^{d} :ξv(x)≥0,∀v∈V} (2.5)
is a convex polytope, with the barycentre ofV as an interior point. We callN_{V} theregion
of nonnegativity for the coordinates given by V (See Fig. 1).

Fig 1. The region of nonnegativityN_{V} for V given by the vertices of a triangle,
square and pentagon.

We write V \w for the sequence (or multiset) obtained by removing the point w from V (once), and Aff(V) for the affine hull of the points inV. We recall:

Proposition 2.6 ([W11_{2}]). The generalised barycentric coordinates satisfy
(a) _{m}^{1} < ξ_{v}(v)≤1.

(b) ξ_{v}(w) =ξ_{w}(v).

(c) ξv(v) = 1if and only if v6∈Aff(V \v), in which case ξv = 0 on Aff(V \v).

(d) P

vξ_{v}(v) =d+ 1.

Expanding the monomial basis for Πn(IR^{d}) in terms of ξ using (2.3), shows that the
polynomials

Bα :=

|α|

α

ξ^{α}, |α|=n (2.7)

span Πn(IR^{d}). By a dimension count, these ^{n+m−1}_{m−1}

polynomials are basis if and only if V consists of d+ 1 affinely independent points. Thus, we refer to {Bα : |α| = n} as the Bernstein frame given by the points V. This is a partition of unity, since applying the multinomial theorem to (2.1) gives

X

|α|=n

Bα = X

|α|=n

n α

ξ^{α} = X

v∈V

ξv

n

= 1. (2.8)

A Bernstein frame is nonnegative onN_{V}, the region of nonnegativity given by (2.5). There
have been studies of the approximation properties of linear operators given by partitions
of unity which may take negative values on the region of interest, see, e.g., [BD85].

Example 2.9. (See Fig. 2) For V ={0,^{1}_{2},1} ⊂IR the generalised barycentric coordinates
are

ξ_{0}(x) = 5

6 −x, ξ^{1}

2(x) = 1

3, ξ_{1}(x) =x− 1
6.

Here, some polynomials in the Bernstein frame have degree < n. This is the case if and
only if the barycentre of V is a point of V. The coordinates for V ={0,^{1}_{3},^{2}_{3},1} are

ξ_{0}(x) = 7
10 − 9

10x, ξ^{1}

3(x) = 2 5 − 3

10x, ξ^{2}

3(x) = 3

10x+ 1

10, ξ_{1}(x) = 9
10x− 1

5.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.2

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

−0.2

−0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

Fig. 2. The Bernstein frame for Π_{2}(IR) for V ={0,^{1}_{2},1} and V ={0,^{1}_{3},^{2}_{3},1}.

We say a polynomial p∈Πn(IR^{d}) is in (Bernstein-B´ezier) B-formif

p= X

|α|=n

c_{α}B_{α}, (2.10)

The mesh function c : α 7→ c_{α} is unique if and only if V consists of d+ 1 points. The
mesh function with minimalℓ_{2}–norm gives a canonical B-form, i.e., what[W11_{1}] calls the
canonical coordinates of p with respect to {Bα}_{|α|=n}.

Many familiar formulas for the Bernstein basis extend to a Bernstein frame. Here are
a couple of examples (also see Section 4). Let e_{v} be the multiindex given by

e_{v}(w) :=

1, v=w;

0,
and define B_{α} := 0 if α6≥0.

Proposition 2.11. The Bernstein frame {Bα}_{|α|=n} can be calculated recursively via
Bα = X

v∈V

ξvB_{α−e}_{v}, B_{0} = 1, (2.12)

and expressed in terms of the Bernstein frame for polynomials of degree n+ 1 via Bα = X

v∈V

α_{v}+ 1

|α|+ 1B_{α+e}_{v}. (2.13)

Proof: We calculate X

v∈V

ξ_{v}B_{α−e}_{v} = X

v∈V

|α| −1
α−e_{v}

ξ^{α} = X

v∈V

α_{v}

|α|

|α|

α

ξ^{α} =B_{α},
and, using P

vξ_{v} = 1, that
Bα =Bα

X

v∈V

ξv = X

v∈V

|α|!

α! ξ^{α}ξv = X

v∈V

αv+ 1

|α|+ 1

|α+ev|!

(α+e_{v})!ξ^{α+e}^{v} = X

v∈V

αv + 1

|α|+ 1B_{α+e}_{v}.

Proposition 2.14 (Differentiation). For u, v, w ∈V, we have
D_{v−w}ξu =ξu(v)−ξu(w) =ξv(u)−ξw(u).

Thus the Bernstein frame satisfies

D_{v−w}B_{α} =|α|X

u∈V

ξ_{u}(v)−ξ_{u}(w)

B_{α−e}_{u}.
Proof: Sinceξ_{u} is affine

(D_{v−w}ξu)(x) = lim

t→0

ξu(x+t(v−w))−ξu(x) t

= lim

t→0

ξ_{u}(x) +tξ_{u}(v)−tξ_{u}(w)−ξ_{u}(x)

t =ξ_{u}(v)−ξ_{u}(w).

By the product and chain rules, we have
D_{v−w}Bα = |α|!

α! D_{v−w} Y

u∈V

ξ_{u}^{α}^{u}ξ^{α−α}^{u}^{e}^{u} = |α|!

α!

X

u∈V

αuξ_{u}^{α}^{u}^{−1} ξu(v)−ξu(w)

ξ^{α−α}^{u}^{e}^{u}

=|α|X

u∈V

ξu(v)−ξu(w)(|α| −1)!

(α−e_{u})!ξ^{α−e}^{u} =|α|X

u∈V

ξu(v)−ξu(w)

B_{α−e}_{u}.

### 3. The generalised Bernstein operator

For a Bernstein frame (2.7) given by points V in IR^{d}, we define a (generalised)
Bernstein operator Bn =Bn,V of degree n≥1 by the usual formula

B_{n}(f) := X

|α|=n

B_{α}f(vα), v_{α} := X

v∈V

α_{v}

|α|v, (3.1)

which is equivalent to

B_{n}(f) = X

v1∈V

· · · X

vn∈V

fv_{1}+· · ·+v_{n}
n

ξ_{v}_{1}· · ·ξ_{v}_{n}.

Fig. 3. The points {vα}_{|α|=n} used in the definition of Bn,Vf, where n= 7 and
V is the vertices of a triangle, square, pentagon and hexagon (respectively).

This maps functionsf which are nonnegative at the points (vα)_{|α|=k} (which are contained
in T = conv(V), the convex hull of the points in V) to polynomials of degree ≤n which
are nonnegative on the convex polytope (region of nonnegativity) N_{V} given by (2.5), so
that

f ≥0 on T =⇒ B_{n}f ≥0 on N_{V}, (3.2)

and reproduces the linear polynomials (cf. Theorem 3.19).

We now show that the generalised Bernstein operator is degree reducing, i.e.,
B_{n}(f)∈Πk(IR^{d}), ∀f ∈Πk (k = 0,1, . . .).

Define the univariate and multivariate (falling) shifted factorials by
[x]^{n}:=x(x−1)· · ·(x−n+ 1), [α]^{β} := Y

v∈V

[αv]^{β}^{v},
and the multivariate Stirling numbers of the second kind by

S(τ, β) := Y

v∈V

S(τv, βv),

where S(τv, β_{v}) are the Stirling numbers of the second kind. We note that

S(τ, β) = 0, β 6≤τ, (3.3)

and define

|α|

α

:= 0, α 6≥0. (3.4)

These are related by

α^{τ} = X

β≤τ

S(τ, β)[α]^{β}. (3.5)

Lemma 3.6. For any τ and n, we have X

|α|=n

α^{τ}
|α|

α

ξ^{α} = X

β≤τ

S(τ, β)[n]^{|β|}ξ^{β}. (3.7)

Proof: Since [|α|]^{|β| |}^{α−β|}_{α−β}

= ^{|α|}_{α}

[α]^{β} (without restriction on α and β), (3.5)
gives

X

β≤τ

S(τ, β)[|α|]^{|β|}

|α−β|

α−β

= |α|

α X

β≤τ

S(τ, β)[α]^{β} =
|α|

α

α^{τ}.
Thus, we calculate

X

|α|=n

α^{τ}
|α|

α

ξ^{α} = X

|α|=n

X

β≤τ

S(τ, β)[|α|]^{|β|}

|α−β| α−β

ξ^{α}

= X

β≤τ

S(τ, β)[n]^{|β|}ξ^{β} X

|α|=n α≥β

|α−β|

α−β

ξ^{α−β} = X

β≤τ

S(τ, β)[n]^{|β|}ξ^{β},

with the last equality given by the multinomial identity.

Theorem 3.8 (Degree reducing). The generalised Bernstein operator B_{n} is degree
reducing. More precisely,

B_{n}(ξ^{β}) = [n]^{|β|}

n^{|β|} ξ^{β} + X

0<|γ|<|β|

[n]^{|γ|}

n^{|β|} a(γ, β)ξ^{γ}, (3.9)
where w_{1}, . . . , wm is the sequence of points in V, and

a(γ, β) := X

|τ1|=βw1

· · · X

|τm|=β_{wm}

βw1

τ_{1}

ξ^{τ}^{1}(w_{1})· · ·
βwm

τ_{m}

ξ^{τ}^{m}(w_{m})S(τ_{1}+· · ·+τ_{m}, γ).

(3.10) Proof: Since each ξw is an affine function, and ξw(v) =ξv(w), we have

ξw(vα) =ξw

X

v∈V

α_{v}

|α|v

= X

v∈V

α_{v}

|α|ξw(v) = X

v∈V

α_{v}

|α|ξv(w),

and the multinomial identity gives
(ξ^{β})(v_{α}) = Y

w∈V

X

v∈V

αv

|α|ξ_{v}(w)βw

= Y

w∈V

X

|τ|=βw

βw

τ

α^{τ}

|α|^{β}^{w}ξ^{τ}(w)

= X

|τ1|=βw1

· · · X

|τm|=βwm

βw1

τ_{1}

· · · βwm

τ_{m}

ξ^{τ}^{1}(w_{1})· · ·ξ^{τ}^{m}(wm)α^{τ}^{1}^{+···+τ}^{m}

|α|^{|β|} .
Thus, by rearranging (3.1) and Lemma 3.6, we have

B_{n}(ξ^{β}) = X

|τ1|=βw1

· · · X

|τm|=β_{wm}

β_{w}_{1}
τ_{1}

ξ^{τ}^{1}(w1)· · ·
β_{w}_{m}

τm

ξ^{τ}^{m}(wm) X

|α|=n

α^{τ}^{1}^{+···+τ}^{m}
n^{|β|}

n α

ξ^{α}

= X

|τ1|=βw1

· · · X

|τm|=β_{wm}

β_{w}_{1}
τ_{1}

ξ^{τ}^{1}(w_{1})· · ·
β_{w}_{m}

τ_{m}

ξ^{τ}^{m}(w_{m})

× X

γ≤τ1+···+τm

[n]^{|γ|}

n^{|β|} S(τ_{1}+· · ·+τm, γ)ξ^{γ}.

Here B_{n}(ξ^{β}) is written as a polynomial in ξ of degree≤ |β|, so thatB_{n} is degree reducing.

The terms of degree |β| can be simplified using the multinomial identity, ξv(wj) =ξwj(v), and (2.4), as follows

X

|τ1|=βw1

· · · X

|τm|=β_{wm}

βw1

τ_{1}

ξ^{τ}^{1}(w_{1})· · ·
βwm

τ_{m}

ξ^{τ}^{m}(w_{m})[n]^{|β|}

n^{|β|} ξ^{τ}^{1}^{+···+τ}^{m}

= [n]^{|β|}

n^{|β|}

Ym

j=1

X

|τj|=β_{wj}

β_{w}_{j}
τj

ξ^{τ}^{j}(wj)ξ^{τ}^{j}

= [n]^{|β|}

n^{|β|}

Ym

j=1

X

v∈V

ξ_{v}(wj)ξv

β_{wj}

= [n]^{|β|}

n^{|β|}

Ym

j=1

X

v∈V

ξwj(v)ξv

β_{wj}

= [n]^{|β|}

n^{|β|}

Ym

j=1

ξw^{β}^{wj}j = [n]^{|β|}

n^{|β|} ξ^{β}.

By collecting the terms of degree < |β|, we obtain (3.9). Here, (3.3) allows us to remove
the restrictionγ ≤τ_{1}+· · ·+τ_{m}, and there are no terms of degree 0 since S(1,0) = 0.

Remark 3.11. If V consists of d+ 1 affinely independent points, then a(γ, β) =

S(β, γ), γ ≤β;

0, γ 6≤β,

and (3.9) simplifies to

Bn(ξ^{β}) = [n]^{|β|}

n^{|β|} ξ^{β}+ X

γ<β

[n]^{|γ|}

n^{|β|} S(β, γ)ξ^{γ}.

This was proved in [CW02](Lemma 2.1) for the case when β_{v}_{0} = 0 for some v_{0} ∈V.
Example 3.12 (Linear reproduction). For |β|= 1, we have

Bn(ξv) =ξv, ∀v∈V, (3.13)

i.e., B_{n} reproduces the linear polynomials Π1(IR^{d}) = span{ξv}v∈V. This is equivalent to

x = X

|α|=n

B_{α}(x)v_{α}, X

|α|=n

B_{α}(x) = 1, x ∈IR^{d}. (3.14)

Example 3.15 (Quadratics). For |β|= 2, we recall S(1,0) = 0,S(2,1) = 1, so that
a(eu,2ew) =ξ_{u}^{2}(v) =ξ^{2}_{v}(u), a(eu, e_{v} +e_{w}) =ξ_{u}(v)ξu(w) = (ξvξ_{w})(u), v 6=w, (3.16)
and we obtain

B_{n}(ξ^{β}) =
1− 1

n

ξ^{β} + 1
n

X

u∈V

ξ^{β}(u)ξ_{u} →ξ^{β}, as n→ ∞, |β|= 2. (3.17)

SinceB_{n} is not a positive operator in general, the application of the Korovkin theory
is more involved (see Theorem 3.31). We easily obtain the following encouraging result.

Corollary 3.18 (Convergence). For all polynomialsf, Bn(f)→f, as n→ ∞.

Proof: It suffices to consider f =ξ^{β}. For n≥ |β|, (3.9) gives
B_{n}(ξ^{β})−ξ^{β} =[n]^{|β|}

n^{|β|} −1

ξ^{β} − X

|γ|<|β|

[n]^{|γ|}

n^{|β|} a(γ, β)ξ^{γ},
where ^{[n]}_{n}|β|^{|β|} −1,^{[n]}_{n}|β|^{|γ|} =O(^{1}_{n}), as n→ ∞.

The remaining eigenstructure of B_{n} is as follows.

Theorem 3.19 (Diagonalisation). The generalised Bernstein operator B_{n} given by
points V ⊂IR^{d} is diagonalisable, with eigenvalues

λ^{(n)}_{k} := [n]^{k}

n^{k} , k= 1, . . . , n, 1 =λ^{(n)}_{1} > λ^{(n)}_{2} >· · ·> λ^{(n)}_{n} >0.

Let P_{k,V}^{(n)} denote the λ^{(n)}_{k} –eigenspace. Then

P_{1,V}^{(n)} = Π1(IR^{d}), ∀n. (3.20)

For k >1, P_{k,V}^{(n)} consists of polynomials of exact degree k, and is spanned by
p^{(n)}_{ξ}β =ξ^{β}+ X

0<|α|<|β|

c(α, β, n)ξ^{α}, |β|=k, (3.21)

where the coefficients can be calculated using (3.10) and the recurrence formula c(α, β, n) := a(α, β)

1− |β|, |α|=|β| −1,
c(α, β, n) := [n]^{|α|}

λ^{(n)}_{|β|} −λ^{(n)}_{|α|}

a(α, β)

n^{|β|} + X

|α|<|γ|<|β|

c(γ, β, n)a(α, γ)
n^{|γ|}

, |α|<|β| −1.

(3.22)
Proof: By Example 3.12, theλ^{(n)}_{1} = 1 eigenspace P_{1,V}^{(n)} contains Π_{1}(IR^{d}). Recall,
from (3.9), that B_{n}(ξ^{β}) has the form

B_{n}(ξ^{β}) =λ^{(n)}_{|β|}ξ^{β}+ X

0<|γ|<|β|

[n]^{|γ|}

n^{|β|} a(γ, β)ξ^{γ}. (3.23)
Motivated by this, we seek λ^{(n)}_{k} –eigenfunctions of the form

f =ξ^{β}+ X

0<|α|<|β|

c(α, β, n)ξ^{α}, |β|=k >1.

We observe that for such an eigenfunction the coefficients c(α, β, n) are not unique– even
when V consists of d+ 1 points. Expanding B_{n}(f) =λ^{(n)}_{k} f using (3.23) gives

Bn(f) =Bn(ξ^{β}) + X

0<|γ|<|β|

c(γ, β, n)Bn(ξ^{γ})

=λ^{(n)}_{|β|}ξ^{β}+ X

0<|α|<|β|

[n]^{|α|}

n^{|β|} a(α, β)ξ^{α} + X

0<|γ|<|β|

c(γ, β, n)

λ^{(n)}_{|γ|}ξ^{γ}+ X

0<|α|<|γ|

[n]^{|α|}

n^{|γ|} a(α, γ)ξ^{α}

=λ^{(n)}_{|β|}ξ^{β}+ X

0<|α|<|β|

λ^{(n)}_{|β|}c(α, β, n)ξ^{α}.

(3.24)
Equating the ξ^{α}, 0<|α|<|β| coefficients gives

λ^{(n)}_{|β|}c(α, β, n) = [n]^{|α|}

n^{|β|} a(α, β) +c(α, β, n)λ^{(n)}_{|α|} + X

|α|<|γ|<|β|

c(γ, β, n)[n]^{|α|}

n^{|γ|} a(α, γ). (3.25)
Since λ^{(n)}_{|α|} > λ^{(n)}_{|β|}, this is equivalent to

c(α, β, n) = 1
λ^{(n)}_{|β|} −λ^{(n)}_{|α|}

[n]^{|α|}

n^{|β|} a(α, β) + X

|α|<|γ|<|β|

c(γ, β, n)[n]^{|α|}

n^{|γ|} a(α, γ)
.
From this we can define suitable c(α, β, n) recursively, as in (3.22), starting from

c(α, β, n) := 1
λ^{(n)}_{|β|} −λ^{(n)}_{|α|}

[n]^{|α|}

n^{|β|} a(α, β) = a(α, β)

1− |β|, |α|=|β| −1.

A simple dimension count shows that the eigenfunction {p_{ξ}^{β}}_{|β|=k}, so defined, span a
space P_{k,V}^{(n)} of dimension ^{k+d−1}_{d−1}

. Again, by dimension counting, we conclude that B_{n} is
diagonalisable, withP_{1,V}^{(n)} = Π_{1}(IR^{d}).

Example 3.26 (Quadratic eigenfunctions). Using (3.16), we have
p^{(n)}_{ξ}β =ξ^{β} −X

u∈V

ξ^{β}(u)ξu, |β|= 2. (3.27)

In general, p^{(n)}_{ξ}α does depend on n (cf. [CW00]).

Despite the fact the coefficients in (3.21) are not unique, we can take their limit as n→ ∞. This indicates that the redundant expansion (3.21) is natural.

Corollary 3.28 (Limits of the eigenfunctions). For 0<|α|<|β|,

n→∞lim c(α, β, n) =c^{∗}(α, β),
where

c^{∗}(α, β) := a(α, β)

1− |β|, |α|=|β| −1,

c^{∗}(α, β) := 2

(|β| − |α|)(−|α| − |β|+ 1) X

|γ|=|α|+1

c^{∗}(γ, β)a(α, γ), |α|<|β| −1. (3.29)

Thus, the eigenfunctions of (3.21) satisfy
p^{(n)}_{ξ}β →p^{∗}_{ξ}β :=ξ^{β} + X

0<|α|<|β|

c^{∗}(α, β)ξ^{α}, as n→ ∞. (3.30)

Proof: Fix β. We use strong induction on j = |β| − |α| = 1, . . . ,|β| to prove the limit exists. For |α| = |β| −1 (j = 1) the limit is clear. Suppose the limit of c(γ, β, n) exists for all γ with |α|<|γ|<|β|. Then taking the limit of (3.22) gives

n→∞lim c(α, β, n) = 2

(|β| − |α|)(|β| − |α| −2|β|+ 1) X

|γ|=|α|+1

c^{∗}(γ, β)a(α, γ).

This follows from the calculations
λ^{(n)}_{|β|} −λ^{(n)}_{|α|} = [n]^{|α|}[n− |α|]^{|β|−|α|}

n^{|β|} − [n]^{|α|}

n^{|α|} = [n]^{|α|}

n^{|β|} ([n− |α|]^{|β|−|α|}−n^{|β|−|α|}),
[n− |α|]^{|β|−|α|}−n^{|β|−|α|} = 1

2(|β| − |α|)(−|α| − |β|+ 1)n^{|β|−|α|−1}+ lower order powers of n.

Since p^{(n)}_{ξ}β , p^{∗}_{ξ}β ⊂Π_{|β|}(IR^{d}), we have the convergence asserted in (3.30).

We now prove the strongest Korovkin theorem that the restricted positivity property (3.2) allows. Since this requires a modification of the usual argument, which is not stated in the literature, we give a self contained proof. This result supercedes Corollary 3.18.

Theorem 3.31 (Korovkin). Let T =T_{V} be the convex hull of V, and N_{V} be the region
of nonnegativity. For f ∈C(T), Bnf →f uniformly on NV.

Proof: Let ε > 0, and M be the maximum of f over T. Since f is uniformly continuous on the compact set T, there is aδ >0 such that |f(s)−f(t)|< ε, ∀ks−tk< δ.

Thus, we obtain the estimate

|f(s)−f(t)| ≤ε+ 2Mks−tk^{2}

δ^{2} , ∀s, t∈T,
For fixed t, we have

−ε− 2M

δ^{2} k · −tk^{2} ≤f−f(t)≤ε+ 2M

δ^{2} k · −tk^{2} onT ,
and so applying B_{n} (which reproduces constants) and using (3.2) gives

−ε− 2M

δ^{2} B_{n} k · −tk^{2}

≤B_{n}f −f(t)≤ε+ 2M

δ^{2} B_{n} k · −tk^{2}

on N_{V}.
This last step is the main difference in argument. For t ∈N_{V}, evaluating at t gives

|Bnf(t)−f(t)| ≤ε+ 2M

δ^{2} B_{n} k · −tk^{2}

(t), ∀t∈N_{V}. (3.32)
We now estimate B_{n}(k · −tk^{2})(t). From (2.1), we obtain

k · −tk^{2} =hX

v

{ξ_{v}v−ξ_{v}t},X

w

{ξ_{w}w−ξ_{w}t}i=X

v

X

w

ξ_{v}ξ_{w}hv−t, w−ti.

Thus, (3.17) gives
Bn k · −tk^{2}

=X

v

X

w

n 1− 1

n

ξvξw+ 1 n

X

u∈V

ξv(u)ξw(u)ξu

o

hv−t, w−ti

= 1− 1

n

k · −tk^{2}+ 1
n

X

v

X

w

X

u∈V

ξv(u)ξw(u)ξuhv−t, w−ti

= 1− 1

n

k · −tk^{2}+ 1
n

X

u∈V

ξuhu−t, u−ti, so that

B_{n} k · −tk^{2}

(t) = 1 n

X

u∈V

ku−tk^{2}ξ_{u}(t) = 1
n

X

u∈V

kuk^{2}ξ_{u}(t)− ktk^{2}

≤ 1
nD^{2},
where D:= diam(T) is the diameter of T. Thus, we obtain the uniform estimate

|Bnf(t)−f(t)| ≤ε+ 1 n

2M D^{2}

δ^{2} , t ∈NV,
and conclude that Bnf →f uniformly on NV.

### 4. Applications to CAGD

We now consider how our generalised Bernstein operator B_{n} = B_{n,V} might be used
in CAGD (computer aided geometric design) to describe polynomials defined on convex
polyhedra which are not simplices. We first observe that the Korovkin theory does not
allow a Bernstein type operator which is positive on the entire region T = conv(V) if the
points of V are not the vertices of a simplex or a cube.

Theorem 4.1 ([MRS94]). LetT be a convex polygon with five or more vertices. There is
no positive linear operator L:C(T)→C(T)which reproduces Π_{1}(IR^{2}) and maps Π_{2}(IR^{2})
to itself other than the identity.

We recall, that our B_{n} reproduces Π1(IR^{2}) and maps Π2(IR^{2}), but has the restricted
positivity property (3.2), i.e.,

f ≥0 on T =⇒ Bnf ≥0 on NV,

where NV is the region of nonnegativity (2.5). In the multivariate case, the description
of other shape preserving properties is involved (cf. [S91]). We mention some which do
generalise easily. These suggest that the target region for representing polynomials on
should perhaps be the convex polyhedron N_{V}, rather than T = conv(V).

Proposition 4.2 (Shape preservation). Let T =TV be the convex hull of V, andNV

be the region of nonnegativity. If f is convex onT, then

B_{n}f ≥B_{n+1}f ≥ · · · ≥f on N_{V}. (4.3)
Proof: Fixx∈N_{V}. By (2.1), and a calculation, we have the convex combinations

x= X

|α|=n

Bα(x)vα, vβ = X

v∈V

β_{v}

|β|v_{β−e}_{v}.
Hence for f convex, Jensen’s inequality gives

Bnf(x) = X

|α|=n

Bα(x)f(vα)≥f(x), X

v∈V

βv

|β|f(v_{β−e}_{v})≥f(vβ). (4.4)
By the degree raising formula (2.13), we have

B_{n}f = X

|α|=n

B_{α}f(v_{α}) = X

v∈V

X

|α|=n

αv + 1

|α|+ 1B_{α+e}_{v}f(v_{α}),
so that

B_{n}f−B_{n+1}f = X

|β|=n+1

c_{β}B_{β}, c_{β} := X

v∈V

β_{v}

|β|f(vβ−ev)−f(vβ).

By (4.4), cβ ≥0, so that Bnf ≥B_{n+1}f on NV.

A polynomial in B-form can be calculated via the de Casteljau algorithm. We present
this in terms of the blossom (polar form) of a polynomial p∈Πn(IR^{d}), which we recall
(cf. [R89], [DLG91], [DMS92]) is the unique symmetric n–affine function ^{ω}p with

ωp(t, . . . , t) =p(t), ∀t ∈IR^{d}. (4.5)
Proposition 4.6 (de Casteljau algorithm). Suppose that p∈Πn(IR^{d}) has the B-form

p= X

|α|=n

c_{α}B_{α}. (4.7)

Then the blossom of p at t_{1}, . . . , t_{n} can be calculated from (cα)_{|α|=n} via
for j = 1 to n do

cα := X

v∈V

ξv(tj)c_{α+e}_{v}, |α|=n−j
end for

with c_{0} =^{ω}p(t_{1}, . . . , tn). In particular, taking t_{1} =· · ·=tn =t gives c_{0} =p(t).

Proof: Suppose p∈Πn(IR^{d}) is in the B-form (4.7), i.e., equivalently
p=X

v1

· · ·X

vn

c_{e}_{v}

1+···+e_{vn}ξ_{v}_{1}· · ·ξ_{v}_{n}. (4.8)
Then its blossom is given by

P(t_{1}, . . . , tn) :=X

v1

· · ·X

vn

cev1+···+e_{vn}ξv1(t_{1})· · ·ξvn(tn), (4.9)
since P clearly defines an n–affine function, which, by (4.8), satisfies (4.5). The n steps of
the algorithm are the averagings given by the n sums in (4.9).

In the terminology of [DMS92], this is a symmetric simplicial algorithm (here the points V need not be the vertices of a simplex).

The coefficients c_{α} in the B-form (4.7) of p are not unique, unlessV is the vertices of
a simplex. Using (2.1) to expand, we have

p(t) =^{ω}p(t, . . . , t) =^{ω}p(X

v1

ξ_{v}_{1}(t)v_{1}, . . . ,X

vn

ξ_{v}_{n}(t)v_{n})

=X

v1

· · ·X

vn

ξ_{v}_{1}(t)· · ·ξ_{v}_{n}(t)^{ω}p(v1, . . . , v_{n}),
and so the coefficients c_{α} in (4.7) can be chosen to be the blossoms

b_{α} =^{ω}p(α·V), α·V := (. . . , v, . . . , v

| {z }

αv times

, . . .), (4.10)

which we call the blossoming coefficients.

Proposition. For any p of the form (4.7), the blossoming coefficients b = (bα)_{|α|=n} are
given by the matrix multiplication

b=Qc, Q= [qαβ]|α|=n,|β|=n, q_{αβ} :=B^{ω}_{β}(α·V).

For n >1, our calculations show that the blossoming coefficients are not the ℓ_{2}–norm
minimising choice (which is given by multiplication by an orthogonal projection matrix).

Remark 4.11. The previous discussion extends to p = (p_{1}, . . . , ps) : IR^{d} → IR^{s}, where
p_{j} ∈Πn(IR^{d}) and c_{α} ∈IR^{s}, by considering coordinates. Since vector valued p are used in
practice, e.g., B´ezier curves in IR^{3}, we henceforth state our results in this setting.

We define the control points of the curve (d = 1), surface (d = 2), etc, given by
t7→p(t)∈IR^{s}, where p has B-form (4.7), to be

nv_{α}
c_{α}

:|α|=no

⊂IR^{d+s}.

Since the B-form is not unique when V has more thand+ 1 points, a given curve (surface, etc) may be given by different sets of control points. This redundancy has advantages from the point of view of design, as a given surface could be arrived at by a number of different choices of control points. Equivalently, each control point has less influence on the the surface, and so moving a control point causes a small change of shape making the fine tuning of a surface easier. Once a suitable surface has been obtained it can be presented in terms of a set of canonical control points – if so desired.

By combining (3.14) and (4.7), we have x

p(x)

= X

|α|=n

cαBα(x), cα :=

vα

c_{α}

. (4.12)

Thus a point x, p(x)

on the curve (surface, etc) can be calculated by the de Casteljau
algorithm applied to the vectors (cα)_{|α|=n} ⊂IR^{d+s}. Many basic properties of B´ezier curves
and surfaces follow from (4.12). We now list these using the terminology of [F02].

Convex hull property: Since P

|α|=nB_{α}(x) = 1 and B_{α}(x) ≥ 0, x ∈ N_{V}, the point
x, p(x)

is an affine combination (also called a barycentric combination in CAGD) of the
control points (cα)_{|α|=n}, which is a convex combination if x∈N_{V}.

Affine invariance: For A : IR^{s} →IR^{ℓ} an affine map, we have
A p(x)

= X

|α|=n

(Ac_{α})B_{α}(x).

Invariance under affine parameter transformations: ForA: IR^{d} →IR^{d} an invertible
affine map, B_{α}(Ax) =A(B_{α}(x)), so that

p(Ax) = X

|α|=n

c_{α}A B_{α}(x)
.

Symmetry: Suppose that A is an affine map which maps (the multiset) V to V. Then
ξv(x) =ξAv(Ax) (see [W11_{2}]), so that

B_{α}(x) =B_{α◦A}^{−1}(Ax),

where α◦A^{−1} :V →ZZ_{+} denotes the multiindex v7→α_{A}^{−1}_{v}. Using this we obtain
p(x) = X

|α|=n

c_{α}B_{α}(x) = X

|α|=n

c_{α}B_{α◦A}^{−1}(Ax) = X

|β|=n

c_{β◦A}B_{β}(Ax).

Invariance under barycentric combinations: The affine combination of two curves (surfaces, etc) is given by corresponding affine combination of the control points, i.e.,

λ X

|α|=n

b_{α}B_{α}(x) + (1−λ) X

|α|=n

c_{α}B_{α}(x) = X

|α|=n

λb_{α} + (1−λ)c_{α} B_{α}(x), λ∈IR.

Endpoint interpolation: If v∈V satisfies v6∈Aff(V \v), then Prop. 2.6 gives (Bnf)(v) =f(v).

The remaining property given in[F02]ispseudolocal control, i.e., the fact thatB_{α}
is peaked at the point vα, and so moving the control point cα has the most influence on
the curve (surface, etc) near the point v_{α}. Numerical calculations indicate some degree of
localisation of the Bα near vα, but we do not make any quantitative statements here.

To summarise, the main features of our Bernstein polynomial approximations on a nonsimplicial convex polytope are:

• A surface may be determined by several different choices of control points.

• Moving the control points has less influence on the shape if there are many of them.

• On N_{V} the surface is a convex combination of the control points.

• Convex functions are monotonely approximated from above by B_{n} on N_{V}.

Typical applications, that require a single polynomial defined on a convex polytope, include the description and construction of finite elements on regular polygons and the design of hexagonal lenses (cf. [D87]).

### 5. Smoothness and multivariate splines

We now consider the possible application of our results to multivariate spline theory, where surfaces are constructed by joining polynomial pieces together smoothly. There is a highly developed theory (and associated software) that involves polynomials defined on simplices, based on the description of the smoothness conditions in terms of the control

points (cf. [B87], [LS07]). We assume that the reader is familiar with this, and we inves-
tigate how it extends to partitions involving nonsimplical cells. There was work in this
direction (see [L89], [CL90]) on mixed grid partitionswhere in the bivariate case the cells
were trianglesandparallelograms, and in the trivariate case they weretetrahedrons,prisms
and parallelopipeds. The B-form developed for nonsimplicial cells was for polynomials of
coordinate degree n (as were the spline spaces), rather than total degree n, though the
positioning of the control points is the same as we propose. Bivariate C^{1}–quadratics on a
polygonal partition were considered in [Wh90].

We will give our results in terms of the blossoming coefficients. The following is adapted from the[B87],[R89]and[L91](cf.[W13]). If (cα) are the blossoming coefficients, then we refer to the (cα) of (4.12) as the blossoming control points of p = P

cαBα. Let

r·x:=x, . . . , x

| {z }

r times

,

Theorem 5.1.. LetW be a set of points inIR^{d}, with L= Aff(W). Then the polynomials
f, g ∈Πn(IR^{d}) have all derivatives of order ≤r onL equal if and only if

ω

f(w1, . . . , w_{n−r}, r·x) =^{ω}g(w1, . . . , w_{n−r}, r·x), w_{1}, . . . , w_{n−r} ∈W, x∈IR^{d}. (5.2)
Let

α·V := (. . . , v, . . . , v

| {z }

αv times

, . . .), α∈ZZ^{+}_{V}.

Lemma 5.3. Let cbe the blossoming coefficients of (4.10) for p∈Πn(IR^{d}). Then

ωp(α·V, r·x) = X

|β|=r

c_{α+β}Bβ(x), |α|=n−r, r = 0, . . . , n.

Proof: Since the blossom^{ω}p is n–affine, we have

ωp(α·V, r·x) =^{ω}p(α·V, X

v1∈V

ξ_{v}_{1}(x)v1, . . . , X

vr∈V

ξ_{v}_{r}(x)vr)

= X

v1∈V

· · · X

vr∈V

ξ_{v}_{1}(x)· · ·ξ_{v}_{r}(x)^{ω}p(α·V, v_{1}, . . . , v_{r})

= X

|β|=r

c_{α+β}B_{β}(x).

Combining Theorem 5.1 and Lemma 5.3 gives smoothness conditions for the C^{r}–
joining of polynomials in terms of their blossoming B–form coefficients.

We now illustrate this for C^{1}–quadratics on a region consisting of a triangle and a
quadrilateral with a common edge.

Example 5.4 (C^{1}–quadratics on triangle and quadrilateral). Without loss of gen-
erality, suppose the vertices of the triangle and quadrilateral are

V ={v1, . . . , v_{3}}={0, e1, a}, a_{2} 6= 0, V˜ ={˜v_{1}, . . . ,v˜_{4}}={0, e1, b,−e2},
so the common edge has endpoints v_{1} = ˜v_{1} = 0. and v_{2} = ˜v_{2} =e_{1} = (1,0). We index the
generalised barycentric coordinatesξ and ˜ξ by both their order and the points themselves,
e.g., ξ_{3} and ξ_{a} =ξ_{(a}_{1}_{,a}_{2}_{)}, whichever is the most convenient. We have

ξ_{1}(x, y) = (1−x) + a_{1}−1

a_{2} y, ξ_{2}(x, y) =x− a_{1}

a_{2}y, ξ_{3}(x, y) = 1
a_{2}y,
ξ˜_{(0,0)}(x, y) = −(b^{2}_{2}+b_{1}b_{2}+b_{1}+ 1)x+ (b^{2}_{1}+b_{1}b_{2}−b_{1}+ 1)y+b^{2}_{1}+b^{2}_{2}+ 1

2(b^{2}_{1}+b^{2}_{2}−b_{1}b_{2}+b_{2}−b_{1}+ 1)

ξ˜_{(1,0)}(x, y) = (2b^{2}_{2}−b_{1}b_{2}+ 2b_{2}−b_{1}+ 2)x+ (b^{2}_{1}−2b_{1}b_{2}−b_{1})y+b^{2}_{1}−b_{1}b_{2}−b_{1}
2(b^{2}_{2}+b_{2}+ 1−b_{1}b_{2}−b_{1}+b^{2}_{1}) ,
ξ˜_{(b}_{1}_{,b}_{2}_{)}(x, y) = (2b1−b_{2}−1)x+ (2b2−b_{1}+ 1)y+b_{2}−b_{1}+ 1

2(b^{2}_{2}+b_{2}+ 1−b_{1}b_{2}−b_{1}+b^{2}_{1}) ,

ξ˜_{(0,−1)}(x, y) = (−b^{2}_{2}+ 2b_{1}b_{2}−b_{2})x+ (−2b^{2}_{1}+b_{1}b_{2}+ 2b_{1}−b_{2}−2)y+b^{2}_{2}−b_{1}b_{2}+b_{2}
2(b^{2}_{2}+b_{2}+ 1−b_{1}b_{2}−b_{1}+b^{2}_{1}) ,
Suppose that f, g ∈Π_{2}(IR^{2}) are quadratics with blossoming coeffient B–forms

f = X

|ga|=2 α∈ZZV

+

c_{α}B_{α}, g= X

|ga|=2 α∈ZZV˜ +

˜
c_{α}B˜_{α}.

In Theorem 5.1, we take n= 2, and

W ={v_{1}, v_{2}}={˜v_{1},v˜_{2}}={0, e_{1}}.

For r = 0, the condition forf and g to have a continuous join on the line L= Aff(W) is
given by the two element sequences from W ={v_{1}, v_{2}}, i.e.,

ω

f(v_{1}, v_{1}) =^{ω}g(v_{1}, v_{1}), f^{ω}(v_{1}, v_{2}) =^{ω}g(v_{1}, v_{2}), ^{ω}f(v_{2}, v_{2}) =^{ω}g(v_{2}, v_{2}).

which by Lemma 5.3 is

c_{(2,0,0)} = ˜c_{(2,0,0,0)}, c_{(1,1,0)} = ˜c_{(1,1,0,0)}, c_{(0,2,0)} = ˜c_{(0,2,0,0)}.

For r = 1, the condition for f and g to have a C^{1}–join on L is given by the one element
sequences from W ={v_{1}, v_{2}}, i.e.,

ωf(v_{1}, x) =^{ω}g(v_{1}, x), f^{ω}(v_{2}, x) =^{ω}g(v_{2}, x).

By Lemma 5.3, these equalities of linear polynomials in x can be written

c_{(2,0,0)}ξ_{1}+c_{(1,1,0)}ξ_{2}+c_{(1,0,1)}ξ_{3} = ˜c_{(2,0,0,0)}ξ˜_{1}+ ˜c_{(1,1,0,0)}ξ˜_{2}+ ˜c_{(1,0,1,0)}ξ˜_{3}+ ˜c_{(1,0,0,1)}ξ˜_{4},
c_{(1,1,0)}ξ_{1}+c_{(0,2,0)}ξ_{2}+c_{(0,1,1)}ξ_{3} = ˜c_{(1,1,0,0)}ξ˜_{1}+ ˜c_{(0,2,0,0)}ξ˜_{2}+ ˜c_{(0,1,1,0)}ξ˜_{3}+ ˜c_{(0,1,0,1)}ξ˜_{4}.
Since the blossoming control points for the quadrilateral are simply those for corresponding
triangles, it follows that these smoothness conditions have the usual geometric interpreta-
tion: that all the control points involved (3 and 4 respectively, see Fig. 4) lie in a common
plane.

Fig. 4. The C^{1}–smoothness conditions of Example 5.4.

### 6. Future work

For applications where it is desirable to have generalised barycentric coordinates which
are nonnegative on the convex hull of V, one could modify the definition of (ξv(x))_{v∈V} for
x∈conv(V) to be the unique minimal ℓ_{2}–norm coefficients a∈IR^{V} satisfying

x= X

v∈V

avv, X

v∈V

av = 1, av ≥0. (6.1)

With this definition ξv is a continuous piecewise linear polynomial, which coincides with
with the original on N_{V}.

### References

[AC94] F. Altomare and M. Campiti, “Korovkin-type approximation theory and its applica- tions”, Walter de Gruyter, Berlin, 1994.

[B87] C. de Boor, B-form basics, in “Geometric Modeling: Algorithms and New Trends”

(G. E. Farin Ed.), pp. 131–148, SIAM Publications, Philadelphia, 1987.

[BD85] C. de Boor and R. DeVore, Partitions of unity and approximation,Proc. Amer. Math.

Soc. 93 (1985), 705–709.

[C03] O. Christensen, “An introduction to frames and Riesz bases”, Birkh¨auser, Boston, 2003.

[CL90] C. K. Chui and M.-J. Lai, Multivariate vertex splines and finite elements,J. Approx.

Theory 60 (1990), 245–343.

[CW00] S. Cooper and S. Waldron, The eigenstructure of the Bernstein operator, J. Approx.

Theory 105 (2000), 133-165.

[CW02] S. Cooper and S. Waldron, The diagonalisation of the multivariate Bernstein operator, J. Approx. Theory 117 (2002), 103–131.

[DMS92] W. Dahmen, C. A. Micchelli, and H.-P. Seidel, Blossoming begetsB-spline bases built better by B-patches, Math. Comp. 59(1992), 97–115.

[DLG91] T. DeRose, M. Lounsbery, and R. Goldman, A tutorial introduction to blossoming, Comput. Graph. Systems Appl. (1991), 267–286.

[D11] Tan Duc Do, Generalised Bernstein operators on polytopes, Honours dissertation, 2011.

[D87] C. F. Dunkl, Orthogonal polynomials on the hexagon, SIAM J. Appl. Math. 47 (1987), 343–351.

[F02] G. Farin, “Curves and surfaces for computer-aided geometric design: A practical guide (5th ed.)”, Academic Press, San Diego, 2002.

[L89] M.-J Lai, On construction of bivariate and trivariate vertex splines on arbitrary mixed grid partitions, Dissertation (Texas A & M Univ.), 1989.

[L91] Ming Jun Lai, A characterization theorem of multivariate splines in blossoming form, Comput. Aided Geom. Design 8 (1991), 513-521.

[LS07] M.-J. Lai and L. L. Schumaker, “Spline functions on triangulations”, Cambridge Uni- versity Press, Cambridge, 2007.

[L53] G. G. Lorentz, “Bernstein Polynomials”, Toronto Press, Toronto, 1953.

[MRS94] F.-J. Mu˜noz-Delgado, V. Ram´irez Gonz´alez, and T. Sauer, Domains for Bernstein polynomials, Appl. Math. Lett. 7 (1994), 7–9.

[S91] T. Sauer, Multivariate Bernstein polynomials and convexity, Comput. Aided Geom.

Design 8 (1991), 465–478.

[R89] L. Ramshaw, Blossoms are polar forms, Comput. Aided Geom. Design 6 (1989), 323- 358.

[W11_{1}] S Waldron, Frames for vector spaces and affine spaces, Linear Algebra Appl. 435
(2011), 77–94.

[W11_{2}] S Waldron, Affine generalised barycentric coordinates, Jaen J. Approx. 3 (2011),
209–226.

[W13] S Waldron, Blossoming and smoothly joining polynomials on affine subspaces, Preprint, 2013.

[Wh90] W. Whiteley, The geometry of bivariate C_{2}^{1}-splines, Preprint, 1990.