Necessary to convert anisotropic diffusion tensor to polar coordinates? Say I want to solve the steady-state diffusion equation
$$0 = \nabla \cdot \left\{D \nabla c\right\},$$
where I know the diffusion tensor $D$ is anisotropic, and can be described in Cartesian coordinates by 
$$D = \begin{bmatrix} D_x & 0 \\ 0 & D_y \end{bmatrix}.$$
However, my boundary conditions are specified in polar coordinates (2D), so I want to solve the equation in those coordinates, $c(r,\theta)$. 
It's straightforward to convert the differential operators into polar coordinates 
$$\nabla \cdot A = \frac{1}{r} \partial_r \left\{r A_r\right\} + \frac{1}{r} \partial_\theta A_\theta, \qquad  \nabla A = \partial_r A \hat{r} + \frac{1}{r}\partial_\theta A \hat{\theta}.$$
Using these relations, and not doing anything fancy to the diffusion tensor,  I end up with 
$$
0= \frac{D_x r \left(\partial_r c+r \partial_{rr}c \right)+D_y \partial_{\theta \theta} c}{r^2}$$
which seems suspicious because the $D_x$ is only on the $\partial_r$ terms and $D_y$ on the $\partial_\theta$, so when you try separation of variables, you end up with the ratio $D_y/D_x$, and if these were equal, would cancel. 
The only place I can see that went awry is that I didn't "convert" the diffusion tensor into this new coordinate system. However, I'm not really sure how to go about doing this or if it even makes sense to do so. In the original PDE, (I think) the gradient is the gradient, regardless of coordinates so $D\nabla c$ means the same thing.
This might be better posed as a more general question about tensors and coordinate systems, but this was the context it arose in. Any insight is much appreciated.
 A: Yes, the components of the diffusion tensor must be written in cylindrical coordinates too. The problem is that $\nabla\cdot$ takes the divergence of a vector expressed in the basis $(\hat r,\hat\theta)$, while the matrix-vector product $D\nabla c$ where $D = \text{diag}(D_x, D_y)$ is written for a vector $\nabla c$ expressed in the basis $(\hat x,\hat y)$.
To perform the change of basis of $D$, let us introduce the change of coordinates
$$
\hat r = \frac{x \hat x + y \hat y}{\sqrt{x^2 + y^2}} \, ,\qquad
\hat \theta = \frac{-y \hat x + x \hat y}{\sqrt{x^2 + y^2}}\, ,
$$
so that
\begin{aligned}
\nabla c &= \partial_r c\, \hat r + \frac{1}{r} \partial_\theta c\, \hat\theta\\
&= \left(\cos\theta\,\partial_r c - \frac{\sin\theta}{r}\partial_\theta c \right) \hat x + \left( \sin\theta\, \partial_r c + \frac{\cos\theta}{r}\,\partial_\theta c\right) \hat y \, ,
\end{aligned}
which can be multiplied now by the matrix $D = \text{diag}(D_x, D_y)$. The inverse change of coordinates
$$
\hat x = \cos\theta\, \hat r - \sin\theta\, \hat \theta \, ,\qquad
\hat y = \sin\theta\, \hat r + \cos\theta\, \hat \theta \, .
$$
makes it possible to express $D\nabla c$ in the basis $(\hat r ,\hat \theta)$, and then to take the polar divergence $\nabla \cdot (D\nabla c)$.

Alternatively, one can derive by hand the 2D anisotropic Laplace equation in polar coordinates. The equation $\nabla\cdot \left(D\nabla c\right) = 0$ writes as follows in Cartesian coordinates:
$$
D_x \partial^2_{xx} c + D_y \partial^2_{yy} c = 0 \, .
$$
With the convention $(x,y) = (r\cos\theta,r\sin\theta )$, one uses the chain rule and the identity $r^2 = x^2 + y^2$ to obtain
\begin{aligned}
\partial^2_{xx} c &= \partial^2_{rr} c\, \frac{x^2}{r^2} + \partial_{r} c\, \frac{y^2}{r^3} - \partial^2_{r\theta} c\, \frac{2xy}{r^3} + \partial_{\theta} c\, \frac{2xy}{r^4} + \partial^2_{\theta\theta}c\, \frac{y^2}{r^4} \, , \\
\partial^2_{yy} c &= \partial^2_{rr} c\, \frac{y^2}{r^2} + \partial_{r} c\, \frac{x^2}{r^3} + \partial^2_{r\theta} c\, \frac{2xy}{r^3} - \partial_{\theta} c\, \frac{2xy}{r^4} + \partial^2_{\theta\theta}c\, \frac{x^2}{r^4} \, .
\end{aligned}
Now, the PDE can be rewritten in polar coordinates:
\begin{aligned}
&\left(1 + \epsilon\cos^2\theta\right) \partial^2_{rr} c + \frac{1 + \epsilon\sin^2\theta}{r} \partial_{r} c - \epsilon \frac{\sin 2\theta}{r} \partial^2_{r\theta} c \\
&+ \epsilon \frac{\sin 2\theta}{r^2} \partial_{\theta} c + \frac{1 + \epsilon\sin^2\theta}{r^2} \partial^2_{\theta\theta} c = 0 \, ,
\end{aligned}
where $\epsilon = (D_x - D_y)/D_y$ is the diffusion contrast.
