rhs_neumann.m

function res = rhs_neumann(x, model, df_info, i, elind, edgeind)
%function res = rhs_neumann(x, model, df_info, i, elind, edgeind)
%

% IM 21.07.2016


res = zeros(length(elind), 1);

for local_edge_ind = 1:3
    
    lcoord = llocal2local(df_info.grid, local_edge_ind, x);
    
    hat_phi_i = evaluate_basis_function(df_info, lcoord, i);
    
    inds = (edgeind == local_edge_ind);
    
    n = cache_function(@(xelind, xlocal_edge_ind, xarg, xmodel)...
                       xmodel.neumann_values(...
                           df_info.grid, xelind, xlocal_edge_ind, ...
                           xarg, xmodel), ...
                       elind(inds), local_edge_ind, x, model);
    
    if ~iscell(n)
        
        res(inds) = n * hat_phi_i;
    else
        
        if ~iscell(res)
            
            res = cell(1, length(n));
            for q = 1:length(n)
                res{q} = zeros(length(elind), df_info.dimrange);
            end
        end
        
        for q = 1:length(n)
            
            res{q}(inds) = n{q} * hat_phi_i;
        end
    end
end

end