function  f = gm_neumann(vtx, facetlist, neuval)
%   f = gm_neumann(vtx, facetlist, neuval)
% Computes the Neumann BC contributions for one face of the brep.
% vtx is a vertex list (zba).
%   f is the return value (integrated Neumann vals).  Entries not
%     on the face in question are set to zeros in f (zba).
%   facetlist is the list of simplex facets on which to compute the function
%   neuval is the nodal evaluations of the Neumann BC's (zba)
%     at nodes
%  This function is called by gmfem.

[numnodes,di] = size(vtx);
[num_fac, scrap] = size(facetlist);

% find the areas of the facets.
if di == 2
  % If the problem is 2D, the boundary facets are 1-d and their
  % areas are their lengths.
  areas = sqrt((vtx(facetlist(:,1),0) - vtx(facetlist(:,0),0)) .^2 + ...
     (vtx(facetlist(:,1),1) - vtx(facetlist(:,0),1)) .^2);
else
  % If the problem is 3D, the boundary facets are 2-D and their
  % area is half the crossproduct of the side vectors.
  v1 = [(vtx(facetlist(:,1),0)-vtx(facetlist(:,0),0)), ...
        (vtx(facetlist(:,1),1)-vtx(facetlist(:,0),1)), ...
        (vtx(facetlist(:,1),2)-vtx(facetlist(:,0),2))];
  v2 = [(vtx(facetlist(:,2),0)-vtx(facetlist(:,0),0)), ...
        (vtx(facetlist(:,2),1)-vtx(facetlist(:,0),1)), ...
        (vtx(facetlist(:,2),2)-vtx(facetlist(:,0),2))];
  areas = sqrt((v2(:,2).*v1(:,1)-v2(:,1).*v1(:,2)).^2 + ...
    (v2(:,2).*v1(:,0)-v2(:,0).*v1(:,2)).^2 + ...            
    (v2(:,1).*v1(:,0)-v2(:,0).*v1(:,1)).^2);
end  


is_on_face = zba(zeros(numnodes,1));
for i = 0 : di - 1
  is_on_face(facetlist(:,i)) = ones(num_fac,1);
end

% if any(is_on_face1 ~= is_on_face)
%   error('Some problem determining whether a vertex is on a face')
% end

% Make an index vector to renumber vertices on this face so that they
% are consecutive 0,1,2,3,4,..

renum = cumsum(is_on_face) - ones(numnodes,1);
renuminv = find(is_on_face);
num_on_face = double(sum(is_on_face));

% Now renumber facetlist so that vertices are consecutive 0,1,2,3,...
rfacetlist = reshape(renum(facetlist),num_fac,di);

% Make a sparse matrix U whose columns we sum to get f.
% There is one row of this matrix per facet, and one
% column per vertex on the face.

is = kron(1:num_fac, ones(1,di));
js = double(rfacetlist + ones(num_fac,di))';
js = js(:)';
ss = zeros(1,di * num_fac);

if di == 2
  ss(1 : 2 : 2 * num_fac) = double((neuval(rfacetlist(:,0)) / 3 + ...
        neuval(rfacetlist(:,1)) / 6) .* areas)';
  ss(2 : 2 : 2 * num_fac) = double((neuval(rfacetlist(:,1)) / 3 + ...
        neuval(rfacetlist(:,0)) / 6) .* areas)';
else
  ss(1 : 3 : 3 * num_fac) = double((neuval(rfacetlist(:,0)) / 12 + ...
        neuval(rfacetlist(:,1)) / 24 +...
        neuval(rfacetlist(:,2)) / 24) .* areas)';

  ss(2 : 3 : 3 * num_fac) = double((neuval(rfacetlist(:,1)) / 12 + ...
        neuval(rfacetlist(:,0)) / 24 +...
        neuval(rfacetlist(:,2)) / 24) .* areas)';

  ss(3 : 3 : 3 * num_fac) = double((neuval(rfacetlist(:,2)) / 12 + ...
        neuval(rfacetlist(:,1)) / 24 +...
        neuval(rfacetlist(:,0)) / 24) .* areas)';
end

U = sparse(is,js,ss, num_fac, num_on_face);


% An 'if' statement to deal with Matlab's irritating overloading of 'sum'
if num_fac == 1
  f0 = full(U);
else
  f0 = full(sum(U));
end


f = zba(zeros(numnodes,1));
f(renuminv) = f0;

% ------------------------------------------------------------------
% Copyright (c) 1998 by Cornell University.  All rights reserved.
% See the accompanying file 'Copyright' for authorship information,
% the terms of the license governing this software, and disclaimers
% concerning this software.
% ------------------------------------------------------------------
% This file is part of the QMG software.  
% Version 2.0 of QMG, release date RELDATE
% ------------------------------------------------------------------