function [elMat,elVec]=ElementContribution(corners,aFunc,bFunc,fFunc)
% [...] =  ElementContribution(...)
%   compute the element stiffness matrix of one element
%
%[elMat,elVec]=ElementContribution(corners,aFunc,bFunc,fFunc)
%  corners  coordinates of the three corners forming the triangular element
%  aFunc bFunc fFunc  function files or scalars for the coefficient functions
%  elMat element stiffness matrix
%  elVec vector contribution to the RHS of the linear equation
 
elMat=zeros(3);
elVec=zeros(3,1);
if is_scalar(aFunc) aV=aFunc*ones(3,1); 
               else aV=feval(aFunc,corners);endif
if is_scalar(bFunc) bV=bFunc*ones(3,1); 
               else bV=feval(bFunc,corners);endif
if is_scalar(fFunc) fV=fFunc*ones(3,1); 
               else fV=feval(fFunc,corners);endif
area=((corners(2,1)-corners(1,1))*(corners(3,2)-corners(1,2))-...
      (corners(2,2)-corners(1,2))*(corners(3,1)-corners(1,1))  )/2;
M=[corners(3,2)-corners(2,2),corners(1,2)-corners(3,2),corners(2,2)-corners(1,2);...
   corners(3,1)-corners(2,1),corners(1,1)-corners(3,1),corners(2,1)-corners(1,1)];

elMat=sum(aV)/(12*area)*M'*M+area/3*diag(bV);
elVec=area/3*fV;
endfunction
