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;