function [mesh, bw, bwo] = renumber_mesh(mesh)
%RENUMBER_MESH Renumber mesh to reduce finite element matrix bandwidths
%   [MESH, BW, BWO] = RENUMBER_MESH(MESH)
%       Renumbers and reorders the nodes of MESH in order to reduce the
%       bandwidth of the finite element matrices. The optional outputs BW
%       and BWO are the bandwidths of the new and the original mesh,
%       respectively.
%
%   Renumbering is performed based on the adjacency of the meah and the
%   reverse Cuthill-McKee algorithm
%
%   See also: adjacency, symrcm

node_IDs = mesh.Nodes(:, 1);
A = adjacency(mesh);
r = symrcm(A);
mesh.Nodes(:, 2:4) = mesh.Nodes(r, 2:4);

[~, re] = ismember((1 : size(mesh.Nodes, 1)).', r);
elements = drop_IDs(mesh);
connections = elements(:, 5:end);
sel = connections ~= 0;
connections(sel) = node_IDs(re(connections(sel)));
mesh.Elements(:, 5:end) = connections;

% Bandwidth calculation
if nargout > 1
    [i, j] = find(A(r,r));
    bw = max(abs(i-j)) + 1;
end

% Bandwidth of the original system
if nargout > 2
    [i, j] = find(A);
    bwo = max(abs(i-j)) + 1;
end


end % of function renumber_mesh
