function b=loadplustraction(T,f,g,t,Free)

% b=loadplustraction(T,f,g,t)
% 
% Input:
%    T    : expanded triangulation
%    f    : f(x,y,t) vectorized
%    g    : g(x,y,t) vector-valued & vectorized
%    t    : current time
%    Free : list of Free nodes
% Output:
%    b    : load + traction vector

nC=size(T.coordinates,1);

% Load vector

f=(f(T.baryc(:,1),T.baryc(:,2),t).*T.detB)/6;
b=accumarray(T.elements(:),[f;f;f],[nC,1]);

% Traction vector

if length(T.neumann)>0
    g=0.5*sum(g(T.midptNeu(:,1),T.midptNeu(:,2),t).*T.normal,2);
    b2=accumarray(T.neumann(:),[g;g],[nC,1]);
else
    b2=zeros(nC,1);
end

b=b+b2;
b=b(Free);

return