CHEMICAL ENGINEERING TRANSACTIONS VOL. 51, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Tichun Wang, Hongyang Zhang, Lei Tian Copyright © 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-43-3; ISSN 2283-9216 Triangular C2 Interpolants by Piecewise Rational Functions Wenbo Yu Jiangxi University of Science and Technology, Nanchang 330013, China ywbo1026@163.com In this paper a C2 interpolation scheme on a triangle is presented. The interpolant assumes given values and derivatives of orders up to 2 at the vertices of the triangle. It is made up of partial interpolants blended with weight functions. Any partial interpolant with respect to a vertex is a piecewise quintic defined on a split of the triangle and interpolates the data at the vertex and on the two sides sharing the vertex, while the corresponding weight function is a rational polynomial in the barycentric coordinates. The resulted interpolant is a piecewise rational polynomials. By using the Bernstein-Bezier representation of polynomials, the interpolant is easy to describe and evaluate. 1. Introduction Local interpolation over triangulations has received widespread discussion as a tool for use in free-form surface design in Computer Aided Geometric Design, finite element computation, and scattered data processing. The general setting for local interpolation of a prescribed smoothness is to setup a model on one single representative triangle. In order to attain global smoothness on the whole triangulation, a successful model must allow the interpolants on any two adjacent triangles satisfying the prescribed smoothness across the common edge. For the sake of computation, the interpolant has to be simple enough and could be a polynomial, a rational polynomial, and even their piecewise versions defined on a split of the model triangle. In practice, the C1 and C2 models are of essential importance. There have been many C1 schemes but few for C 2 ones; see Farin (1986), Farin (1990), Goodman and Said(1991), Peters (1990), and Strang and Fix (1973). Since C2 schemes supersede C1 ones in better approximation and visual effect (esthetic feeling), endeavor has been paid in designing C2 schemes (Alfred, 1984; Wang, 1992). In general, a Cr scheme of polynomial requires vertex data of order 2r, and the degree of the polynomial can reach 4r +1. So for a C2 scheme, the derivatives assigned at vertices must be of orders up to 4 and the degree of the polynomial has to be 9 (Zenisek, 1970). Higher degrees mean more complexity and instability. Ideas to overcome this are to use rational polynomials or to split the macrotriangle into microtriangles on which a C2 spline is to be defined; see Alfred (1984) and Zhan (1996) for instance. The drawback lying in the splitting trick is that the splitting causes more microtriangles. For this, a rational polynomial interpolant can be adopted. Alfeld and Barnhill (1984) developed a transfinite C2 scheme the discrete version of which results in rational interpolant. Liu and Zhu (1995) characterized C2 rational schemes and presented certain C2 discrete triangular interpolants. Similar schemes can be found in Zenisek (1970). The rational approach is able to prevent the triangle T from being split. Nonetheless, this approach exhibits a possible flaw that the degrees of the denominator and numerator could still be very high. For instance, for the scheme with quintic precision in Whelan (1986), the rational interpolant has degrees (9, 4) of the pairs of the denominator and numerator. By combining the rational formation with the splitting technique, in the present paper, we’ll develop a C 2 scheme in a quite natural way. The interpolant is a weighted sum of three partial interpolants, each of which is a C 2 piecewise quintic on a split of the triangle T with respect to a vertex. And the weight functions are rational polynomials of degrees (2, 2). Then the degrees of the resulted interplant are (7, 2). The paper is organized as follows. In §2, as a preliminary, we formulate bivariate quintics in Bernstein-Bezier forms. Then in §3, the C2 triangular interpolation scheme is described and the explicit formulation of the interpolant is displayed in §4. At last, a conclusion is made in §5. DOI: 10.3303/CET1651203 Please cite this article as: Yu W.B., 2016, Triangular c2 interpolants by piecewise rational functions, Chemical Engineering Transactions, 51, 1213-1218 DOI:10.3303/CET1651203 1213 2. A formulation of quintic polynomials Take T = A0A1A2 as a triangle. Let ei = Ai+1Ai+2 denote the opposite edge of Ai and ni the outer normal vector of ei, i = 0,1,2, counting modulo 3. For α = (p, q) with |α|= p+q, denote by Dα the partial differential operator ∂p+q / ∂xp ∂y q . Suppose f C2(T). For a C2 scheme, an interpolant gC2(T) of f must at least satisfies Dα( g－f )(Ai) = 0, 0 ≤ |α| ≤ 2, iZ3. Let Ai R2, i  Z3, form a triangle T = A0A1A2  R2. For A = (x, y)R2, denote by u = (u0, u1, u2) its barycentric coordinates (see Farin (1986) for the detail) with respect to A0, A1, and A2 (or to T ), i.e. A=u0A0+u1A1+u2A2, with , u0 +u1 +u2 = 1. An index λ = (λ0, λ1, λ2), λi Z+, has length |λ| = λ0 + λ1 + λ2. A polynomial p in Pk, the set of all the bivariate polynomials of degrees up to k, can be expressed in Bernstein-Bezier form (short for B-form) p(u) :=p(x, y)=|λ|=k Bk, λ(u), 210 210 210 , !!! ! =)u( λλλλk uuu λλλ k B where bλ is the B-ordinates at index λ with the Bernstein basis function Bk, λ(u). Let e R2 be a vector e=θ0A0 + θ1A1 + θ2A2, with θ0 + θ1 + θ2 = 0, and θ = (θ0, θ1, θ2). Take ε0 = (1, 0, 0), ε1 = (0, 1, 0), and ε2 = (0, 0, 1). Define recursively ;=|| ,=)0( kλbb λλ .|| ,)!( ! 3 )( rkλbu rk k b Zi ελi r λ i    ∑ ∈ Then the r-th derivative of p along direction e is formulated as ∂ r | e      rkλ λrk r λ Bθb rk k p || , )( )u()( )!( ! )u( . Take ei = Ai+1Ai+2 as the vector from Ai+1 to Ai+2 opposite to Ai, i  Z3. Then ni = ei+2  hi ei with hi = (ei •ei+2)/(ei •ei) is an outer normal vector to ei, i  Z3. Let .3/)2( ,3/)2( ,2/)( 21221121   iiiiiiiii AAVAAVAAV Taking Λ= {λ = (λ0, λ1, λ2); |λ| = 5}, Λi = {λ Λ; λi ≥ 3}, iZi ΛΛ 3V   and cΛ = {(1, 2, 2), (2, 1, 2), (2, 2, 1)}, and ,≠,,∈,,),(∂∂),(∂),( ikjZkjiApFGApFApE iAAAAiijkiAAijii kijiji 3 === we have (Farin, 1990) Lemma 2.1 Let p be a quintic of the form   5|| ,5 )()(   uu Bbp . Define     )()( ,5 uu Bbq . Then    ijkZkji ijkijk jiZji ijij Zi ii GFEq ,,,,,, 333 )()()()( uuuu  where ),61510()u( 23 iiii uuuΦ  , ),34()u( 3 jiuuuΦ ijiij  , , 2 1 )u( 23 jiuuΦ jiijj  }.2,1,0{},,{ ,)u( 3  kjiuuuΦ kjiijk This lemma implies that the B-coordinates of p with indices Λv are determined by the data assigned at the vertices of T. Further, given 3 normal derivatives on one side of T, the B-ordinates with Λc can be determined. Precisely, we have, for instance. Lemma 2.2 On the basis of Lemma 2.1, b122, b212 and b221 are determined in addition by )( 0n0 Vp , ),( 01 2 n0 Vp and )( 02 2 n0 Vp , )( 12 8 0n122 0 Vpb  )1(44([ 6 1 0104113131140 hbbbb  )464( 014023032041050 bbbbb  )]464( 0050140230320410 bbbbbh  , H )]()(2[ 40 3 02 2 n01 2 n221 00 VpVpb  )2( 021 )2( 003 )2( 030 3 1 6 5 bbb  , )]()(2[ 40 3 01 2 n02 2 n212 00 VpVpb  )2( 012 )2( 030 )2( 003 3 1 6 5 bbb  . 1214 Lemma 2.3 With the assumptions that p 0n  and p 2 n0  on edge e0 reduce to polynomials of degrees 3 and 1, respectively, the formulas in Lemma 2.2 are replaced with )44[( 6 1 104113131140122 bbbbb  )464)(1( 0140230320410500 bbbbbh  )],464( 0050140230320410 bbbbbh  ,3/)2( )2(021 )2( 003 )2( 030221 bbbb  .3/)2( )2( 012 )2( 030 )2( 003212 bbbb  In the sequel, we’ll use c , Λλbiλ  , to specify the B-ordinates λb obtained with respect to the data on the edge ei, i ∈ Z3. 3. The C 2 interpolation scheme In this section, we’re presenting a C 2 interpolation scheme. Given a triangulation τ of a polygonal region 2RΩ  , and integers k and r, k ≥ r ≥ 0, the spline space of piecewise Cr polynomials of degree k is defined as (Wang, 2001) σPsΩCsΩτSτS kσ rr k r k ,|);({:),()(  being any triangle in τ }. 3.1 Partial interpolants For the triangle T = A0A1A2, take .2 ,1 , ,10 ,)1( 32121   mZittAtAtB iiiimiimim For any i  Z3, by connecting Ai and Bim, m = 1, 2, we form a split iΔ of T with subtriangles Ti1 = AiAi+1Bi1, Ti2 = AiBi1Bi2, and Ti3 = AiBi2Ai+2. See Figure 1. By using cofactor conforming approach developed by Wang(2001), it’s easy to verify that dim )(25 iΔS 33. Figure 1: The partial split Δ0 Figure 2: The stencil of {bλ ; l , l = 1, 2, 3} Now consider the following partial interpolation problem for i  Z3. Interpolation Problem (P) For f  C2(T), find )(25 ii ΔSs  such that ; ,2||0 ,0))(( 3ZjαAfsD ji α  (1) 3 ,2 ,1 , , ,0))(( ,0))(( 3 2 nn  mZlilVfsVfs lmili ll (2) Suppose si is an interpolant of (P), denote sil = si|Til, for i  Z3 and i = 1, 2, 3. Note that the number of conditions given in problem (P) is 27, which is 6 less than the dimension of the underlying interpolation space )(25 iΔS . So, for the interpolant of (P), if it exists, to be uniquely determined, additional restrictions must be imposed on the interpolation. For the sake of simplicity, without loss of generality, the discussion in this subsection is only on Δ0. Suppose s0 is an interpolant of (P) and    5|| ,5;0 )u( λ lλlλl Bbs (3) where ul = ( u0;l, u1;l, u2;l ) are the barycentric coordinates of point (x, y)  R2 with respect to T0l. Note that for any A = (x, y)  R2, , 23;2023;103;0022;2012;102;0011;211;101;0221100 AuBuAuBuBuAuBuAuAuAuAuAuA  and then u0 = u0;1 = u0;2 = u0;3. Base on this observation, each s0l can be written as 1215             5 0 ;2;1; 5 00 ),( 5 )u( k lllk k l uugu k s where gk; l ( u1; l, u2; l ) is a (univariate) polymomial of degree k, which can be seen as the k-th layer of s0l. This allows s0 another formulation             5 0 21 5 00 ),( 5 )u( k k k uugu k s where gk is a univariate piecewise polynomials of degree k, k = 0, 1, ..., 5, with successive pieces gk; 1, gk; 2, and gk; 3, defined on the partition 0 < t01 < t02 < 1. The supports of the B-ordinates are depicted in Figure 2. From Lemma 2.1, the B-ordinates marked with  are determined by conditions (1), and from Lemma 2.2, the conditions (2) in addition determine those marked with * in Figure 2. In particular, the gk with k = 0, 1, and 2, are all overall polynomials determined by conditions (2) and it is easy to see that, with C2 relations, g3 is uniquely solved by the previously determined B-ordinates. Refer to Schumaker (1981) for univariate splines. For g4 and g5, we make the following convention that (i) g5 reduces to an overall quintic; (ii) g4 is under C3 relations and some two successive pieces of g4 reduce to a quartic. Under restriction (i), g5 is obviously determined by conditions (1). The rule (ii) means that g4 is a C3 cubics of two pieces with 6 parameters, just the same as the number of * and  on g4 . Then we conclude that Theorem 3.1 Under the restrictions (i) and (ii), there uniquely exists an interpolant )( 0 2 5 ΔSs  satisfying the conditions (1) and (2). 3.2 The macro interpolants Similar to s0, the partial interpolants )( 1 2 51 ΔSs  and )( 2 2 52 ΔSs  of (P) are obtained uniquely under (i) and (ii). From the discussion above, we can also see that the partial interpolants 3, Zisi  , hold si | ei = si+1 | ei = si+2 | ei , |1n  i r s i ei = |2n  i r s i ei , r = 1, 2. Define weight functions 3 2 2 2 1 2 0 2 ),/()u( Ziuuuuw ii  . Obviously, these functions have the properties that  Cwi , 1u)(3 Zi iw , ,0|)u( ieiw and 3n ,0|)u( Ziw ii ei  . )u(01s )u(02s )u(03s Figure 3: The stencil of s0 in the uniform BCS u. Now it’s ready to the main result of this paper. Theorem 3.2 Suppose 3 2 5 ),( ZiΔSs ii  , are interpolants of (P) for f C 2(T) under the restrictions of (i) and (ii). Then the blending function    3 )u()u()u( Zi iiT swf (4) satisfies that 3 ,2 || 0 ,0))(( ZiαAffD iT α  (5) 0))((n  iT Vffi , 2 ,1 , ,0))(( 3 2 n  mZiVff imTi . Still, 3n ,e| Zif iT r i  , are polynomials of degrees 5 − r, r = 0, 1, 2. In addition, fT reproduces polynomials of degree 5.       1216 Theorem 3.3 Suppose 3 2 5 ),( ZiΔSs ii  , satisfies ; ,2||0 ,0))(( 3ZjαAfsD ji α  lisl e|n and lisl e| 2 n are of degrees 3 and 1, respectively, 3 , Zlil  . In addition, suppose 3 , Zisi  , are under the restrictions of (i) and (ii), and fT is defined as in (4). Then fT reproduces polynomials of degree 3 and satisfies (5). On ie , Tfin  and ,2n Tfi , 3Zi  are polynomials of degrees 3 and 1, respectively. The usage of rational functions can be dated back to the well-known BBG scheme (Barnhill et al, 1973) and Nielson’s side-vertex scheme (Nielson, 1979). The interpolation schemes in Liu and Zhu (1995) and Xu, et al (2000) can also be regarded as side-vertex schemes, since a partial interpolant on a triangle in those schemes interpolates data at a vertex and on its opposite side. In our method, a partial interpolant interpolates data on two sides sharing a vertex and then our interpolation can be called a side-side scheme. The weight functions used in side-vertex schemes and in side-side schemes are different. For instance, the weight functions used in Liu and Zhu (1995) are 3 2 3 2 2 2 2 2 1 2 1 2 0 2 2 2 1 ),/()u( Ziuuuuuuuuw iii   . Obviously, the macro interpolant as a whole is a piecewise rational function with degrees (7, 2) for the pairs of the denominator and numerator. By contrast, for the scheme with quintic precision in Liu and Zhu (1995), the rational interpolant has degrees (9, 4). 4. Evaluation of the interpolant Suppose f  C2(T). Let q(u) =   V )u(,5Λλ λλBb be the polynomial in Lemma 2.1, with the symbol p replaced with f. Then q is nothing but the polynomial interpolating function f at the vertices of T such that 3 ,2 || 0 ,0))(( ZiαAfqD i α  . Now fix one i  Z3. By assigning 2 ,1 ,0))(( ,0))(( 2nn  mVfpfVfp imi ii , the B-ordinates C , Λλb i λ  , can be determined via Lemma 2.2. Continuing the preceding, we will give an expression of the partial interpolant 3 2 5 ),( ZiΔSs ii  , and hence the macro interpolant fT . In order to express si, owing to the symmetry, we only illustrate the formulation of s0, and the others can be obtained by rotation with the index i. Note that a B-form polynomial in a barycentric coordinates system (BCS) can be expressed in any other BCS with an affine transform between the two systems. For the interpolant ),( 0 2 50 ΔSs  we prefer to write 3 ,2 ,1 ),,( 5 )u(|)u( 21; 5 0 5 0,5 5|| ;000              luugu k BbTss lk k k λ λ lλll (6) in the uniform BCS u with T rather than that with lT0 as in (3). This form has the advantages that , , , 22211;212 2 1221;221 2 2121;122 bbbbbb  ;1 , , 2 1223;212 1 2123;221 1 2213;122 bbbbbb  ,; λlλ bb  for ;3 ,2 ,1 ,2 ,10  lλ ,1; λλ bb  for ;1Λλ ,3; λλ bb  for .2Λλ Since the segments A0B01 and A0B02 are on the lines ,)1(: : 20110111 ututγΓ  ,)1(: : 20210212 ututγΓ  respectively, we can write ),(5),(10)u()u( 2140213 2 0 2 ,1 ,5|| ,50 0 uuguuuguBbs λλ λλ    (7) where g3 and g4 are piecewise C2 cubics and quartics in u1 and u2, and          ,0 , ,0 , ,0 , 233 2132 131 3 γg γγg γg g          ,0 , ,0 , ,0 , 243 2142 141 4 γg γγg γg g 4.1 The representation of g3 The unknown B-ordinates for g3 include b203;1, b230;2, b221;2, b212;2, b203;2, and b230;3. We only provide the formula of b203;1 and the other can be derived similarly. 1217 The C2 smoothness of g3 implies that there are two numbers c1 and c2, such that ),,()],([),( 2132 3 21112131 uuguuγcuug  ).,()],([),( 2133 3 21222132 uuguuγcuug  Identifying the coefficients of u1 and u2 in these equations, we get ))(1( )()1)(( 01020101 021;2123;212021;2213;221 1 tttt tbbtbb c    , ))(1( )()1)(( 02010202 011;2123;212011;2213;221 2 tttt tbbtbb c    .)1()1( 3022 3 0113;2031;203  tctcbb 4.2 The representation of g4 Under (ii), g4 is of C3, and two pieces of g4, e.g. g41 and g42 are identical, i.e. ),,(),( 21422141 uuguug  ),,()],([),( 2143 4 21222142 uuguuγduug  where d2 is the C3 cofactor of g4 on 2Γ . Then 2 02 2 02 1;1223;122 2 )1( tt bb d    , ,)1( 3020223;1131;113 ttdbb  .)1( 4 0223;1041;104 tdbb  The other B-ordinates b140;2, b131;2, b122;2, b113;2, b104;2, b140;3, and b131;3 are formulated similarly and omitted here. 5. Conclusion We have shown in this paper a C2 local interpolation method and the detailed representation of the interpolant. In literature, Alfeld (1984) developed a C2 scheme on the so called twice HCT split for only C2 data. Whelan (1986) displayed a C2 interpolant of piecewise nonics with C4 data, which indeed a special case of the schemes in Zenisek(1970). Both are comparatively complicated because of the large number of split pieces of the macrotriangle of the former and the high degree of polynomials and high order of derivatives of the latter. We can see that the scheme in this paper overcomes these drawbacks to a large extent. With the explicit formulation of the interpolant with Bernstern-Bezier technique, it’s also easy to implement the evaluation of this interpolation scheme. Acknowledgments This paper was supported by Science and Technology Research Project of Jiangxi Provincial Education Department (No. GJJ150696). References Alfeld P., 1984, A bivariate C2 Clough-Tocher scheme, Computer Aided Geometric Design 1, 257-267. Alfeld P., Barnhill R. E., 1984, A transfinite C2 interpolant over triangles, Rocky Mountain J. Math. 14, 17-39. Barnhill R. E., Birkhoff G., Gordon W. J., 1973, Smooth interpolation in triangles, J. Approx. Theory 8, 114-128. Farin G., 1986, Triangular Berstein-Bezier patches, Computer Aided Geometric Design 3, 83-127. Farin G., 1990, Curves and surfaces for computer aided geometric design, Academic Press. Goodman, T.N.T., Said H. B., 1991, A C1 triangular interpolant suitable for scattered data interpolation, Communcation in Applied Numerical Methods 7, 479-485. Liu X., Zhu Y., 1995, A characterization of certain C2 discrete triangular interpolants, CAGD 12, 329-348. Nielson G. M., 1979, The side-vertex method for interpolation in triangles, J. Approx. Theory, 25,318–336. Peters J., 1990, Local smooth surface interpolation: a classification, CAGD 7,191-195. Schumaker L., 1981, Spline Functions: Basic Theory, Wiley, New York. Strang G., Fix G., 1973, An Analysis of the Finite Element Method, Prentice-Hall, New York. Wang R., 2001, Multivariate Spline Functions and Their Applications, Kluwer Academic Publishers. Wang T., 1992, A C2-quintic spline interpolation scheme on triangulation, CAGD 9, 379-386. Whelan T., 1986, A representation of a C2 interpolant over triangles, CAGD 3, 53–66. Xu G., Chu C., Xue W., 2000, Triangular Cm interpolation by rational functions, Approximation Theory and its Applications 16, no. 4, 35–54. Zenisek A., 1970, Interpolation polynomials on the triangle, Numer, Math. 15, 383-296. Zhan Y., 1996, A C2 interpolation scheme over triangulations, SPIE 2644, 232–237. 1218