Using the Laplacian definition $\nabla^2f(x,y,z) = \frac{\partial^2f(x,y,z)}{\partial x^2}+\frac{\partial^2f(x,y,z)}{\partial y^2}+\frac{\partial^2f(x,y,z)}{\partial z^2}$ it is easy to get a numerical estimation if I have the $f$ values at the points of a regular orthogonal grid, with spacings $d_x, d_y, d_z$:
\begin{equation} \nabla^2f(x,y,z) \approx \frac{f(x+d_x,y,z)+f(x-d_x,y,z)}{d_x^2} + \frac{f(x,y+d_y,z)+f(x,y-d_y,z)}{d_y^2} + \frac{f(x,y,z+d_z)+f(x,y,z-d_z)}{d_z^2} - 2\left(\frac{1}{d_x^2}+\frac{1}{d_y^2}+\frac{1}{d_z^2}\right)f(x,y,z) \end{equation}
Even if the axes of the grid are not aligned with the Cartesian axes, I just "rename" $x,y,z$.
But this doesn't work if I have a non-orthogonal grid (or curvilinear, but I'm interested in the linear case now). Is there an "easy" formula like the above I could use? I guess I would need to use 18 points instead of 6.
OK, so if I got it right here is the solution.
We start from a regular grid, each point of the grid obeys to $\mathbf{r}_0 + n_a\mathbf{a} + n_b\mathbf{b} + n_c\mathbf{c}$, where $\mathbf{r}_0$ is some point of the grid, $n_a,n_b,n_c$ are integers and $\mathbf{a},\mathbf{b},\mathbf{c}$ are three linearly independent vectors (not necessarily unit vectors, nor orthogonal). We have the values of the scalar function $f$ at each point of the grid.
According to https://en.wikipedia.org/wiki/Tensors_in_curvilinear_coordinates#Scalar_field_2 the Laplacian is:
\begin{equation} \nabla^2 f = \sum_{i=a,b,c}\frac{1}{\sqrt{g}}\frac{\partial}{\partial n_i}\left(\sum_{j=a,b,c} g^{ji}\frac{\partial f}{\partial n_j}\sqrt{g}\right) \end{equation}
$[g_{ij}]$ is a symmetric matrix with elements $\mathbf{i}\cdot\mathbf{j}$, $[g^{ij}]$ is its inverse and $g^{ij} = g^{ji}$ is an element of the latter, while $g = \det{[g_{ij}]}$. Since $\mathbf{a},\mathbf{b},\mathbf{c}$ are constant, $g$ and $g^{ij}$ are constant, and:
\begin{equation} \nabla^2 f = \frac{\sqrt{g}}{\sqrt{g}}\sum_{i=a,b,c}\frac{\partial}{\partial n_i}\left(\sum_{j=a,b,c} g^{ji}\frac{\partial f}{\partial n_j}\right) = \sum_{i,j=a,b,c}g^{ij}\frac{\partial^2 f}{\partial n_i \partial n_j} \end{equation}
Now evaluating the second derivatives numerically should be easy. In hopefully evident notation:
\begin{gather} \frac{\partial^2 f}{\partial n_a^2} \approx f(+1,0,0)+f(-1,0,0)-2f(0,0,0) \\ \frac{\partial^2 f}{\partial n_a \partial n_b} \approx \frac{f(+1,+1,0)+f(-1,-1,0)+f(+1,-1,0)-f(-1,+1,0)}{4} \end{gather}
So finally the Laplacian should be:
\begin{equation} \nabla^2 f \approx -2(g^{aa}+g^{bb}+g^{cc})f(0,0,0) + g^{aa}[f(+1,0,0)+f(-1,0,0)] + g^{bb}[f(0,+1,0)+f(0,-1,0)] + g^{cc}[f0,0,+1)+f(0,0,-1)] + \frac{g^{ab}}{2}[f(+1,+1,0)+f(-1,-1,0)-f(+1,-1,0)-f(-1,+1,0)] + \frac{g^{ac}}{2}[f(+1,0,+1)+f(-1,0,-1)-f(+1,0,-1)-f(-1,0,+1)] + \frac{g^{bc}}{2}[f(0,+1,+1)+f(0,-1,-1)-f(0,+1,-1)-f(0,-1,+1)] \end{equation}
If $\mathbf{a},\mathbf{b},\mathbf{c}$ are orthogonal, $[g_{ij}]$ is diagonal, $g^{ii}$ is simply $1/|\mathbf{i}|^2$ and $g^{ij}=0$ for $i\neq j$, so this reduces to the formula in the question. Numerical experiments seem to indicate this is correct.