function [Nodeu, maxbcid, bc] = translate_Node(Nodeu, bc_u, Nbcids)
%function [Nodeu, maxbcid, bc] = translate_Node(Nodeu, bc_u, Nbcids)

% We initially set all boundary IDs to negative numbers so as not to
% lose information.
nextid = -1;

% 1: Deal with multi-valued Dirichlet BCs ('dmv')
id2type = find(bc_u == 1);
ids = false(size(Nodeu));
for i=1:length(id2type)
    ids = ids | (Nodeu == id2type(i));
end;
Nodeu(ids) = nextid: -1:(nextid + 1 - sum(ids(:)));
nextid = nextid - sum(ids(:));

% 2: Deal with single-valued Dirichlet BCs ('d')
id2type = find(bc_u == 2);
for i=1:length(id2type)
    ids = (Nodeu == id2type(i));
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end;
bc(1) = -nextid - 1;

% 3: Deal with open BCs ('o')
id2type = find(bc_u == 3);
for i=1:length(id2type)
    ids = (Nodeu == id2type(i));
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end;
bc(2) = -nextid - 1;

% 4: Deal with multi-valued Neumann BCs ('nmv')
id2type = find(bc_u == 4);
ids = false(size(Nodeu));
for i=1:length(id2type)
    ids = ids | (Nodeu == id2type(i));
end;
Nodeu(ids) = nextid: -1:(nextid + 1 - sum(ids(:)));
nextid = nextid - sum(ids(:));

% 5: Deal with single-valued Neumann BCs ('n')
id2type = find(bc_u == 5);
for i=1:length(id2type)
    ids = (Nodeu == id2type(i));
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end;
bc(3) = -nextid - 1;

% Finally, we set the interior node IDs and set boundary IDs back to
% positive numbers.
ids = (Nodeu > 0);
Nodeu(ids) = Nodeu(ids) - Nbcids + (-nextid - 1);
Nodeu(~ids) = -Nodeu(~ids);

maxbcid = -nextid - 1;