Nx = 32;
Ny = 32;
Nz = 32;
Lx = 2*pi;
Ly = 2*pi;
Lz = 2*pi;
x = (0:Nx-1)/Nx*2*pi; % x coordinate in Fourier, equally spaced grid
y = (0:Ny-1)/Ny*2*pi;
z = (0:Nz-1)/Nz*2*pi;
dx = Lx/Nx;
dy = Ly/Ny;
dz = Lz/Nz;
kx = makeK(Lx,Nx)';
ky = makeK(Ly,Ny)';
kz = makeK(Lz,Nz)';
%$$$$$$$$$$$$$$$$$$$$
ksqu = (sin( kx * dx/2)/(dx/2)).^2 + (sin( ky * dy/2)/(dy/2)).^2 + (sin( kz * dz/2)/(dz/2)).^2 ;
kx = sin(kx * dx) / dx;
ky = sin(ky * dy) / dy;
kz = sin(kz * dz) / dz;
%-----------
%make mesh
[X,Y,Z] = meshgrid(x,y,z);
A = 2*pi / Lx;
B = 2*pi / Ly;
C = 2*pi / Lz;
%%
u = sin(C*Z).^2 .* sin(A*X) .* cos(B*Y);
uh = fft(fft(fft(u, [], 1), [], 2), [], 3);
duxk = derivk(uh,kx);
dux = ifft(ifft(ifft(duxk, [], 1), [], 2), [], 3);
duyk = derivk(uh,reshape(ky, [Ny,1,1]));
duy = ifft(ifft(ifft(duyk, [], 1), [], 2), [], 3);
duzk = derivk(uh,reshape(kz, [1,1,Nz]));
duz = ifft(ifft(ifft(duzk, [], 1), [], 2), [], 3);
ExactDux = sin(C*Z).^2 .* A.* cos(A*X) .* cos(B*Y) ;
DEx = ExactDux - dux;
ExactDuy = -B .* sin(A*X) .*sin(B*Y) .*sin(C*Z).^2;
DEy = ExactDuy - duy;
ExactDuz = 2*C .* cos(B*Y) .* cos(C*Z) .*sin(A*X) .* sin(C*Z);
DEz = ExactDuz - duz;
%%
figure
%============================================================
% isosurface receive real data only. use real or abs here
isosurface(X,Y,Z,real(dux)); %ERROR here
title('\partial u/\partial x');
figure
isosurface(X,Y,Z,ExactDux);
title('Exact \partial u/\partial x');
function data_deriv = derivk(fk,k)
% Takes derivative of a 2D matrix using Fourier transforms
data_deriv = 1i * k .* fk;
end
function k = makeK(L, n)
% Makes a k vector for Matlab fft using the length and number of points
% Tested and works
k = (2 * pi / L) * [0:n/2-1 , -n/2:-1]';
end
