Is it possible to construct a mapping from a square $S: (\xi, \eta) \in [0, 1] \times [0,1]$ to a quadrilateral $Q$ with vertices $\vec p_1, \vec p_2, \vec p_3, \vec p_4$ that has constant jacobian $$ \det \mathbf J = \operatorname{const}\\ \mathbf J(\xi, \eta) = \begin{pmatrix} x_\xi & x_\eta\\ y_\xi & y_\eta \end{pmatrix}. $$ I'm only interested in mapping that is affine at the edges of the square.
In other words, given $x_1, \dots, x_4, y_1, \dots, y_4$ can the following problem have a solution for some value $J$? $$ x_\xi y_\eta - x_\eta y_\xi = J = \operatorname{const}\\ x\big|_{\xi = 0} = (1 - \eta)x_1 + \eta x_4, \quad y\big|_{\xi = 0} = (1 - \eta)y_1 + \eta y_4\\ x\big|_{\xi = 1} = (1 - \eta)x_2 + \eta x_3, \quad y\big|_{\xi = 1} = (1 - \eta)y_2 + \eta y_3\\ x\big|_{\eta = 0} = (1 - \xi)x_1 + \xi x_2, \quad y\big|_{\eta = 0} = (1 - \xi)y_1 + \xi y_2\\ x\big|_{\eta = 1} = (1 - \xi)x_4 + \xi x_3, \quad y\big|_{\eta = 1} = (1 - \xi)y_4 + \xi y_3. $$
By differentiating w.r.t. $\xi$ and $\eta$ one can obtain a system of two PDEs of second order. This makes me think that the problem may be posed correctly.
One of the simplest mapping from $Q$ to $S$ is a bilinear one, but its jacobian is constant iff $Q$ is a parallelogram.
If you force the map to take linearly interpolated values at the boundary as described, I think the answer is no.
Imagine taking the unit square and pushing the top right corner downwards halfway to the bottom right vertex, while keeping all other vertices the same. If the map is differentiable and you zoom in far enough, at any location you look the map must locally look like it maps squares to parallelograms. (Of course, since the map is nonlinear the shape and size of these parallelograms may change as you move around)
Now, if you consider a tiny square in the bottom left corner of the domain of size $dx$-by-$dy$, in order to satisfy the boundary conditions, the parallelogram it gets mapped to must also be a square of size $dx$-by-$dy$ (to first order). On the other hand, if you consider a tiny square in the bottom right, in order to satisfy the boundary conditions it must get mapped to a rectangle of size $dx$-by-$dy/2$, which has half the area.
So, as stated the goal is not possible to achieve in general. However, this sort of counterexample relies heavily on the fixed nature of the mapping of the boundary. Perhaps if the exact natures of the boundary-to-boundary map are left unspecified, one may be able to construct a volume preserving map. (e.g., suppose the in the above example the spacing between field lines on the bottom boundary increases as one goes left to right, so that the bottom right rectangle has dimensions $2dx$-by-$dy/2$, or something like that)