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

nextid = -1;
id2type = find(bc_u == 1);
ids = false(size(Nodeu));
for i=1:length(id2type)
    ids = ids + Nodeu == id2type(i);
end
%We first add all the boundary id's as negative numbers, such as to not
%lose information
Nodeu(ids) = nextid: -1:(nextid + 1 - sum(ids(:)));
nextid = nextid - sum(ids(:));

%Now the single-valued Dirichlet BC's
id2type = find(bc_u == 2);
for i=1:length(id2type)
    ids = ids + Nodeu == id2type(i);
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end
bc(1) = -nextid - 1;

%Now do the open BC's
id2type = find(bc_u == 3);
for i=1:length(id2type)
    ids = ids + Nodeu == id2type(i);
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end
bc(2) = -nextid - 1;

% Now the multi-valued Neumann Bcs
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(:));

%Finally, the single-valued Neumann Bcs
id2type = find(bc_u == 5);
for i=1:length(id2type)
    ids = ids + Nodeu == id2type(i);
    Nodeu(ids) = nextid;
    nextid = nextid - 1;
end

% Now we have to correct the interior -- pad them to be the correct numbers
ids = Nodeu > 0;
Nodeu(ids) = Nodeu(ids) - Nbcids - nextid - 1;
Nodeu(~ids) = - Nodeu(~ids);

maxbcid = -nextid - 1;
bc(3) = -nextid - 1;