function [A, DOF] = nedelec_a(model, fun)
%NEDELEC_A Loading vector from a Neumann BC using Nedelec finite elements
%   [A, DOF] = NEDELEC_A_FUN(MODEL, FUN) Computes the loading vector A 
%       of the finite element model defined by MODEL using linear Nedelec
%       shape functions.
%
% See also: nedelec_mk, nedelec_f_fun, nedelec_shapefun

% 2020

%% Get boundary-edges and their directions 

bou = get_boundary(model);

Elements = drop_IDs(bou); % drop IDs

edges = get_edges(model.Elements);

bedges = get_edges(bou.Elements);

[is1, where1] = ismember(bedges, edges, 'rows');
[is2, where2] = ismember(bedges, fliplr(edges), 'rows');

[dofs, dirs] = get_edge_dofs(edges);

bdofs = zeros(size(bedges, 1), 1);
bdirs = zeros(size(bedges, 1), 1);
bdofs(is1) = dofs(where1(is1));
bdirs(is1) = dirs(where1(is1));
bdofs(~is1 & is2) = dofs(where2(~is1 & is2));
bdirs(~is1 & is2) = -1*dirs(where2(~is1 & is2));


%% 

i_edge = 0;
types = [12 23 24]';
numedges = [1 3 4]';

[~,type_idx] = ismember(model.Elements(:,2),types);
nEdges = numedges(type_idx);
max_dof = max(dofs);  

nDim = floor(Elements(:,2) / 10);       % number of dimensions of elements
if(length(unique(nDim))~=1)
    error('nedelec_mk:inhomogeneous_mesh','inhomogeneous mesh');
end

% dofs = bsxfun(@plus,(dofs-1)*nDim(1),1:nDim(1));

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



% N = Start(end); % total number of matrix entries (non overlapped)
A  = zeros(max_dof,1);  % preallocating A values


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


% Note: 1D and 2D elements only
for type = [12 23 24]
    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+1)).';  % corner vertices
        % Assemble element matrix
        dir_idx = bdirs(i_edge+(1:nEdges),:);
%         dir_idx = repmat(dir_idx,1,nDim(1));
%         dir_idx = reshape(dir_idx',[],1);
        
        a = elem_a(type, vert, w, N, dN, L, dL, dir_idx, fun);
        dof_idx = bdofs(i_edge+(1:nEdges),:);
        dof_idx = reshape(dof_idx',[],1);
        
        A(dof_idx) = a(:);
        iElem = iElem+1;
        progbar(1, nElem, iElem);
        
        i_edge = i_edge + nEdges;
    end

end


%% Postprocessing

% Degree of Freedom vector
% DOF = boundary.Nodes(:,1);


end

%%

function [A] = elem_a(type, xnode, w, N, dNxi, L, dLxi, dir, fun)
%ELEM_A  Acoustic excitation matrices
%   A = ELEM_A(TYPE, XNODE, MAT) Computes the element excitation
%   matrix A of an acoustic element.
%   TYPE  : Element type selector [12 23 24]
%   XNODE : nDim x nNodes matrix containing the coordinates of each node.
%   MAT   : [rho c] vector containing the mass density and sound velocity
%   A     : nNodes x nNodes excitation matrix
%

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

nGauss = size(w,1);

switch type
    case 12
        %% line element (boundary of 2D model)
        dim = 2;
        dr = diff(xnode,1,2);                        
        J = dr / 2;
        t = dr / norm(dr);              % Tangential vector
        skal = J.' * J;
        A = 0;
        for iGauss = 1 : nGauss
            Rt = f_x(iGauss) * t;
            A = A + w(iGauss) * 1/2 * 1/skal * J.' * Rt * norm(dr)/2 * dir;
        end
        
        if(0)
        %j = repmat(sqrt(dot(dr,dr,1))/2,nGauss,1);  % Jacobian
        j = repmat(2./sqrt(dot(dr,dr,1)),nGauss,1);  % Jacobian
        q = kron(w.*j,ones(dim,1));
        %n = [dr(2); -dr(1)];
        Nx = repmat(dr/2,nGauss,1);
        %A = diag(dir)*Nx.'*diag(q)*reshape(repmat(f_x, 1, 2).', [], 1);
        A = diag(dir)*Nx.'*diag(q)*reshape(f_x.', [], 1);
        end
%     case 23
%         %% tria element (boundary of 3D model)
%         df = cross(diff(xnode(:,[1 2]), 1, 2), diff(xnode(:,[1 3]), 1, 2));
%         A = [
%             2 -1 -1
%             -1 2 3
%             -1 3 6
%             ]*sqrt(dot(df,df))/24*(material(1)^2*material(2)^2);
%     case 24
%         %% quad element (boundary of 3D model)
%         load(fullfile(fileparts(mfilename('fullpath')),...
%             'data', 'gauss_quad_3'), 'w', 'N', 'dNxi');
%         ngauss = length(w);
%         J = xnode * dNxi.';
%         df = zeros(3,ngauss);
%         for n = 1 : ngauss
%             df(:,n) = cross(J(:,2*n-1), J(:,2*n));
%         end
%         j = sqrt(dot(df,df,1)).';
%         A = (N.'*diag(w.*j)*N)*(material(1)^2*material(2)^2);
end
end



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