function b = applyNeumannBoundaryConditions(mesh, pGroup, g, eps )

if isa(g, 'function_handle')
  x = mesh.P(1,:);
  y = mesh.P(2,:);
  temp = g(x,y);
  g = temp;
elseif length(g) == 1
  g = g*ones(1,mesh.NuDoF);
elseif length(g) ~= mesh.NuDoF
  error('g does not fit in size or type');
end
if length(eps) == 1
  temp=eps;
  eps = ones(mesh.NuTri,1);
  eps = temp*eps;
  clear 'temp';
elseif length(eps) == mesh.NuDoF
  temp = (eps(tri(:,1)) + eps(tri(:,2)) + eps(tri(:,3)))./3;
  eps = temp;
  clear 'temp';
end

b = zeros (mesh.NuDoF,1);

accessRow = mesh.BoundaryEdges(mesh.GBoundaryEdges == pGroup);
firstNode  = mesh.edges(mesh.BoundaryEdges(:),1);
secondNode = mesh.edges(mesh.BoundaryEdges(:),2);

h(mesh.BoundaryEdges(:)) = (mesh.P(1,firstNode)-mesh.P(1,secondNode) ).^2 + ( mesh.P(2,firstNode)-mesh.P(2,secondNode) ).^2;
h = sqrt(h);

for i=accessRow'
  myeps = eps(mesh.edge2triangle(i));
  firstNode  = mesh.edges(i,1);
  secondNode = mesh.edges(i,2);
  b(firstNode) = b(firstNode) + myeps/6*h(i)*(2*g(firstNode) + g(secondNode)); 
  b(secondNode) = b(secondNode) + myeps/6*h(i)*(2*g(secondNode) + g(firstNode)); 
end

