function A = assembConvection (k, mesh)

gradN(1,:) = [-1, -1]; gradN(2,:) = [1, 0]'; gradN(3,:) = [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

c = (abs(mesh.Bdet)./mesh.Bdet)/6;
for y=1:3
  MVProduct1=  squeeze(mesh.B(2,2,:))*gradN(y,1)-squeeze(mesh.B(2,1,:))*gradN(y,2);
  MVProduct2= -squeeze(mesh.B(1,2,:))*gradN(y,1)+squeeze(mesh.B(1,1,:))*gradN(y,2);
  K =  c.*(k(1,:)'.*MVProduct1 + k(2,:)'.*MVProduct2);
  for x=1:3
    tempMatrix = sparse(mesh.T(:,x), mesh.T(:,y),K,NuDoF,NuDoF);
    A = A + tempMatrix;
  end
end