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

% Advection term's spatial discretization using central differences

if nargin<6    % implicit TS assembles the coeff matrix 'mat' (LHS)
  vector = [];
  
  % neighbor contributions
  mat = 1/2*sparse(I, J, EdgeVelocities.*S, NumInsideNodes,NumInsideNodes);
  
  % main diagonal contributions
  mat = sparse(1:NumInsideNodes, 1:NumInsideNodes, sum(mat, 2), NumInsideNodes, NumInsideNodes) + mat;
  
else            % explicit TS assembles the 'vector = mat*sol' (RHS)
  mat = [];
  
  % neighbor contributions
  A = 1/2*sparse(I, J, EdgeVelocities.*S, NumInsideNodes,NumInsideNodes);
  
  % main diagonal contributions
  A = sparse(1:NumInsideNodes, 1:NumInsideNodes, -sum(A, 2), NumInsideNodes, NumInsideNodes) - A;
  vector = A*Sol;
  
end

end  % end of function 'AdvectionTerm'