Discretization of the Anisotropic Diffusion Operator for Finite Difference Method I have to derive and apply a Finite Difference scheme to solve a steady state, anisotropic, diffusion equation. So I have to find a discretization of the following equation

$$\nabla \cdot ( \mathbf{K}(x,y) \mathbf{\nabla} u(x,y) ) + f = 0$$

where 

$$\mathbf{K}(x,y) = \left(\begin{matrix} A(x,y) &amp;&amp; B(x,y) \\ C(x,y) &amp;&amp; D(x,y) \end{matrix}\right)$$

are the components of the space-dependent diffusion/permeability tensor.

Using some standard approach found for example on "Morton K. W., Mayers D.F. - Numerical Solution of Partial Differential Equations, An Introduction (Cambridge University Press, 2005" I translated the local, differential, equation, into an integral form and I wrote a flux conservation law for each rectangular cell of my domain.

If $\mathbf{K}$ is diagonal all is fine, I end un with something like a 5-point stencil for the laplacian with adjusted "weigths".

But if $B \ne 0$ (usually $B = C$) I end-up with an issue, because, I think, in the contribution to the total flux from a generic edge of the cell, appear the derivative of $u$ along the direction of the edge.

I'm not sure how to estimate this derivative, because the edge is at half way between grid nodes...

$$
0 \approx \Delta x_i \Delta y_j f_{ij} + \Delta y_j \left(A u_x + B u_y\right)_{x_{i-1/2},y_j}^{x_{i+1/2},y_j} + ...
$$

$$
0 \approx \Delta x_i \Delta y_j f_{ij} + \Delta y_j \left\{ A_{i+1/2,j} \times 2 \frac{ u_{i+1,j} - u_{i,j} }{\Delta x_i + \Delta x_{i+1}} - A_{i-1/2,j} \times 2 \frac{ u_{i,j} - u_{i-1,j} }{\Delta x_{i-1} + \Delta x_{i}} + B_{i+1/2,j} \times ? \right\} + ...
$$

The only ideas that comes to my mind are:


Use some near grid nodes, f.i:


$$
(u_y)_{i+1/2,j} \approx (u_y)_{i+1,j} \approx 2 \times \frac{u_{i+1,j+1} - u_{i+1,j-1}}{\Delta y_{j-1} + 2 \Delta y_j + \Delta y_{j+1}}
$$


Use some weighted mean / interpolation between nearest nodes: 


$$
(u_y)_{i+1/2,j} \approx \frac{\Delta x_{i+1} (u_y)_{i,j} + \Delta x_{i} (u_y)_{i+1,j}}{\Delta x_{i} + \Delta x_{i+1}} 
$$


Use a non-compact, bigger, 13-point (?) stencil and write an equation every $2\times2$ cells of the discretization "mesh"


Unfortunately I didn't find any reference on this topic.
Maybe this would help https://www.mia.uni-saarland.de/weickert/Papers/book.pdf#page=107

It gives the stencil for a laplacian of the form $\text{div}(D\nabla)$ where : 
$D = \left(\begin{matrix} a &amp; b \\ b &amp; c \end{matrix}\right)$
Look this:
https://www.sciencedirect.com/science/article/pii/S0021999114003155

This causes diffusion, heat conduction, and viscous momentum loss, to effectively be aligned with the magnetic field lines. This alignment leads to different values for the respective diffusive coefficients in the magnetic field direction and in the perpendicular direction, to the extent that heat diffusion coefficients can be up to 1012 times larger in the parallel direction than in the perpendicular direction. This anisotropy puts stringent requirements on the numerical methods used to approximate the MHD-equations since any misalignment of the grid may cause the perpendicular diffusion to be polluted by the numerical error in approximating the parallel diffusion. Currently the common approach is to apply magnetic field-aligned coordinates, an approach that automatically takes care of the directionality of the diffusive coefficients. This approach runs into problems at x-points and at points where there is magnetic re-connection, since this causes local non-alignment. It is therefore useful to consider numerical schemes that are tolerant to the misalignment of the grid with the magnetic field lines, both to improve existing methods and to help open the possibility of applying regular non-aligned grids. To investigate this, in this paper several discretization schemes are developed and applied to the anisotropic heat diffusion equation on a non-aligned grid.