219 lines
7.7 KiB
Matlab
219 lines
7.7 KiB
Matlab
function [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin)
|
|
% [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin)
|
|
% Determines the global index of grid points which lie within the bounds
|
|
% xb,yb,zb specified in the parameters, i.e. the grid points for which
|
|
% xb(1) <= x <= xb(2),
|
|
% yb(1) <= y <= yb(2),
|
|
% zb(1) <= z <= zb(2).
|
|
% If the domain is periodic, the specified region can lie "outside" of the domain,
|
|
% i.e. it is automatically translated to coordinates in the inner domain.
|
|
% If xb(1)==xb(2), a slice located at the nearest neighbor point will be created.
|
|
% Input
|
|
% fdir directory with files
|
|
% iseq sequence index
|
|
% field field to be indexed
|
|
% {'u','v','w','p','s1','s2',...}
|
|
% xb,yb,zb arrays which describe interval, e.g. [0,1] for all points between 0 and 1
|
|
% ? fuzzy Include the first point outside of the region? (default: 0)
|
|
% ? verbosity verbose output? (default: 0)
|
|
% ? debug debug output? (default: 0)
|
|
% Output
|
|
% ii,jj,kk array with field indices
|
|
% xq,yq,zq corresponding grid, monotonically increasing, even if periodic boundary is crossed
|
|
|
|
% Parse input
|
|
par = inputParser;
|
|
addParamValue(par,'verbosity',0,@isnumeric);
|
|
addParamValue(par,'debug',0,@isnumeric);
|
|
addParamValue(par,'fuzzy',0,@isnumeric);
|
|
parse(par,varargin{:});
|
|
iverb = par.Results.verbosity;
|
|
idebug= par.Results.debug;
|
|
ifuzzy= par.Results.fuzzy;
|
|
|
|
% Parse bounds
|
|
xmin = xb(1);
|
|
xmax = xb(2);
|
|
ymin = yb(1);
|
|
ymax = yb(2);
|
|
zmin = zb(1);
|
|
zmax = zb(2);
|
|
|
|
% Parse field
|
|
switch field(1)
|
|
case 'u'; ifield=1; fbase='uvwp'; cmesh='u';
|
|
case 'v'; ifield=2; fbase='uvwp'; cmesh='v';
|
|
case 'w'; ifield=3; fbase='uvwp'; cmesh='w';
|
|
case 'p'; ifield=4; fbase='uvwp'; cmesh='p';
|
|
case 's'; ifield=str2double(field(2:end)); fbase='scal'; cmesh='p';
|
|
otherwise; error('Invalid field: %s',field);
|
|
end
|
|
|
|
% Read auxillary files
|
|
fparam = sprintf('%s/parameters_%04d.asc',fdir,iseq);
|
|
fgrid = sprintf('%s/grid_%04d.bin',fdir,iseq);
|
|
fproc = sprintf('%s/proc_%04d.bin',fdir,iseq);
|
|
params = read_parameters_ucf(fparam,'tarmode',0,'verbosity',iverb);
|
|
[x.u,y.u,z.u,x.v,y.v,z.v,x.w,y.w,z.w,x.p,y.p,z.p] = read_grid_ucf(fgrid,'tarmode',0,'verbosity',iverb,'debug',idebug);
|
|
x = x.(cmesh);
|
|
y = y.(cmesh);
|
|
z = z.(cmesh);
|
|
[ibeg.u,iend.u,jbeg.u,jend.u,kbeg.u,kend.u,...
|
|
ibeg.v,iend.v,jbeg.v,jend.v,kbeg.v,kend.v,...
|
|
ibeg.w,iend.w,jbeg.w,jend.w,kbeg.w,kend.w,...
|
|
ibeg.p,iend.p,jbeg.p,jend.p,kbeg.p,kend.p] = read_procgrid_ucf(fproc,'tarmode',0,'verbosity',iverb,'debug',idebug);
|
|
ibeg = ibeg.(cmesh);
|
|
jbeg = jbeg.(cmesh);
|
|
kbeg = kbeg.(cmesh);
|
|
iend = iend.(cmesh);
|
|
jend = jend.(cmesh);
|
|
kend = kend.(cmesh);
|
|
nx = numel(x);
|
|
ny = numel(y);
|
|
nz = numel(z);
|
|
|
|
% Parse domain information
|
|
a = params.geometry.a;
|
|
b = params.geometry.b;
|
|
c = params.geometry.c;
|
|
d = params.geometry.d;
|
|
e = params.geometry.e;
|
|
f = params.geometry.f;
|
|
x_periodic = params.geometry.xperiodic;
|
|
y_periodic = params.geometry.yperiodic;
|
|
z_periodic = params.geometry.zperiodic;
|
|
|
|
% Check if bounds are valid
|
|
if xmax<xmin
|
|
error('xmax must be greater than xmin.');
|
|
end
|
|
if ymax<ymin
|
|
error('ymax must be greater than ymin.');
|
|
end
|
|
if zmax<zmin
|
|
error('zmax must be greater than zmin.');
|
|
end
|
|
|
|
% Check if min and max are the same point (if so, treat it as slice)
|
|
feps = 1e-12;
|
|
islice = false;
|
|
jslice = false;
|
|
kslice = false;
|
|
if abs(xmax-xmin)<feps
|
|
islice = true;
|
|
end
|
|
if abs(ymax-ymin)<feps
|
|
jslice = true;
|
|
end
|
|
if abs(zmax-zmin)<feps
|
|
kslice = true;
|
|
end
|
|
|
|
% Check values outside of domain, circular shift if periodic, otherwise error
|
|
% Keep the number of shifts and initial offset to reconstruct grid
|
|
if x_periodic
|
|
nxloop = floor((xmax-xmin)/(b-a));
|
|
nxoffset = floor(xmin/(b-a));
|
|
xmin = mod(xmin-a,b-a)+a;
|
|
xmax = mod(xmax-a,b-a)+a;
|
|
elseif xmin<a || xmin>b || xmax<a || xmax>b
|
|
error('Interval must be located inside domain.');
|
|
else
|
|
nxloop = 0;
|
|
nxoffset = 0;
|
|
end
|
|
if y_periodic
|
|
nyloop = floor((ymax-ymin)/(d-c));
|
|
nyoffset = floor(ymin/(d-c));
|
|
ymin = mod(ymin-c,d-c)+c;
|
|
ymax = mod(ymax-c,d-c)+c;
|
|
elseif ymin<c || ymin>d || ymax<c || ymax>d
|
|
error('Interval must be located inside domain.');
|
|
else
|
|
nyloop = 0;
|
|
nyoffset = 0;
|
|
end
|
|
if z_periodic
|
|
nzloop = floor((zmax-zmin)/(f-e));
|
|
nzoffset = floor(zmin/(f-e));
|
|
zmin = mod(zmin-e,f-e)+e;
|
|
zmax = mod(zmax-e,f-e)+e;
|
|
elseif zmin<e || zmin>f || zmax<e || zmax>f
|
|
error('Interval must be located inside domain.');
|
|
else
|
|
nzloop = 0;
|
|
nzoffset = 0;
|
|
end
|
|
|
|
% Find gridpoints closest to the bounds
|
|
if ifuzzy
|
|
[~,imin] = find((x-xmin)<0,1,'last');
|
|
[~,imax] = find((x-xmax)>0,1,'first');
|
|
[~,jmin] = find((y-ymin)<0,1,'last');
|
|
[~,jmax] = find((y-ymax)>0,1,'first');
|
|
[~,kmin] = find((z-zmin)<0,1,'last');
|
|
[~,kmax] = find((z-zmax)>0,1,'first');
|
|
else
|
|
[~,imin] = find((x-xmin)>=0,1,'first');
|
|
[~,imax] = find((x-xmax)<=0,1,'last');
|
|
[~,jmin] = find((y-ymin)>=0,1,'first');
|
|
[~,jmax] = find((y-ymax)<=0,1,'last');
|
|
[~,kmin] = find((z-zmin)>=0,1,'first');
|
|
[~,kmax] = find((z-zmax)<=0,1,'last');
|
|
end
|
|
if isempty(imax); imax=1; end
|
|
if isempty(jmax); jmax=1; end
|
|
if isempty(kmax); kmax=1; end
|
|
|
|
% If a slice is requested, correct previous results
|
|
if islice
|
|
[~,imin] = min(abs(x-xmin));
|
|
imax = imin;
|
|
end
|
|
if jslice
|
|
[~,jmin] = min(abs(y-ymin));
|
|
jmax = jmin;
|
|
end
|
|
if kslice
|
|
[~,kmin] = min(abs(z-zmin));
|
|
kmax = kmin;
|
|
end
|
|
|
|
% Construct the points which are inside of the interval and reconstruct the grid
|
|
if imin<=imax
|
|
ii = [repmat(mod(imin+(0:nx-1)-1,nx)+1,[1,nxloop]),imin:imax];
|
|
shift = [reshape((repmat([zeros(1,nx-imin+1),ones(1,imin-1)],nxloop,1)+(1:nxloop)')',1,[]),...
|
|
ones(1,imax-imin+1)+nxloop]-1;
|
|
xq = x(ii)+(shift+nxoffset)*(b-a);
|
|
else
|
|
ii = [imin:nx,repmat(1:nx,[1,nxloop]),1:imax];
|
|
shift = [ones(1,nx-imin+1),...
|
|
reshape((repmat(ones(1,nx),nxloop,1)+(1:nxloop)')',1,[]),...
|
|
ones(1,imax)+nxloop+1]-1;
|
|
xq = x(ii)+(shift+nxoffset)*(b-a);
|
|
end
|
|
if jmin<=jmax
|
|
jj = [repmat(mod(jmin+(0:ny-1)-1,ny)+1,[1,nyloop]),jmin:jmax];
|
|
shift = [reshape((repmat([zeros(1,ny-jmin+1),ones(1,jmin-1)],nyloop,1)+(1:nyloop)')',1,[]),...
|
|
ones(1,jmax-jmin+1)+nyloop]-1;
|
|
yq = y(jj)+(shift+nyoffset)*(d-c);
|
|
else
|
|
jj = [jmin:ny,repmat(1:ny,[1,nyloop]),1:jmax];
|
|
shift = [ones(1,ny-jmin+1),...
|
|
reshape((repmat(ones(1,ny),nyloop,1)+(1:nyloop)')',1,[]),...
|
|
ones(1,jmax)+nyloop+1]-1;
|
|
yq = y(jj)+(shift+nyoffset)*(d-c);
|
|
end
|
|
if kmin<=kmax
|
|
kk = [repmat(mod(kmin+(0:nz-1)-1,nz)+1,[1,nzloop]),kmin:kmax];
|
|
shift = [reshape((repmat([zeros(1,nz-kmin+1),ones(1,kmin-1)],nzloop,1)+(1:nzloop)')',1,[]),...
|
|
ones(1,kmax-kmin+1)+nzloop]-1;
|
|
zq = z(kk)+(shift+nzoffset)*(f-e);
|
|
else
|
|
kk = [kmin:nz,repmat(1:nz,[1,nzloop]),1:kmax];
|
|
shift = [ones(1,nz-kmin+1),...
|
|
reshape((repmat(ones(1,nz),nzloop,1)+(1:nzloop)')',1,[]),...
|
|
ones(1,kmax)+nzloop+1]-1;
|
|
zq = z(kk)+(shift+nzoffset)*(f-e);
|
|
end
|
|
end |