function [S,M,Mlump,T,Dir,Free]=FEMprepare2wave(T)

% [S,M,Mlump,T,Dir,Free]=FEMprepare2wave(T)
%
% Input:
%     T    : basic triangulation
% Output:
%     S    : stiffness matrix (sparse)
%     M    : mass matrix (sparse)
%     Mlump: lumped mass matrix (vector)
%     T    : expanded triangulation (.baryc, .midpointNeu,....)
%     Dir  : Dirichlet nodes
%     Free : Free (non-Dirichlet) nodes

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);
Mlump=sum(M);

Dir=unique(T.dirichlet);
Free=(1:nC)'; Free(Dir)=[];

return