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;