The essence of this question is that I don't know how to convert an equation from Cartesian coordinates into curvilinear coordinates, and would like to know how, preferably using the language of co/contravariance of physical quantities and vector notation.
In Daniel Appelö's A stable finite difference method for the elastic wave equation on complex geometries with free surfaces, he formulates the elastic wave equation in component form, $$\begin{align}\rho\partial_t^2u&=\partial_x\color{red}{\left((2\mu+\lambda)\partial_xu+\lambda\partial_yv\right)}+\partial_y\color{blue}{\left(\mu(\partial_xv+\partial_yu)\right)}\\ \rho\partial_t^2v&=\partial_x\color{blue}{\left(\mu(\partial_xv+\partial_yu)\right)}+\partial_y\color{green}{\left(\lambda\partial_xu+(2\mu+\lambda)\partial_yv\right)}\end{align}$$ where $(u,v)$ is the material displacement and $(x,y)$ are the physical coordinates.
He then defines curvilinear coordinates $(q(x,y),r(x,y))$, from which we obtain $$\begin{align}\partial_x&=q_x\partial_q+r_x\partial_r\\ \partial_y&=q_y\partial_q+r_y\partial_r.\end{align}$$
He then reformulates the wave equation in the curvilinear coordinate system, and I'm having trouble seeing how he gets it. My attempt is as follows: a direct replacement of all $\partial_x$ and $\partial_y$ gives
$$\begin{align}\rho\partial_t^2u=&(q_x\partial_q+r_x\partial_r)\color{red}{\left((2\mu+\lambda)(q_x\partial_q+r_x\partial_r)u+\lambda(q_y\partial_q+r_y\partial_r)v\right)}+\\ &(q_y\partial_q+r_y\partial_r)\color{blue}{\left(\mu((q_x\partial_q+r_x\partial_r)v+(q_y\partial_q+r_y\partial_r)u)\right)}\\ \rho\partial_t^2v=&(q_x\partial_q+r_x\partial_r)\color{blue}{\left(\mu((q_x\partial_q+r_x\partial_r)v+(q_y\partial_q+r_y\partial_r)u)\right)}+\\ &(q_y\partial_q+r_y\partial_r)\color{green}{\left(\lambda(q_x\partial_q+r_x\partial_r)u+(2\mu+\lambda)(q_y\partial_q+r_y\partial_r)v\right)}\end{align} $$
but am not sure how to proceed. Appelö ends up getting (c.f. equation 2.6 and 2.7), in condensed notation,
$$\begin{align}J\rho\partial_t^2u&=\partial_q\left[Jq_x\color{red}{\sigma_{11}}+Jq_y\color{blue}{\sigma_{12}}\right]+\partial_r\left[Jr_x\color{red}{\sigma_{11}}+Jr_y\color{blue}{\sigma_{12}}\right]\\ J\rho\partial_t^2v&=\partial_q\left[Jq_x\color{blue}{\sigma_{21}}+Jq_y\color{green}{\sigma_{22}}\right]+\partial_r\left[Jr_x\color{blue}{\sigma_{21}}+Jr_y\color{green}{\sigma_{22}}\right] \end{align}$$
where $J=x_qy_r-x_ry_q$ is apparently the determinant of the Jacobian of the transform.
How does one proceed to derive a result like this? Is it simply a matter of applying product and chain rule dozens of times and cleaning up the mess? Or is there a simple mechanized way to do it using tensor calculus? I suspect there may be an easier way, since the original equation (in $x,y$ coordinates) can be written in vector form as
$$\rho\ddot{\mathbf{r}}=\sigma\cdot\nabla$$
where $$\sigma=\pmatrix{\color{red}{\sigma_{11}}&\color{blue}{\sigma_{12}}\\\color{blue}{\sigma_{21}}&\color{green}{\sigma_{22}}}$$
is the (symmetric) material stress tensor. I've heard that the stress tensor transforms double-contravariantly, and $\nabla$ transforms covariantly, but I'm not sure how to put together all the pieces and get Appelö's result.
Can anyone provide help on the matter?
Edit 1
So I noticed that Equation 2.6 and 2.7 (the transformed equations of motion) can further be condensed to the following:
$$J\rho\ddot{\mathbf{r}}=\left[J\pmatrix{\color{red}{\sigma_{11}}&\color{blue}{\sigma_{12}}\\\color{blue}{\sigma_{21}}&\color{green}{\sigma_{22}}}\pmatrix{q_x&r_x\\q_y&r_y}\right]\cdot\nabla_{(q,r)}.$$
This is tantalizingly close to the original equations of motion $\rho\ddot{\mathbf{r}}=\sigma\cdot\nabla$ with the covariant substitution
$$\nabla=\pmatrix{q_x&r_x\\q_y&r_y}\nabla_{(q,r)}$$
except I'm missing a pair of $J$'s. I feel like I'm missing some conceptual aspect that would put this all together.