- Email: [email protected]
Local lagrange interpolation with cubic C1 splines on tetrahedral partitions
Journal of Approximation Theory 160 (2009) 89 – 102 www.elsevier.com/locate/jat
Local lagrange interpolation with cubic C 1 splines on tetrahedral partitions Gero Hecklina , Günther Nürnbergera , Larry L. Schumakerb,∗ , Frank Zeilfeldera a Institute for Mathematics, University of Mannheim, 68 131 Mannheim, Germany b Department of Mathematics, Vanderbilt University, Nashville, TN 37240, USA
Received 24 August 2007; accepted 31 January 2008 Communicated by Carl de Boor Available online 17 October 2009 Dedicated to Professor Paul L. Butzer on the occasion of his 80th birthday
Abstract We describe an algorithm for constructing a Lagrange interpolation pair based on C 1 cubic splines defined on tetrahedral partitions. In particular, given a set of points V ∈ R3 , we construct a set P containing V and a spline space S13 () based on a tetrahedral partition whose set of vertices include V such that interpolation at the points of P is well-defined and unique. Earlier results are extended in two ways: (1) here we allow arbitrary sets V, and (2) the method provides optimal approximation order of smooth functions. © 2008 Elsevier Inc. All rights reserved. MSC: 41A15; 41A05; 65D05; 65D07; 65D17; 41A63 Keywords: Trivariate splines; Lagrange interpolation
1. Introduction m be a set of points in R3 . In this paper we are interested in the following problem. Let V:={i }i=1
Problem 1.1. Find a tetrahedral partition whose set of vertices includes V, an N-dimensional N space S of C 1 cubic splines defined on , and a set of additional points {i }i=m+1 such that for ∗ Corresponding author.
E-mail addresses: [email protected] (G. Hecklin), [email protected] (G. Nürnberger), [email protected] (L.L. Schumaker), [email protected] (F. Zeilfelder). 0021-9045/$ - see front matter © 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.jat.2008.01.012
90
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
N , there is a unique spline s ∈ S1 () satisfying every choice of real numbers {ri }i=1 3
s(i ) = ri , i = 1, . . . , N . We call
N P:={i }i=1
and
S13 ()
(1.1) a Lagrange interpolation pair.
It is easy to solve this problem using C 0 splines, see Remark 1. However, the situation is much more complicated if we want to use C 1 splines. In [1] we solved the problem in the special case where the points V lie on a rectangular grid and the partition is a Freudenthal partition, see Remark 3. In [7] we presented a C 1 cubic spline method that works for arbitrary tetrahedral partitions, but with suboptimal approximation power. The aim of this paper is to describe a method which works for arbitrary tetrahedral partitions while giving optimal approximation order four. To achieve this aim, we start with an arbitrary initial tetrahedral partition 0 with vertex set V, and then define an appropriate refinement and an associated set of points P such that P and the spline space S13 () form a Lagrange interpolation pair. The construction of both and P is based on special orderings of the tetrahedra in 0 arising from two different decompositions of 0 , but using different and more sophisticated techniques than in [7]. In addition, here we make use of the partial Worsey–Farin splits introduced in [1], rather than the simpler Worsey–Farin splits used in [7]. The paper is organized as follows. In the next section, we recall some basic notation and a few key concepts from the Bernstein–Bézier theory of splines, and present two useful results on bivariate splines. We discuss partial Worsey–Farin splits in Section 3. In Section 4 we give two algorithms for decomposing tetrahedral partitions. The construction of our Lagrange-interpolation pair follows from the algorithms of Sections 5 and 6, which show how to construct and P, respectively. The main result of the paper can be found in Section 7. In Section 8, we show that the resulting interpolation method is local and stable, and use these facts to establish error bounds for how well the interpolating spline approximates a given smooth function. We conclude with some remarks. 2. Preliminaries We recall that for any given tetrahedral partition , the associated space of C 1 cubic splines is S13 ():={s ∈ C 1 () : s|T ∈ P3 , for all T ∈ }, where P3 is the 20-dimensional space of trivariate cubic polynomials. Throughout the paper we use standard Bernstein–Bézier techniques for dealing with trivariate splines. Here we recall only some of the most critical concepts and notation. For a detailed treatment and more on the notation, see the book [2] or our earlier papers [1,7]. Given a tetrahedron T = [v1 , v2 , v3 , v4 ] with vertices v1 , v2 , v3 , and v4 , we write iv1 + jv2 + kv3 + v4 T DT := i jkl := 3 i+ j+k+=3 for the associated set of domain points, and write D for the union of the sets DT over all T ∈ . The ball (of radius 1) around v1 is the set D T (v1 ):={iTjk : i 2}. Associated with an edge
e:=v1 , v2 of T, the tube (of radius 1) around e is the set t T (e):={iTjk : k, 1}, see [2, p. 439], where balls and tubes of arbitrary radius are defined. If is a tetrahedral partition, then the balls and tubes are defined as D(v) = ∪{D T (v) : T has a vertex at v}, and t(e) = {t T (e) : T
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
91
contains the edge e}. Associated with a tetrahedron T in a partition , we set star 0 (T ):=T , and for 1, define star (T ) to be the union of the set of all tetrahedra in which touch a tetrahedron in star −1 (T ). We recall that a spline s is uniquely defined by its set {c }∈D of B-coefficients, and that a spline s belongs to C 1 () if and only these coefficients satisfy an appropriate set of simple linear side conditions, see [2, p. 454]. We also make use of the concepts of minimal determining sets and nodal minimal determining sets. Recall that a set M of domain points of a spline space S is called a minimal determining set for S provided it is the smallest set of such points such that the corresponding coefficients {c }∈M can be set independently, and all other coefficients of s can be consistently determined from smoothness conditions, i.e., in such a way that all smoothness n conditions are satisfied, see [2, p. 485]. Suppose N:={i }i=1 is a set of linear functionals of 1 2 3 the form i :=i || m i ai D , where D :=Dx D y Dz , and i denotes point-evaluation at the point i in . Then N is called a nodal determining set for S if s ∈ S and s = 0 for all ∈ N implies s ≡ 0. If there is no smaller such set, then N is called a nodal minimal determining set, see [2, p. 490]. For later use, we now present two results on bivariate interpolation. Let F:=v1 , v2 , v3 be an arbitrary triangle. The first result is a minor extension of Lemma 5.1 in [1]. Lemma 2.1. Suppose that we are given all of the coefficients {ciFjk }i+ j+k=3 of a bivariate cubic F . Then for any given real number r, there exists a unique c F so that polynomial p except for c111 111 F F is a stable process in the sense that p(111 ) = r . Moreover, the computation of c111 F | C |r | + max |ciFjk | , |c111 (i, j,k) (1,1,1)
where C = 29 . Proof. The interpolation condition gives F F F F c111 B111 (111 )=r− ciFjk BiFjk (111 ), (i, j,k) (1,1,1)
where BiFjk are the bivariate Bernstein-basis polynomials associated with the triangle F. The result F ( F ) = 2/9 and the fact that the B F form a partition of then follows from the fact that B111 111 i jk unity.
Now suppose FC T is the well-known Clough–Tocher split of F into three subtriangles Fi :=v F , vi , vi+1 , i = 1, 2, 3, where v F is a point in the interior of F. The following result is similar to Lemma 5.2 in [1]. Lemma 2.2. Suppose that we are given all of the B-coefficients of a C 1 cubic bivariate spline s F1 F1 F1 F1 , c210 , c201 , c111 . Then for any given real number r, there exists a defined on FC T except for c300 F1 unique choice of these coefficients so that s(111 ) = r . The computation of these coefficients is stable in the sense that there is an absolute constant C such that
|ciFjk1 | C |r | + max |c | , (i jk) ∈ {(300), (210), (201), (111)}, c ∈K
where K is the set of known coefficients of s.
92
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
F1 Proof. The interpolation condition s(111 ) = r leads to the equation F1 F1 F1 F1 + 3c210 + 3c201 + 6c111 = R1 , c300
where R1 is a combination of r and the known B-coefficients of s. Suppose b1 , b2 , b3 are the barycentric coordinates of v F relative to F. Writing down the C 1 smoothness conditions across the interior edges of F1 , we are led to the linear system ⎛ F ⎞ 1 ⎞ c300 ⎞ ⎛ ⎛ 1 3 3 6 R1 ⎜ F ⎟ 1 ⎟ ⎜ 1 −b1 −b2 0 ⎟ ⎜ c210 R2 ⎟ ⎟ ⎜ ⎟⎜ ⎟, ⎜ (2.1) ⎜ ⎟=⎜ ⎝0 1 0 −b2 ⎠ ⎜ c F1 ⎟ ⎝ R3 ⎠ 201 ⎠ ⎝ 0 0 1 −b1 R4 F1 c111 where R2 , R3 , R4 are combinations of r and the coefficients in K. The determinant of this matrix is D:= − 6 − 3b1 − 3b2 − 2b1 b2 . It is easy to see that as we vary v F in T, (b1 , b2 ) runs over the set {(b1 , b2 ) : 0 b1 , b2 1, b1 + b2 1}. This gives |D| 6 for any choice of v F , which implies both the existence of a solution and the stability of its computation. 3. Partial Worsey–Farin splits The following definition is taken from [1]. Definition 3.1. Let T be a tetrahedron, and let vT be a point in its interior. Given an integer 0 m 4, let F1 , . . . , Fm be distinct faces of T, and for each i = 1, . . . , m, let v Fi be a point in the interior of Fi . Then we define the corresponding m-th order partial Worsey–Farin split m WF of T to be the tetrahedral partition obtained by the following steps: (1) connect vT to each of the four vertices of T, (2) connect vT to the points v Fi for i = 1, . . . , m, (3) connect v Fi to the three vertices of Fi for i = 1, . . . , m. The m-th order partial Worsey–Farin split of a tetrahedron results in 4 + 2m subtetrahedra, see Fig. 1. The splits 0W F and 4W F are the well-known Alfeld and Worsey–Farin splits, respectively, m see [2, Section 16.7]. We now recall an important fact about the space S13 (m W F ), where W F is an m-th order partial Worsey–Farin split of a tetrahedron T :=v1 , v2 , v3 , v4 . Theorem 3.2. ([1, Theorem 6.3]). Fix 0 m 4. Let Mm be the union of the following sets of domain points in DmW F : (1) for each i = 1, . . . , 4, D(vi )∩Ti for some tetrahedron Ti ∈ m W F containing vi ; F ; (2) for each face F of T that is not split, the point 111 Fi 3 (3) for each face F of T that has been subjected to a Clough–Tocher split, the points {111 }i=1 , where F1 , F2 , F3 are the subfaces of F. Then Mm is a minimal determining set for S13 (m W F ).
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
93
Fig. 1. Partial Worsey–Farin splits subdividing a given tetrahedron into 4, 6, 8, 10, or 12 subtetrahedra.
4. Decomposition of tetrahedral partitions In this section we describe two algorithms for decomposing a given tetrahedral partition 0 into classes of tetrahedra. These decompositions will be used to define our Lagrange interpolation pair. Algorithm 4.1. For i = 0 to 4, repeat until no longer possible: choose a tetrahedron T such that exactly i vertices of T belong to tetrahedra that have been chosen earlier. We call these vertices the marked vertices of T. Put T into the class Ai . This algorithm decomposes the partition 0 into five classes A0 , . . . , A4 . It also produces an ordering T1 , . . . , Tn of the tetrahedra of 0 . Clearly, if T is a tetrahedron in the class A0 , then all of its vertices are unmarked. Lemma 4.2. Suppose v is a marked vertex of a tetrahedron Tn of class A j with 1 j 4. Then there exists a tetrahedron Tm of class Ai with i < j such that v is a vertex of Tm . Proof. Let v be a marked vertex of Tn , and let Tm be the first of the tetrahedra T1 , . . . , Tn−1 to contain v. Then Tm must be in some class Ai with 0 i j. We claim that it cannot be in class A j , since at the time Tm was chosen, Tn would have had at most j − 1 vertices in common with the tetrahedra T1 , . . . , Tm−1 , and thus would have been chosen before Tm unless Tm is in Ai for some i < j. The following algorithm produces a different decomposition based on shared edges. Algorithm 4.3. For i = 0 to 6, repeat until no longer possible: choose a tetrahedron T such that exactly i edges of T belong to tetrahedra that have been chosen earlier. We call these edges marked edges of T. Put T into the class Bi .
94
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
This algorithm partitions 0 into seven classes B0 , . . . , B6 . It also produces an ordering 1 , . . . , T n of the tetrahedra of 0 . It is easy to see that A0 ⊆ B0 , and if T is a tetrahedron T in the class B0 , then all of its edges are unmarked. n of class B j with 1 j 6. Then Lemma 4.4. Suppose e is a marked edge of a tetrahedron T m . there exists a tetrahedron Tm of class Bi with i < j such that v is an edge of T n , and let T m be the first of the tetrahedra T 1 , . . . , T n−1 to Proof. Let e be a marked edge of T contain e. Then Tm must be in some class Bi with 0 i j. We claim that it cannot be in class m was chosen, T n would have had at most j − 1 edges in common with the B j , since at the time T m . tetrahedra T1 , . . . , Tm−1 , and thus would have been chosen before T 5. Construction of As a first step towards creating a Lagrange interpolating pair solving Problem 1.1, in this section we describe an algorithm for constructing a suitable tetrahedral partition . As a starting point, let 0 be an arbitrary tetrahedral partition with vertices at the points V. We shall construct by splitting certain of the tetrahedra of 0 . The following algorithm proceeds in two steps. In the first step we apply Clough–Tocher splits to some of the triangular faces of 0 . The choice of which faces to split is controlled by the n of the tetrahedra in 0 produced by Algorithm 4.3. The Clough–Tocher split 1 , . . . , T ordering T of a face involves inserting a point v F in the interior of F and connecting it to each of the three vertices of F. In the second step we apply partial Worsey–Farin splits to the tetrahedra of 0 that have one or more split faces. Algorithm 5.1 (Construct ). i that is not shared with any tetrahedron T j with j < i, (1) For i = 1, . . . , n, if F is a face of T apply a Clough–Tocher split to F when either two or three of the edges of F are marked i . edges of T (2) For each tetrahedron T ∈ 0 , let m be the number of its faces that have been split in step 1. If m > 0, apply an m-th order partial Worsey–Farin split to T using its incenter. In Algorithm 5.1 we did not specify how to choose the split points v F to be used in creating the Clough–Tocher splits of faces. For our purposes, these points must be chosen in a special way: (1) if F is a boundary face of 0 , choose v F to be the barycenter of F, (2) if F is an interior face of 0 , choose v F to be the intersection of F with the line connecting the incenters of the two tetrahedra sharing F. It is easy to see that none of the tetrahedra in the classes B0 or B1 is split. Tetrahedra in class B2 are either not split, or are subjected to a partial-Worsey–Farin split of order m = 1. The types of splits that may be applied to tetrahedra in the various classes are shown in Table 1, where m-WF stands for the m-th order Worsey–Farin split. In this table the symbol “–” indicates that the corresponding split does not occur. The symbol “◦” identifies cases where a tetrahedron in a class B j may be given the split in the indicated column because the tetrahedron shares one or more faces with tetrahedra in lower classes. The symbol “×” identifies cases where a tetrahedron in B j has no neighbor in a lower class.
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
95
Table 1 Possible splits for the different classes . Class
No split
1-WF
2-WF
3-WF
4-WF
B0 B1 B2 B3 B4 B5 B6
× × × ◦ – – –
– – × ×o – – –
– – – × ◦ ◦ ◦
– – – × ×o ◦ ◦
– – – – × ×◦ ×◦
Fig. 2. Points inserted on faces (black dots) by Algorithm 6.1.
6. Construction of P In this section we present an algorithm for constructing a set P of interpolation points which together with the spline space S13 () based on the tetrahedral partition of the previous section will form a Lagrange interpolation pair. The algorithm makes use of both of the ordered lists 1 , . . . , T n created by Algorithms 4.1 and 4.3. T1 , . . . , Tn and T Algorithm 6.1 (Construction of P). For i = 1, . . . , n, (1) for each edge u, v of Ti ; (a) choose the points v and (u + 2v)/3 if v is not a marked vertex of Ti ; (b) choose the points u and (2u + v)/3 if u is not a marked vertex of Ti . i that is not a face of a tetrahedron T j with j < i; (2) for each face F of T i , (a) choose the barycenter of F if F has no marked edges of T (b) choose the barycenter of the subtriangle of F with no marked edge if F has exactly two i . marked edges of T We emphasize that the interpolation points chosen by this algorithm belong to D , and all lie on edges and faces of 0 . In particular, P contains all of the points in V, i.e., the vertices of 0 . Interpolation points are chosen on faces of tetrahedra in 0 when the face is not split and contains no marked edges, or when the face is split but contains exactly two marked edges, see Fig. 2, where the thicker lines indicate marked edges, i.e. edges in common with tetrahedra appearing n . 1 , . . . , T earlier in the list T 7. The main result In this section we prove the main result of this paper, Theorem 7.3 below, which states that P and S13 () form a Lagrange interpolating pair. We begin by showing that P is a minimal determining set for S13 ().
96
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
Theorem 7.1. The set M:=P produced by Algorithm 6.1 is a minimal determining set for S13 (). Proof. Suppose that for each ∈ M, we fix the corresponding coefficient of a spline s ∈ S03 (). We now show that all other coefficients of s are consistently determined, i.e., in such a way that all C 1 smoothness conditions are satisfied so that s lies in S13 (). Step 1: It is easy to see by the construction of P that for each vertex v of 0 , there is a tetrahedron T in 0 with vertex at v such that P contains the four domain points in D T (v). Then the coefficients associated with the remaining domain points in D(v) are consistently determined by the C 1 smoothness at v. The consistency is insured by the fact that these balls do not overlap. Step 2: To show that all remaining coefficients are consistently determined, we consider the n of Algorithm 4.3, and apply Theorem 3.2. We begin with T 1 . 1 , . . . , T tetrahedra in the order T We already have determined the coefficients associated with the balls around the vertices of T1 . F 1 is in the class B0 , none of its faces is split and P contains the point 111 for each face Since T F, which means that the corresponding coefficient has been fixed. Theorem 3.2 then implies that 1 are consistently determined. The all remaining coefficients corresponding to domain points in T 1 C conditions imply that all coefficients associated with domain points in the tubes around the 1 are also consistently determined. Now suppose we have determined the coefficients edges of T 1 , . . . , T j−1 , and consider T j . Using corresponding to all of the domain points in the tetrahedra T 1 C smoothness, we get all coefficients of s corresponding to domain points in the balls around the j . To apply Theorem 3.2, we j , as well as in the tube around each marked edge of T vertices of T now show that we have already determined the coefficients corresponding to the domain points j . If F is shared with a tetrahedron T i in items (2) and (3) of the theorem. Consider a face F of T with i < j, then all coefficients corresponding to domain points on F are already known. Suppose now that F is not shared with such a tetrahedron. There are four cases: F . We have already (a) F has no marked edges. Then F is not split, and P contains the point 111 fixed the corresponding coefficient. (b) F has one marked edge e. Then F is not split, and we already know the coefficient correF since it lies in the tube around e. sponding to 111 (c) F has two marked edges. Then F has been subjected to a CT split. In this case the coefficients corresponding to the barycenters of the three subtriangles are known since the corresponding domain points are either in a tube around a marked edge, or the point is in P. (d) F has three marked edges. Then F has been subjected to a CT split. In this case the coefficients corresponding to the barycenters of the three subtriangles are known since the corresponding domain points are in tubes around the marked edges.
We have shown that the coefficients of s corresponding to domain points in the minimal determining set of Theorem 3.2 are determined, and we conclude from the theorem that the coefficients j are consistently determined. Continuing of s corresponding to all remaining domain points in T with the remaining tetrahedra, we get all coefficients of s. There is one subtle point concerning consistency that remains to be discussed. It concerns the C 1 smoothness conditions across faces of 0 . Suppose F is a face of 0 that is shared by two tetrahedra Tm and Tn , and suppose F has been subjected to a Clough–Tocher split. Then a necessary and sufficient condition that all C 1 smoothness conditions across F are satisfied is that the split point v F lie on the straight line between the interior split points vTm and vTn , see the proof of Theorem 18.11 in [2]. But we have made sure that this is the case in our construction of .
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
97
The fact that P is a minimal determining set for S13 () implies that the dimension of S13 () is equal to the number N of points in the set P, see [2, Section 17.3]. Clearly, N = 4n V + n F , where n V = #V is the number of vertices of 0 , and n F is the number of faces where an interpolation point was added to P in step 2 of Algorithm 6.1. To state our next theorem, we recall that denotes point-evaluation at the point . Theorem 7.2. Let P be the set of domain points chosen by Algorithm 6.1. Then the set N:={ }∈P is a nodal minimal determining set for S13 (). Proof. We have already observed that the dimension of S13 () is equal to #P = #N. Thus, to prove N is a nodal minimal determining set for S13 (), it suffices to show that if we fix the values of s() for all ∈ P, then all B-coefficients of s are determined. Step 1: We show that all coefficients of s associated with domain points on the edges of 0 and in the balls around the vertices of 0 are determined. To this end, we examine the tetrahedra of 0 in the order T1 , . . . , Tn defined by Algorithm 4.1. Suppose we have established the assertion for the edges and vertices of all of the tetrahedra T1 , . . . , T j−1 , and consider T j . We have already determined the coefficients of s in the balls around marked vertices of T j . Let e:=u, v be an edge of T j . Suppose first that u, v are both unmarked vertices of T j . Then P contains all four domain points on e, call them 1 , . . . , 4 . We are given the values of the univariate cubic polynomial s|e at these four points. This leads to a linear system of four equations for the four B-coefficients of s associated with 1 , . . . , 4 . It is easy to see that this system is nonsingular with determinant equal to 12 81 , independent of the length or orientation of the edge e. Now suppose that one vertex of e is marked, say u. Then we know the coefficients associated with D(u), but not those associated with D(v). But now we can use the interpolation conditions at the two points (u + 2v)/3 and v to get a linear system of 2 equations for the B-coefficients associated with these two domain points. This system is nonsingular with determinant equal to 12 81 , independent of the length or orientation of the edge e. We have shown how to compute the coefficients of s associated with all domain points on the edges of T j . We then use these to determine the coefficients associated with domain points in the ball D(v) around each unmarked vertex of T j . Step 2: We now show that all remaining coefficients of s can be computed from the interpolation 1 , . . . , T n defined by Algorithm conditions. To this end, we consider the tetrahedra in the order T 4.3. Suppose we have computed all coefficients of s associated with domain points in the tetrahedra 1 , . . . , T j−1 , and consider T j . It suffices to show how to compute the coefficients associated with T the minimal determining set Mm of Theorem 3.2. We already have the coefficients associated j and in the tubes around its marked edges. with domain points in the balls around the vertices of T We now deal with the coefficients corresponding to the remaining points in Mm . These lie on j . Let F be a face of T j . If F is shared with a tetrahedron T i with i < j, then we have faces of T already computed the coefficients of s associated with all domain points on F. Thus, it suffices to j that are not shared with a tetrahedron T i with i < j. consider faces F of T Suppose that F is not split. We need to show how to compute the coefficient associated with F F j , then P contains the point 111 111 . There are two subcases. If F contains no marked edges of T , and we can use the interpolation condition at this point to compute the corresponding coefficient, F j , then 111 lies in the tube around that edge, see Lemma 2.1. If F contains one marked edge of T and the corresponding coefficient is already known. Now suppose F has been subjected to a Clough–Tocher split. We have to show how to compute Fi the coefficients associated with the domain points 111 , where F1 , F2 , F3 are the subfaces of F.
98
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
j , then the coefficient corresponding to Fi is already known If Fi contains a marked edge of T 111 since this point lies in the tube around the marked edge. Thus, if all three edges of F are marked j , we are done. Now suppose only two of the edges of F are marked edges of T j , say edges of T F1 those contained in F2 , F3 . Then P contains the barycenter 111 of F1 , and we can apply Lemma F1 . We now apply Theorem 3.2 to determine the 2.2 to compute the corresponding coefficient c111 j . coefficients of s corresponding to all remaining domain points in T We are ready to establish the main result of the paper. N together with the spline space S1 () form a Lagrange interTheorem 7.3. The set P:={i }i=1 3 polation pair.
Proof. The proof of Theorem 7.2 shows that for any set of real numbers r1 , . . . , r N , there is a unique spline s ∈ S13 () satisfying s(i ) = ri for i = 1, . . . , N . 8. Bounds on the error of the interpolant N and S1 () are the Lagrange interpolation pair constructed above Suppose that P:={i }i=1 3 from an initial tetrahedral partition 0 of a domain ⊂ R3 . Then for every f ∈ C(), there is a unique spline s f ∈ S13 () such that
s f (i ) = f (i ), i = 1, . . . , N . This interpolation process defines a linear projector mapping C() onto S13 (). Our goal in this section is to provide a bound on f − s for smooth functions, where the error is measured in the maximum norm on . To accomplish this goal, we will apply Theorem 17.22 of [2]. To do so, we need to show that nodal minimal determining set P of Theorem 7.2 is local and stable, see Definition 17.21 of [2]. For each ∈ D , let :={ ∈ P : c depends on s()}, where {c }∈D are the B-coefficients of s f . Theorem 8.1. The nodal minimal determining set P for S13 () is local in the sense that for all ∈ D , ⊆ star (T ),
(8.1)
with = 10, where T is a tetrahedron in 0 that contains . It is also stable in the sense that there exists a constant C depending only on the smallest solid and face angles in 0 such that for all ∈ D , |c | C max | f ()|. ∈
(8.2)
Proof. For the definition of solid and face angles of a tetrahedron, see [2, p. 462]. We consider three cases. Case 1: ( lies on an edge of a tetrahedron T in 0 ). Let e:=u, v be the edge of T containing . If neither u nor v is a marked vertex of T , then P contains all four domain points on e, and the
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
99
A0 A1
A2
A3
A4
Fig. 3. Longest chain of influence of an interpolation value along the edges of 0 . Interpolation points are shown as gray and black dots, and the direction of propagation (of the value f at the black dot) into neighboring tetrahedra of higher classes Ai is indicated by the arrows. The value f at the black dot can influence the B-coefficient indicated as an open circle, which is associated with a domain point in a tetrahedron from class A4 .
corresponding coefficients depend only on the values of s at these four points, see the proof of Theorem 7.2. For these coefficients, (8.1) holds with = 0. This argument applies to all edges of tetrahedra T of class A0 . If T lies in a higher class, the situation is more complicated, and there is a certain propagation of influence. We claim that if T is of class A j , then (8.1) holds with = j. We have already established this for j = 0, and can now proceed by induction. Suppose we have established the claim for classes A0 , . . . , A j−1 , and consider T ∈ A j . If neither end of e is marked, then as before we get (8.1) with = 0. Now suppose u is marked. Let ∈ D(u). By Lemma 4.2, there exists a tetrahedron Tn in class A j−1 sharing the vertex u, and we can compute c as a combination of coefficients associated with D Tn (u). By the induction hypothesis, these can be computed from data at points lying in star j−1 (Tn ), and we conclude that (8.1) holds with = j. It remains to consider the case ∈ D(v), where u is marked but v is not. In this case c is computed by solving a 2 × 2 linear system to find the coefficients associated with the domain points (u + 2v)/2 and v. But the right-hand side of this system includes coefficients associated with domain points in the ball D(u). It follows that ⊆ star j (T ). We conclude that in all cases, if lies on any edge of a tetrahedron T , then (8.1) holds with = 4. The worst case is illustrated in Fig. 3. We now verify the stability condition (8.2) when lies on an edge of T . By the proof of Theorem 7.2, coefficients with domain points on edges of 0 are only computed in one of three ways: (1) by solving a 2 × 2 linear system whose determinant is 49 ; (2) by solving a 4 × 4 linear system whose determinant is 12 81 ; (3) by employing the C 1 smoothness at a vertex to compute a coefficient from previously computed coefficients. The constant of stability in (3) depends on the smallest solid and face angles of 0 . Case 2: ( lies on a face F of a tetrahedron T in 0 , but not on an edge). We claim that if T is of class B j , then (8.1) holds with = 4+ j. This is immediate for j = 0, since when T is of class B0 , then c is computed from Lemma 2.1, which shows that c depends on the value of s(), but also on the coefficients corresponding to domain points in the balls around the vertices of F. We showed in Case 1 that such coefficients depend only on data at points in star4 (T ). Suppose now that we have established the claim for all tetrahedra in the classes B0 , . . . , B j−1 , and consider T ∈ B j . Now c will be computed from a C 1 smoothness condition around a marked edge e, or by applying Lemma 2.1 or 2.2. Thus, c may depend on the value of s at a point in F, but more importantly can also depend on coefficients associated with points in another tetrahedron sharing the marked edge e, which by Lemma 4.4 is in a lower class. Then, applying the induction
100
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
hypothesis, we conclude that ⊆ star4+ j (T ). The worst case is for j = 6, and we conclude that (8.1) holds with = 10 for all . We now verify the stability condition (8.2) when lies on a face F of a tetrahedron T in 0 , but not on an edge. By the proof of Theorem 7.2, coefficients with domain points in the interior of faces of 0 are only computed in one of three ways: (1) by enforcing a C 1 smoothness condition around an edge, (2) by Lemma 2.1, (3) by Lemma 2.2. All three of these computations are stable, where the constant of stability in (1) depends on the smallest solid and face angles of 0 . Case 3: ( lies in the interior of a tetrahedron T in 0 ). In this case c is computed from known coefficients on the faces of T , cf. Theorem 3.2, by enforcing C 1 smoothness conditions. Thus, (8.1) holds with = 10. These computations are stable with a constant depending on the largest solid and face angles in 0 . Given a compact set B ⊆ and an integer m > 0, let | f |m,B :=
D f B
(8.3)
||=m
be the usual Sobolev seminorm, where . B denotes the infinity norm on B. For any tetrahedron T in 0 , let |T | be the length of its longest edge. The following result follows immediately from Theorem 17.22 of [2]. Theorem 8.2. Given T ∈ 0 , let T :=star10 (T ) in 0 . Then for every f ∈ C m+1 (T ) with 0 m 3,
D ( f − s f ) T K |T |m+1−|| | f |m+1,T ,
(8.4)
for all || m. If T is convex, the constant K depends only on the smallest solid and face angles of T as a subpartition of . If T is not convex, K also depends on the Lipschitz constant of the boundary of T . The global version of this theorem also holds with T and T replaced by , and with |T | replaced by the mesh size |0 | of 0 , i.e., the length of the longest edge in 0 . 9. Remarks Remark 1. Suppose V is an arbitrary set of points in R3 , and that is a tetrahedral partition with vertices V. Set P:=V, and let S:=S01 () be the space of continuous linear splines. Then clearly P and S form a Lagrange interpolation pair. It is also straightforward to create a Lagrange interpolation pair using C 0 splines of higher degree, provided we add an appropriate set of additional interpolation points. Remark 2. There is a fairly extensive history of Lagrange interpolation with bivariate splines, although even in this case, it is quite difficult to construct Lagrange interpolating pairs using C 1 splines. See [3–6,8,9].
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
101
Remark 3. Much less is known about Lagrange interpolation with C 1 trivariate splines. The first local Lagrange method was constructed in [10] using quintic splines, but only for the case where the points of V lie on a grid. With this same restriction on V, we recently constructed a Lagrange interpolation pair using cubic splines, see [1]. In [7] we constructed Lagrange interpolation pairs for an arbitrary initial point set V using both quadratic and cubic splines. However, as pointed out in the introduction, that cubic spline method did not provide optimal approximation order. Remark 4. For given V and initial tetrahedral partition 0 , there are generally many different associated Lagrange interpolation pairs. They arise by different choices of the decompositions in Section 4. Remark 5. Our construction shows that the minimal determining set M of Theorem 7.1 is stable and 1-local in the sense of Definition 17.11 of [2]. But then Theorem 17.14 of [2] applies to show that S13 () approximates functions in the Sobolev spaces Wqm+1 () with 0 m 3 and 1q ∞ to order |0 |m+1 . In particular, with m = 3, this shows that the space has full approximation power in all of the q-norms. Remark 6. With a little additional work, it can be shown that (8.1) in Theorem 8.1 holds with = 9. This is based on the fact that if T is in class B1 , then none of its faces is split, and so coefficients associated with a tube around one edge do not interact via C 1 smoothness conditions with coefficients associated with a tube around its opposite edge. This implies that if T is in class B2 , then (8.1) actually holds with = 5, and induction then gives a worst case of = 9. We note that for most tetrahedral partitions, will be smaller than 9. Remark 7. It should be emphasized that care is needed in the choice of the interior split points vT in Step 2 of Algorithm 5.1 for constructing the refined partition from 0 . To insure C 1 smoothness across interior faces of 0 , we need to choose the split point on the face to lie on the line joining the interior split points of the neighboring tetrahedra. To accomplish this, we have chosen the incenters of these tetrahedra, rather than the more natural barycenters, since then the line joining them always intersects the common face F at an interior point, see Lemma 16.24 in [2]. Remark 8. It was shown in Corollary 6.2 of [1] that a C 1 cubic spline defined on a partialWorsey–Farin split of a tetrahedron is actually C 2 at the interior split point vT . It follows that the splines in our space S13 () enjoy this property at all such vertices of . Remark 9. Given a set of n vertices V, an associated initial triangulation 0 can be computed with standard algorithms. Once we have 0 , the operation count for constructing our Lagrange interpolation pair is O(n). Finding the coefficients of an interpolating spline is also of linear complexity in n. Remark 10. In step 2 of Algorithm 5.1 we apply a partial-Worsey–Farin split to T whenever the number of CT-split faces m is positive. We could have also applied a 0-th order split (also called the Alfeld split) in the case when m = 0, but this would lead to a refinement with more tetrahedra, so instead we do not split the tetrahedron at all when m = 0. Remark 11. Trivariate local Lagrange interpolation methods are useful for the construction and reconstruction of volumetric models, as well as for scattered data fitting problems. A major
102
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
advantage is that they only require function values, but no derivatives. For example, the method here could be used in a two-stage process, where in the first stage we construct a C 0 linear spline that interpolates at the vertices of a very fine tetrahedral partition. In the second stage, we construct a C 1 cubic spline on a coarser tetrahedral subpartition , where the data needed for interpolation are taken directly from the linear spline. Because of the optimal approximation order of the interpolation methods, it is natural to choose the partitions such that |1 |2 is about |0 |4 . For numerical results based on this idea in the bivariate setting, see [9]. Remark 12. As shown in Theorem 17.17 of [2], there is a natural stable local basis { }∈M for S13 () associated with the minimal determining set M of Theorem 7.1. Since M is 1-local, each has support at most star(T ), where T is a tetrahedron of 0 containing . Remark 13. There is a different stable local basis { }∈M for S13 () associated with the nodal minimal determining set N of Theorem 7.2. Since N is just the set of point-evaluation functionals at the points of P, the are Lagrange basis splines in the sense that () = , for all , ∈ P. Theorem 8.1 implies that each has support at most star10 (T ), where T is a tetrahedron of 0 containing . For typical tetrahedral partitions, most if not all of these basis functions will have much smaller supports. Remark 14. In the theory of trivariate splines, nodal determining sets are usually constructed using point-evaluation of both function values and derivatives. Here we have restricted ourselves to the use of function values only in order to get a Lagrange interpolation method rather than a Hermite interpolation method. References [1] G. Hecklin, G. Nürnberger, L.L. Schumaker, F. Zeilfelder, A local Lagrange interpolation method based on C 1 cubic splines on Freudenthal partitions, Math. Comp. 77 (2008) 1017–1036. [2] M.-J. Lai, L.L. Schumaker, Splines on Triangulations, Cambridge University Press, Cambridge, 2007. [3] G. Nürnberger, V. Rayevskaya, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation with C 2 splines of degree seven on triangulations, in: M. Neamtu, E. Saff (Eds.), Advances in Constructive Approximation, Nashboro Press, Brentwood, TN, 2004, pp. 345–370. [4] G. Nürnberger, V. Rayevskaya, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation with bivariate splines of arbitrary smoothness, Constr. Approx. 23 (2006) 33–59. [5] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation by bivariate C 1 cubic splines, in: T. Lyche, L.L. Schumaker (Eds.), Mathematical Methods for Curves and Surfaces III, Oslo, 2000, Vanderbilt University Press, Nashville, 2000, pp. 393–404. [6] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Lagrange interpolation by C 1 cubic splines on triangulated quadrangulations, Adv. Comp. Math. 21 (2004) 357–380. [7] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Two Lagrange interpolation methods based on C 1 splines on tetrahedral partitions, in: C.K. Chui, M. Neamtu, L.L. Schumaker (Eds.), Approximation Theory XI: Gatlinburg 2004, Nashboro Press, Brentwood, 2005, pp. 327–344. [8] G. Nürnberger, F. Zeilfelder, Developments in bivariate spline interpolation, J. Comp. Appl. Math. 121 (2000) 125–152. [9] G. Nürnberger, F. Zeilfelder, Lagrange interpolation by bivariate C 1 splines with optimal approximation order, Adv. Comp. Math. 21 (2004) 381–419. [10] L.L. Schumaker, T. Sorokina, C 1 quintic splines on type-4 tetrahedral partitions, Adv. Comp. Math. 21 (2004) 421–444.
Local lagrange interpolation with cubic C 1 splines on tetrahedral partitions Gero Hecklina , Günther Nürnbergera , Larry L. Schumakerb,∗ , Frank Zeilfeldera a Institute for Mathematics, University of Mannheim, 68 131 Mannheim, Germany b Department of Mathematics, Vanderbilt University, Nashville, TN 37240, USA
Received 24 August 2007; accepted 31 January 2008 Communicated by Carl de Boor Available online 17 October 2009 Dedicated to Professor Paul L. Butzer on the occasion of his 80th birthday
Abstract We describe an algorithm for constructing a Lagrange interpolation pair based on C 1 cubic splines defined on tetrahedral partitions. In particular, given a set of points V ∈ R3 , we construct a set P containing V and a spline space S13 () based on a tetrahedral partition whose set of vertices include V such that interpolation at the points of P is well-defined and unique. Earlier results are extended in two ways: (1) here we allow arbitrary sets V, and (2) the method provides optimal approximation order of smooth functions. © 2008 Elsevier Inc. All rights reserved. MSC: 41A15; 41A05; 65D05; 65D07; 65D17; 41A63 Keywords: Trivariate splines; Lagrange interpolation
1. Introduction m be a set of points in R3 . In this paper we are interested in the following problem. Let V:={i }i=1
Problem 1.1. Find a tetrahedral partition whose set of vertices includes V, an N-dimensional N space S of C 1 cubic splines defined on , and a set of additional points {i }i=m+1 such that for ∗ Corresponding author.
E-mail addresses: [email protected] (G. Hecklin), [email protected] (G. Nürnberger), [email protected] (L.L. Schumaker), [email protected] (F. Zeilfelder). 0021-9045/$ - see front matter © 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.jat.2008.01.012
90
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
N , there is a unique spline s ∈ S1 () satisfying every choice of real numbers {ri }i=1 3
s(i ) = ri , i = 1, . . . , N . We call
N P:={i }i=1
and
S13 ()
(1.1) a Lagrange interpolation pair.
It is easy to solve this problem using C 0 splines, see Remark 1. However, the situation is much more complicated if we want to use C 1 splines. In [1] we solved the problem in the special case where the points V lie on a rectangular grid and the partition is a Freudenthal partition, see Remark 3. In [7] we presented a C 1 cubic spline method that works for arbitrary tetrahedral partitions, but with suboptimal approximation power. The aim of this paper is to describe a method which works for arbitrary tetrahedral partitions while giving optimal approximation order four. To achieve this aim, we start with an arbitrary initial tetrahedral partition 0 with vertex set V, and then define an appropriate refinement and an associated set of points P such that P and the spline space S13 () form a Lagrange interpolation pair. The construction of both and P is based on special orderings of the tetrahedra in 0 arising from two different decompositions of 0 , but using different and more sophisticated techniques than in [7]. In addition, here we make use of the partial Worsey–Farin splits introduced in [1], rather than the simpler Worsey–Farin splits used in [7]. The paper is organized as follows. In the next section, we recall some basic notation and a few key concepts from the Bernstein–Bézier theory of splines, and present two useful results on bivariate splines. We discuss partial Worsey–Farin splits in Section 3. In Section 4 we give two algorithms for decomposing tetrahedral partitions. The construction of our Lagrange-interpolation pair follows from the algorithms of Sections 5 and 6, which show how to construct and P, respectively. The main result of the paper can be found in Section 7. In Section 8, we show that the resulting interpolation method is local and stable, and use these facts to establish error bounds for how well the interpolating spline approximates a given smooth function. We conclude with some remarks. 2. Preliminaries We recall that for any given tetrahedral partition , the associated space of C 1 cubic splines is S13 ():={s ∈ C 1 () : s|T ∈ P3 , for all T ∈ }, where P3 is the 20-dimensional space of trivariate cubic polynomials. Throughout the paper we use standard Bernstein–Bézier techniques for dealing with trivariate splines. Here we recall only some of the most critical concepts and notation. For a detailed treatment and more on the notation, see the book [2] or our earlier papers [1,7]. Given a tetrahedron T = [v1 , v2 , v3 , v4 ] with vertices v1 , v2 , v3 , and v4 , we write iv1 + jv2 + kv3 + v4 T DT := i jkl := 3 i+ j+k+=3 for the associated set of domain points, and write D for the union of the sets DT over all T ∈ . The ball (of radius 1) around v1 is the set D T (v1 ):={iTjk : i 2}. Associated with an edge
e:=v1 , v2 of T, the tube (of radius 1) around e is the set t T (e):={iTjk : k, 1}, see [2, p. 439], where balls and tubes of arbitrary radius are defined. If is a tetrahedral partition, then the balls and tubes are defined as D(v) = ∪{D T (v) : T has a vertex at v}, and t(e) = {t T (e) : T
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
91
contains the edge e}. Associated with a tetrahedron T in a partition , we set star 0 (T ):=T , and for 1, define star (T ) to be the union of the set of all tetrahedra in which touch a tetrahedron in star −1 (T ). We recall that a spline s is uniquely defined by its set {c }∈D of B-coefficients, and that a spline s belongs to C 1 () if and only these coefficients satisfy an appropriate set of simple linear side conditions, see [2, p. 454]. We also make use of the concepts of minimal determining sets and nodal minimal determining sets. Recall that a set M of domain points of a spline space S is called a minimal determining set for S provided it is the smallest set of such points such that the corresponding coefficients {c }∈M can be set independently, and all other coefficients of s can be consistently determined from smoothness conditions, i.e., in such a way that all smoothness n conditions are satisfied, see [2, p. 485]. Suppose N:={i }i=1 is a set of linear functionals of 1 2 3 the form i :=i || m i ai D , where D :=Dx D y Dz , and i denotes point-evaluation at the point i in . Then N is called a nodal determining set for S if s ∈ S and s = 0 for all ∈ N implies s ≡ 0. If there is no smaller such set, then N is called a nodal minimal determining set, see [2, p. 490]. For later use, we now present two results on bivariate interpolation. Let F:=v1 , v2 , v3 be an arbitrary triangle. The first result is a minor extension of Lemma 5.1 in [1]. Lemma 2.1. Suppose that we are given all of the coefficients {ciFjk }i+ j+k=3 of a bivariate cubic F . Then for any given real number r, there exists a unique c F so that polynomial p except for c111 111 F F is a stable process in the sense that p(111 ) = r . Moreover, the computation of c111 F | C |r | + max |ciFjk | , |c111 (i, j,k) (1,1,1)
where C = 29 . Proof. The interpolation condition gives F F F F c111 B111 (111 )=r− ciFjk BiFjk (111 ), (i, j,k) (1,1,1)
where BiFjk are the bivariate Bernstein-basis polynomials associated with the triangle F. The result F ( F ) = 2/9 and the fact that the B F form a partition of then follows from the fact that B111 111 i jk unity.
Now suppose FC T is the well-known Clough–Tocher split of F into three subtriangles Fi :=v F , vi , vi+1 , i = 1, 2, 3, where v F is a point in the interior of F. The following result is similar to Lemma 5.2 in [1]. Lemma 2.2. Suppose that we are given all of the B-coefficients of a C 1 cubic bivariate spline s F1 F1 F1 F1 , c210 , c201 , c111 . Then for any given real number r, there exists a defined on FC T except for c300 F1 unique choice of these coefficients so that s(111 ) = r . The computation of these coefficients is stable in the sense that there is an absolute constant C such that
|ciFjk1 | C |r | + max |c | , (i jk) ∈ {(300), (210), (201), (111)}, c ∈K
where K is the set of known coefficients of s.
92
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
F1 Proof. The interpolation condition s(111 ) = r leads to the equation F1 F1 F1 F1 + 3c210 + 3c201 + 6c111 = R1 , c300
where R1 is a combination of r and the known B-coefficients of s. Suppose b1 , b2 , b3 are the barycentric coordinates of v F relative to F. Writing down the C 1 smoothness conditions across the interior edges of F1 , we are led to the linear system ⎛ F ⎞ 1 ⎞ c300 ⎞ ⎛ ⎛ 1 3 3 6 R1 ⎜ F ⎟ 1 ⎟ ⎜ 1 −b1 −b2 0 ⎟ ⎜ c210 R2 ⎟ ⎟ ⎜ ⎟⎜ ⎟, ⎜ (2.1) ⎜ ⎟=⎜ ⎝0 1 0 −b2 ⎠ ⎜ c F1 ⎟ ⎝ R3 ⎠ 201 ⎠ ⎝ 0 0 1 −b1 R4 F1 c111 where R2 , R3 , R4 are combinations of r and the coefficients in K. The determinant of this matrix is D:= − 6 − 3b1 − 3b2 − 2b1 b2 . It is easy to see that as we vary v F in T, (b1 , b2 ) runs over the set {(b1 , b2 ) : 0 b1 , b2 1, b1 + b2 1}. This gives |D| 6 for any choice of v F , which implies both the existence of a solution and the stability of its computation. 3. Partial Worsey–Farin splits The following definition is taken from [1]. Definition 3.1. Let T be a tetrahedron, and let vT be a point in its interior. Given an integer 0 m 4, let F1 , . . . , Fm be distinct faces of T, and for each i = 1, . . . , m, let v Fi be a point in the interior of Fi . Then we define the corresponding m-th order partial Worsey–Farin split m WF of T to be the tetrahedral partition obtained by the following steps: (1) connect vT to each of the four vertices of T, (2) connect vT to the points v Fi for i = 1, . . . , m, (3) connect v Fi to the three vertices of Fi for i = 1, . . . , m. The m-th order partial Worsey–Farin split of a tetrahedron results in 4 + 2m subtetrahedra, see Fig. 1. The splits 0W F and 4W F are the well-known Alfeld and Worsey–Farin splits, respectively, m see [2, Section 16.7]. We now recall an important fact about the space S13 (m W F ), where W F is an m-th order partial Worsey–Farin split of a tetrahedron T :=v1 , v2 , v3 , v4 . Theorem 3.2. ([1, Theorem 6.3]). Fix 0 m 4. Let Mm be the union of the following sets of domain points in DmW F : (1) for each i = 1, . . . , 4, D(vi )∩Ti for some tetrahedron Ti ∈ m W F containing vi ; F ; (2) for each face F of T that is not split, the point 111 Fi 3 (3) for each face F of T that has been subjected to a Clough–Tocher split, the points {111 }i=1 , where F1 , F2 , F3 are the subfaces of F. Then Mm is a minimal determining set for S13 (m W F ).
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
93
Fig. 1. Partial Worsey–Farin splits subdividing a given tetrahedron into 4, 6, 8, 10, or 12 subtetrahedra.
4. Decomposition of tetrahedral partitions In this section we describe two algorithms for decomposing a given tetrahedral partition 0 into classes of tetrahedra. These decompositions will be used to define our Lagrange interpolation pair. Algorithm 4.1. For i = 0 to 4, repeat until no longer possible: choose a tetrahedron T such that exactly i vertices of T belong to tetrahedra that have been chosen earlier. We call these vertices the marked vertices of T. Put T into the class Ai . This algorithm decomposes the partition 0 into five classes A0 , . . . , A4 . It also produces an ordering T1 , . . . , Tn of the tetrahedra of 0 . Clearly, if T is a tetrahedron in the class A0 , then all of its vertices are unmarked. Lemma 4.2. Suppose v is a marked vertex of a tetrahedron Tn of class A j with 1 j 4. Then there exists a tetrahedron Tm of class Ai with i < j such that v is a vertex of Tm . Proof. Let v be a marked vertex of Tn , and let Tm be the first of the tetrahedra T1 , . . . , Tn−1 to contain v. Then Tm must be in some class Ai with 0 i j. We claim that it cannot be in class A j , since at the time Tm was chosen, Tn would have had at most j − 1 vertices in common with the tetrahedra T1 , . . . , Tm−1 , and thus would have been chosen before Tm unless Tm is in Ai for some i < j. The following algorithm produces a different decomposition based on shared edges. Algorithm 4.3. For i = 0 to 6, repeat until no longer possible: choose a tetrahedron T such that exactly i edges of T belong to tetrahedra that have been chosen earlier. We call these edges marked edges of T. Put T into the class Bi .
94
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
This algorithm partitions 0 into seven classes B0 , . . . , B6 . It also produces an ordering 1 , . . . , T n of the tetrahedra of 0 . It is easy to see that A0 ⊆ B0 , and if T is a tetrahedron T in the class B0 , then all of its edges are unmarked. n of class B j with 1 j 6. Then Lemma 4.4. Suppose e is a marked edge of a tetrahedron T m . there exists a tetrahedron Tm of class Bi with i < j such that v is an edge of T n , and let T m be the first of the tetrahedra T 1 , . . . , T n−1 to Proof. Let e be a marked edge of T contain e. Then Tm must be in some class Bi with 0 i j. We claim that it cannot be in class m was chosen, T n would have had at most j − 1 edges in common with the B j , since at the time T m . tetrahedra T1 , . . . , Tm−1 , and thus would have been chosen before T 5. Construction of As a first step towards creating a Lagrange interpolating pair solving Problem 1.1, in this section we describe an algorithm for constructing a suitable tetrahedral partition . As a starting point, let 0 be an arbitrary tetrahedral partition with vertices at the points V. We shall construct by splitting certain of the tetrahedra of 0 . The following algorithm proceeds in two steps. In the first step we apply Clough–Tocher splits to some of the triangular faces of 0 . The choice of which faces to split is controlled by the n of the tetrahedra in 0 produced by Algorithm 4.3. The Clough–Tocher split 1 , . . . , T ordering T of a face involves inserting a point v F in the interior of F and connecting it to each of the three vertices of F. In the second step we apply partial Worsey–Farin splits to the tetrahedra of 0 that have one or more split faces. Algorithm 5.1 (Construct ). i that is not shared with any tetrahedron T j with j < i, (1) For i = 1, . . . , n, if F is a face of T apply a Clough–Tocher split to F when either two or three of the edges of F are marked i . edges of T (2) For each tetrahedron T ∈ 0 , let m be the number of its faces that have been split in step 1. If m > 0, apply an m-th order partial Worsey–Farin split to T using its incenter. In Algorithm 5.1 we did not specify how to choose the split points v F to be used in creating the Clough–Tocher splits of faces. For our purposes, these points must be chosen in a special way: (1) if F is a boundary face of 0 , choose v F to be the barycenter of F, (2) if F is an interior face of 0 , choose v F to be the intersection of F with the line connecting the incenters of the two tetrahedra sharing F. It is easy to see that none of the tetrahedra in the classes B0 or B1 is split. Tetrahedra in class B2 are either not split, or are subjected to a partial-Worsey–Farin split of order m = 1. The types of splits that may be applied to tetrahedra in the various classes are shown in Table 1, where m-WF stands for the m-th order Worsey–Farin split. In this table the symbol “–” indicates that the corresponding split does not occur. The symbol “◦” identifies cases where a tetrahedron in a class B j may be given the split in the indicated column because the tetrahedron shares one or more faces with tetrahedra in lower classes. The symbol “×” identifies cases where a tetrahedron in B j has no neighbor in a lower class.
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
95
Table 1 Possible splits for the different classes . Class
No split
1-WF
2-WF
3-WF
4-WF
B0 B1 B2 B3 B4 B5 B6
× × × ◦ – – –
– – × ×o – – –
– – – × ◦ ◦ ◦
– – – × ×o ◦ ◦
– – – – × ×◦ ×◦
Fig. 2. Points inserted on faces (black dots) by Algorithm 6.1.
6. Construction of P In this section we present an algorithm for constructing a set P of interpolation points which together with the spline space S13 () based on the tetrahedral partition of the previous section will form a Lagrange interpolation pair. The algorithm makes use of both of the ordered lists 1 , . . . , T n created by Algorithms 4.1 and 4.3. T1 , . . . , Tn and T Algorithm 6.1 (Construction of P). For i = 1, . . . , n, (1) for each edge u, v of Ti ; (a) choose the points v and (u + 2v)/3 if v is not a marked vertex of Ti ; (b) choose the points u and (2u + v)/3 if u is not a marked vertex of Ti . i that is not a face of a tetrahedron T j with j < i; (2) for each face F of T i , (a) choose the barycenter of F if F has no marked edges of T (b) choose the barycenter of the subtriangle of F with no marked edge if F has exactly two i . marked edges of T We emphasize that the interpolation points chosen by this algorithm belong to D , and all lie on edges and faces of 0 . In particular, P contains all of the points in V, i.e., the vertices of 0 . Interpolation points are chosen on faces of tetrahedra in 0 when the face is not split and contains no marked edges, or when the face is split but contains exactly two marked edges, see Fig. 2, where the thicker lines indicate marked edges, i.e. edges in common with tetrahedra appearing n . 1 , . . . , T earlier in the list T 7. The main result In this section we prove the main result of this paper, Theorem 7.3 below, which states that P and S13 () form a Lagrange interpolating pair. We begin by showing that P is a minimal determining set for S13 ().
96
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
Theorem 7.1. The set M:=P produced by Algorithm 6.1 is a minimal determining set for S13 (). Proof. Suppose that for each ∈ M, we fix the corresponding coefficient of a spline s ∈ S03 (). We now show that all other coefficients of s are consistently determined, i.e., in such a way that all C 1 smoothness conditions are satisfied so that s lies in S13 (). Step 1: It is easy to see by the construction of P that for each vertex v of 0 , there is a tetrahedron T in 0 with vertex at v such that P contains the four domain points in D T (v). Then the coefficients associated with the remaining domain points in D(v) are consistently determined by the C 1 smoothness at v. The consistency is insured by the fact that these balls do not overlap. Step 2: To show that all remaining coefficients are consistently determined, we consider the n of Algorithm 4.3, and apply Theorem 3.2. We begin with T 1 . 1 , . . . , T tetrahedra in the order T We already have determined the coefficients associated with the balls around the vertices of T1 . F 1 is in the class B0 , none of its faces is split and P contains the point 111 for each face Since T F, which means that the corresponding coefficient has been fixed. Theorem 3.2 then implies that 1 are consistently determined. The all remaining coefficients corresponding to domain points in T 1 C conditions imply that all coefficients associated with domain points in the tubes around the 1 are also consistently determined. Now suppose we have determined the coefficients edges of T 1 , . . . , T j−1 , and consider T j . Using corresponding to all of the domain points in the tetrahedra T 1 C smoothness, we get all coefficients of s corresponding to domain points in the balls around the j . To apply Theorem 3.2, we j , as well as in the tube around each marked edge of T vertices of T now show that we have already determined the coefficients corresponding to the domain points j . If F is shared with a tetrahedron T i in items (2) and (3) of the theorem. Consider a face F of T with i < j, then all coefficients corresponding to domain points on F are already known. Suppose now that F is not shared with such a tetrahedron. There are four cases: F . We have already (a) F has no marked edges. Then F is not split, and P contains the point 111 fixed the corresponding coefficient. (b) F has one marked edge e. Then F is not split, and we already know the coefficient correF since it lies in the tube around e. sponding to 111 (c) F has two marked edges. Then F has been subjected to a CT split. In this case the coefficients corresponding to the barycenters of the three subtriangles are known since the corresponding domain points are either in a tube around a marked edge, or the point is in P. (d) F has three marked edges. Then F has been subjected to a CT split. In this case the coefficients corresponding to the barycenters of the three subtriangles are known since the corresponding domain points are in tubes around the marked edges.
We have shown that the coefficients of s corresponding to domain points in the minimal determining set of Theorem 3.2 are determined, and we conclude from the theorem that the coefficients j are consistently determined. Continuing of s corresponding to all remaining domain points in T with the remaining tetrahedra, we get all coefficients of s. There is one subtle point concerning consistency that remains to be discussed. It concerns the C 1 smoothness conditions across faces of 0 . Suppose F is a face of 0 that is shared by two tetrahedra Tm and Tn , and suppose F has been subjected to a Clough–Tocher split. Then a necessary and sufficient condition that all C 1 smoothness conditions across F are satisfied is that the split point v F lie on the straight line between the interior split points vTm and vTn , see the proof of Theorem 18.11 in [2]. But we have made sure that this is the case in our construction of .
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
97
The fact that P is a minimal determining set for S13 () implies that the dimension of S13 () is equal to the number N of points in the set P, see [2, Section 17.3]. Clearly, N = 4n V + n F , where n V = #V is the number of vertices of 0 , and n F is the number of faces where an interpolation point was added to P in step 2 of Algorithm 6.1. To state our next theorem, we recall that denotes point-evaluation at the point . Theorem 7.2. Let P be the set of domain points chosen by Algorithm 6.1. Then the set N:={ }∈P is a nodal minimal determining set for S13 (). Proof. We have already observed that the dimension of S13 () is equal to #P = #N. Thus, to prove N is a nodal minimal determining set for S13 (), it suffices to show that if we fix the values of s() for all ∈ P, then all B-coefficients of s are determined. Step 1: We show that all coefficients of s associated with domain points on the edges of 0 and in the balls around the vertices of 0 are determined. To this end, we examine the tetrahedra of 0 in the order T1 , . . . , Tn defined by Algorithm 4.1. Suppose we have established the assertion for the edges and vertices of all of the tetrahedra T1 , . . . , T j−1 , and consider T j . We have already determined the coefficients of s in the balls around marked vertices of T j . Let e:=u, v be an edge of T j . Suppose first that u, v are both unmarked vertices of T j . Then P contains all four domain points on e, call them 1 , . . . , 4 . We are given the values of the univariate cubic polynomial s|e at these four points. This leads to a linear system of four equations for the four B-coefficients of s associated with 1 , . . . , 4 . It is easy to see that this system is nonsingular with determinant equal to 12 81 , independent of the length or orientation of the edge e. Now suppose that one vertex of e is marked, say u. Then we know the coefficients associated with D(u), but not those associated with D(v). But now we can use the interpolation conditions at the two points (u + 2v)/3 and v to get a linear system of 2 equations for the B-coefficients associated with these two domain points. This system is nonsingular with determinant equal to 12 81 , independent of the length or orientation of the edge e. We have shown how to compute the coefficients of s associated with all domain points on the edges of T j . We then use these to determine the coefficients associated with domain points in the ball D(v) around each unmarked vertex of T j . Step 2: We now show that all remaining coefficients of s can be computed from the interpolation 1 , . . . , T n defined by Algorithm conditions. To this end, we consider the tetrahedra in the order T 4.3. Suppose we have computed all coefficients of s associated with domain points in the tetrahedra 1 , . . . , T j−1 , and consider T j . It suffices to show how to compute the coefficients associated with T the minimal determining set Mm of Theorem 3.2. We already have the coefficients associated j and in the tubes around its marked edges. with domain points in the balls around the vertices of T We now deal with the coefficients corresponding to the remaining points in Mm . These lie on j . Let F be a face of T j . If F is shared with a tetrahedron T i with i < j, then we have faces of T already computed the coefficients of s associated with all domain points on F. Thus, it suffices to j that are not shared with a tetrahedron T i with i < j. consider faces F of T Suppose that F is not split. We need to show how to compute the coefficient associated with F F j , then P contains the point 111 111 . There are two subcases. If F contains no marked edges of T , and we can use the interpolation condition at this point to compute the corresponding coefficient, F j , then 111 lies in the tube around that edge, see Lemma 2.1. If F contains one marked edge of T and the corresponding coefficient is already known. Now suppose F has been subjected to a Clough–Tocher split. We have to show how to compute Fi the coefficients associated with the domain points 111 , where F1 , F2 , F3 are the subfaces of F.
98
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
j , then the coefficient corresponding to Fi is already known If Fi contains a marked edge of T 111 since this point lies in the tube around the marked edge. Thus, if all three edges of F are marked j , we are done. Now suppose only two of the edges of F are marked edges of T j , say edges of T F1 those contained in F2 , F3 . Then P contains the barycenter 111 of F1 , and we can apply Lemma F1 . We now apply Theorem 3.2 to determine the 2.2 to compute the corresponding coefficient c111 j . coefficients of s corresponding to all remaining domain points in T We are ready to establish the main result of the paper. N together with the spline space S1 () form a Lagrange interTheorem 7.3. The set P:={i }i=1 3 polation pair.
Proof. The proof of Theorem 7.2 shows that for any set of real numbers r1 , . . . , r N , there is a unique spline s ∈ S13 () satisfying s(i ) = ri for i = 1, . . . , N . 8. Bounds on the error of the interpolant N and S1 () are the Lagrange interpolation pair constructed above Suppose that P:={i }i=1 3 from an initial tetrahedral partition 0 of a domain ⊂ R3 . Then for every f ∈ C(), there is a unique spline s f ∈ S13 () such that
s f (i ) = f (i ), i = 1, . . . , N . This interpolation process defines a linear projector mapping C() onto S13 (). Our goal in this section is to provide a bound on f − s for smooth functions, where the error is measured in the maximum norm on . To accomplish this goal, we will apply Theorem 17.22 of [2]. To do so, we need to show that nodal minimal determining set P of Theorem 7.2 is local and stable, see Definition 17.21 of [2]. For each ∈ D , let :={ ∈ P : c depends on s()}, where {c }∈D are the B-coefficients of s f . Theorem 8.1. The nodal minimal determining set P for S13 () is local in the sense that for all ∈ D , ⊆ star (T ),
(8.1)
with = 10, where T is a tetrahedron in 0 that contains . It is also stable in the sense that there exists a constant C depending only on the smallest solid and face angles in 0 such that for all ∈ D , |c | C max | f ()|. ∈
(8.2)
Proof. For the definition of solid and face angles of a tetrahedron, see [2, p. 462]. We consider three cases. Case 1: ( lies on an edge of a tetrahedron T in 0 ). Let e:=u, v be the edge of T containing . If neither u nor v is a marked vertex of T , then P contains all four domain points on e, and the
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
99
A0 A1
A2
A3
A4
Fig. 3. Longest chain of influence of an interpolation value along the edges of 0 . Interpolation points are shown as gray and black dots, and the direction of propagation (of the value f at the black dot) into neighboring tetrahedra of higher classes Ai is indicated by the arrows. The value f at the black dot can influence the B-coefficient indicated as an open circle, which is associated with a domain point in a tetrahedron from class A4 .
corresponding coefficients depend only on the values of s at these four points, see the proof of Theorem 7.2. For these coefficients, (8.1) holds with = 0. This argument applies to all edges of tetrahedra T of class A0 . If T lies in a higher class, the situation is more complicated, and there is a certain propagation of influence. We claim that if T is of class A j , then (8.1) holds with = j. We have already established this for j = 0, and can now proceed by induction. Suppose we have established the claim for classes A0 , . . . , A j−1 , and consider T ∈ A j . If neither end of e is marked, then as before we get (8.1) with = 0. Now suppose u is marked. Let ∈ D(u). By Lemma 4.2, there exists a tetrahedron Tn in class A j−1 sharing the vertex u, and we can compute c as a combination of coefficients associated with D Tn (u). By the induction hypothesis, these can be computed from data at points lying in star j−1 (Tn ), and we conclude that (8.1) holds with = j. It remains to consider the case ∈ D(v), where u is marked but v is not. In this case c is computed by solving a 2 × 2 linear system to find the coefficients associated with the domain points (u + 2v)/2 and v. But the right-hand side of this system includes coefficients associated with domain points in the ball D(u). It follows that ⊆ star j (T ). We conclude that in all cases, if lies on any edge of a tetrahedron T , then (8.1) holds with = 4. The worst case is illustrated in Fig. 3. We now verify the stability condition (8.2) when lies on an edge of T . By the proof of Theorem 7.2, coefficients with domain points on edges of 0 are only computed in one of three ways: (1) by solving a 2 × 2 linear system whose determinant is 49 ; (2) by solving a 4 × 4 linear system whose determinant is 12 81 ; (3) by employing the C 1 smoothness at a vertex to compute a coefficient from previously computed coefficients. The constant of stability in (3) depends on the smallest solid and face angles of 0 . Case 2: ( lies on a face F of a tetrahedron T in 0 , but not on an edge). We claim that if T is of class B j , then (8.1) holds with = 4+ j. This is immediate for j = 0, since when T is of class B0 , then c is computed from Lemma 2.1, which shows that c depends on the value of s(), but also on the coefficients corresponding to domain points in the balls around the vertices of F. We showed in Case 1 that such coefficients depend only on data at points in star4 (T ). Suppose now that we have established the claim for all tetrahedra in the classes B0 , . . . , B j−1 , and consider T ∈ B j . Now c will be computed from a C 1 smoothness condition around a marked edge e, or by applying Lemma 2.1 or 2.2. Thus, c may depend on the value of s at a point in F, but more importantly can also depend on coefficients associated with points in another tetrahedron sharing the marked edge e, which by Lemma 4.4 is in a lower class. Then, applying the induction
100
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
hypothesis, we conclude that ⊆ star4+ j (T ). The worst case is for j = 6, and we conclude that (8.1) holds with = 10 for all . We now verify the stability condition (8.2) when lies on a face F of a tetrahedron T in 0 , but not on an edge. By the proof of Theorem 7.2, coefficients with domain points in the interior of faces of 0 are only computed in one of three ways: (1) by enforcing a C 1 smoothness condition around an edge, (2) by Lemma 2.1, (3) by Lemma 2.2. All three of these computations are stable, where the constant of stability in (1) depends on the smallest solid and face angles of 0 . Case 3: ( lies in the interior of a tetrahedron T in 0 ). In this case c is computed from known coefficients on the faces of T , cf. Theorem 3.2, by enforcing C 1 smoothness conditions. Thus, (8.1) holds with = 10. These computations are stable with a constant depending on the largest solid and face angles in 0 . Given a compact set B ⊆ and an integer m > 0, let | f |m,B :=
D f B
(8.3)
||=m
be the usual Sobolev seminorm, where . B denotes the infinity norm on B. For any tetrahedron T in 0 , let |T | be the length of its longest edge. The following result follows immediately from Theorem 17.22 of [2]. Theorem 8.2. Given T ∈ 0 , let T :=star10 (T ) in 0 . Then for every f ∈ C m+1 (T ) with 0 m 3,
D ( f − s f ) T K |T |m+1−|| | f |m+1,T ,
(8.4)
for all || m. If T is convex, the constant K depends only on the smallest solid and face angles of T as a subpartition of . If T is not convex, K also depends on the Lipschitz constant of the boundary of T . The global version of this theorem also holds with T and T replaced by , and with |T | replaced by the mesh size |0 | of 0 , i.e., the length of the longest edge in 0 . 9. Remarks Remark 1. Suppose V is an arbitrary set of points in R3 , and that is a tetrahedral partition with vertices V. Set P:=V, and let S:=S01 () be the space of continuous linear splines. Then clearly P and S form a Lagrange interpolation pair. It is also straightforward to create a Lagrange interpolation pair using C 0 splines of higher degree, provided we add an appropriate set of additional interpolation points. Remark 2. There is a fairly extensive history of Lagrange interpolation with bivariate splines, although even in this case, it is quite difficult to construct Lagrange interpolating pairs using C 1 splines. See [3–6,8,9].
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
101
Remark 3. Much less is known about Lagrange interpolation with C 1 trivariate splines. The first local Lagrange method was constructed in [10] using quintic splines, but only for the case where the points of V lie on a grid. With this same restriction on V, we recently constructed a Lagrange interpolation pair using cubic splines, see [1]. In [7] we constructed Lagrange interpolation pairs for an arbitrary initial point set V using both quadratic and cubic splines. However, as pointed out in the introduction, that cubic spline method did not provide optimal approximation order. Remark 4. For given V and initial tetrahedral partition 0 , there are generally many different associated Lagrange interpolation pairs. They arise by different choices of the decompositions in Section 4. Remark 5. Our construction shows that the minimal determining set M of Theorem 7.1 is stable and 1-local in the sense of Definition 17.11 of [2]. But then Theorem 17.14 of [2] applies to show that S13 () approximates functions in the Sobolev spaces Wqm+1 () with 0 m 3 and 1q ∞ to order |0 |m+1 . In particular, with m = 3, this shows that the space has full approximation power in all of the q-norms. Remark 6. With a little additional work, it can be shown that (8.1) in Theorem 8.1 holds with = 9. This is based on the fact that if T is in class B1 , then none of its faces is split, and so coefficients associated with a tube around one edge do not interact via C 1 smoothness conditions with coefficients associated with a tube around its opposite edge. This implies that if T is in class B2 , then (8.1) actually holds with = 5, and induction then gives a worst case of = 9. We note that for most tetrahedral partitions, will be smaller than 9. Remark 7. It should be emphasized that care is needed in the choice of the interior split points vT in Step 2 of Algorithm 5.1 for constructing the refined partition from 0 . To insure C 1 smoothness across interior faces of 0 , we need to choose the split point on the face to lie on the line joining the interior split points of the neighboring tetrahedra. To accomplish this, we have chosen the incenters of these tetrahedra, rather than the more natural barycenters, since then the line joining them always intersects the common face F at an interior point, see Lemma 16.24 in [2]. Remark 8. It was shown in Corollary 6.2 of [1] that a C 1 cubic spline defined on a partialWorsey–Farin split of a tetrahedron is actually C 2 at the interior split point vT . It follows that the splines in our space S13 () enjoy this property at all such vertices of . Remark 9. Given a set of n vertices V, an associated initial triangulation 0 can be computed with standard algorithms. Once we have 0 , the operation count for constructing our Lagrange interpolation pair is O(n). Finding the coefficients of an interpolating spline is also of linear complexity in n. Remark 10. In step 2 of Algorithm 5.1 we apply a partial-Worsey–Farin split to T whenever the number of CT-split faces m is positive. We could have also applied a 0-th order split (also called the Alfeld split) in the case when m = 0, but this would lead to a refinement with more tetrahedra, so instead we do not split the tetrahedron at all when m = 0. Remark 11. Trivariate local Lagrange interpolation methods are useful for the construction and reconstruction of volumetric models, as well as for scattered data fitting problems. A major
102
G. Hecklin et al. / Journal of Approximation Theory 160 (2009) 89 – 102
advantage is that they only require function values, but no derivatives. For example, the method here could be used in a two-stage process, where in the first stage we construct a C 0 linear spline that interpolates at the vertices of a very fine tetrahedral partition. In the second stage, we construct a C 1 cubic spline on a coarser tetrahedral subpartition , where the data needed for interpolation are taken directly from the linear spline. Because of the optimal approximation order of the interpolation methods, it is natural to choose the partitions such that |1 |2 is about |0 |4 . For numerical results based on this idea in the bivariate setting, see [9]. Remark 12. As shown in Theorem 17.17 of [2], there is a natural stable local basis { }∈M for S13 () associated with the minimal determining set M of Theorem 7.1. Since M is 1-local, each has support at most star(T ), where T is a tetrahedron of 0 containing . Remark 13. There is a different stable local basis { }∈M for S13 () associated with the nodal minimal determining set N of Theorem 7.2. Since N is just the set of point-evaluation functionals at the points of P, the are Lagrange basis splines in the sense that () = , for all , ∈ P. Theorem 8.1 implies that each has support at most star10 (T ), where T is a tetrahedron of 0 containing . For typical tetrahedral partitions, most if not all of these basis functions will have much smaller supports. Remark 14. In the theory of trivariate splines, nodal determining sets are usually constructed using point-evaluation of both function values and derivatives. Here we have restricted ourselves to the use of function values only in order to get a Lagrange interpolation method rather than a Hermite interpolation method. References [1] G. Hecklin, G. Nürnberger, L.L. Schumaker, F. Zeilfelder, A local Lagrange interpolation method based on C 1 cubic splines on Freudenthal partitions, Math. Comp. 77 (2008) 1017–1036. [2] M.-J. Lai, L.L. Schumaker, Splines on Triangulations, Cambridge University Press, Cambridge, 2007. [3] G. Nürnberger, V. Rayevskaya, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation with C 2 splines of degree seven on triangulations, in: M. Neamtu, E. Saff (Eds.), Advances in Constructive Approximation, Nashboro Press, Brentwood, TN, 2004, pp. 345–370. [4] G. Nürnberger, V. Rayevskaya, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation with bivariate splines of arbitrary smoothness, Constr. Approx. 23 (2006) 33–59. [5] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Local Lagrange interpolation by bivariate C 1 cubic splines, in: T. Lyche, L.L. Schumaker (Eds.), Mathematical Methods for Curves and Surfaces III, Oslo, 2000, Vanderbilt University Press, Nashville, 2000, pp. 393–404. [6] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Lagrange interpolation by C 1 cubic splines on triangulated quadrangulations, Adv. Comp. Math. 21 (2004) 357–380. [7] G. Nürnberger, L.L. Schumaker, F. Zeilfelder, Two Lagrange interpolation methods based on C 1 splines on tetrahedral partitions, in: C.K. Chui, M. Neamtu, L.L. Schumaker (Eds.), Approximation Theory XI: Gatlinburg 2004, Nashboro Press, Brentwood, 2005, pp. 327–344. [8] G. Nürnberger, F. Zeilfelder, Developments in bivariate spline interpolation, J. Comp. Appl. Math. 121 (2000) 125–152. [9] G. Nürnberger, F. Zeilfelder, Lagrange interpolation by bivariate C 1 splines with optimal approximation order, Adv. Comp. Math. 21 (2004) 381–419. [10] L.L. Schumaker, T. Sorokina, C 1 quintic splines on type-4 tetrahedral partitions, Adv. Comp. Math. 21 (2004) 421–444.