function [mat,vector] = AdvectionTermUpwind(EdgeVelocities, I, J, S, NumInsideNodes, Sol)

% Advection term's spatial discretization using first-order upwind method

if nargin<6    % implicit TS assembles the coeff matrix 'mat' (LHS)
  vector = [];
  
  % EdgeVelocities : the edge velocity .dot-product. with the edge normal,
  % therefore, if it has a positive sign, the flux goes out of the CV
  % through the edge, otherwise the flux comes inside the CV.
  indI = EdgeVelocities > 0;
  indJ = ~indI;
  
  % for outgoing fluxes, assign the flux contribution to the column
  % pertained to the CV itself in the coefficient matrix 
  mat = sparse(I(indI), I(indI), EdgeVelocities(indI).*S(indI), NumInsideNodes,NumInsideNodes);
  
  % for incomming fluxes, assign the flux contribution to the column
  % corresponding the neghibor node in the coefficient matrix
  mat = mat + sparse(I(indJ), J(indJ), EdgeVelocities(indJ).*S(indJ), NumInsideNodes,NumInsideNodes);
  
else            % explicit TS assembles the 'vector = mat*sol' (RHS)
  mat = [];
  
  indI = EdgeVelocities > 0;
  indJ = ~indI;
  A = sparse(I(indI), I(indI), EdgeVelocities(indI).*S(indI), NumInsideNodes,NumInsideNodes);
  
  A = A + sparse(I(indJ), J(indJ), EdgeVelocities(indJ).*S(indJ), NumInsideNodes,NumInsideNodes);
  vector = -A*Sol;
  
end

end  % end of function 'AdvectionTermUpwind'