function F = nedelec_f_fun(model, fun)
%NEDELEC_F_FUN Nedelec-type volume excitation vector for a known function
%   F = NEDELEC_F_FUN(MODEL, FUN) Assembles the excitation vector F for 

% Check input
if ~isa(fun, 'function_handle')
    error('nedelec_f:invalid_argument', ...
        'Input parameter FUN must be a function handle');
end

% Prepare elements and edges
Elements = drop_IDs(model);             % Drop the user IDs
edges = get_edges(Elements);            % Extract edges of the mesh
[dofs, dirs] = get_edge_dofs(edges);    % Assign DOFs to each edge
max_dof = max(dofs);                    % Maximum DOF number

% Possible element types
% Note: it is important to have the same order of types as in get_edges
types = [12 23 24 34 36 38].';      
numedges = [1 3 4 6 9 12].';        % Number of edges for each element type
% Find the type index for each element
[~, type_idx] = ismember(model.Elements(:,2), types);
% Number of edges for each element
nEdges = numedges(type_idx);

nDim = floor(Elements(:,2) / 10);       % number of dimensions of elements
if(length(unique(nDim))~=1)
    error('nedelec_f:inhomogeneous_mesh', ...
        'Model contains elements with different number of dimensions');
end

nDof = nEdges;
nVal = (nDof);                  % number of matrix entries per elem
Start = [0; cumsum(nVal)];      % starting indices of elements

N = Start(end); % total number of matrix entries (non overlapped)
%I  = zeros(N,1);  % preallocating row indices
%J  = zeros(N,1);  % preallocating column indices
F  = zeros(max_dof, 1);  % preallocating mass values

% Matrix assembly 
nElem = size(Elements,1);               % number of elements
iElem = 0;                              % 

i_edge = 0; 
for type = [12 23 24 34 36 38]
    typeind = Elements(:,2) == type;
    if ~any(typeind)
        continue;
    end
    dim = floor(type/10);               % dimension number
    nnodes = mod(type,10);
    elements = Elements(typeind,:);     % the elements of this type
    start = Start(typeind);             % start indices of elements
    [w, N, dN, L, dL] = elem_quad(type);       % Gaussian quadarture for this type
    [~,type_idx] = ismember(type,types);
    nEdges = numedges(type_idx);
    
    for ielem = 1 : size(elements,1)
        elem = elements(ielem,:);
        nodes = elem(4+(1:nnodes));             % elem nodes
        vert = model.Nodes(nodes,1+(1:dim)).';  % corner vertices
        % Assemble element matrix
        dir_idx = dirs(i_edge+(1:nEdges),:);
        
        f = elem_f(type, vert, w, N, dN, L, dL, dir_idx, fun);
        dof_idx = dofs(i_edge+(1:nEdges),:);
        dof_idx = reshape(dof_idx',[],1);
        
        i = dof_idx;         % row indices
        F(i) = F(i) + f;
        
        iElem = iElem+1;
        progbar(1, nElem, iElem);
        
        i_edge = i_edge + nEdges;
    end

end

end % of function nedelec_f_fun


function [w, N, dNxi, L, dLxi] = elem_quad(type)

% Compute Gaussian quadrature weights and shape function samples
nnod = mod(type, 10);

switch type
    case 12
        [xi, w] = gaussquad1(2);
    case {23, 24}
        [xi, w] = gaussquad2(2, nnod);
    case 34
        [xi, w] = gaussquad3(2, nnod);
    case {36, 38}
        [xi, w] = gaussquad3(3, nnod);
end

[L, dLxi] = shapefun(xi, type);
dLxi = reshape(shiftdim(dLxi,2),numel(dLxi)/nnod,nnod);

[N, dNxi] = nedelec_shapefun(xi, type);
dNxi = reshape(shiftdim(dNxi,2),numel(dNxi)/nnod,nnod);

end


function F = elem_f(type, xnode, w, N, dNxi, L, dLxi, dir, fun) %#ok<INUSL>

dim = floor(type/10);
ngauss = length(w);

x = L * xnode.';                % Physical coords of the quadrature points
f_x = fun(x);

% Coordinate transform
J = dLxi * xnode.';             % Jacobian of the coordiate transform

% Preallocate
Nx = zeros(size(N));            % Nedelec shape function transformed
j = zeros(ngauss,1);            % Jacobian determinant per quadrature point

for n = 1 : ngauss
    ind = dim*(n-1)+(1:dim);
    Jac = J(ind,:);                     % Jacobi transposed
    
    Nx(ind,:) = Jac \ N(ind,:);
    
    j(n) = abs(det(Jac));
    
end

% Vector assembly
q = kron(w.*j,ones(dim,1));
F = diag(dir)*Nx.'*diag(q)*reshape(f_x.', [], 1);

end