function [iel] = findloc(xp, yp, xy, tri2nod)

% This routine does a brute force search through all elements until 
% it finds the triangle that contains the point, (xp, yp).  It
% returns the element number if successful.  It returns 0 if it
% does not find an element containing it.

[d, Ne] = size(tri2nod);

% Loop over elements until element is found.  The check 
% is based on taking the cross-product of each edge with
% the vector from the first node of the edge to the desired
% point. If the point is outside the triangle, at least on of the
% three cross products will be negative.  This assumes that the
% nodes are number counterclockwise.

xyp = [xp; yp];

for i = 1:Ne,

  xy1  = xy(:,tri2nod(1,i));
  xy2  = xy(:,tri2nod(2,i));
  xy3  = xy(:,tri2nod(3,i));

  xyp1 = xyp - xy1;
  xyp2 = xyp - xy2;
  xyp3 = xyp - xy3;
  
  xy21 = xy2 - xy1;
  xy32 = xy3 - xy2;
  xy13 = xy1 - xy3;
  
  c1 = cross([xy21; 0],[xyp1; 0]);
  c2 = cross([xy32; 0],[xyp2; 0]);
  c3 = cross([xy13; 0],[xyp3; 0]);

  if ( (c1(3) >= 0) & (c2(3) >= 0) & (c3(3) >= 0) ),
    iel = i;
    return;
  end
end

% If it get's here, no element was found which contains the point.

iel = 0;
return;