Transfinite Interpolation - Algebraic way of transposing computational domain to physical domain

Transfinite Interpolation was first described by William Gordon in 1973, it's used for transcribing the computational domain to the physical domain. I've been reading Handbook of Grid Generation by J.Thompson. Transfinite interpolation can be really confusing to people who are new to CFD and Grid work.

There is limited resources online on how Transfinite Interpolation really works. The goal of this post is to explain TFI in detail:

General Concept

From J.Thompson. Shows computational domain transformed into the physical domain

TFI comes from the Boolean Sum Formulation - J.Thompson Handbook of Grid Generation
\begin{align}
U(\xi,\eta,\zeta) = \Sigma_{i=1}^L \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \frac{\partial^n X(\xi_i, \eta, \zeta)}{\partial \xi^n}\\
V(\xi,\eta,\zeta) = \Sigma_{j=1}^M \Sigma_{m=0}^Q \beta_{j}^{m}(\eta) \frac{\partial^m X(\xi, \eta_j, \zeta)}{\partial \eta^m}\\
W(\xi,\eta,\zeta) = \Sigma_{k=1}^N \Sigma_{l=0}^R \gamma_{k}^{l}(\zeta) \frac{\partial^l X(\xi, \eta, \zeta_k)}{\partial \zeta^l}\\
\end{align}

$$X(\xi,\eta,\zeta)  =
\begin{bmatrix}
x(\xi,\eta,\zeta) \\
y(\xi,\eta,\zeta) \\
z(\xi,\eta,\zeta) \\
\end{bmatrix}$$

Tensor Products - not defined correctly in Thompson, but defined correctly here.
\begin{align}
UV = VU = \Sigma_{i=1}^L \Sigma_{j=1}^M \Sigma_{m=0}^Q \Sigma_{n=0}^R \alpha_{i}^{n}(\xi) \beta_{j}^{m}(\eta) \frac{\partial^{nm} X(\xi_i,\eta_j,\zeta)}{\partial \eta^{NM} \partial \xi^n}\\
UW = WU = \Sigma_{i=1}^L \Sigma_{k=1}^N \Sigma_{l=0}^R \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \gamma_{k}^{l}(\zeta) \frac{\partial^{ln} X(\xi_i,\eta,\zeta_k)}{\partial \zeta^{l} \partial \xi^{n}}\\
VW = WV = \Sigma_{j=1}^M \Sigma_{k=1}^N \Sigma_{m=0}^Q \Sigma_{l=0}^R \beta_{j}^{m}(\eta) \gamma_{k}^k (\zeta) \frac{\partial^{lm} X(\xi,\eta_j,\zeta_k)}{\partial \zeta^l \partial \eta^m}
UVW = \Sigma_{i=1}^L \ Sigma_{j=1}^M \Sigma_{k=1}^N \Sigma_{l=0}^R \ Sigma_{m=0}^Q \Sigma_{n=0}^P \alpha_{i}^{n}(\xi) \beta_{j}^{m}(\eta) \gamma_{k}^{l}(\zeta) \frac{\partial^{lmn} X(\xi_i,\eta_j,\zeta_k) }{\partial \zeta^{l} \partial \eta^m \xi^{n}}
\end{align}

$$ X(\xi_i,\eta_j,\zeta_k)  = U+V+W - UV-UW-VW-UVW$$

Looks confusing, lots of variables. The coefficients $\alpha$, $\beta$, and $\gamma$ depend on what method we use, they are called blending functions. Thompson goes into detail with them however, for the purpose of this post we will consider the simple case, linear interpolation.

Linear Interpolation the blending functions are simple, P=Q=R=0, L=M=N=2.
We are going to consider the 2D case
\begin{align}
\alpha_{1}^0 = 1-\xi \\
\alpha_{2}^0 = \xi \\
\beta_{1}^0 = 1-\eta \\
\beta_{2}^0 = \eta\\
U(\xi_i,\eta_j) = (1-\xi_i) X(0,\eta_j) + \xi_i X(1,\eta_j) \\
V(\xi_i,\eta_j) = (1-\eta_j)X(\xi_i,0) + \eta_j X(\xi_i,1)\\
UV(\xi_i,\eta_j) = (1-\xi_i)(1-eta_j) X(0,0) + \xi_i (1-\eta_j) X(1,0) + (1-\xi_i)(\eta_j) X(0,1) + \xi_i \eta_j X(1,1)
\end{align}

How is it used?
Take the simple duct example below. In order for us to map a computational domain onto the duct we have to convert this geometry from a function of x(x,y) y(x,y) to $x(\xi,\eta)$ and $y(\xi,\eta)$. Then we can use the above equations to map the computational domain.

In my code I do a simple mapping, nothing complicated. There may be a few bugs if you change geometries. $\xi$ and $\eta$ are passed into to a function and that outputs (x,y) coordinates. So where you see $X(\xi,\eta)$ the output is a vector $(x,y)$ then that vector is multiplied by $(1-\xi)$ or $\xi$ to get U,V, UV.
$$X(\xi,\eta) = U+V-UV$$

Grid mapped onto duct

Link to Matlab Code: GitHub