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

3 views (last 30 days)
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');

Accepted Answer

Chunru
Chunru on 4 Oct 2022
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
  4 Comments
Chunru
Chunru on 4 Oct 2022
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.

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!