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

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,y\rangle}dy = \frac1{|\det\Phi|}\mathscr F(f)(\Phi^{-T}s)$$

Here a substitution $$x = \Phi^{-1}y$$ 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$$.