NAT'L INST. OF STAND & TECH R.I.C llll Hill III!Ill WIST A111D5 t.SOt.30 PUBLICATIONS NISTIR 6246 Surface Energy Formula for a Hsieh-Clough-Tocher Element Marjorie McClain Christoph Witzgall U.S. DEPARTMENTOF COMMERCE Technology Administration National institute ofStandards and Technology Gaithersburg, MD 20899 Mathematical and Computational Sciences Division January 1999 NET QC 100 .U56 NO. 6246 1999 NISTIR 6246 Surface Energy Formula for a i Hsieh-Clough-Tocher Element Marjorie McClain Christoph Witzgall U.S. DEPARTMENTOFCOMMERCE Technology Administration National Institute of Standards and Technology Gaithersburg, MD 20899 Mathematical and Computational Sciences Division January 1999 U.S. DEPARTMENT OF COMMERCE William M. Daley, Secretary TECHNOLOGY ADMINISTRATION Gary R. Bachula, Acting Under Secretary for Technology NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY iContribution ofthe National Institute of Standards and Technology (NIST) Raymond G. Kammer, Director andnotsubjecttocopyrightintheUnitedStates.Researchsupportedinpart by: NSFGrant ECS-8709795, ONR Grant N00014-89-J1537, ONR Contract U.S. DEPARTMENT OF ENERGY N-0014-87-F0063. Washington. D.C. 20858 ) - Abstract Given a triangulation of a planar region, the reduced Hsieh-Clough-Tocher (rHCT) triangular element enables the construction of a smooth piecewise-cubic surface. In prepa- ration for the use of energy minimization to select interpolants from that family of surfaces, a formula for the bending energy of the generic rHCT element is determined, treating the element as a thin and almost flat plate. The aim is to find an approach to interpolation which displays low sensitivity to the choice of triangulation. The problem arose in terrain modeling. Keywords: Delaunay triangulation, finite elements, Hsieh-Clough-Tocher element, interpo- lation, spline, surface, terrain modeling, triangulation Introduction 1. Lawson [9],[10] pioneered the use of (^-compatible finite elements for the task of interpolating a continuously differentiable (“C1”) function z — f{x,y) from a finite set of spatial points = (xi,Vi,Zi), i 1 so that Zi = z(xi,yi). An important application is to terrain modeling from large sets of elevation data, where the function values Zi are elevations at specified data locations (xl,yl). For his work, Lawson chose an element generally ascribed to Clough and Tocher but [4], frequently referred to in the literature (e.g. Ciarlet [2]) as the reduced Hsieh-Clough-Tocher (rHCT) element For a description of that and related elements, we recommend the text by Bernadou and Boisserie [1]. Farin [7] formulated Clough-Tocher interpolation in terms of Bernstein-Bezier polynomials and also provided a modification of the rHTC element which trades C2-continuity at the barycenter - where subtriangles with different cubics join (see Sections 2 and 4 below) - for improved approximation to C2-continuity at the edges of the triangular element. For general information, the reader may want to consult Ciarlet and [3] Zienkiewicz [18], to name just two representatives of a large body of literature on the Finite Element Method. The rHCT element is a triangular surface patch, that is, it is defined over a footprint triangle in the x, y-plane. Its use 'in the context of interpolation thus presupposes a “trian = gulation” of the data locations (x^,?/*), i l,...,n, that is, a set of triangles whose interiors do not meet, whose vertices are the data locations, and whose union covers the convex hull of the latter. Two different triangles of the triangulation thus have either a single edge or a single vertex in common, or do not meet at all. In the context of terrain modeling, a set of elevation data locations triangulated in this fashion is now frequently called a “triangulated irregular network (TIN)”. Any set of planar data locations Xi,yi can be triangulated in (. many different ways and this affects the surface to be constructed. A frequently used generic method is the Delaunay triangulation [5]. 1 The rHCT elements - by their construction - join smoothly at triangle boundaries, thus ensuring a continously differentiable fit. The specification ofan rHCT element over a particular triangle requires that the function value Zi and the gradient (zjX, Ziy) be given at each triangle vertex, that is, at each data location (x^yi). The resulting function is to agree with those prescribed values and slopes. Unless such gradient information is also provided, it must be estimated from given function values Zi in order to complete the specification of rHCT elements for each triangle of the TIN. Thus one distinction between alternate methods for surface interpolation based on the rHCT element is what method for estimating gradients at data locations is chosen. Lawson’s approach [10] is to select at least six neighbors of a data point, assign them weights which decrease with distance from the data point, and then consider 6-parameter quadratic bivariate functions which pass through the data point. From among those functions select the one which minimizes the weighted least squares error. The tangent to that quadratic function at the data point provides the gradient estimate. Franke and Stead compare [8] [15] various procedures for triangulation based interpolation, in particular, aspects of gradient estimation. In the context of terrain modeling, with elevation data given along digitized contour lines, Mandel, Bernal, and Witzgall [11] arrived at gradient estimates by minimizing the elastic energy of a mechanical surrogate structure for the surface consisting of thin beams along the edges of the triangular patches and joined tangentially to thin plates, one at each vertex. The orientation of the plates can be adjusted to minimize the elastic energy of the beams, thus defining gradients at triangulation vertices, that is, at data locations. That procedure was chosen because of the distribution of data points along lines, which rendered local fitting procedures less attractive. It also eliminated the somewhat arbitrary choice of neighboring data points and their weights. The two rHCT methods mentioned above are less sensitive to changes in the underlying triangulation than the still most frequently employed linear TIN methods. In those linear methods, each footprint triangle in the TIN gives rise to a planar triangular facet spanned by the elevations at the triangle vertices. The resulting piecewise linear surface is, indeed, extremely sensitive to the choice of triangulation. With the goal of further reducing the sensitivity of the interpolating surface to the choice of triangulation, we propose to examine an alternate method for estimating the gradients at data locations. The idea is to replace the surrogate mechanical structure consisting of thin beams used in [11] by a surrogate mechanical structure consisting ofthe actual rHCT elements, joined together smoothly, and to minimize the total elastic energy ofthe resulting interpolating surface considered as consisting of almost flat thin plates. More precisely, we determine that 2 interpolating surface z — z(x,y) which consists of rHCT elements and minimizes d2z a2 d2 (i.i) + 2 . + . dx dy dx2 dx dy dy'4 with respect to the choice of gradients at triangulation vertices. A critical first step in this direction is to develop closed formulas for that energy operator as applied to the generic rHCT element, and that is the purpose ofthis report. Those formulas are unwieldy, and the availability of symbolic computation packages was instrumental in their derivation. This work relied on Mathematica (see for instance [16]). Powell-Sabin splines [14] also provide a unique piecewise polynomial Ci-function which meets prescribed values and gradients at the vertices of a triangulation. The concept of thin plate energy minimization has also been considered for various other approaches to bivariate interpolation (e.g. Powell [13]). Mansfield [12] minimizes the full strain energy functional. The reader may want to consult Dierckx [6]. We also revisit the definition of the rHCT element and verify its main properties. This exposition as well as the derivation ofthe energy formulas is conducted in terms of barycentric coordinates. In this fashion, the inherent symmetries of the element are preserved. In particular, the roles of the vertices of the element and their associated quantities are interchangeable. More precisely, many formulas involving indexed quantities follow from each other by — — — (1.2) cyclic substitution : 1 2 ^ 3 > 1, where each index value is replaced by its cyclic successor. We will not always write all three instances of formulas that arise from each other by cyclic substitution of indices. Instead, we will write one instance of the formula, and indicate that the remaining instances can be generated by cyclic substitution. rHCT 2. Definition of the Element In this section we define various quantities which are used to calculate the generic element. Suppose we are given three triangle vertices, (£i,2/i), (£2,2/2), (£3,2/3), with associated function values, i.e., elevations, *i, z2, z3, and partial derivatives with respect to x and y, i.e. coordinate slopes, Zlx, ^lyi ^2xi %2i/, xi Z3y. 3 . We wish to define the elevation z — z(x, y at any location (x, y) in the footprint triangle of ) rHCT the element, so that dz dz Zj Z,{Xi, yi), Zix (^Xx5 yi)i Ziy (*^z: Vi)> ^ 1,2,3. For that purpose, it is convenient to express functions over triangles using ubarycentric coordinates” For any point (x, y) in the plane, its barycentric coordinates Ai, A2, A3 are defined by the relations Ai + A2 4- A3 = 1 (2.1) X\X\ + A2X2 + A3X3 = x A12/1 + A22/2 + A3y3 = y . Following Zienkiewicz we use the abbreviations [17], Xij := Xi - xj: yio := yx - yjx % := z{ - Zj, ij = 1, 2,3, i / j. The following determinants will play a role: 1 1 1 (2.2) 'xy II IQa := det Xi x2 23 — 242/23 + ^22/31 + x3Vl2, V\ 2/2 2/3 1 1 1 zy II IQ := det Z\ 22 23 Z\y23 + ^22/31 + 232/12: V\ 2/2 2/3 = -Dxz := det 1 1 1 Z\X^3 + 22^31 + 23X12. 2l 22 23 Xi x2 £3 D The determinant xy is double the area of the footprint triangle. It is assumed that the area of that triangle does not vanish. The footprint triangle is required, moreover, to have positive orientation, that is, D > xy 0. 4 ) } Solving the linear system of equations (2.1) yields (2.3) Xi = {y23X + x32y + x2ys-y2X3)/DXy X2 — {y3ix + x\sy + x^yi — yzX\)IDxy A3 — (ynx + x2iy + x\y2 ~ yix2)/Dxy. Note that the above formulas arise from each other by cyclic substitution (1.2). Indeed, barycentric coordinates treat all vertices of a triangle symmetrically. Their signs indicate whether a location lies inside or outside the footprint triangle: a negative barycentric coor- dinate means the location is outside the triangle; if all three coordinates are positive, the location lies in the interior of the triangle. Edges of the footprint triangle are characterized by the vanishing of one of the barycentric coordinates. At vertex i, A* = 1 and X0 = 0 for j / i. The barycenter or centroid x fx\ + x2 + ^3 Vi + 2/2 + yz ( o, 2/o V 3 5 3 of the triangle is characterized by Ax = A2 = A3 = 1/3. The barycenter plays a key role in the definition of the rHCT element. It is used to define a “barycentric partition” of the triangle into four subsets, the barycenter itself and three subtriangles (see Figure 1): = = = Bo '= {(Ai, A2, A3) : Ai A2 A3 1/3} B\ := {(Ai,A2,A3) : 0 < Ax < A2, Ai < A3} B2 := {(Ai, A2, A3) : 0 < A2 < A3, A2 < Ai B3 {(Ai, A2, A3) : 0 < A3 < Al5 A3 < A2} . = Note that in each of the subtriangles Bi,i 1, 2, 3, the corresponding barycentric coordinate is dominated by the remaining ones; e.g., (2.4) Ax < min{A2, A3} in B\. The rHCT function is defined in terms of three correction functions [10] Pi, P2, P3, which are piecewise cubic polynomials in the barycentric coordinates, each defined in a piece- wise manner with respect to the subtriangles of the barycentric partition. 5 6