ucftools/matlab/read_uvwp_complete_ucf.m

function [u,v,w,p] = read_uvwp_complete_ucf(file,varargin)
    % [u,v,w,p] = read_uvwp_complete_ucf(file,varargin)
    % Reads u,v,w,p from all processor chunks and combines them to a single field.
    % Input
    %   file              file name (if tar-mode: ustar handle)
    % ? step              index of timestep to be read (default: 1)
    % ? verbosity         verbose output? (default: no)
    % ? debug             debug output? (default: no)
    % ? tarmode           read from tar-file? (default: 0)
    % Output
    %   u,v,w,p           complete fields with dim(nx,ny,nz)

    % Parse optional input arguments
    par = inputParser;
    addParamValue(par,'step',1,@isnumeric);
    addParamValue(par,'verbosity',0,@isnumeric);
    addParamValue(par,'debug',0,@isnumeric);
    addParamValue(par,'tarmode',0,@isnumeric); % deprecated
    parse(par,varargin{:});
    istep = par.Results.step;

    % Define sets to be read
    sets = {'u','v','w','p'};
    nset = numel(sets);

    % Open first file
    obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
    switch class(file)
    case 'char'
      [fdir,fbase,fext] = fileparts(file);
      fname = sprintf('%s/%s.%05d',fdir,fbase,0);
      obj.open(fname);
    case {'ustar','ucfmulti'}
      fname = sprintf('uvwp.%05d',0);
      obj.opentar(file.pointer(fname));
    otherwise
      error('Input file type not supported: %s',class(file));
    end
    
    if ~obj.validateType('field')
      error('read error: no field data.');
    end

    % Read raw data
    q  = cell(1,nset);
    ib  = cell(1,nset);
    jb  = cell(1,nset);
    kb  = cell(1,nset);
    nxl = cell(1,nset);
    nyl = cell(1,nset);
    nzl = cell(1,nset);
    nx  = cell(1,nset);
    ny  = cell(1,nset);
    nz  = cell(1,nset);
    for iset=1:nset
      [data,params] = obj.readSet(istep,iset);
      params = cast(params,'double');
      ighost   = params(1);
      ib{iset} = params(2);
      jb{iset} = params(3);
      kb{iset} = params(4);
      nxl{iset}= params(5);
      nyl{iset}= params(6);
      nzl{iset}= params(7);
      nx{iset} = params(8);
      ny{iset} = params(9);
      nz{iset} = params(10);
      q{iset}  = reshape(data,nxl{iset}+2*ighost,nyl{iset}+2*ighost,nzl{iset}+2*ighost);
      if ighost
        q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
      end
    end

    % Reassign content from loop and create global arrays
    iset = find(strcmp(sets,'u'));
    u = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
    u(ib{iset}:ib{iset}+nxl{iset}-1,...
      jb{iset}:jb{iset}+nyl{iset}-1,...
      kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

    iset = find(strcmp(sets,'v'));
    v = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
    v(ib{iset}:ib{iset}+nxl{iset}-1,...
      jb{iset}:jb{iset}+nyl{iset}-1,...
      kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

    iset = find(strcmp(sets,'w'));
    w = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
    w(ib{iset}:ib{iset}+nxl{iset}-1,...
      jb{iset}:jb{iset}+nyl{iset}-1,...
      kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

    iset = find(strcmp(sets,'p'));
    p = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
    p(ib{iset}:ib{iset}+nxl{iset}-1,...
      jb{iset}:jb{iset}+nyl{iset}-1,...
      kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

    % Close the first file
    obj.close();

    % Now loop consecutively through files
    ifile = 1;
    switch class(file)
    case 'char'
      loopCondition = @(x) exist(x,'file');
      fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
    case {'ustar','ucfmulti'}
      loopCondition = @(x) file.isSubfile(x);
      fname = sprintf('uvwp.%05d',ifile);
    end
    while loopCondition(fname)
      % Open file
      obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
      switch class(file)
      case 'char'
        obj.open(fname);
      case {'ustar','ucfmulti'}
        obj.opentar(file.pointer(fname));
      end
      
      if ~obj.validateType('field')
        error('read error: no field data.');
      end

      % Read raw data
      q  = cell(1,nset);
      ib  = cell(1,nset);
      jb  = cell(1,nset);
      kb  = cell(1,nset);
      nxl = cell(1,nset);
      nyl = cell(1,nset);
      nzl = cell(1,nset);
      for iset=1:nset
        [data,params] = obj.readSet(istep,iset);
        params = cast(params,'double');
        ighost   = params(1);
        ib{iset} = params(2);
        jb{iset} = params(3);
        kb{iset} = params(4);
        nxl{iset}= params(5);
        nyl{iset}= params(6);
        nzl{iset}= params(7);
        q{iset}  = reshape(data,nxl{iset}+2*ighost,nyl{iset}+2*ighost,nzl{iset}+2*ighost);
        if ighost
          q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
        end
      end

      % Reassign content from loop and create global arrays
      iset = find(strcmp(sets,'u'));
      u(ib{iset}:ib{iset}+nxl{iset}-1,...
        jb{iset}:jb{iset}+nyl{iset}-1,...
        kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

      iset = find(strcmp(sets,'v'));
      v(ib{iset}:ib{iset}+nxl{iset}-1,...
        jb{iset}:jb{iset}+nyl{iset}-1,...
        kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

      iset = find(strcmp(sets,'w'));
      w(ib{iset}:ib{iset}+nxl{iset}-1,...
        jb{iset}:jb{iset}+nyl{iset}-1,...
        kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

      iset = find(strcmp(sets,'p'));
      p(ib{iset}:ib{iset}+nxl{iset}-1,...
        jb{iset}:jb{iset}+nyl{iset}-1,...
        kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};

      % Close file
      obj.close();

      % Move on to next file
      ifile = ifile+1;
      switch class(file)
      case 'char'
        fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
      case {'ustar','ucfmulti'}
        fname = sprintf('uvwp.%05d',ifile);
      end
    end
end