There are eight nodes in the Treil maxillo-facial frame [4, 5]. Only one shape function will be evaluated in the (x,y,z) coordinates system following a suitable method, such as Taig's isoparametric formulation [6] or Wachspress' rational polynomial interpolation [7, 8]. (Mathematica coding of these two popular and effective procedures are included in the appendix.) This initial shape function in terms of (x,y,z) will be indicated by the special symbol s[0], and the associated coordinates will be denoted by (x0,y0,z0). Here the string "0" reflects the starting choice for the proposed algorithm. For the eight-node maxillo-facial frame the remaining seven interpolants, Array[s,7], will be calculated to reproduce the set of exact fields denoted by f as a list of functions in (x,y,z). In particular, here: f={1,x,y,z,x*y,y*z,z*x}. Its algebraic relations are crucial to this presentation, as explained below.

Consider a pre-assigned field u[x,y,z] to be an element of the list f. With the prescribed nodal values u[x[i],y[i],z[i]], the set of shape functions s[i] in terms of (x,y,z) should lead to

The aforementioned relation will be exactly satisfied in terms of the coordinate variables (x,y,z). Hence the name "exact interpolated fields" for f.

The following code generates the coefficient matrix mat and the right-hand side vec for LinearSolve:

The anatomical points are symmetrically located on the left and right side of the face. The canonical maxillo-facial frame is then symmetric about x-axis. Hence the following is a specialization for which a symmetric frame is generated according to the rule.

Now the coefficient matrix and the right-hand side vector, matr, vecr, respectively, become:

Before attempting to solve for the unknown shape functions, the consistency of the assumption of the interpolated field should be established.

In general, the nonzero determinant signifies that the assumption:

could be acceptable and the matrix matr is indeed invertible.

For given numerical data for the nodal coordinates, the shape functions will be obtained from the expression below, after removing the comments (* and *).

Hence the following evaluation is not attempted since the coefficient matrix is not invertible.