function [ A, b ] = assembPDE( eps, k, r, f, mesh)

A = assembDiffusion (eps, mesh.T, mesh.P);
if ~isempty(k)
  A = A + assembConvection (k, mesh);
end
if ~isempty(r)
  A = A + assembReaction (r, mesh.T, mesh.P);
end
b = assembRightSide(mesh, f);

if ~isempty(k)
  kNorm     = max(sqrt(k(1,:).^2 + k(2,:).^2));
  epsMin = min(eps);
  if (kNorm / epsMin) > 1
    if  length(k) == mesh.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
    
    deltaT = computeSDCoefficient(eps, k, mesh);
    b = b + assembRightSideStreamlineDiffusion(mesh, f, k, deltaT);
    A  = A + assembStreamlineDiffusion (k, mesh, deltaT);
    if ~isempty(r)
      A  = A + ( assembConvection ([deltaT .* r; deltaT .* r] .* k, mesh) )';
    end
  end
end
