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;