function [ M,K ] = nedelec_mk( model )
%NEDELEC_MK Summary of this function goes here
%   Detailed explanation goes here

% drop User IDs
Elements = drop_IDs(model);

Edges = get_edges(Elements);
[dofs, dir] = get_edge_dofs(Edges);
max_dof = max(dofs);

% Preallocating space for the matrices
% The elements of the sparse matrices are collected in arrays. The matrix
% values of each element are stored in non overlapping parts of the arrays.
% do not calculate infinite elements

i_edge = 0;
types = [12 23 24 34 36 38]';
numedges = [1 3 4 6 9 12]';

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

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)
I  = zeros(N,1);  % preallocating row indices
J  = zeros(N,1);  % preallocating column indices
M  = zeros(N,1);  % preallocating mass values
K  = zeros(N,1);  % preallocating stiffness values

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


% Note: 3D elements only
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 = dir(i_edge+(1:nEdges),:);
%         dir_idx = repmat(dir_idx,1,nDim(1));
%         dir_idx = reshape(dir_idx',[],1);
        
        [m, k] = elem_mk(type, vert, w, N, dN, dL, dir_idx);
        dof_idx = dofs(i_edge+(1:nEdges),:);
        dof_idx = reshape(dof_idx',[],1);
        
        i = repmat(dof_idx,1,nEdges);         % row indices
        j = i.';                                % column indices
        ind = start(ielem)+(1:(nEdges)^2);        % index positions
        I(ind) = i(:);
        J(ind) = j(:);
        M(ind) = m(:);
        K(ind) = k(:);
        iElem = iElem+1;
        progbar(1, nElem, iElem);
        
        i_edge = i_edge + nEdges;
    end

end


M = sparse(I,J,M,max_dof,max_dof);
K = sparse(I,J,K,max_dof,max_dof);

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


function [M, K] = elem_mk(type, xnode, w, N, dNxi, dLxi, dir)

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

% Coordinate transform
J = dLxi * xnode.';                 % Jacobian of the coordiate transform
%dLx = zeros(size(dLxi));            % Gradient in xyz space
%dNx = zeros(size(dNxi)); 
Nx = zeros(size(N));

if dim==2
    rotNx = zeros(ngauss,size(N,2));
elseif dim==3
    rotNx = zeros(ngauss*3,size(N,2));
end

j = zeros(ngauss,1);

for n = 1 : ngauss
    ind = dim*(n-1)+(1:dim);
    Jac = J(ind,:);                     % Jacobi transposed
    
    Nx(ind,:) = Jac \ N(ind,:);
    
%    dLx(ind,:) = (Jac) \ dLxi(ind,:);

    % Pre-expand shape function derivatives
%     for i = 1 : dim
%        row_inds = (n-1)*dim*dim + (i-1)*dim + (1 : dim);
%        dNx(row_inds,:) = Jac \ dNxi(row_inds,:);
%     end
    
    j(n) = abs(det(Jac));
    
    if dim==2
        rotNx(n,:) = (dNxi((n-1)*dim*dim+3,:) - dNxi((n-1)*dim*dim+2,:)) / j(n);
    elseif dim==3
        rotNx((n-1)*3+(1:3),:) = Jac.' * [ ...
                dNxi((n-1)*dim*dim+8,:) - dNxi((n-1)*dim*dim+6,:) ;
                dNxi((n-1)*dim*dim+3,:) - dNxi((n-1)*dim*dim+7,:) ;
                dNxi((n-1)*dim*dim+4,:) - dNxi((n-1)*dim*dim+2,:)
                        ] / j(n);
    end
        
        
    
end

% Numerical integration with Gaussian quadrature
% Use the same shape functions in all dimensions

% Mass matrix assembply
q = kron(w.*j,ones(dim,1));
M = diag(dir)*Nx.'*diag(q)*Nx*diag(dir);

% Stiffness matrix assembly

if dim==2
    q = w.*j;
elseif dim==3
    q = kron(w.*j, ones(dim, 1));
end

K = diag(dir)*rotNx.'*diag(q)*rotNx*diag(dir);




end