Mapping from a square $\left[-\frac{1}{2},\frac{1}{2}\right]\times\left[-\frac{1}{2},\frac{1}{2}\right]$ with local coordinate system $\,(\xi,\eta)\,$ to an arbitrary quadrilateral with global coordinates $\,(x,y)\,$ and bi-linear interpolation can be visualized as: The bi-linear transformation is an isoparametric one. This means that the global coordinates $\,(x,y)\,$ transform like any other function $T$ at the quadrilateral: $$\begin{cases} T = N_1 T_1 + N_2 T_2 + N_3 T_3 + N_4 T_4 \\ x = N_1 x_1 + N_2 x_2 + N_3 x_3 + N_4 x_4 \\ y = N_1 y_1 + N_2 y_2 + N_3 y_3 + N_4 y_4 \end{cases}$$ Where: $$\begin{cases} N_1 = (\frac{1}{2}-\xi)(\frac{1}{2}-\eta)\\ N_2 = (\frac{1}{2}+\xi)(\frac{1}{2}-\eta)\\ N_3 = (\frac{1}{2}-\xi)(\frac{1}{2}+\eta)\\ N_4 = (\frac{1}{2}+\xi)(\frac{1}{2}+\eta)\\ \end{cases}$$ From an answer by Steven Stadnicki the conditions on convexity (and non-self-intersection, which is just a special-case) are required to guarantee that the mapping is $1:1$ .
More in general, the Inverse Function Theorem gives sufficient conditions for a function to be invertible in a neighborhood of a point in its domain. The theorem can be generalized to any continuously differentiable, vector-valued function whose Jacobian determinant is nonzero at a point in its domain. For the bi-linear mapping at hand the Jacobian determinant can be easily calculated: $$\Delta = \frac{dx}{d\xi}\frac{dy}{d\eta} - \frac{dx}{d\eta}\frac{dy}{d\xi} > 0$$ Where: $$\frac{dx}{d\xi} = \left(\frac{1}{2}-\eta\right)(x_2-x_1) + \left(\frac{1}{2}+\eta\right)(x_4-x_3) \\ \frac{dy}{d\xi} = \left(\frac{1}{2}-\eta\right)(y_2-y_1) + \left(\frac{1}{2}+\eta\right)(y_4-y_3) \\ \frac{dx}{d\eta} = \left(\frac{1}{2}-\xi\right)(x_3-x_1) + \left(\frac{1}{2}+\xi\right)(x_4-x_2) \\ \frac{dy}{d\eta} = \left(\frac{1}{2}-\xi\right)(y_3-y_1) + \left(\frac{1}{2}+\xi\right)(y_4-y_2)$$ Substitution gives: $$\Delta = N_1 \Delta_1 + N_2 \Delta_2 + N_3 \Delta_3 + N_4 \Delta_4$$ Where: $$\Delta_1 = (x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)\\ \Delta_2 = (x_2 - x_1)(y_4 - y_2) - (x_4 - x_2)(y_2 - y_1)\\ \Delta_3 = (x_4 - x_3)(y_3 - y_1) - (x_3 - x_1)(y_4 - y_3)\\ \Delta_4 = (x_4 - x_3)(y_4 - y_2) - (x_4 - x_2)(y_4 - y_3)$$ Remember the triangle result: $$\Delta = (x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)$$ Now replace $(1,2,3)$ by $(2,4,1)$ , $(3,1,4)$ , $(4,3,2)$ and respectively attach labels $1$ , $2$ , $3$ and $4$ to these triangles. Then the Jacobians $\Delta_k$ at a quadrilateral are identical to the Jacobians $\Delta_k$ at the triangles labeled $k = 1,2,3,4$ , according to the above. Therefore the Jacobian is expressed by the shape functions into its values $\Delta_k$ at (the triangles beloning to) the nodal points of the quad, as it should be.

Numerical experiments have been performed to see where the Jacobian determinant is $\color{blue}{\mbox{zero}}$, $\color{green}{\mbox{positive}}$ or $\color{red}{\mbox{negative}}$. It turns out that there are also places where the Jacobian is positive and negative, i.e. has more than one value (!) These places have been colored $\color{gray}{\mbox{gray}}$ . Five different cases are distinguished: (1) convex area negative, (2) convex area positive, (3) concave area negative, (4) concave area positive, (5) self-intersecting with positive and negative area: The local coordinates are such that $\;-1/2<\xi<+1/2,-1/2<\eta<+1/2\;$ and thus extend beyond the quadrilateral for the three degenerated cases (3), (4) and (5) : this is the $\color{gray}{\mbox{gray}}$ area. Careful inspection leads to the following

Conjecture. The Jacobian determinant is $\color{blue}{\mbox{zero}}$ at the curved boundary of the $\color{gray}{\mbox{gray}}$ area.

Two issues:

• Can we prove this conjecture?
• And what is the shape of that curved boundary?

• UPDATE. Look up @ Wikipedia: Jacobian determinant. A nonlinear map $f: \mathbb{R}^2\to\mathbb{R}^2$ sends a small square to a distorted parallelogram close to the image of the square under the best linear approximation of $f$ near the point. That's what the caption says. $$\vec{r}_1 = \vec{r}(\xi+d\xi,\eta)-\vec{r}(\xi,\eta) = \begin{bmatrix} \large \frac{\partial x}{\partial \xi} \\ \large \frac{\partial y}{\partial \xi} \end{bmatrix} d\xi \\ \vec{r}_2 = \vec{r}(\xi,\eta+d\eta)-\vec{r}(\xi,\eta) = \begin{bmatrix} \large \frac{\partial x}{\partial \eta} \\ \large \frac{\partial y}{\partial \eta} \end{bmatrix} d\eta$$ The area of the distorted parallelogram spanned by the vectors $\vec{r}_1$ and $\vec{r}_2$ is the determinant which is well known from multivariable integration: $$\begin{vmatrix} \large \frac{\partial x}{\partial \xi} & \large \frac{\partial x}{\partial \eta} \\ \large \frac{\partial y}{\partial \xi} & \large \frac{\partial y}{\partial \eta} \end{vmatrix} d\xi\,d\eta = \left(\frac{\partial x}{\partial \xi}\frac{\partial y}{\partial \eta} -\frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \xi}\right) d\xi\,d\eta$$ More pictures, with Jacobian determinant $\color{red}{\mbox{negative}}$ and $\color{green}{\mbox{positive}}$ : It is seen that the Jacobian determinant always switches sign at the parabolic boundary, which is inevitable because one of the answers below shows that there is only one such boundary, where the determinant is zero. The only way to change sign is to change orientation, because continuity must be ensured, i.e. the quads must remain connected with their neighbor infinitesimal quadrilaterals. This may not be a die hard proof, but it makes me to believe that the Conjecture is indeed true.

• Why was this downvoted? it is a perfectly good question – Rab Sep 15 '17 at 17:27
• @RabMakh: Thank you for this kind remark. – Han de Bruijn Sep 15 '17 at 18:02

A possible clue may be that the interpolation of the Jacobian is not really bi-linear. To see this, take a look at its Finite Difference representation, which is a different way of telling the same as with the Finite Element representation mentioned in the question: $$\Delta = N_1\Delta_1+N_2\Delta_2+N_3\Delta_3+N_4\Delta_4 = A_\Delta + B_\Delta\xi + C_\Delta\eta + D_\Delta\xi\eta$$ Where: $$A_\Delta = \Delta(0,0) = \frac{1}{4}\left(\Delta_1+\Delta_2+\Delta_3+\Delta_4\right)\\ B_\Delta = \frac{\partial \Delta}{\partial \xi}(0,0) = \frac{1}{2}\left(-\Delta_1+\Delta_2-\Delta_3+\Delta_4\right)\\ C_\Delta = \frac{\partial \Delta}{\partial \eta}(0,0) = \frac{1}{2}\left(-\Delta_1-\Delta_2+\Delta_3+\Delta_4\right)\\ D_\Delta = \frac{\partial^2 \Delta}{\partial \xi \partial \eta}(0,0) = \Delta_1-\Delta_2-\Delta_3+\Delta_4$$ The last term $D_\Delta$ is zero. This is seen by considering the fact that the area of the quadrilateral is equal to the sum of the areas of the triangles (1) and (4), but also to the sum of the areas of the triangles (2) and (3): Therefore $\Delta_1/2 + \Delta_4/2 = \Delta_2/2 + \Delta_3/2$ , which proves the above claim.
We conclude that $\Delta$ is not a bi-linear but only a linear expression in the $\Delta_k$'s. At places where the Jacobian determinant is zero we thus have an equation of the form: $$\Delta = A_\Delta + B_\Delta\xi + C_\Delta\eta = 0 \quad \Longrightarrow \quad C_\Delta\eta = -(A_\Delta + B_\Delta\xi)$$ Now have a look at the Finite Difference representation of the global coordinates: $$\begin{cases} x = A_x + B_x\xi + C_x\eta + D_x\xi\eta\\ y = A_x + B_y\xi + C_y\eta + D_y\xi\eta \end{cases}$$ It follows that, in general, $x$ and $y$ are quadratic functions of the parameter $\xi$ : $$\begin{cases} C_\Delta x &=& C_\Delta A_x + C_\Delta B_x\xi - C_x(A_\Delta + B_\Delta\xi) - D_x\xi(A_\Delta + B_\Delta\xi) \\&=& (C_\Delta A_x - C_x A_\Delta) + (C_\Delta B_x - C_x B_\Delta - D_x A_\Delta)\xi - D_x B_\Delta\xi^2\\ C_\Delta y &=& (C_\Delta A_y - C_y A_\Delta) + (C_\Delta B_y - C_y B_\Delta - D_y A_\Delta)\xi - D_y B_\Delta\xi^2 \end{cases}$$ At Mathematics Stack Exchange we have the following reference to prove that the resulting curve must be a parabola:
Apart from special cases eventually, this answeres the second issue of the question. Free software source confirming this numerically is found at my website . The range of the local coordinates (displayed as $\color{gray}{\mbox{gray}}$ pixels) has been enlarged beyond $[-1/2,+1/2]\times[-1/2,+1/2]$ and (part of) the parabolic curve is drawn in $\color{blue}{\mbox{blue}}$ : 