Skip navigation links
(Wind-US Documentation Home Page) (Wind-US User's Guide) (GMAN User's Guide) (MADCAP User's Guide) (CFPOST User's Guide) (Wind-US Utilities) (Common File User's Guide) (Wind-US Installation Guide) (Wind-US Developer's Reference) (Guidelines Documents)

(b4wind User's Guide) (Converting Files) (Computing Initial Conditions) (Interpolation) (Computing the Reynolds Number) (Preparing NPARC or XAIR Namelists) (Installation and Execution) (Conclusions) (Computing the Mach Number From the Area Ratio) (Trilinear Interpolation)

Trilinear Interpolation

This section describes the trilinear interpolation algorithm using pseudo-Fortran code.

First, number the eight vertices of a hexahedron.

   x1 = x(i  ,j  ,k  )
   x2 = x(i+1,j  ,k  )
   x3 = x(i  ,j+1,k  )
   x4 = x(i+1,j+1,k  )
   x5 = x(i  ,j  ,k+1)
   x6 = x(i+1,j  ,k+1)
   x7 = x(i  ,j+1,k+1)
   x8 = x(i+1,j+1,k+1)

Defined as below, notice that xa = x1 when a = −1 and that xa = x2 when a = 1 and similar observations for the rest of the equations.

   xa = [(1-a)*x1 + (1+a)*x2] / 2
   xb = [(1-a)*x3 + (1+a)*x4] / 2
   xc = [(1-a)*x5 + (1+a)*x6] / 2
   xd = [(1-a)*x7 + (1+a)*x8] / 2
   xe = [(1-b)*xa + (1+b)*xb] / 2
   xf = [(1-b)*xc + (1+b)*xd] / 2
   x  = [(1-g)*xe + (1+g)*xf] / 2

Take three steps of substitutions.

   x = {(1-g) * [(1-b)*xa + (1+b)*xb] +
        (1+g) * [(1-b)*xc + (1+b)*xd]} / 4

   x = {(1-g) * [(1-b) * ((1-a)*x1 + (1+a)*x2)  +
                 (1+b) * ((1-a)*x3 + (1+a)*x4)] +
        (1+g) * [(1-b) * ((1-a)*x5 + (1+a)*x6)  +
                 (1+b) * ((1-a)*x7 + (1+a)*x8)]} / 8

   x = {(1-g)*(1-b)*(1-a)*x1 + (1-g)*(1-b)*(1+a)*x2 +
        (1-g)*(1+b)*(1-a)*x3 + (1-g)*(1+b)*(1+a)*x4 +
        (1+g)*(1-b)*(1-a)*x5 + (1+g)*(1-b)*(1+a)*x6 +
        (1+g)*(1+b)*(1-a)*x7 + (1+g)*(1+b)*(1+a)*x8} / 8

Regrouping

   x =       (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 +
           a*(x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8 +
           b*(x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8 +
           g*(x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8 +
         a*b*(x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8 +
         a*g*(x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8 +
         b*g*(x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8 +
       a*b*g*(x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Define

   f0 = (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1) / 8 - x
   f1 = (x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1) / 8
   f2 = (x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1) / 8
   f3 = (x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1) / 8
   f4 = (x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1) / 8
   f5 = (x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1) / 8
   f6 = (x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1) / 8
   f7 = (x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1) / 8

Then define

   f = f0     + f1*a   + f2*b   + f3*g     +
       f4*a*b + f5*a*g + f6*b*g + f7*a*b*g

Exactly similar relations can be derived for y and z.

   g = g0     + g1*a   + g2*b   + g3*g     +
       g4*a*b + g5*a*g + g6*b*g + g7*a*b*g

   h = h0     + h1*a   + h2*b   + h3*g     +
       h4*a*b + h5*a*g + h6*b*g + h7*a*b*g
where the g's and h's are defined exactly similarly to the f's. Note that in the three equations above everything is known except for a, b and g. From the expression for x, and the similar equations for y and z, it is seen that
   f = 0
   g = 0
   h = 0

So the three equations can be solved for the unknowns a, b and g. A solution can be obtained from Newton's method by letting

  df = -f
  dg = -g
  dh = -h
where
   df = f1*da     + f2*db     + f3*dg     +
        f4*b*da   + f4*a*db
        f5*g*da   + f5*a*dg
        f6*g*db   + f6*b*dg
        f7*b*g*da + f7*a*g*db + f7*a*b*dg

Regrouping

   df = (f1 + f4*b + f5*g + f7*b*g)*da +
        (f2 + f4*a + f6*g + f7*a*g)*db +
        (f3 + f5*a + f6*b + f7*a*b)*dg
      = -f

Similarly

   dg = (g1 + g4*b + g5*g + g7*b*g)*da +
        (g2 + g4*a + g6*g + g7*a*g)*db +
        (g3 + g5*a + g6*b + g7*a*b)*dg
      = -g

   dh = (h1 + h4*b + h5*g + h7*b*g)*da +
        (h2 + h4*a + h6*g + h7*a*g)*db +
        (h3 + h5*a + h6*b + h7*a*b)*dg
      = -h

Using an initial guess of a = b = g = 0 the three equations can be solved for da, db and dg. Then one has

   anew = aold + da
   bnew = bold + db
   gnew = gold + dg

If the matrix of the coefficients are zero it is because the hexahedron is degenerate. The new a, b and g can be used to compute new coefficients and solve for new deltas and iteration should continue util the deltas are negligible. The interpolating program has a hard-wired limit of 20 iterations. If after 20 iterations or if the absolute values of any of the three parameters exceed 5.0, then the method fails and the starting stencil is returned as the "best" one.

After a, b and g have been computed the interplated value of any variable known at the vertices can be computed by

   v  = [(1-g)*(1-b)*(1-a)*v1 + (1-g)*(1-b)*(1+a)*v2 +
         (1-g)*(1+b)*(1-a)*v3 + (1-g)*(1+b)*(1+a)*v4 +
         (1+g)*(1-b)*(1-a)*v5 + (1+g)*(1-b)*(1+a)*v6 +
         (1+g)*(1+b)*(1-a)*v7 + (1+g)*(1+b)*(1+a)*v8] / 8

This equation was copied from the equation for X and can easily be verified by letting a, b and g take on their extreme values of ±1. If (x,y,z) is inside the hexahedron, then the parameters a, b and g will be between ±1. If one of them is outside those limits then the given point is outside the hexahedron and the above formula becomes an extrapolation, not an interpolation.