source: issm/trunk-jpl/src/m/mesh/NodeInElement.m@ 13730

function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity)
% NODEINELEMENT - find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
%
%    Usage:
%      node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
%
%  See also Nodeconnectivity
%
epsilon=10^-10;

%compute some quantities that will speed up the process
x3x1=x(elements(:,1))-x(elements(:,3));
y3y1=y(elements(:,1))-y(elements(:,3));
x3x2=x(elements(:,2))-x(elements(:,3));
y3y2=y(elements(:,2))-y(elements(:,3));
x3=x(elements(:,3));
y3=y(elements(:,3));
delta=x(elements(:,2)).*y(elements(:,3))-y(elements(:,2)).*x(elements(:,3))-x(elements(:,1)).*y(elements(:,3))+y(elements(:,1)).*x(elements(:,3))+x(elements(:,1)).*y(elements(:,2))-y(elements(:,1)).*x(elements(:,2));

%max connectivity:
max_connectivity=max(nodeconnectivity(:,end));
node_in_element=zeros(length(newx),max_connectivity+1); %last column is the number of elements to which the row node is connected.

for i=1:length(newx),
        x0=newx(i);
        y0=newy(i);

        %first area coordinate
        area_1=(y3y2.*(x0-x3)-x3x2.*(y0-y3))./delta;
        %second area coordinate
        area_2=(x3x1.*(y0-y3)-y3y1.*(x0-x3))./delta;
        %third area coordinate
        area_3=1-area_1-area_2;

        %get elements for which all area coordinates are positive (meaning (x0,y0) belongs to these elements
        pos=find((area_1>=0-epsilon) & (area_2>=0-epsilon) & (area_3>=0-epsilon));

        num_elements=length(pos);

        node_in_element(i,1:num_elements)=pos;
        node_in_element(i,end)=num_elements;
end