- Email: [email protected]
C1 trivariate polynomial interpolation
Computer Aided Geometric Design 4 (1987) 237-244 North-Holland
237
C a trivariate polynomial interpolation K.L. RESCORLA Department of Mathematics, Eastern Michigan University, Ypsilanti, MI 48197, U.S.A. Received July 1985 Revised July 1987 Abstract. A trivariate polynomial interpolating to function and derivative values through order four at various points on a tetrahedron is constructed. This interpolant is represented in terms of a monomial basis which triangularizes the Vandermonde matrix. Then, given a domain in R 3 which has been tessellated into tetrahedra, the piecewise polynomial over this domain has C 1 continuity. Keywords. Polynomial, interpolation, trivariate, tetrahedron, continuously differentiable.
1. Introduction
In this work a trivariate polynomial of degree nine interpolating to data on a tetrahedron is constructed. It is represented in terms of a monomial basis which triangularizes the Vandermonde matrix. This polynomial can be used to define a continuously differentiable piecewise function on an arbitrary triangulated domaim in R 3. There are numerous applications for such an interpolation method. [Webb '80] gives an example where pollen percentages are modeled as a function of latitude, longitude and time. [Stead '83] discusses an application in geology whereby soil characteristics are given as a function of latitude longitude and depth. Applications for smooth trivariate interpolation methods also arise in mining. See [Clark '77a, b] and [Pauncz & Johnson '76]. Finally, R. Franke, in a recent communication, discussed the application of modeling gas motion in engine cylinders. The following notation is used. Let X220 be the vector space of trivariate polynomials of degree nine. The dimension of X220 is 220. Let L 1. . . . . L220 denote given linear functionals on X. Let b 1, b 2, b3, b4 be barycentric coordinates with respect to the tetrahedron with vertices V1, V2, 1,'3, V4 and denote the edges of the tetrahedron by eji=Vj-Vi,
i,j=1,2,3,4,
i4:j.
The goal is to construct functions x l , . . . , x220 in X220 satisfying
Lix j = 3ij
whenever i < j , j = 1 , . . . , 2 2 0 .
(1)
The properties (1) are called the semicardinalproperties. They imply the Vandermonde matrix is lower triangular with diagonal entries equal to one. The following statements are simple linear algebra results under the assumption that x l , . . . , x220 satisfying (1) exist. 1. Both sets ( L 1. . . . , L 2 2 0 } and ( x I . . . . , X220} are linearly independent: This implies the x~ form a basis of X. Such a basis is called a semicardinal basis. 0167-8396/87/$3.50 © 1987, Elsevier Science Publishers B.V. (North-Holland)
K.L. Rescorla / C’ trivariate interpolation
238
2. The interpolation problem of finding x E X,,, such that prespecified parameters) for i = 1,. , . ,220 is uniquely solvable.
Lix = d,
(these
dj are are
3. The solution x may be written in the form x = UlX1 + - * * +a,,,x,,,
(2)
where the a, are given by the forward substitution i-l a, = d,,
aj=di-
c aj(L,x,)
for i=2
,...,
220.
j=l
(3)
2. The linear functionals The goal in this section is to define the 220 linear functionals L,, . . . , L,,, on X,,, in such a way that the resulting piecewise interpolant will have C’ continuity. Let nikj denote the direction in the face of the tetrahedron determined by the vertices F, I$ and V, and normal to edge ekj. The direction niki is defined by the equation
nikj = V; - ajkjj?$ - aikjkVk
(4)
where the scalars a;kjj and aikjk are chosen to satisfy the equations nikj
-
ekj
=
aikjj
0,
+
aikjk
=
1. (9
Let denote the inward normal direction to the face determined by I$, V,, and V,. It is defined by the equation n;jk[
nijkl
=
y
-
ajjt$
-
a;,Vk
-
ailv[ (6)
where the scalars are determined nijk[
’
ekj
=
0,
*
by the equations elk --0,
nijkl’
(7)
a,j+ajk+ai[=l.
Let
f
be a primitive function;
then the given data are of the following forms.
f (O?
(8) r+s+t a f
r+s+ t= I
a(e~i)‘a(eki)“a(ej;)’
[
1
af anijk (v/+v,./z’ a”f a(nijk)‘a(n,k>”
f(&+ af [
a(nijkl)
‘k+
1
(9)
(10)
1 (zv,+
&)/3),
1,...,4,
v,’
v,),3’
r +
’ =
2y
(11) 02)
03)
K.L. Rescorla / C 1 trivariate interpolation
239
r + s = 1, 2,
(14)
or+,+1f 3(etj)ro(ekj)s3 (nijkt) (v~+vk+v,)/3 f ( ( V , + Vj + Vk + V t ) / 4 ) ,
(15)
[
(16)
3~] (v,+ ~ + vk + v~)/4
It is now possible to explicitly define the 220 linear functionals. They fall naturally into thirteen categories. Although the ordering within any one category is arbitrary, it is assumed to be fixed. 1. Let L l f , . . . , L 4 f denote the four functions values, f(V,), at the four vertices. 2. Let L 5 f , . . . , L16f denote the twelve derivatives of the form (9) with r + s + t = 1. 3. Let L17f . . . . . L40f denote the twenty-four derivatives of the form (9) with r + s + t = 2. 4. Let L4a f ..... Lso f denote the forty derivatives Qf the form (9) with r + s + t = 3. 5. Let L81 f .... , L14of denote the sixty derivatives of the form (9) with r + s + t = 4. 6. Let L l n l f , . . . , La52f denote the twelve derivatives of the form (10). 7. Let Lls3f .... , L188f denote the thirty-six derivatives of the form (11). 8. Let La89f,... , L192f denote the four function values of the form (12). 9. Let L 1 9 3 f , . . . , Z 1 9 6 f denote the four derivatives of the form (13). 10. Let Z 1 9 7 f , . . . , L2o4f denote the eight derivatives of the form (14) with r + s = 1. 11. Let L2osf,..., Ln6 f denote the twelve derivatives of the form (14) with r + s = 2. 12. Let L217f denote the function value (15) at the centroid of the tetrahedron. 13. Let L2a8f,... , L220f denote the three derivatives of the form (16). This completes the definition of the linear functionals.
3. The semicardinal basis T h e goal in this section is to construct a semicardinal basis corresponding to the 220 linear functionals defined in the previous section. The following l e m m a is helpful.
Lemma 1. Let g be a differentiable real-valued function on R 3 and consider the function bqg for some fixed i. I f 1. D is a directional derivative operator involving only directions in which b i is constant, 2. V ~ R 3 with bi(V ) = O, 3. q > 1, then for any direction e,
240
K.L. Rescorla / C / trivariate interpolation
Proof
+b~-~e
)]
[ _1[ Obi Og)] D~ qg-~e + b i ~ e , j v
v
= [ bq
(byl)
= 0
(by 2 a n d 3).
[]
The semicardinal basis functions are n o w defined. A brief discussion as to w h y these functions satisfy the semicardinal properties (1) is given. F o r a m o r e formal a n d detailed p r o o f see [Rescorla '85]. 1. For p = 1 . . . . . 4, Lpf is of the f o r m f(V/). C o r r e s p o n d i n g to this Lgxj= 8~j for i < j , j = 1 , . . . , 4 .
Lpf is of the f o r m (9). C o r r e s p o n d i n g to this Lp define
2. For p = 5 . . . . ,140, m
Lp define Xp = bi. T h e n
r
s
l
blbkb) m = 1 + m a x ( r, s, t } , Xp = r!s!t! ' r + s + t = l ..... 4. be
(18)
N o t e that this function has a zero of order at least r + s + t + 1 at the vertices Vj, V/, , a n d Vt. It has a zero of order r + s + t at V~. These facts imply the semicardinal properties L~xj = 6 , / f o r i
Lpf is of the f o r m (10). C o r r e s p o n d i n g to this Lp define
3. F o r p = 1 4 1 , . . . , 152, 8
4
4
xp = 2 bibj bk.
(19)
N o t e that xp has a zero of order at least five at each vertex so that all derivatives t h r o u g h order four vanish at all vertices. Also Xp has a zero of order at least four along all edges except edge eka on which it has a zero of o r d e r one. This implies all first-order derivatives vanish along all edges except edge ekj. These properties imply the semicardinal properties L~xj = ~ij for i < j , j= 1,...,152. 4. Next, consider Lp define
x, =
Lpf, for p = 1 5 3 , . . . , 158. Each Lpf is of the f o r m (11). C o r r e s p o n d i n g to this 37bib,bab3k(bj- ½) 8r!s!
(20)
N o t e that this xp has a zero of o r d e r at least five at each vertex. It has a zero of order at least two along all edges so that all first-order derivatives vanish along all edges. In fact, it has a zero of order at least four along all edges except edge ekj SO that all derivatives of order two vanish along all edges except edge ekj. O n edge ekj at the p o i n t (Vj + 2Vk)/3 there is a zero of order three. These conditions are sufficient to i m p l y the semicardinal properties Lixj = 6ij for i < j , j = 1,...,188. 5. For p = 189 . . . . ,192, Xp
Lpf is of the f o r m (12). C o r r e s p o n d i n g to this Lp define
m- " 1 9 A 3 ~ 3 A 3
.,, uj ~,kVl .
(21)
This function has a zero of order at least six at each vertex and a zero of order at least three along each edge. It has the value o n e at (Vj + Vk + V+)/3 a n d zero at the centroids of the other three faces. These Conditions i m p l y the semicardinal properties Lixj = 8,j for i < j , j = 1 , . . . , 192.
K.L. Rescorla / C 1 trivariate interpolation
241
6. N o w consider L p f for p = 1 9 3 , . . . , 196. Each Lp is of the f o r m (13). C o r r e s p o n d i n g to this Lp define Xp b y 2 2 Xp = 3 6 bib)2 bkb I•
(22)
This function has a zero of order at least five at each vertex, and it has a zero of order at least three along each edge. It has a zero of order two on each face except the face d e t e r m i n e d by the vertices Vj, Vk, and V~ on which it has a zero of order one. These conditions are sufficient to imply the semicardinal properties L~xj = ~ / for i < j , j = 1 . . . . ,196. 7. F o r p = 1 9 7 , . . . , 204, L p f is of the f o r m (14) with r + s = 1. C o r r e s p o n d i n g to this Lp define
Xp = 36bibZb~b2( b, - ½)s( b , - ½)'.
(23)
This function has a zero of order at least five at each vertex, a zero of o r d e r at least three along each edge a n d a zero of order at least two at the centroid of each face. Also, L e m m a 1 implies that Lixp = 0 for i = 1 9 7 , . . . , 204, i 4: p. Hence, the semicardinal properties Lixj = 6~j for i < j , j - - 1 , . . . , 204 are satisfied. 8. F o r p = 2 0 5 , . . . , 216, L p f is of the form (14) with r + s = 2. C o r r e s p o n d i n g to this Lp define
36bibZb~b~ ( bk _ ½)s ( b t _ ½)r Xp = r!s!
(24)
This function has a zero of order at least five at each vertex, a zero of o r d e r at least three along each edge and a zero of order at least two at the centroid of each face. L e m m a I together with the fact that there is a zero of order three at the centroid of the face d e t e r m i n e d by Vj, Vk, a n d Vt imply that Lixp = 0 for i = 197,...,216, i =/=p. These conditions imply the semicardinal properties Lix j = 6ij for i < j , j = 1 , . . . , 216. 9. Corresponding to L217f = f((V1 + V2 + V3 + V4)/4 ) define X217 :
4
8222
bi
bj bkb t2 .
(25)
This function has a zero of order six at each vertex, a zero of order four along each edge a n d a zero of order two on each face. L e m m a 1 implies that Lix217 = 0 for i = 1 9 7 , . . . , 216. Hence, the semicardinal properties Lix j = 6ij for i < j , j = 1 , . . . , 2 1 7 are satisfied. 10. F o r p = 2 1 8 , . . . , 2 2 0 , L p f is of the form (16). C o r r e s p o n d i n g to this
X p = 4 8 b2 2i2b j b k b2, ( b j - ] ) .
Lp define (26)
The semicardinal basis has been defined.
4. Precision and smoothness A s s u m e that arbitrarily spaced points in R 3 have been given, that these points have b e e n tessellated into tetrahedra and that data, as defined in Section 2, at vertices, edges, faces a n d centroids of the tetrahedra have been taken from a primitive function f . D e n o t e the peicewise interpolant to these data by Pf. 1. f f f ~ X22o, then P f = f. 2. P/has C 1 continuity whenever f does.
Theorem 2.
K.L. Rescorla /
242
C1
trivariate interpolation
Proof. If f is in )(220, then both f and Pf are solutions, in 3(220, to the interpolation problem. Since the solution, in X220, is unique, it follows that Pf = f. To show 2, consider the interpolants over two adjacent tetrahedra. The restrictions of these polynomials to the common face is a bivariate degree-nine polynomial and is uniquely determined by the following fifty-five linearly independent conditions: three of the form (8), forty-two of the form (9), three of the form (10), six of the form (11) and one of the form (12). Hence, the piecewise interpolant has C O continuity. Since the restrictions to the common face agree, derivatives of all orders and in any direction parallel to the face agree on the face. Note the derivative in the inward direction normal to the face restricted to the face is a bivariate polynomial of degree eight and is uniquely determined by the following forty-five linearly independent conditions: thirty of the form (9), three of the form (10), six of the form (11), one of the form (13) and five of the form (14). Hence, the piecewise interpolant has C 1 continuity as desired. The second-order derivative normal to a face restricted to the face is a bivariate polynomial of degree seven and requires thirty-six linearly independent functionals to be uniquely determined. However, the given data specify only twenty-three conditions on this derivative. So, in general, the piecewise interpolant will not be twice continuously differentiable. []
5. Examples To test the interpolant, P, the unit cube [0,1] 3 w a s tessellated into tetrahedra as shown in Fig. 1. Data as defined in Section 2 were taken from the triogonometric primitive
fa(x, Y, z ) = cos(3.1nx) c o s ( y - 0.5) sin(3.1n(z - 0.5)).
(27)
This primitive was used by Stead [Barnhill & Stead '84] to compare trivariate interpolants. The 0.2 contour of fl is shown in Fig. 2, and the corresponding contour of Pfl is shown in Fig. 3. These two figures are identical within the accuracy of the graphics devices which produced them. This is due to the large amount of data from fl to which Pf interpolates. The contouring was accomplished by subdividing the unit cube into subcubes. The function was evaluated at each corner of each subcube allowing linear interpolants to be contructed along each edge of each subcube. These linear interpolants were then used to define piecewise linear approximations, in the faces of the subcubes, to the true contours.
Fig. 1. Tessellation of the unit cube into tetrahedra.
K.L. Rescorla / C l trivariate interpolation
Fig. 2. The 0.2 contour of fl-
243
Fig. 3. The corresponding contour of Pfl.
The precision of the implementation of P was verified to be )(220 by choosing f2 in X220 with 220 nonzero randomly generated coefficients and observing that Pf2 ( x, y, z) = f2 ( x, y, z) for arbitrary (x, y, z) in the domain. Also, it was verified that the implementation of Pfl is C 1 continuous by MICROSCOPE, a multivariate function analyzer [Alfeld & Harris '84].
6. Conclusions A tetrahedral interpolant, P, with the following properties has been constructed. 1. 2. 3. 4.
The required data are C 4. The smoothness is C 1. The precision set is all polynomials of degree less than or equal to nine. The interpolant is a trivariate degree-nine piecewise polynomial. It is represented in terms of a monomial basis which triangularized the Vandermonde matrix. Hence, the coefficients are given by a simple recursion.
Previous to this work, several tetrahedral interpolants were known. [Alfeld '84a] constructs a C 1 scheme for C 1 data by forming convex combinations of projectors. [Barnhill & Little '84] describe a C 1 scheme defined as a boolean sum of projectors. Due to the nature of boolean sums, C 1 data and some auxiliary data are required. [Alfeld '85b] gives an algorithm for constructing an interpolant on simplices of arbitrary dimension having arbitrary smoothness, a C 1 tetrahedral scheme being a special case. All three of these schemes when discretized, are piecewise rational. However, polynomial intrpolants are desirable for a number of reasons. To begin, they are easily evaluated. Secondly, [Alfeld '85a] describes a technique for generating unknown derivative data via functional minimization. This well-motivated method requires differentiation and integration of the interpolant. These operations are greatly simplified when the interpolant is polynomial rather than rational. Thirdly, the fact that polynomials are easy to integrate makes them useful in the finite element method. Previous to this work, [Alfeld '84b] constructed a C 1 piecewise polynomial interpolant. This is a Clough-Tocher scheme: That is, each tetrahedron is subdivided into four subtetrahedra and a trivariate polynomial is defined on each subtetrahedron. Alfeld's scheme has several
244
K.L. Rescorla /
C1
trivariate interpolation
advantages over P. To begin, it requires only C 2 data versus C 4. It is degree-five rather than degree-nine. Finally, it is represented in terms of the Bernstein basis functions which are particularly efficient to evaluate as demonstrated in [Farin '83]. However, on the side of P, there is simplicity. Subdividing the tetrahedra for the Clough-Tocher scheme complicates certain data structures, and evaluation requires the minor additional task of determining in which subtetrahedron the point of evaluation lies. The coefficients of P are defined by a simple forward substitution which is easy to program, whereas the coefficients for the Bernstein basis functions are given recursively by rather lengthy formulas which must be explicitly written. Once programmed, though, the Clough-Tocher scheme is more efficient to evaluate than P due to its lower degree (and hence, fewer terms) and Bernstein basis functions. One application where P could be used and the Clough-Tocher scheme could not is in the generation of high-order (higher than two) derivative data.
Acknowledgements This research was supported in part by the Department of Energy with Contract No. DE-AC02-82ER-12046 to the University of Utah. Thanks is given to Robert Barnhill for providing the stimulating atmosphere which spawned this problem and to Peter Alfeld for his support.
References Alfeld, P. (1984a), A discrete C 1 interpolant for tetrahedral data, Rocky Mountain J. Math. 14, 5-16. Alfeld, P. (1984b), A trivariate Ctough-Tocher scheme for tetrahedral data, Computer Aided Geometric Design 1, 169-181. Alfeld, P. (1985a), Derivative generation by functional minimization, Computer Aided Geometric Design 2, 281-296. Alfeld, P. (1985b), Multivariate perpendicular interpolation, SIAM J. Numer. Anal. 22 (1), 95-106. Alfeld, P. and Harris, W. (1984), MICROSCOPE: a software system for multivariate analysis, MRC Technical Summary Report #2701. Barnhill, R.E. and Little, F.F. (1984), Three- and four-dimensional surfaces, Rocky Mountain J. Math. 14, 77-102. Barnhill, R.E. and Stead, S.E. (1984), Multistage trivariate surfaces, Rocky Mountain J. Math. 14, 103-118. Clark, I. (1977a), Practical kriging in three dimensions, Computers and Geosciences 3, 173-180. Clark, I. (1977b), SNARK - a four-dimensional trend surface computer program, Computers and Geosciences 3, 283-308. Farin, G. (1983), Smooth interpolation to scattered 3D data, in: R.E. Bamhill and W. Boehm, eds., Surfaces in Computer Aided Geometric Design, North-Holland, Amsterdam. Paun¢e, I. and Johnson, K.R. (1976), Three-dimensional interpolation, in: D.F. Merriam, ed., Recent Advances in Geomathematics, 135-155. Rescorla, K. (1985), Multivariate interpolation, Ph.D. Thesis, Department of Mathematics, University of Utah. Stead, S.E. (1983), Smooth multistage multivariate approximation, Ph.D. Thesis, Brown University. Webb, T. (1980), The reconstruction of climatic sequences from botanical data, J. Interdisciplinary History, X:4, 749-772.
237
C a trivariate polynomial interpolation K.L. RESCORLA Department of Mathematics, Eastern Michigan University, Ypsilanti, MI 48197, U.S.A. Received July 1985 Revised July 1987 Abstract. A trivariate polynomial interpolating to function and derivative values through order four at various points on a tetrahedron is constructed. This interpolant is represented in terms of a monomial basis which triangularizes the Vandermonde matrix. Then, given a domain in R 3 which has been tessellated into tetrahedra, the piecewise polynomial over this domain has C 1 continuity. Keywords. Polynomial, interpolation, trivariate, tetrahedron, continuously differentiable.
1. Introduction
In this work a trivariate polynomial of degree nine interpolating to data on a tetrahedron is constructed. It is represented in terms of a monomial basis which triangularizes the Vandermonde matrix. This polynomial can be used to define a continuously differentiable piecewise function on an arbitrary triangulated domaim in R 3. There are numerous applications for such an interpolation method. [Webb '80] gives an example where pollen percentages are modeled as a function of latitude, longitude and time. [Stead '83] discusses an application in geology whereby soil characteristics are given as a function of latitude longitude and depth. Applications for smooth trivariate interpolation methods also arise in mining. See [Clark '77a, b] and [Pauncz & Johnson '76]. Finally, R. Franke, in a recent communication, discussed the application of modeling gas motion in engine cylinders. The following notation is used. Let X220 be the vector space of trivariate polynomials of degree nine. The dimension of X220 is 220. Let L 1. . . . . L220 denote given linear functionals on X. Let b 1, b 2, b3, b4 be barycentric coordinates with respect to the tetrahedron with vertices V1, V2, 1,'3, V4 and denote the edges of the tetrahedron by eji=Vj-Vi,
i,j=1,2,3,4,
i4:j.
The goal is to construct functions x l , . . . , x220 in X220 satisfying
Lix j = 3ij
whenever i < j , j = 1 , . . . , 2 2 0 .
(1)
The properties (1) are called the semicardinalproperties. They imply the Vandermonde matrix is lower triangular with diagonal entries equal to one. The following statements are simple linear algebra results under the assumption that x l , . . . , x220 satisfying (1) exist. 1. Both sets ( L 1. . . . , L 2 2 0 } and ( x I . . . . , X220} are linearly independent: This implies the x~ form a basis of X. Such a basis is called a semicardinal basis. 0167-8396/87/$3.50 © 1987, Elsevier Science Publishers B.V. (North-Holland)
K.L. Rescorla / C’ trivariate interpolation
238
2. The interpolation problem of finding x E X,,, such that prespecified parameters) for i = 1,. , . ,220 is uniquely solvable.
Lix = d,
(these
dj are are
3. The solution x may be written in the form x = UlX1 + - * * +a,,,x,,,
(2)
where the a, are given by the forward substitution i-l a, = d,,
aj=di-
c aj(L,x,)
for i=2
,...,
220.
j=l
(3)
2. The linear functionals The goal in this section is to define the 220 linear functionals L,, . . . , L,,, on X,,, in such a way that the resulting piecewise interpolant will have C’ continuity. Let nikj denote the direction in the face of the tetrahedron determined by the vertices F, I$ and V, and normal to edge ekj. The direction niki is defined by the equation
nikj = V; - ajkjj?$ - aikjkVk
(4)
where the scalars a;kjj and aikjk are chosen to satisfy the equations nikj
-
ekj
=
aikjj
0,
+
aikjk
=
1. (9
Let denote the inward normal direction to the face determined by I$, V,, and V,. It is defined by the equation n;jk[
nijkl
=
y
-
ajjt$
-
a;,Vk
-
ailv[ (6)
where the scalars are determined nijk[
’
ekj
=
0,
*
by the equations elk --0,
nijkl’
(7)
a,j+ajk+ai[=l.
Let
f
be a primitive function;
then the given data are of the following forms.
f (O?
(8) r+s+t a f
r+s+ t= I
a(e~i)‘a(eki)“a(ej;)’
[
1
af anijk (v/+v,./z’ a”f a(nijk)‘a(n,k>”
f(&+ af [
a(nijkl)
‘k+
1
(9)
(10)
1 (zv,+
&)/3),
1,...,4,
v,’
v,),3’
r +
’ =
2y
(11) 02)
03)
K.L. Rescorla / C 1 trivariate interpolation
239
r + s = 1, 2,
(14)
or+,+1f 3(etj)ro(ekj)s3 (nijkt) (v~+vk+v,)/3 f ( ( V , + Vj + Vk + V t ) / 4 ) ,
(15)
[
(16)
3~] (v,+ ~ + vk + v~)/4
It is now possible to explicitly define the 220 linear functionals. They fall naturally into thirteen categories. Although the ordering within any one category is arbitrary, it is assumed to be fixed. 1. Let L l f , . . . , L 4 f denote the four functions values, f(V,), at the four vertices. 2. Let L 5 f , . . . , L16f denote the twelve derivatives of the form (9) with r + s + t = 1. 3. Let L17f . . . . . L40f denote the twenty-four derivatives of the form (9) with r + s + t = 2. 4. Let L4a f ..... Lso f denote the forty derivatives Qf the form (9) with r + s + t = 3. 5. Let L81 f .... , L14of denote the sixty derivatives of the form (9) with r + s + t = 4. 6. Let L l n l f , . . . , La52f denote the twelve derivatives of the form (10). 7. Let Lls3f .... , L188f denote the thirty-six derivatives of the form (11). 8. Let La89f,... , L192f denote the four function values of the form (12). 9. Let L 1 9 3 f , . . . , Z 1 9 6 f denote the four derivatives of the form (13). 10. Let Z 1 9 7 f , . . . , L2o4f denote the eight derivatives of the form (14) with r + s = 1. 11. Let L2osf,..., Ln6 f denote the twelve derivatives of the form (14) with r + s = 2. 12. Let L217f denote the function value (15) at the centroid of the tetrahedron. 13. Let L2a8f,... , L220f denote the three derivatives of the form (16). This completes the definition of the linear functionals.
3. The semicardinal basis T h e goal in this section is to construct a semicardinal basis corresponding to the 220 linear functionals defined in the previous section. The following l e m m a is helpful.
Lemma 1. Let g be a differentiable real-valued function on R 3 and consider the function bqg for some fixed i. I f 1. D is a directional derivative operator involving only directions in which b i is constant, 2. V ~ R 3 with bi(V ) = O, 3. q > 1, then for any direction e,
240
K.L. Rescorla / C / trivariate interpolation
Proof
+b~-~e
)]
[ _1[ Obi Og)] D~ qg-~e + b i ~ e , j v
v
= [ bq
(byl)
= 0
(by 2 a n d 3).
[]
The semicardinal basis functions are n o w defined. A brief discussion as to w h y these functions satisfy the semicardinal properties (1) is given. F o r a m o r e formal a n d detailed p r o o f see [Rescorla '85]. 1. For p = 1 . . . . . 4, Lpf is of the f o r m f(V/). C o r r e s p o n d i n g to this Lgxj= 8~j for i < j , j = 1 , . . . , 4 .
Lpf is of the f o r m (9). C o r r e s p o n d i n g to this Lp define
2. For p = 5 . . . . ,140, m
Lp define Xp = bi. T h e n
r
s
l
blbkb) m = 1 + m a x ( r, s, t } , Xp = r!s!t! ' r + s + t = l ..... 4. be
(18)
N o t e that this function has a zero of order at least r + s + t + 1 at the vertices Vj, V/, , a n d Vt. It has a zero of order r + s + t at V~. These facts imply the semicardinal properties L~xj = 6 , / f o r i
Lpf is of the f o r m (10). C o r r e s p o n d i n g to this Lp define
3. F o r p = 1 4 1 , . . . , 152, 8
4
4
xp = 2 bibj bk.
(19)
N o t e that xp has a zero of order at least five at each vertex so that all derivatives t h r o u g h order four vanish at all vertices. Also Xp has a zero of order at least four along all edges except edge eka on which it has a zero of o r d e r one. This implies all first-order derivatives vanish along all edges except edge ekj. These properties imply the semicardinal properties L~xj = ~ij for i < j , j= 1,...,152. 4. Next, consider Lp define
x, =
Lpf, for p = 1 5 3 , . . . , 158. Each Lpf is of the f o r m (11). C o r r e s p o n d i n g to this 37bib,bab3k(bj- ½) 8r!s!
(20)
N o t e that this xp has a zero of o r d e r at least five at each vertex. It has a zero of order at least two along all edges so that all first-order derivatives vanish along all edges. In fact, it has a zero of order at least four along all edges except edge ekj SO that all derivatives of order two vanish along all edges except edge ekj. O n edge ekj at the p o i n t (Vj + 2Vk)/3 there is a zero of order three. These conditions are sufficient to i m p l y the semicardinal properties Lixj = 6ij for i < j , j = 1,...,188. 5. For p = 189 . . . . ,192, Xp
Lpf is of the f o r m (12). C o r r e s p o n d i n g to this Lp define
m- " 1 9 A 3 ~ 3 A 3
.,, uj ~,kVl .
(21)
This function has a zero of order at least six at each vertex and a zero of order at least three along each edge. It has the value o n e at (Vj + Vk + V+)/3 a n d zero at the centroids of the other three faces. These Conditions i m p l y the semicardinal properties Lixj = 8,j for i < j , j = 1 , . . . , 192.
K.L. Rescorla / C 1 trivariate interpolation
241
6. N o w consider L p f for p = 1 9 3 , . . . , 196. Each Lp is of the f o r m (13). C o r r e s p o n d i n g to this Lp define Xp b y 2 2 Xp = 3 6 bib)2 bkb I•
(22)
This function has a zero of order at least five at each vertex, and it has a zero of order at least three along each edge. It has a zero of order two on each face except the face d e t e r m i n e d by the vertices Vj, Vk, and V~ on which it has a zero of order one. These conditions are sufficient to imply the semicardinal properties L~xj = ~ / for i < j , j = 1 . . . . ,196. 7. F o r p = 1 9 7 , . . . , 204, L p f is of the f o r m (14) with r + s = 1. C o r r e s p o n d i n g to this Lp define
Xp = 36bibZb~b2( b, - ½)s( b , - ½)'.
(23)
This function has a zero of order at least five at each vertex, a zero of o r d e r at least three along each edge a n d a zero of order at least two at the centroid of each face. Also, L e m m a 1 implies that Lixp = 0 for i = 1 9 7 , . . . , 204, i 4: p. Hence, the semicardinal properties Lixj = 6~j for i < j , j - - 1 , . . . , 204 are satisfied. 8. F o r p = 2 0 5 , . . . , 216, L p f is of the form (14) with r + s = 2. C o r r e s p o n d i n g to this Lp define
36bibZb~b~ ( bk _ ½)s ( b t _ ½)r Xp = r!s!
(24)
This function has a zero of order at least five at each vertex, a zero of o r d e r at least three along each edge and a zero of order at least two at the centroid of each face. L e m m a I together with the fact that there is a zero of order three at the centroid of the face d e t e r m i n e d by Vj, Vk, a n d Vt imply that Lixp = 0 for i = 197,...,216, i =/=p. These conditions imply the semicardinal properties Lix j = 6ij for i < j , j = 1 , . . . , 216. 9. Corresponding to L217f = f((V1 + V2 + V3 + V4)/4 ) define X217 :
4
8222
bi
bj bkb t2 .
(25)
This function has a zero of order six at each vertex, a zero of order four along each edge a n d a zero of order two on each face. L e m m a 1 implies that Lix217 = 0 for i = 1 9 7 , . . . , 216. Hence, the semicardinal properties Lix j = 6ij for i < j , j = 1 , . . . , 2 1 7 are satisfied. 10. F o r p = 2 1 8 , . . . , 2 2 0 , L p f is of the form (16). C o r r e s p o n d i n g to this
X p = 4 8 b2 2i2b j b k b2, ( b j - ] ) .
Lp define (26)
The semicardinal basis has been defined.
4. Precision and smoothness A s s u m e that arbitrarily spaced points in R 3 have been given, that these points have b e e n tessellated into tetrahedra and that data, as defined in Section 2, at vertices, edges, faces a n d centroids of the tetrahedra have been taken from a primitive function f . D e n o t e the peicewise interpolant to these data by Pf. 1. f f f ~ X22o, then P f = f. 2. P/has C 1 continuity whenever f does.
Theorem 2.
K.L. Rescorla /
242
C1
trivariate interpolation
Proof. If f is in )(220, then both f and Pf are solutions, in 3(220, to the interpolation problem. Since the solution, in X220, is unique, it follows that Pf = f. To show 2, consider the interpolants over two adjacent tetrahedra. The restrictions of these polynomials to the common face is a bivariate degree-nine polynomial and is uniquely determined by the following fifty-five linearly independent conditions: three of the form (8), forty-two of the form (9), three of the form (10), six of the form (11) and one of the form (12). Hence, the piecewise interpolant has C O continuity. Since the restrictions to the common face agree, derivatives of all orders and in any direction parallel to the face agree on the face. Note the derivative in the inward direction normal to the face restricted to the face is a bivariate polynomial of degree eight and is uniquely determined by the following forty-five linearly independent conditions: thirty of the form (9), three of the form (10), six of the form (11), one of the form (13) and five of the form (14). Hence, the piecewise interpolant has C 1 continuity as desired. The second-order derivative normal to a face restricted to the face is a bivariate polynomial of degree seven and requires thirty-six linearly independent functionals to be uniquely determined. However, the given data specify only twenty-three conditions on this derivative. So, in general, the piecewise interpolant will not be twice continuously differentiable. []
5. Examples To test the interpolant, P, the unit cube [0,1] 3 w a s tessellated into tetrahedra as shown in Fig. 1. Data as defined in Section 2 were taken from the triogonometric primitive
fa(x, Y, z ) = cos(3.1nx) c o s ( y - 0.5) sin(3.1n(z - 0.5)).
(27)
This primitive was used by Stead [Barnhill & Stead '84] to compare trivariate interpolants. The 0.2 contour of fl is shown in Fig. 2, and the corresponding contour of Pfl is shown in Fig. 3. These two figures are identical within the accuracy of the graphics devices which produced them. This is due to the large amount of data from fl to which Pf interpolates. The contouring was accomplished by subdividing the unit cube into subcubes. The function was evaluated at each corner of each subcube allowing linear interpolants to be contructed along each edge of each subcube. These linear interpolants were then used to define piecewise linear approximations, in the faces of the subcubes, to the true contours.
Fig. 1. Tessellation of the unit cube into tetrahedra.
K.L. Rescorla / C l trivariate interpolation
Fig. 2. The 0.2 contour of fl-
243
Fig. 3. The corresponding contour of Pfl.
The precision of the implementation of P was verified to be )(220 by choosing f2 in X220 with 220 nonzero randomly generated coefficients and observing that Pf2 ( x, y, z) = f2 ( x, y, z) for arbitrary (x, y, z) in the domain. Also, it was verified that the implementation of Pfl is C 1 continuous by MICROSCOPE, a multivariate function analyzer [Alfeld & Harris '84].
6. Conclusions A tetrahedral interpolant, P, with the following properties has been constructed. 1. 2. 3. 4.
The required data are C 4. The smoothness is C 1. The precision set is all polynomials of degree less than or equal to nine. The interpolant is a trivariate degree-nine piecewise polynomial. It is represented in terms of a monomial basis which triangularized the Vandermonde matrix. Hence, the coefficients are given by a simple recursion.
Previous to this work, several tetrahedral interpolants were known. [Alfeld '84a] constructs a C 1 scheme for C 1 data by forming convex combinations of projectors. [Barnhill & Little '84] describe a C 1 scheme defined as a boolean sum of projectors. Due to the nature of boolean sums, C 1 data and some auxiliary data are required. [Alfeld '85b] gives an algorithm for constructing an interpolant on simplices of arbitrary dimension having arbitrary smoothness, a C 1 tetrahedral scheme being a special case. All three of these schemes when discretized, are piecewise rational. However, polynomial intrpolants are desirable for a number of reasons. To begin, they are easily evaluated. Secondly, [Alfeld '85a] describes a technique for generating unknown derivative data via functional minimization. This well-motivated method requires differentiation and integration of the interpolant. These operations are greatly simplified when the interpolant is polynomial rather than rational. Thirdly, the fact that polynomials are easy to integrate makes them useful in the finite element method. Previous to this work, [Alfeld '84b] constructed a C 1 piecewise polynomial interpolant. This is a Clough-Tocher scheme: That is, each tetrahedron is subdivided into four subtetrahedra and a trivariate polynomial is defined on each subtetrahedron. Alfeld's scheme has several
244
K.L. Rescorla /
C1
trivariate interpolation
advantages over P. To begin, it requires only C 2 data versus C 4. It is degree-five rather than degree-nine. Finally, it is represented in terms of the Bernstein basis functions which are particularly efficient to evaluate as demonstrated in [Farin '83]. However, on the side of P, there is simplicity. Subdividing the tetrahedra for the Clough-Tocher scheme complicates certain data structures, and evaluation requires the minor additional task of determining in which subtetrahedron the point of evaluation lies. The coefficients of P are defined by a simple forward substitution which is easy to program, whereas the coefficients for the Bernstein basis functions are given recursively by rather lengthy formulas which must be explicitly written. Once programmed, though, the Clough-Tocher scheme is more efficient to evaluate than P due to its lower degree (and hence, fewer terms) and Bernstein basis functions. One application where P could be used and the Clough-Tocher scheme could not is in the generation of high-order (higher than two) derivative data.
Acknowledgements This research was supported in part by the Department of Energy with Contract No. DE-AC02-82ER-12046 to the University of Utah. Thanks is given to Robert Barnhill for providing the stimulating atmosphere which spawned this problem and to Peter Alfeld for his support.
References Alfeld, P. (1984a), A discrete C 1 interpolant for tetrahedral data, Rocky Mountain J. Math. 14, 5-16. Alfeld, P. (1984b), A trivariate Ctough-Tocher scheme for tetrahedral data, Computer Aided Geometric Design 1, 169-181. Alfeld, P. (1985a), Derivative generation by functional minimization, Computer Aided Geometric Design 2, 281-296. Alfeld, P. (1985b), Multivariate perpendicular interpolation, SIAM J. Numer. Anal. 22 (1), 95-106. Alfeld, P. and Harris, W. (1984), MICROSCOPE: a software system for multivariate analysis, MRC Technical Summary Report #2701. Barnhill, R.E. and Little, F.F. (1984), Three- and four-dimensional surfaces, Rocky Mountain J. Math. 14, 77-102. Barnhill, R.E. and Stead, S.E. (1984), Multistage trivariate surfaces, Rocky Mountain J. Math. 14, 103-118. Clark, I. (1977a), Practical kriging in three dimensions, Computers and Geosciences 3, 173-180. Clark, I. (1977b), SNARK - a four-dimensional trend surface computer program, Computers and Geosciences 3, 283-308. Farin, G. (1983), Smooth interpolation to scattered 3D data, in: R.E. Bamhill and W. Boehm, eds., Surfaces in Computer Aided Geometric Design, North-Holland, Amsterdam. Paun¢e, I. and Johnson, K.R. (1976), Three-dimensional interpolation, in: D.F. Merriam, ed., Recent Advances in Geomathematics, 135-155. Rescorla, K. (1985), Multivariate interpolation, Ph.D. Thesis, Department of Mathematics, University of Utah. Stead, S.E. (1983), Smooth multistage multivariate approximation, Ph.D. Thesis, Brown University. Webb, T. (1980), The reconstruction of climatic sequences from botanical data, J. Interdisciplinary History, X:4, 749-772.