function [b] = assembRightSide(mesh, f)

NuDoF = length(mesh.P);

if length(f) == 1
    temp=f;
    f = ones(NuDoF,1);
    f = temp*f;
    clear 'temp';
elseif size(f,1)<size(f,2)
  f=f';
end

b=zeros(NuDoF, 1); 
S4= [2, 1, 1; 1, 2, 1; 1, 1, 2]/24;
for i = 1:length(mesh.T) 
   mydet = abs(mesh.Bdet(i));
   b(mesh.T(i,:)) =  b(mesh.T(i,:)) + mydet * S4*f(mesh.T(i,:));
end 
%b = assembLoop(mesh.Bdet, mesh.T,f);
