Given a vector field $C$ on $[0,1]\times [0,1]$ which is divergence free, i.e., $divC=0$. Now for any scalar field $f$ on $[0,1]\times [0,1]$, consider the operator $\mathcal Lf=C\cdot \nabla f$. Then the adjoint of $\mathcal L$, $\mathcal L^* f=-div(fC)=-\nabla f\cdot C-fdivC=-C\cdot \nabla f$.
I would like to discretize the above operator $\mathcal L$ making it to be matrix form. Suppose the space is partitioned evenly with size $h$. Every function is periodic to the boudary. To discretize the gradient operator, I use: $$ \nabla f(x) = \frac{f(x+he_1)-f(x-he_1)}{2h}i+\frac{f(x+he_2)-f(x-he_2)}{2h}j. $$
Now here comes the problem. By the calculation in the first part, $L^Tf=-Lf$, where $L$ is the discretized matrix of $\mathcal L$. But using the origianl $C$, $L$ may not still be skew-symmetric. My attempt shows that $C$ must satisfy: $$ C(x_1,x_2)=C_1(x_1,x_2)i+C_2(x_1,x_2)j=C_1(x_2)i+C_2(x_1)j $$ that is $C_1$ is a constant with respect to $x_1$ and $C_2$ similary.
These $C$'s are rather limited. Is there any way to consider the larger class of $C$ and $L$ still skew-symmetric? (Maybe changing how I discretize the gradient operator?)
Thanks.