function [elems2,nodes2,bathy2,meshcd2,bndele2] = ...
     reorder(elems,nodes,bathy,meshcd,bndele)
%  [meshele2,meshnod2,bathy2,meshcd2,bndele2] = ...
%     reorder(meshele,meshnod,bathy,meshcd,bndele);
% Reorders the mesh to achieve a narrow bandwidth of S
%  Same as ordrmesh.m, but it includes meshcd 
sb = size(bndele)
if sb(1) ==0
bndele = ones(3,3);
end

polys = elems;    
polys(:,1)=elems(:,4);   %   triangles close on themselves

[m,n]=size(polys);
c = ones(4*m,1);
c(:) = polys;

xy=nodes(:,2:3);
i1 = polys(:,1);
j1 = polys(:,2);
i2 = polys(:,2);
j2 = polys(:,3);
i3 = polys(:,3);
j3 = polys(:,4);
i4 = polys(:,4);
j4 = polys(:,1);
S = sparse([i1;i2;i3],[j1;j2;j3],1);

figure(4)

Sr=symrcm(S);
%spy(S(Sr,Sr)) 
size(Sr) 
R = randperm(length(S));
SR = symrcm(S(R,R));
Sr = R(SR);
% spy(S(Sr,Sr))

l = zeros(length(nodes),1);
for i = 1:length(nodes)
v = Sr(i);
l(v) = i;
end
%Change values inside these arrays to new node numbers
elems2=elems;
elems2(:,2) = l(elems(:,2));
elems2(:,3) = l(elems(:,3));
elems2(:,4) = l(elems(:,4));
[ju,l]=sort(elems2(:,2));
elems2=elems2(l,:);
bndele2=bndele;
bndele2(:,2) = l(bndele(:,2));
bndele2(:,3) = l(bndele(:,3));

%Change the order of these arrays 
nodes2=nodes;
nodes2(:,2) = nodes(Sr,2);
nodes2(:,3) = nodes(Sr,3);
bathy2=bathy;
bathy2(:,2) = bathy(Sr,2);
meshcd2=meshcd;
meshcd2(:,2) = meshcd(Sr,2);


 
%   Plot before and after pictures
%polys = elems2;    
%polys(:,1)=elems2(:,4);   %   triangles close on themselves

%xy=nodes2(:,2:3);

%X=polys;
%X(:)=xy(polys,1);
%Y=polys;
%Y(:)=xy(polys,2);
%X=X';
%Y=Y';

%C = [ 1:length(X)]/length(X);

%fill(X,Y,C)
%title('Newly reordered')
%figure(gcf)
