%function  [isonbnd,index]=onboundary(indp,e)
%
% indp = indexes of points (=of columns of matrix 'p')
% boudary information matrix 'e'
% 
% Returns the maximum number of points of each column of indp 
% that lie on boundary.
%
% E.g. To test whether some edge of the triangle tr=[p0;p1;p2] 
% lies on boundary we say isonbnd(tr,e)==2. 
% 
% in 'index' we return the column indexes in 'e' matrix for corresponding 
% nodes in 't'; the nodes not on boundary have zero index.


function [isonbnd,index]=onboundary(indp,e)


isonbnd = zeros(1,size(indp,2)); 
index = sparse([],[],[],size(indp,1),size(indp,2),20); 

e1 = meshgrid(e(1,:),ones(size(indp,1),1));
e2 = meshgrid(e(2,:),ones(size(indp,1),1));

   % assume the edges are positively oriented!
 
for i=1:size(indp,2); 
     meshgrid(indp(:,i),ones(1,size(e,2)))';
     (ans == e1) + (ans == e2);
     [y,j] = max(sum(ans)); 
     index(:,i) = j*ans(:,j);
     isonbnd(i) = y; 
end



