Why does the Jacobian of the Bilinear transform appear in the Piola transform? I'm doing work with finite element methods and I often do mappings from the reference square $[0,1]^2$ to the physical quadrilateral by a bilinear/trilinear mapping (we'll stick with 2D for now).
Suppose this mapping maps our quadrilateral vertices by


*

*$(0,0)$ to $(x_1,y_1)$

*$(1,0)$ to $(x_2,y_2)$

*$(0,1)$ to $(x_3,y_3)$

*$(1,1)$ to $(x_4,y_4)$
I worked on my own and found that correct bilinear transformation is
$$
  \left[\begin{array}{cc}
    x \\ y
  \end{array}\right]
  =
  \left[\begin{array}{cc}
    x_1 \\ y_1
  \end{array}\right]
  +
  \left[\begin{array}{cc}
    x_2-x_1 \\ y_2-y_1
  \end{array}\right]
  \hat{x}
  +
  \left[\begin{array}{cc}
    x_3-x_1 \\ y_3-y_1
  \end{array}\right]
  \hat{y}
  +
  \left[\begin{array}{cc}
    x_4-x_3-x_2+x_1 \\ y_4-y_3-y_2+y_1
  \end{array}\right]
  \hat{x}\hat{y}
$$
However, in most cases, we would like an additional transformation that maps normal vectors on the reference square to normal vectors on the physical quadrilateral. I've read that the correct choice is the Piola transformation which is defined by
$$
  \vec{n} = \frac{1}{\text{det}J(\hat{x},\hat{y})} J(\hat{x},\hat{y})\hat{\vec{n}}
$$
where $J(\hat{x},\hat{y})$ is the jacobian of the bilinear transformation above. Why is this true? What is the proof for this?
I started working on my own by trying to map the vectors $(1,0)^T$ and $(0,1)^T$ to the corresponding normals on the physical cell, but I was already erroneous because I assumed the transform was a constant matrix (which is only valid if we map to rectangle without rotating, but this transform is also automatically diagonal in that case), so I have no clue where this comes from or what the intuition for it should be. Any proofs or links to proofs are greatly appreciated.
 A: When we consider the normal vector we only care what happens locally, and this is exactly what the Jacobian lets us do: It tells us how a transformation behaves locally. So it is just a matter of plugging the right definitions into the transform and seing what we get out of it: 
First let's write down the Jacobian.
$$J = \begin{bmatrix}
(x_2-x_1) + (x_4-x_3-x_2+x_1)\hat y  &
(x_3-x_1) + (x_4-x_3-x_2+x_1)\hat x \\
(y_2-y_1) + (y_4-y_3-y_2+y_1)\hat y  &
(y_3-y_1) + (y_4-y_3-y_2+y_1)\hat x \\\end{bmatrix}$$
For simplicity I will write $p_i = (x_i, y_i)^t$.
Let us consider the first edge between $p_1$ and $p_2$ where $\hat y=0$.
In this case $\hat n=(0,-1)$. So for any point on this edge $\hat n$ is definde by the difference $\hat n = \hat q - \hat p$. between $\hat p = (\hat x, 0)$ and $\hat q = (\hat x, -1)$. To find the normal $n$ in the physical space we can just transform $\hat p, \hat q$ via the transform you've given and find $p,q$ which lets us define $n=q-p$ which will be the transformed normal vector (ignoring the scaling factor). We can use this definition of the transformed normal vector because the bilinear map preserves lines along the coordinate axes.
Let us do this: So $\hat p$ is transformed to 
$$q = p_1 + (p_2-p_1)\hat x.$$
Similarly $\hat q$ is transformed to
$$
\begin{align}
p &= p_1 + (p_2-p_1)\hat x - (p_3-p_1) - (p_4 - p_3 - p_2 + p_1)\hat x \\
\end{align}$$
Therefore 
$$
\begin{align}
n &= q - p \\
&= p_1 + (p_2-p_1)\hat x - (p_3-p_1) - (p_4 - p_3 - p_2 + p_1)\hat x - p_1 - (p_2-p_1) \hat x\\
&= p_1-p_3 + (-p_4+p_3+p_2-p_1)\hat x \\
&= -[(p_3-p_1) + (p_4 - p_3  - p_2 + p_1)\hat x] \\
&= J (0, -1)^t \\
&= J \hat n
\end{align}$$
So we see that we find indeed the transformed normal vector by multiplying it with the Jacobian. (Ignoring the scaling, we can always rescale the vector afterwards.) And I'd expect to see similar results if you repeat these computations for the other three edges. 
