Thompson's equations are set of coupled differential equations used to generate body-fitted meshes to be used in numerical simulations such as computational fluid dynamics (CFD). I have been unable to derive these equations for myself and I could do with some help.
Mesh generation via Thompson's equations
We assume there is an object with a non-trivial boundary situated in Euclidean space $\mathbb{R}^2$ with Cartesian coordinates $(x,y)$. The task of generating a mesh is to find a pair of coordinate transformation functions $\xi(x,y)$ and $\eta(x,y)$ such that the $(\xi , \eta)$ plane is square and the boundary conditions become simple. For example, we could take $\xi$ as a radial-like coordinate such that $\xi = 0$ on the boundary and $\eta$ as a periodic coordinate.
Thompson achieves this by defining the coordinate transformation functions as solutions of two Poisson equations
$$\xi_{xx} + \xi_{yy} = P, \quad \eta_{xx} + \eta_{yy} = Q$$
where $P$ and $Q$ are a pair of functions used to control the coordinate system locally. Now these equations cannot be solved directly because the boundary conditions in the $(x,y)$ plane are non-trivial, so Thompson instead solves for the inverse functions $x(\xi, \eta)$ and $y(\xi, \eta)$ by swapping the independent and dependent variables and arrives at the pair of equations
$$ \begin{aligned} \alpha x_{\xi \xi } - 2 \beta x_{\xi \eta} + \gamma x_{\eta \eta} & = - J^2 (P x_\xi + Q x_\eta) \\ \alpha y_{\xi \xi } - 2 \beta y_{\xi \eta} + \gamma y_{\eta \eta} & = - J^2 (P y_\xi + Q y_\eta) \end{aligned} $$ where
$$ \begin{aligned} \alpha & = x_\eta^2 + y_\eta^2 \\ \beta & = x_\xi x_\eta + y_\xi y_\eta \\ \gamma & = x_\xi^2 + y_\xi^2 \\ J & = x_\xi y_\eta - x_\eta y_\xi \end{aligned} $$
I have tried to derive these equations and cannot. I have looked at many references, such as the original paper, a related paper by the author (page 15), and textbooks such as this (Equations 2 to 4) and this (page 196), and some lecture slides (Page 22) but they do not give the derivation explicitly.
Attempt 1
I can swap the independent and dependent variables with the help of the Jacobian:
$$ \begin{pmatrix} \xi_x & \xi_y \\ \eta_x & \eta_y \end{pmatrix} = \begin{pmatrix} x_\xi & x_\eta \\ y_\xi & y_\eta \end{pmatrix}^{-1} = \frac{1}{J} \begin{pmatrix} y_\eta & -x_\eta \\ -y_\xi & x_\xi \end{pmatrix} $$
This gives me
$$ \begin{aligned} \xi_{xx} & = \partial_x \xi_x \\ & = (\xi_x \partial_\xi + \eta_x \partial_\eta ) \xi_x \\ & = \frac{1}{J} y_\eta \partial_\xi \left(\frac{1}{J} y_\eta \right) - \frac{1}{J} y_\xi \partial_\eta \left( \frac{1}{J} y_\eta \right) \\ & = \frac{y_\eta}{J} \left( - \frac{J_\xi}{J^2} y_\eta + \frac{1}{J} y_{\xi \eta} \right) - \frac{y_\xi}{J} \left( - \frac{J_\eta}{J^2} y_\eta + \frac{1}{J} y_{\eta \eta} \right) \end{aligned} $$ and \begin{aligned} \xi_{yy} & = \ldots \\ & = \frac{x_\eta}{J} \left( - \frac{J_\xi}{J^2} x_\eta + \frac{1}{J} x_{\xi \eta} \right) - \frac{x_\xi}{J} \left( - \frac{J_\eta}{J^2} x_\eta + \frac{1}{J} x_{\eta \eta} \right) \end{aligned}
Inserting these into the Poisson equation and multiply through by $J^2$ I get
$$ \begin{aligned} \beta \frac{J_\eta}{J} - \alpha \frac{J_\xi}{J} + y_\eta y_{\xi \eta} - y_\xi y_{\eta \eta} + x_\eta x_{\xi \eta} - x_\xi x_{\eta \eta} & = J^2 P \end{aligned}$$
and similarly for the Poisson equation for $\eta$ we arrive at
$$ \begin{aligned} \beta \frac{J_\xi}{J} - \gamma \frac{J_\eta}{J} + y_\eta y_{\xi \xi} - y_\xi y_{\eta \xi} - x_\eta x_{\xi \xi} + x_\xi x_{\xi \eta} & = J^2 Q \end{aligned}$$
From this point onwards it's not clear to me how to bring this into the form of Thompson's equations. I have tried expanding out the Jacobians and their derivatives but the algebra becomes incredibly tedious.
Attempt 2
Instead I refer to the Laplace-Beltrami operator. If I calculate the Laplacian of the inverse functions $x(\xi , \eta)$ and $y(\xi , \eta)$ in the coordinate system $(\xi , \eta)$ then we have
$$ \begin{aligned} \nabla^2 x & = \frac{1}{\sqrt{|g|}} \partial_\mu ( \sqrt{|g|} g^{\mu \nu} \partial_\nu x ) \\ & = \frac{1}{\sqrt{|g|}} \partial_\mu ( \sqrt{|g|} g^{\mu \nu}) \partial_\nu x + g^{\mu \nu} \partial_\mu \partial_\nu x \end{aligned} $$
where $g_{\mu \nu} $ is the metric tensor and $|g| = \det g_{\mu \nu}$. If we set $\nabla^2 x = 0$ then we find that
$$ \begin{aligned} g^{\mu \nu} \partial_\mu \partial_\nu x & = - \frac{1}{\sqrt{|g|}} \partial_\mu ( \sqrt{|g|} g^{\mu \nu}) \partial_\nu x \\ \Rightarrow \alpha x_{\xi \xi} - 2\beta x_{\xi \eta} + \gamma x_{\eta \eta} & = - J \partial_\mu ( J g^{\mu \nu}) \partial_\nu x \end{aligned} $$
where I used the fact that the matrix elements of $g^{\mu \nu}$ are the given by inverse of the matrix $g_{\mu \nu}$, and used the fact that the metric in the $(\xi, \eta)$ frame is related to the metric in the Cartesian $(x,y)$ frame is
$$
g_{\mu \nu} = \frac{\partial x^i}{\partial \xi^\mu } \frac{\partial x^j}{\partial x^\nu} \delta_{ij}
$$
which relates the metric to the Jacobian matrix. This looks very close, and takes the form of Thompson's equations, but it is not clear where $P$ and $Q$ come from as I started from the equation $\nabla^2 x = 0$ and not $\nabla^2 \xi = P$ etc.
Any help would be appreciated.