function [bndelenew] = findbndy(meshele,meshnod,plottingon)
% [bndelenew] = findbndy(meshele,meshnod);
% [bndelenew] = findbndy(meshele,meshnod,plottingon);
% This finds the boundary elements of meshele
%  I.E.  all segments with only one element
disp(['mesh/meshtools/findbndy:  Finding the boundary, Watch for red triangles'])

if nargin==2
plottingon = 0;
else
plottingon = 1;
end

%  Make sure the smallest index comes first for all meshelements
for i = 1:length(meshele)
if (meshele(i,2)<meshele(i,3) & meshele(i,2)<meshele(i,4) )
%  Correct order
elseif (meshele(i,3)<meshele(i,2) & meshele(i,3)<meshele(i,4) )
me2 = meshele(i,2);
meshele(i,2) = meshele(i,3);
meshele(i,3) = meshele(i,4);
meshele(i,4) = me2;
elseif (meshele(i,4)<meshele(i,3) & meshele(i,4)<meshele(i,2) )
me2 = meshele(i,2);
meshele(i,2) = meshele(i,4);
meshele(i,4) = meshele(i,3);
meshele(i,3) = me2;
end
end

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

[m,n]=size(polys);
c = ones(4*m,1);
c(:) = polys;
l = min(c):max(c);    % this will be harder    
mc = min(c)-1 ;

xy=meshnod(l,2:3);
i1 = polys(:,1)-mc;
j1 = polys(:,2)-mc;
i2 = polys(:,2)-mc;
j2 = polys(:,3)-mc;
i3 = polys(:,3)-mc;
j3 = polys(:,4)-mc;
i4 = polys(:,4)-mc;
j4 = polys(:,1)-mc;
S = sparse([i1;i2;i3],[j1;j2;j3],1);
%spy(S,'g')

St = S';
Sd = S-St;
%hold on
%spy(Sd,'r')
[m,n]=size(Sd);

disp(['start loop on non zeros'])

[i,j,junk]=find(Sd==1);
bndele=[ [1:length(i)]' i j];

if plottingon == 1
figure(2)
clf
gplot(S,xy,'w')
hold on
plot([meshnod(i,2) meshnod(j,2)]',[meshnod(i,3) meshnod(j,3)]','g')
hold off
axis('equal')
orient landscape
grid
end



bndelenew = bndele;
n = length(bndele);
% ld = find(meshnod(:,3)<10000);
l =  find(meshnod(:,2)==min(meshnod(:,2))) ;
%start with the smallest x 

if length(l)>1
l = l(1);
end

ilast = find(bndele(:,2)==l);
%ilast = 1
bndelenew(1,:)= bndele(ilast,:);

N = n
for j = 2:n
l=find(bndele(:,2)==bndele(ilast,3));

if length(l)==1
bndelenew(j,:)= bndele(l,:);
bndele(ilast,:) = -bndele(ilast,:);
ilast = l;
else
N = N-1;
    if plottingon == 1
    disp(['length(l)~=1  ',int2str(ilast),'  ', int2str(bndele(ilast,3))...
     ,'  N =',int2str(N)])
    end
end

end
bndelenew = bndelenew(1:N,:);
bndelenew(:,1) = [1:N]';
n = N;

if plottingon == 1


%disp('size bndelenew')
%size(bndelenew)

figure(3)
clf
ax = [ -2000 2000 21500 23500 ; 6000 8000 19500 21500; ...
  -2000 2000 0 4000;  18000 22000 20000 23000];
plot([meshnod(bndelenew(:,2),2) meshnod(bndelenew(:,3),2)]',...
  [meshnod(bndelenew(:,2),3) meshnod(bndelenew(:,3),3)]','g')

for i = 1:n
       text(meshnod(bndelenew(i,2),2),meshnod(bndelenew(i,2),3),[int2str(i)])
end
title('bndele numbers     Use zoom to examine')
axis('equal')
grid

figure(4)
clf
plot([meshnod(bndelenew(:,2),2) meshnod(bndelenew(:,3),2)]',...
  [meshnod(bndelenew(:,2),3) meshnod(bndelenew(:,3),3)]','g')

for i = 1:n
        text(meshnod(bndelenew(i,2),2),meshnod(bndelenew(i,2),3),...
         [int2str(bndelenew(i,2))])
end
title('meshnod numbers     Use zoom to examine')
axis('equal')
grid
end    % if plottingon == 1



bndelenew = [bndelenew zeros(n,1) ones(n,1)];
% new one needs to be changed
% Set shoreline to value of  one
% Set open water to value of five
%  Do clever find statements to get indexs of ocean in lopen
l1 = find(meshnod(bndelenew(:,2),3) ==min(meshnod(bndelenew(:,2),3)));
l2 = find(meshnod(bndelenew(:,2),2) ==max(meshnod(bndelenew(:,2),2)));
lopen = l1:l2;
bndelenew(lopen,5)=ones(length(lopen),1)*5;    % open sea


lplus = find(bndele(:,1)>0);
 if plottingon == 1 | (length(lplus) > 1)
figure(5)
l = find(bndelenew(:,5)==5);
nl = find(bndelenew(:,5)~=5);
size(bndelenew)
size(meshnod)
if ~isempty(l)
  plot([meshnod(bndelenew(l,2),2) meshnod(bndelenew(l,3),2)]',...
  [meshnod(bndelenew(l,2),3) meshnod(bndelenew(l,3),3)]','b')
hold on
end
if ~isempty(nl)
  plot([meshnod(bndelenew(nl,2),2) meshnod(bndelenew(nl,3),2)]',...
  [meshnod(bndelenew(nl,2),3) meshnod(bndelenew(nl,3),3)]','g')
hold off
end

% Bad ones from bndele
lplus = find(bndele(:,1)>0);
if length(lplus) > 1
  bndele(lplus,:)
hold on
  plot([meshnod(bndele(lplus,2),2) meshnod(bndele(lplus,3),2)]',...
  [meshnod(bndele(lplus,2),3) meshnod(bndele(lplus,3),3)]','r')
  plot([meshnod(bndele(lplus,2),2) meshnod(bndele(lplus,3),2)]',...
  [meshnod(bndele(lplus,2),3) meshnod(bndele(lplus,3),3)]','r*')
title(' Interior triangles in Red','Color',[1 0 0])
hold off
end

end %  if plottingon == 1 figure 5
