function [a, b, c] = elefunc(meshele,meshnod)
% [a, b, c] = elefunc(meshele,meshnod);
% Returns the element functions for meshnod and meshele
% a, b, c are ( length(meshele) , 3) arrays
% f[iele,i](x,y) = a(iele,i)*x + b(iele,i)*y + c(iele,i)
% f[iele,i](x(j),y(j)) = 1 i==j,  =0 i~=j
%  
% An interpolation inside the triangle (Z(1:3) at nodes):
%            e.g.  Z = bathy(meshele(iele,2:4));
% Zinner(x,y) = Z(:).*f[iele,:](x,y)
% x,y scalars:
% Zinner = (a(iele,:)*x + b(iele,:)*y + c(iele,:))*Z(:); 
% x,y vectors, :
% Zinner = (x(:)*a(iele,:) + y(:)*b(iele,:) + ones(length(x),1)*c(iele,:)...
%             )*Z(:); 



% Reserve space for a, b, c
a = meshele(:,1:3);
b = a;
c = a;

for iele= 1:length(meshele)
  x = meshnod(meshele(iele,2:4),2);
  y = meshnod(meshele(iele,2:4),3);
  D = [x y ones(3,1)];
  d= det(D);
  for i = 1:3
    a(iele,i) = (y(2)-y(3))/d;
    b(iele,i) = (x(3)-x(2))/d;
    c(iele,i) = (x(2)*y(3)-x(3)*y(2))/d;
    % rotate the x,y for next function
    x=x([2 3 1]);
    y=y([2 3 1]);
  end
end



%  Force to columns
%x = x(:);
%y = y(:);
%for i = 1:3
%f = a(iele,i)*x + b(iele,i)*y + c(iele,i)*ones(length(x),1);
%end
