function edgeVec=ElementContributionEdge(corners,gNFunc)
% [...] =  ElementContributionEdge(...)
%   compute the vector contribution due to one line segement of the Neumann boundary
%
% edgeVec=ElementContributionEdge(corners,gNFunc)
%  corners  coordinates of the two points forming the line segment
%  gNFunc  function file or scalar for the coefficient function
%  edgeVec vector contribution to the RHS of the linear equation

% this is a Gauss integration
w = 0.5/sqrt(3);
wm = 0.5 - w;
wp = 0.5 + w;

if is_scalar(gNFunc)
    g = gNFunc*ones(2,1);
else
    p1 = wp*corners(1,:) + wm*corners(2,:);
    p2 = wm*corners(1,:) + wp*corners(2,:);
    g = feval(gNFunc,[p1;p2]);
endif

edgeVec=0.5*norm(corners(2,:)-corners(1,:))*[wp*g(1)+wm*g(2);wm*g(1)+wp*g(2)];
endfunction
