Here's a vectorial approach that will work in $\mathbb R^3$ for any pair of quadrilaterals that are not even co-planar. Let the vertices of each quad be represented by $x-$, $y-$, and $z-$ coordinates such that the position vectors at these points are
$$
\vec{P_{01}} = [x_{01},y_{01},z_{01}] ,\hspace 0.25 cm
\vec{P_{02}} = [x_{02},y_{02},z_{02}] ,\hspace 0.25 cm
\vec{P_{03}} = [x_{03},y_{03},z_{03}] ,\hspace 0.25 cm
\vec{P_{04}} = [x_{04},y_{04},z_{04}]
$$
for the first quadrilateral (index $0$), and likewise for the second quadrilateral (index 1)
$$
\vec{P_{11}} = [x_{11},y_{11},z_{11}] ,\hspace 0.25 cm
\vec{P_{12}} = [x_{12},y_{12},z_{12}] ,\hspace 0.25 cm
\vec{P_{13}} =[x_{13},y_{13},z_{13}] ,\hspace 0.25 cm
\vec{P_{14}} = [x_{14},y_{14},z_{14}]
$$
The position vector at any given point within the first quadrilateral (index $0$) can be represented by the bilinear parametrisation
$$
\vec{P_{0}} =\vec{A_{0}}+\vec{B_{0}}\xi +\vec{C_{0}}\eta+\vec{D_{0}}\xi\eta
$$
where
$$
\vec{A_{0}}=\vec{P_{01}}\\
\vec{B_{0}}=\vec{P_{02}}-\vec{P_{01}}\\
\vec{C_{0}}=\vec{P_{03}}-\vec{P_{01}}\\
\vec{D_{0}}=\vec{P_{04}}-\vec{P_{02}}-\vec{P_{03}}+\vec{P_{01}}\\
$$
The objective is to find the normalised $\xi$ and $\eta$ coordinates corresponding to this point. This is achieved by first defining two vectors $\vec u$ orthogonal to $\vec{D_{0}}$ and $\vec v$ orthogonal to $\vec{C_{0}}$ as follows
$$
\vec u=[\vec{D_{0}} \times \vec{C_{0}}] \times\vec{D_{0}}
$$
and
$$
\vec v=[\vec{D_{0}} \times \vec{C_{0}}] \times\vec{C_{0}}
$$
The inner products of the position vector $\vec{P_{0}}$ with $\vec u$ and $\vec v$ yield the following set of equations for $\xi$ and $\eta$
$$
a=b\xi+c\eta\hspace 2 cm \text{[EQ. 1]}\\
f=g\xi+h\xi\eta\hspace 2 cm \text{[EQ. 2]}
$$
where
$$
a=(\vec{P_{0}}-\vec{A_{0}})\cdot \vec u\\
b=\vec{B_{0}}\cdot \vec u\\
c=\vec{C_{0}}\cdot \vec u
$$
and
$$
f=(\vec{P_{0}}-\vec{A_{0}})\cdot \vec v\\
g=\vec{B_{0}}\cdot \vec v\\
h=\vec{D_{0}}\cdot \vec v
$$
From EQ. 1, we have
$$
\eta=\frac{a-b\xi}{c}\hspace 2 cm \text{[EQ. 3]}
$$
Substituting EQ. 3 in EQ. 2 yields a quadratic in $\xi$ and only the root lying in the interval $[0,1]$ is relevant. Once this root is calculated, $\eta$ can be calculated from Eq. 3.
The corresponding point $\vec{P_{1}}$ in the second quadrilateral (with index 1) can now be evaluated as
$$
\vec{P_{1}} =\vec{A_{1}}+\vec{B_{1}}\xi +\vec{C_{1}}\eta+\vec{D_{1}}\xi\eta
$$
where
$$
\vec{A_{1}}=\vec{P_{11}}\\
\vec{B_{1}}=\vec{P_{12}}-\vec{P_{11}}\\
\vec{C_{1}}=\vec{P_{13}}-\vec{P_{11}}\\
\vec{D_{1}}=\vec{P_{14}}-\vec{P_{12}}-\vec{P_{13}}+\vec{P_{11}}\\
$$
If the first quadrilateral (index $0$) is a rectangle, then $\vec D_0$ (and hence $\vec u$ and $\vec v$) will be the null vector. In this case, set $\vec u=\vec B_0$ and $\vec v=\vec C_0$ whereupon EQs. 1 and 2 will be replaced by a set of simultaneous linear equations in $\xi$ and $\eta$.
Hope this helps.