function [w, N, dNx, x, e] = fe_quad(nodes, elements, type, order)
%FE_QUAD Quadrature for finite element system matrix assembly
%
%   Output:
%       W       - Quadrature weights
%       N       - Shape function values
%       dNx     - Shape function derivatives w.r.t. x
%       x       - Quadrature base points in x y z
%       e       - Element indices


[w, N, dNxi] = elem_quad(mod(type, 100), order);
n_gauss = size(N, 1);

n_nod = mod(type, 10);
dim = floor(mod(type, 100) / 10);

elem_sel = (elements(:, 2) == type);
node_inds = elements(elem_sel, 4 + (1 : n_nod));
n_elems = size(node_inds, 1);

if n_elems == 0
    w = zeros(0, 1);
    dNx = zeros(0, size(N, 2));
    x = zeros(0, dim);
    if nargout > 4
        e = 0;
    end
else

    x_nodes = nodes(reshape(node_inds.', [], 1), 1 + (1:dim));
    jacobi = block_product(dNxi,  x_nodes); 
    det_jac = block_det(jacobi);
    dNx = multi_block_inverse(jacobi, repmat(dNxi, n_elems, 1));
    w = abs(det_jac) .* repmat(w, n_elems, 1);
    if nargout > 3
        x = block_product(N, x_nodes);
        if nargout > 4
            e = kron(find(elem_sel), ones(n_gauss, 1));
        end
    end
end

end % of function fe_quad


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

% Compute Gaussian quadrature weights and shape function samples
nnod = mod(type, 10);
switch type
    case 12
        [xi, w] = gaussquad1(order);
    case {23, 24}
        [xi, w] = gaussquad2(order, nnod);
    case 34
        [xi, w] = gaussquad3(order, nnod);
    case {36, 38}
        [xi, w] = gaussquad3(order, nnod);
end
[N, dNxi] = shapefun(xi, type);
dNxi = reshape(shiftdim(dNxi,2), numel(dNxi)/nnod, nnod);

end % of function elem_quad