Error: Edge vector must be monotonically non-decreasin with isosurface

I have the following code where I am taking 3D FFT for 3D matrix and comparing its derivatives to the "exact" values, but I am getting the error:

Edge vector must be monotonically non-decreasing.

Error in hist (line 152)
    nn = histc(y,edgesc,1);

Error in isovalue (line 16)
[n, ctrs] = hist(data(1:r:end),100);

Error in isosurface (line 87)
  value = isovalue(data);

Error in Fourier3D_nonlinear (line 108)
isosurface(X,Y,Z,dux); 

The full code:

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(X,Y,Z,dux); %ERROR here
title('\partial u/\partial x');
figure
isosurface(X,Y,Z,ExactDux); 
title('Exact \partial u/\partial x');

Answer 1

[Chunru]

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

[Jamie Al]

Thanks a lot. But, I wonder why this instead won't work:

dux = real(ifft(ifft(ifft(duxk, [], 1), [], 2), [], 3));

[Jamie Al]

Also, why the y and z derivative plots look slightly different? Does this have to do with the way I am taking 3D FFT?

figure(3)
isosurface(X,Y,Z,real(duy)); 
title('\partial u/\partial y');
figure(4)
isosurface(X,Y,Z,ExactDuy); 
title('Exact \partial u/\partial y');
%%
figure(5)
isosurface(X,Y,Z,real(duz)); 
title('\partial u/\partial z');
figure(6)
isosurface(X,Y,Z,ExactDuz); 
title('Exact \partial u/\partial z');

[Chunru]

dux = real(ifft(ifft(ifft(duxk, [], 1), [], 2), [], 3));

The above works for me. Your ifft has no problem.

For your second question, I am not sure what is the problem. But my guess is something to do with the exact solution formula. Double check it.

[Jamie Al]

Thanks!