120 lines
2.1 KiB
Matlab
120 lines
2.1 KiB
Matlab
function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
|
|
xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
|
|
generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
|
|
% [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
|
|
% xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,dx,dy,dz]=...
|
|
% generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,xperiodic,yperiodic,zperiodic)
|
|
% Generates a staggered grid.
|
|
% Input
|
|
% a,b,c,d,e,f domain bounds
|
|
% nxp,nyp,nzp number of points for pressure grid
|
|
% x/y/zperiodic domain periodicity
|
|
|
|
if xperiodic==1
|
|
nxu=nxp;
|
|
nxv=nxp;
|
|
nxw=nxp;
|
|
else
|
|
nxu=nxp+1;
|
|
nxv=nxp;
|
|
nxw=nxp;
|
|
end
|
|
if yperiodic==1
|
|
nyu=nyp;
|
|
nyv=nyp;
|
|
nyw=nyp;
|
|
else
|
|
nyu=nyp;
|
|
nyv=nyp+1;
|
|
nyw=nyp;
|
|
end
|
|
if zperiodic==1
|
|
nzu=nzp;
|
|
nzv=nzp;
|
|
nzw=nzp;
|
|
else
|
|
nzu=nzp;
|
|
nzv=nzp;
|
|
nzw=nzp+1;
|
|
end
|
|
%grid step:
|
|
if xperiodic==1
|
|
dx=(b-a)/(nxp);
|
|
else
|
|
dx=(b-a)/(nxp-1);
|
|
end
|
|
if yperiodic==1
|
|
dy=(d-c)/(nyp);
|
|
else
|
|
dy=(d-c)/(nyp-1);
|
|
end
|
|
if zperiodic==1
|
|
dz=(f-e)/(nzp);
|
|
else
|
|
dz=(f-e)/(nzp-1);
|
|
end
|
|
%
|
|
if xperiodic==1
|
|
xu=a+((1:nxu)-1)*dx;
|
|
else
|
|
xu=a+((1:nxu)-3/2)*dx;
|
|
end
|
|
if yperiodic==1
|
|
yu=c+((1:nyu)-1/2)*dy;
|
|
else
|
|
yu=c+((1:nyu)-1)*dy;
|
|
end
|
|
if zperiodic==1
|
|
zu=e+((1:nzu)-1/2)*dz;
|
|
else
|
|
zu=e+((1:nzu)-1)*dz;
|
|
end
|
|
%
|
|
if xperiodic==1
|
|
xv=a+((1:nxv)-1/2)*dx;
|
|
else
|
|
xv=a+((1:nxv)-1)*dx;
|
|
end
|
|
if yperiodic==1
|
|
yv=c+((1:nyv)-1)*dy;
|
|
else
|
|
yv=c+((1:nyv)-3/2)*dy;
|
|
end
|
|
if zperiodic==1
|
|
zv=e+((1:nzv)-1/2)*dz;
|
|
else
|
|
zv=e+((1:nzv)-1)*dz;
|
|
end
|
|
%
|
|
if xperiodic==1
|
|
xw=a+((1:nxw)-1/2)*dx;
|
|
else
|
|
xw=a+((1:nxw)-1)*dx;
|
|
end
|
|
if yperiodic==1
|
|
yw=c+((1:nyw)-1/2)*dy;
|
|
else
|
|
yw=c+((1:nyw)-1)*dy;
|
|
end
|
|
if zperiodic==1
|
|
zw=e+((1:nzw)-1)*dz;
|
|
else
|
|
zw=e+((1:nzw)-3/2)*dz;
|
|
end
|
|
%
|
|
if xperiodic==1
|
|
xp=a+((1:nxp)-1/2)*dx;
|
|
else
|
|
xp=a+((1:nxp)-1)*dx;
|
|
end
|
|
if yperiodic==1
|
|
yp=c+((1:nyp)-1/2)*dy;
|
|
else
|
|
yp=c+((1:nyp)-1)*dy;
|
|
end
|
|
if zperiodic==1
|
|
zp=e+((1:nzp)-1/2)*dz;
|
|
else
|
|
zp=e+((1:nzp)-1)*dz;
|
|
end
|
|
end |