function [S,M,b1,b2]=FEMmatrices(T,f,g)

% [S,M,load,trct]=FEMmatrices(T,f,g)
%
% Input:
%     T    : basic triangulation
%     f    : vectorized function of two variables
%     g    : vectorized function of two variables
%            with values in R^2
% Output:
%     S    : stiffness matrix (sparse)
%     M    : mass matrix (sparse)
%     load : load vector (associated to f)
%     trct : traction vector (associated to g)

T=expandGrid(T);
nC=size(T.coordinates,1);
nE=size(T.elements,1);

% Matrices in the reference element

K11=0.5*[1 -1 0;-1 1 0;0 0 0];
K22=0.5*[1 0 -1;0 0 0;-1 0 1];
K12=0.5*[1 0 -1;-1 0 1;0 0 0];
M=1/24*[2 1 1;1 2 1;1 1 2];

% Assembly: 

[J,I]=meshgrid([1 2 3]);
I=T.elements(:,I)';I=reshape(I,3,3*nE);
J=T.elements(:,J)';J=reshape(J,3,3*nE);
S=kron(T.c11',K11)+kron(T.c22',K22)+kron(T.c12',K12+K12');
S=sparse(I,J,S);
M=kron(T.detB',M);
M=sparse(I,J,M);

% Load vector

f=(f(T.baryc(:,1),T.baryc(:,2)).*T.detB)/6;
b1=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.normal,2);
    b2=accumarray(T.neumann(:),[g;g],[nC,1]);
else
    b2=zeros(nC,1);
end

return