31.2 Piecewise Cubic Interpolation

Piecewise cubic interpolation is a technique that can be used to create functional $C^1$ surfaces (1st derivative continuous) throughout a grid domain by using functional and derivative information at each vertex of a grid cell. The domain function $F$ and its derivatives $dF$ are assumed to be known or at least calculable at each vertex. In order to show the essential features of the piecewise cubic interpolation method and to keep the formulas and matrix sizes manageable, we present the method for 2D domain geometries.

Consider a particular rectangular 2D cell. We wish to calculate a domain function $F$ as a piecewise cubic polynomial in terms of $[0,1]$ rescaled variables $x,y$ inside each cell:

$$F(x,y) = \sum^3_{i,j=0} a_{ij}x^iy^j.$$ (31.1)

The $[0,1]$ rescaled variable $x$ is obtained from the domain global variable $X$ and the cell's lower $X_{\ell}$ and upper $X_u$ dimension as

$$x = \frac{X-X_{\ell}}{X_u-X_{\ell}},$$ (31.2)

with a similar expression for $y$. From the expansion in 31.1, there are 16 expansion coefficients which need to be determined (fitted) to an appropriate set of data for specific points on the cell, which, in order to assure continuous representation of $F$ over the entire domain, have to sit on the cell's boundaries. The most obvious choice for these points are the cell vertices, each vertex belonging to four cells. The data needed for each vertex are $F$, $dF/dx$, $dF/dy$ and $d^2F/dxdy$. This leads to a linearly independent and rotationally invariant set of 16 data values, from which the 16 expansion coefficients can be uniquely determined. It can also be shown that for other dimensions only the inclusion of all mixed simple higher order derivatives leads to a linearly independent and rotationally invariant set of data values. For example, for rectangular 3D geometries the data needed at each of the 8 vertices of a cubic cell is $F$, $dF/dx$, $dF/dy$, $dF/dz$, $d^2F/dxdy$, $d^2F/dxdz$, $d^2F/dydyz$ and $d^3F/dxdydz$.

Let us now organize the 16 expansion coefficients $a_{ij}$ into a vector ${\bf a}$, such that each element $a(k)$ of ${\bf a}$ is associated with the following expansion coefficient:

$$a(k) = a_{ij} \;\;\;;\;\;\; k = 1 + i + 4j.$$ (31.3)

Likewise, we stack the data values of the four vertices $v=1,2,3,4$ into a vector ${\bf b}$ as follows:

$${\bf b} = \left\{ \begin{array}{l} F_v \\ (dF/dx)_v \\ (dF/dy)_v \\ ... .../dx)_3 (dF/dx)_4\end{array}\right. \;\;\;;\;\;\;\mbox{etc...}$$ (31.4)

Let us associate the four vertex indices with the following four $[0,1]$ rescaled variables of the cell:

$$\begin{array}{lcl} v = 1 & \longrightarrow & (x,y) = (0,0) \\... ... (0,1) \\ v = 4 & \longrightarrow & (x,y) = (1,1) \end{array}.$$ (31.5)

Substituting these values for $x$ and $y$ into equation 31.1 and its derivatives, we can establish a connection between the expansion coefficient vector and the data vector in the form

$${\bf b} = {\bf Ba},$$ (31.6)

where the 16 x 16 matrix ${\bf B}$ has the following structure

$${\bf B} = \left( \begin{array}{cccccccccccccccc} 1 & 0 & 0 & 0 & 0 & 0 ... ...& 1 & 2 & 3 & 0 & 2 & 4 & 6 & 0 & 3 & 6 & 9 \end{array}\right),$$ (31.7)

containing only positive integers and many zeros (175 out of 256 elements). Since we are ultimately interested in the expansion coefficients for the cell, we invert equation 31.6 and get

$${\bf a} = {\bf B}^{-1}{\bf b}.$$ (31.8)

The inverse ${\bf B}^{-1}$ is also an integer matrix and again contains many zero elements, although not as many as ${\bf B}$ itself (156 out of 256 elements):

$${\bf B}^{-1} = \left( \begin{array}{rrrrrrrrrrrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 ... ...& -2 & -2 & 2 & -2 & 2 & -2 & 1 & 1 & 1 & 1 \end{array}\right).$$ (31.9)

The matrix ${\bf B}^{-1}$ is universal for all 2D cells. Each cell has its specific data vector ${\bf b}$ from which its expansion coefficient vector ${\bf a}$ can be established via 31.8. Rather than storing ${\bf B}^{-1}$ in matrix form and performing a direct matrix times vector operation for each cell, ${\bf B}^{-1}$ is used implicitly when converting ${\bf b}$ to ${\bf a}$. This avoids the redundant zero multiplications and is thus much more efficient. In fact, by using cleverly arranged reusable intermediates, this ${\bf b}$ to ${\bf a}$ transformation can be coded using only additions and subtractions and no multiplications at all.

Once the expansion coefficient vector ${\bf a}$ for a cell has been established, equation 31.1 can be used to obtain the function value $F$ at any point inside the cell, provided the $[0,1]$ rescaled coordinates of that point are given. The function values are given by Eq.(31.1), which can be computed efficiently using a double Horner scheme:

$${\tilde a}_j = ((a_{3j}x+a_{2j})x+a_{1j})x+a_{0j}$$ (31.10)

$$F(x,y) = (({\tilde a}_3y+{\tilde a}_2)y+{\tilde a}_1)y+{\tilde a}_0$$ (31.11)

The evaluation of the global coordinate differentials is done using the chain rule and leads to:

$$\frac{d^{m+n}F(X,Y)}{dX^mdY^n} = (X_u-X_{\ell})^{-m}(Y_u-Y_{\ell})^{-n}\frac{d^{m+n}F(x,y)}{dx^mdy^n}$$ (31.12)

$$\frac{d^{m+n}F(x,y)}{dx^mdy^n} = \sum^3_{i=m}\sum^3_{j=n} (i)_m(j)_na_{ij}x^{i-m}y^{j-n},$$ (31.13)

where the Pochhammer symbol $(i)_m$ is defined as $(i)_m=i(i-1)...(i-m+1)$.