function [NodeP, nbcP, Nodeu, nbcu, Nodev, nbcv, NodeRho, nbcrho] = set_bcs(...
    bcid2type_P, NodeP, bcid2type_u, Nodeu, bcid2type_v, Nodev, bcid2type_Rho, NodeRho)
% function [NodeP, nbcP, Nodeu, nbcu,  Nodev, nbcv, NodeRho, nbcrho] =...
%     set_bcs(bcid2type_P, NodeP,...
%     bcid2type_u, Nodeu, bcid2type_v, Nodev, bcid2type_Rho, NodeRho)
%
% This function re-numbers the Node-matrices to be consistent with the
% solvers boundary-condition rules. The different boundary conditions are:
% bcid2type = {'D'} - Dirichlet boundary with a single value on that
%                     segment
% bcid2type = {'Dmv'} - Dirichlet boundary with multiple values on that
%                       segment (spatially variable)
% bcid2type = {'O'} - Open boundary condition
% bcid2type = {'N'} - Neumann boundary with single value on that segment
% bcid2type = {'Nmv'} - Neumann boundary with multiple values on that
%                       segment (spatially variable)
%
% INPUTS:
%       bcid2type_P: Cell array containing the boundary condition types
%                      for Pressure
%       NodeP:       Node numbering matrix for pressure
%       bcid2type_u: Cell array containing the boundary condition types
%                      for u-velocity
%       Nodeu:       Node numbering matrix for u-velocity
%       bcid2type_v: Cell array containing the boundary condition types
%                      for v-velocity
%       Nodev:       Node numbering matrix for v-velocity
%       bcid2type_rho: Cell array containing the boundary condition types
%                      for density (Optional)
%       NodeRho:       Node numbering matrix for density (Optional)
% OUTPUTS:
%       NodeP:    Re-numbered Node numbering matrix for pressure
%       nbcP:     Array giving maximum id of each pressure boundary type.
%                  Used as an input to init_bcs.m
%       Nodeu:    Re-numbered Node numbering matrix for u-velocity
%       nbcu:     Array giving maximum id of each u-velocity boundary type
%                  Used as an input to init_bcs.m
%       Nodev:    Re-numbered Node numbering matrix for v-velocity
%       nbcv:     Array giving maximum id of each v-velocity boundary type
%                  Used as an input to init_bcs.m
%       NodeRho:  Re-numbered Node numbering matrix for density (Optional)
%       nbcrho:   Array giving maximum id of each density boundary type
%                  Used as an input to init_bcs.m (Optional)
%
% See also: translate_Node.m, bc_maxrule.m, init_bcs.m
%
% Written by:   Matt Ueckermann

%%

Nbcids =  length(bcid2type_P);
if (nargin < 7)
    NodeRho = NodeP;
    bcid2type_Rho = bcid2type_P;
end;

%% First we translate the bcid2type arrays to some nice numbers

bc_u = zeros(length(bcid2type_u), 1);
bc_v = zeros(length(bcid2type_v), 1);
bc_P = zeros(length(bcid2type_P), 1);
bc_rho = zeros(length(bcid2type_Rho), 1);

for i=1 : length(bcid2type_u)
    bc_u(i) = bc_str2num(bcid2type_u{i});
    bc_v(i) = bc_str2num(bcid2type_v{i});
    bc_P(i) = bc_str2num(bcid2type_P{i});
    bc_rho(i) = bc_str2num(bcid2type_Rho{i});
end

%% Do the translation to follow the order rules

[Nodeu, maxbcid_u, nbcu] = translate_Node(Nodeu, bc_u, Nbcids);
[Nodev, maxbcid_v, nbcv] = translate_Node(Nodev, bc_v, Nbcids);
[NodeP, maxbcid_P, nbcP] = translate_Node(NodeP, bc_P, Nbcids);
[NodeRho, maxbcid_rho, nbcrho] = translate_Node(NodeRho, bc_rho, Nbcids);
maxbcid = max([maxbcid_u, maxbcid_v, maxbcid_P, maxbcid_rho]);

%% Do the translation to follow the same number of BCs rule

Nodeu = bc_maxrule(Nodeu, maxbcid_u, maxbcid);
Nodev = bc_maxrule(Nodev, maxbcid_v, maxbcid);
NodeP = bc_maxrule(NodeP, maxbcid_P, maxbcid);
NodeRho = bc_maxrule(NodeRho, maxbcid_rho, maxbcid);
nbcP(3) = maxbcid;
nbcu(3) = maxbcid;
nbcv(3) = maxbcid;
nbcrho(3) = maxbcid;

