function [S, M] = assemble(p, t)
[weight,point]=quad_3;  % get 3 point quadrature rule
nq=length(weight);
yloc=zeros(3,nq);
gyloc=zeros(2,3,nq);
for r=1:nq
  [yloc(:,r),gyloc(:,:,r)]=phiRef(point(:,r));
end
nt = length(t); nv=length(p);
S=sparse(nv,nv); M=sparse(nv,nv);
for it=1:nt
    indices=t(1:3,it)';
    % The coordinates of the vertices of triangle 'it'
    v=p(:,t(1:3,it));    
    B=[v(:,2)-v(:,1),v(:,3)-v(:,1)]; detB=abs(det(B));
	for r=1:nq
        gphi(:,:,r)=(inv(B))'*gyloc(:,:,r);
	end
    % Stiffness matrix
    Sloc=zeros(3,3);
    for r=1:nq
        Sloc=Sloc+((gphi(:,:,r)'*gphi(:,:,r)))*weight(r);
    end
    Sloc=Sloc*detB/2;
    S(indices, indices) = S(indices, indices)+Sloc;
    % Mass matrix
    Mloc=zeros(3,3);
    for r=1:nq
        Mloc=Mloc+((yloc(:,r)*yloc(:,r)'))*weight(r);
    end
    Mloc=Mloc*detB/2;
    M(indices, indices) = M(indices, indices)+Mloc;
end