function [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, bnd, periodic)
% function [NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask, bnd, periodicity)
% This function adds boundary or periodic cells at the four sides of the domain and
% creates the node numbering reference matrices for P, u and v grids.
%
% INPUTS:
%       Mask:   The pressure grid mask matrix (i.e. the tracer grid mask matrix)
%               NOTE: Internal boundaries should be numbered from 1 to #internal BCs!
%       bnd:    Vector containing true or false values of which boundaries to add:
%               bnd = [left right bottom top]
%       periodic: Logical variable: true if there are periodic BCs (optional)
%   
% OUTPUTS:
%       NodeP:  Reference node numbering matrix for pressure
%       NodeU:  Reference node numbering matrix for u-velocity 
%       NodeV:  Reference node numbering matrix for v-velocity 
%       idsP:   Location of active Pressure cells in interior
%       idsu:   Location of active u-velocity cells in interior
%       idsv:   Location of active v-velocity cells in interior
%
% See also: maskuv.m
%
% Written by:   Matt Ueckermann


%% Construct node numbering matrices with boundaries at four sides added

% Get mask size
[Ny, Nx] = size(Mask);

% Get number of different internal bcs
bndnum = max(Mask(:));
Nbcs = max(Mask(:)) + sum(bnd);

% Create u, v masks from P mask.
[Masku, Maskv] = maskuv(Mask);

% Find internal active nodes id's.
idsP = find(Mask == 0);
NodeP = Mask;
idsu = find(Masku == 0);
Nodeu = Masku;
idsv = find(Maskv == 0);
Nodev = Maskv;

% Number only the internal active nodes, starting from #bcs+1
NodeP(idsP) = (1:length(idsP)) + Nbcs;
Nodeu(idsu) = (1:length(idsu)) + Nbcs;
Nodev(idsv) = (1:length(idsv)) + Nbcs;

% Add boundaries
% Left
if bnd(1)
    bndnum = bndnum + 1;
    NodeP = [bndnum*ones(1,Ny  ); NodeP']';
    Nodeu = [bndnum*ones(1,Ny  ); Nodeu']';
    Nodev = [bndnum*ones(1,Ny-1); Nodev']';
    Nx = Nx + 1;
end;
% Right
if bnd(2)
    bndnum = bndnum + 1;
    NodeP = [NodeP'; bndnum*ones(1,Ny  )]';
    Nodeu = [Nodeu'; bndnum*ones(1,Ny  )]';
    Nodev = [Nodev'; bndnum*ones(1,Ny-1)]';
    Nx = Nx + 1;
end;
% Bottom
if bnd(3)
    bndnum = bndnum + 1;
    NodeP = [NodeP; bndnum*ones(1,Nx  )];
    Nodeu = [Nodeu; bndnum*ones(1,Nx-1)];
    Nodev = [Nodev; bndnum*ones(1,Nx  )];
end;
% Top
if bnd(4)
    bndnum = bndnum + 1;
    NodeP = [bndnum*ones(1,Nx  ); NodeP];
    Nodeu = [bndnum*ones(1,Nx-1); Nodeu];
    Nodev = [bndnum*ones(1,Nx  ); Nodev];
end;
Nx = Nx - sum(bnd(1:2));


%% Check for periodicity

% Add a boundary to handle periodicity correctly.
if (nargin == 3)&&(periodic)
    if  size(NodeP,1) == Ny
        NodeP = [NodeP(end,:) ; NodeP ; NodeP(1,:)];
        Nodeu = [Nodeu(end,:) ; Nodeu ; Nodeu(1,:)];
        Nodev = [Nodev(end,:) ; Nodev ; Nodev(1,:)];
    end;
    if size(NodeP,2) == Nx
        NodeP = [NodeP(:,end)' ; NodeP' ; NodeP(:,1)']';
        Nodeu = [Nodeu(:,end)' ; Nodeu' ; Nodeu(:,1)']';
        Nodev = [Nodev(:,end)' ; Nodev' ; Nodev(:,1)']';
    end;
end;
