function [A,b] = assembleNeumann(x,e,DOF,bounde,boundDOF,bilinIntegrators,linIntegrators,bilinBoundIntegrators,linBoundIntegrators,coefficients)

% [A,b] =
% assembleNeumann(x,e,DOF,bounde,boundDOF,bilinIntegrators,linIntegrators,bilinBoundIntegrators,linBoundIntegrators,coefficients)
%          assembles the Neumann finite-element system matrix and the right-hand side vector.
%
%          x is a nodal discretization
%          e is a list of nodal indices of elements
%          DOF is a list of global degrees of freedom relating to elements
%          bounde is a list of nodal indices of boundary or interface elements
%          boundDOF is a list of global degrees of freedom relating to boundary or interface elements
%          bilinIntegrators is a list of names of the bilinear integrators involved
%          linIntegrators is a list of names of the linear integrators involved
%          bilinBoundIntegrators is a list of names of the bilinear integrators involved
%          linBoundIntegrators is a list of names of the linear integrators involved
%          coefficients is a list of vectors of elementwise coefficient values
%            for the bilinear and linear integrators, respectively


% dimension
d = size(x,1);

% number of degrees of freedom
n = max(max(abs(DOF)));

% number of degrees of freedom per element
nlocal = size(DOF,1);

% number of elements
ne = size(e,2);

% number of DOFs per boundary element
nblocal = size(boundDOF,1);

% number of boundary or interface elements
nbe = size(bounde,2);

% numbers of integrators
nb = max(size(bilinIntegrators));
nl = max(size(linIntegrators));

% numbers of boundary or interface integrators
nbbound = max(size(bilinBoundIntegrators));
nlbound = max(size(linBoundIntegrators));

% allocation of the resulting matrix data
matrixRowIdx = zeros(ne*nlocal*nlocal+nbe*nblocal*nblocal,1);
matrixColumnIdx = zeros(ne*nlocal*nlocal+nbe*nblocal*nblocal,1);
matrixData = zeros(ne*nlocal*nlocal+nbe*nblocal*nblocal,1);
matrixIdx = 1;

b = zeros(n,1);

for i=1:ne
    
    if mod(i,10000)==0
        log = sprintf('assembleNeumann: %d out of %d elements',i,ne);
        disp(log)
    end

    ei = e(:,i);
    xe = x(:,ei);
        
    Ae = zeros(nlocal,nlocal);
    be = zeros(nlocal,1);
    
    for j=1:nb
        Ae = Ae + feval(str2func(bilinIntegrators{j}),coefficients{j}(:,i),xe);
    end
    
    for j=1:nl
        be = be + feval(str2func(linIntegrators{j}),coefficients{nb+j}(:,i),xe);
    end
    
    DOFe = DOF(:,i);
    
    for j=1:nlocal
        for k=1:nlocal
            matrixRowIdx(matrixIdx) = abs(DOFe(j));
            matrixColumnIdx(matrixIdx) = abs(DOFe(k));
            matrixData(matrixIdx) = sign(DOFe(j))*sign(DOFe(k))*Ae(j,k);
            matrixIdx = matrixIdx+1;
        end
        b(abs(DOFe(j))) = b(abs(DOFe(j))) + sign(DOFe(j))*be(j);
    end
    
end


for i=1:nbe
    
    if mod(i,10000)==0
        log = sprintf('assembleNeumann: %d bound. elements',i);
        disp(log)
    end

    ei = bounde(:,i);
    xe = x(:,ei);

    Ae = zeros(nblocal,nblocal);
    be = zeros(nblocal,1);
    
    for j=1:nbbound
        Ae = Ae + feval(str2func(bilinBoundIntegrators{j}),coefficients{nb+nl+j}(:,i),xe);
    end
    
    for j=1:nlbound
        be = be + feval(str2func(linBoundIntegrators{j}),coefficients{nb+nl+nbbound+j}(:,i),xe);
    end

    DOFe = boundDOF(:,i);
    
    for j=1:nblocal
        for k=1:nblocal
            matrixRowIdx(matrixIdx) = abs(DOFe(j));
            matrixColumnIdx(matrixIdx) = abs(DOFe(k));
            matrixData(matrixIdx) = sign(DOFe(j))*sign(DOFe(k))*Ae(j,k);
            matrixIdx = matrixIdx+1;
        end
        b(abs(DOFe(j))) = b(abs(DOFe(j))) + sign(DOFe(j))*be(j);
    end

end

A = sparse(matrixRowIdx,matrixColumnIdx,matrixData,n,n);
