function A = assembStreamlineDiffusion (k, mesh, deltaT)

S1 =0.5* [ 1 -1 0; -1  1 0; 0 0 0];
S2 =0.5* [ 2 -1 -1; -1  0  1; -1  1  0];
S3 =0.5* [ 1  0 -1;  0  0  0; -1  0  1];

NuDoF = length(mesh.P);
A = sparse(NuDoF,NuDoF);

if  length(k) == NuDoF
  temp1 = (k(1,mesh.T(:,1)) + k(1,mesh.T(:,2)) + k(1,mesh.T(:,3)))./3;
  temp2 = (k(2,mesh.T(:,1)) + k(2,mesh.T(:,2)) + k(2,mesh.T(:,3)))./3;
  clear 'k';
  k(1,:) = temp1;
  k(2,:) = temp2;
end

term1 = squeeze(mesh.B(2,2,:))'.*k(1,:) - squeeze(mesh.B(1,2,:))'.*k(2,:);
term2 = squeeze(mesh.B(1,1,:))'.*k(2,:) - squeeze(mesh.B(2,1,:))'.*k(1,:);
factor = deltaT./abs(mesh.Bdet)';
gamma1 = factor.*term1.^2;
gamma2 = factor.*term1.*term2;
gamma3 = factor.*term2.^2;

for x=1:3
  for y=1:3
    temp = gamma1*S1(x,y)+gamma2*S2(x,y)+gamma3*S3(x,y);
    tempMatrix = sparse(mesh.T(:,x), mesh.T(:,y),temp,NuDoF,NuDoF);
    A = A + tempMatrix;
  end
end