FFT on non-orthogonal lattice ( for computing convolutions and solving PDEs )

I saw many examples of application of FFT for computation of convolutions and solving PDEs ( like Poisson equation ). It is very strightforward and efficient if I work with rectangular (orthogonal) lattice.

However, I don't see how exatly modify the transformation if the lattice is sampled along non-orthogonal lattice vectors

eg. like this lattice: I guess, I can still use common algorithms or libraries ( like FFTW) because they care just about the n-dimensional array of numbers (they are agnostic to position of sampling points in real space ) ... but after that I have to transfrom the results by the lattice vectors somehow (or it's reciprocal vectors)( ? am I right? ).

Especially I'm interested in three application:

1. solving posisson equation
2. covolution of two functions
3. gradient ofthe function.

That is not very difficult actually. Let $\Lambda\subset\Bbb R^n$ be the lattice on which you sampled your points, and let $\Phi:\Bbb R^n\to\Bbb R^n$ be a linear transformation under which $\Bbb Z^n\to\Lambda$ isomorphically (i.e. the columns of the matrix of $\Phi$ generate $\Lambda$). Let $f:\Bbb R^n\to\Bbb C$ be the function you are interested in, and of which you are interested in its values in $\Lambda$. From now on, we will alternatingly identify $f$ with a function on $\Lambda$ taking the value $f(\lambda)$ at $\lambda\in\Lambda$, and with a function (or rather distribution)

$$\sum_{\lambda\in\Lambda}f(\lambda)\delta(\lambda - x)$$

on $\Bbb R^n$, i.e. zero outside $\Lambda$, and a delta function at lattice points.

We can obtain the properties of $f$ in terms of its Fourier transform from the Fourier transform of $f\circ\Phi$, which is a function on $\Bbb Z^n$:

$$\mathscr F(f\circ\Phi)(s) \equiv \int f\circ\Phi(x)e^{-2\pi i\langle s,x\rangle}dx = \frac1{|\det\Phi|}\int f(y)e^{-2\pi i\langle\Phi^{-T}s,x\rangle}dy = \frac1{|\det\Phi|}\mathscr F(f)(\Phi^{-T}s)$$

Here a substitution $y = \Phi^{-1}x$ was made, and it was used that $\langle s,\Phi^{-1}y\rangle = \langle\Phi^{-T}s,y\rangle$. Note that your conventions w.r.t. the sign of the exponent and factors of $\sqrt{2\pi}$ may differ.

Now what does the formula mean that we just wrote down? $f\circ\Phi$ is just our set of samples of $f$ arranged in an $n$-dimensional array, whose values have to be interpreted as lying on a standard grid. Its Fourier transform gives an array of the same shape, but the interpretation is different: the value at index $s$ in reality is the value at $\Phi^{-T}s$ in $\Bbb R^n$, multiplied by $|\det\Phi|$.

In particular, if the function $f$ is supported on $\Lambda\subset\Bbb R^n$, then its Fourier transform is supported on the image of $\Phi^{-T}$, the span of the columns of the inverse transposed matrix of $\Phi$.