function [ meshele]=gridbndy(meshnodold,bndele)
%  figure(1); [ meshele]=gridbndy(meshnod,bndele);
% Given a meshnod type array this will build meshele array
%  bndele is a boundary to be used to contain the new mesh
%  Make a big circle surrounding the data
R = 1*max( max(meshnodold(:,2))-min(meshnodold(:,2)), ...
           max(meshnodold(:,3))-min(meshnodold(:,3)) );
xo = mean(meshnodold(:,2)); yo = mean(meshnodold(:,2));
theta =0:(2*pi/20):(2*pi);
 noderingx = R*cos(theta)+xo;
noderingy = R*sin(theta)+yo;
[nring,m]=size(noderingx');
 ringnod = [ ones(nring,1) noderingx' noderingy' ];
meshnod= [ meshnodold; ringnod];

%plot(meshnodold(:,2),meshnodold(:,3))
%L = axis;
bathy = meshnod(:,1:2);
meshnodextra=meshnod;

meshnodextra(:,1) = [1:length(meshnodextra)]';
[meshele] = makeeles(meshnodextra);


n = length(meshnodold);
l = find(meshele(:,2)<=n & meshele(:,3)<=n ...
                                    & meshele(:,4)<=n );
meshele =meshele(l,:);


  [ax] = ploteles(meshele,meshnodold);
% Now eliminate elements which have boundary segments in the wrong direction

for i = 1:length(bndele)
l2 = find(meshele(:,2)==bndele(i,2) ) ;
l3 = find(meshele(:,3)==bndele(i,2) ) ;
l4 = find(meshele(:,4)==bndele(i,2) ) ;
lf = [ l2; l3; l4];
l2 = find(meshele(:,2)==bndele(i,3) ) ;
l3 = find(meshele(:,3)==bndele(i,3) ) ;
l4 = find(meshele(:,4)==bndele(i,3) ) ;
ls = [ l2; l3; l4]  ;  
l = intersec(lf,ls)  ;  % elements with both bndele(i,2) and bndele(i,3)

   if length(l)==2
%  This boundary segment has two elements, one is bad
%   disp([ones(2,1)*i l' meshele(l,:)])

    for j = 1:2
      lj = l(j);
      mele = meshele(lj,:); mele(1,1) = mele(1,4);
      l2 = find(mele(1:3)==bndele(i,2));
          if mele(l2+1)~=bndele(i,3)
%       backwards no good meshele(lj)
        meshele(lj,1) = meshele(lj,1)*NaN;
      end
    end
   end

end  

l = find(~isnan(meshele(:,1)));
meshele= meshele(l,:);
figure(gcf+1)
  [ax] = ploteles(meshele,meshnodold);
