function [A, DOF] = model_a_cyl(model)
% MODEL_A_CYL Excitation matrix assembly
%


% Argument check
narginchk(1, 1);
nargoutchk(1, 1);

% Preprocessing
% Drop non-acoustical elements
AcMat = model.Materials(model.Materials(:,2) == 1, 1);
model.Elements = model.Elements(ismember(model.Elements(:,3), AcMat),:);

% drop User IDs
Elements = drop_IDs(model);

% Preallocate space for the acoustic matrices
nNode = mod(Elements(:, 2), 10);        % number of nodes per element
nVal = nNode.^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

A  = zeros(N,1);            % preallocating mass values

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

% Go through the elements type-by-type
for type = 12
    typeind = Elements(:,2) == type;
    if ~any(typeind)
        continue;
    end
    dim = floor(type/10);               % dimension number
    nnodes = mod(type,10);              % number of nodes of an element
    elements = Elements(typeind,:);     % the elements of this type
    nelem = size(elements,1);           % number of elements of this type
    start = Start(typeind);             % start indices of elements
    [w, N, dN] = elem_quad(type);       % Gaussian quadarture for this type

    % Go through elements of the given type
    for ielem = 1 : nelem
        elem = elements(ielem,:);
        nodes = elem(4+(1:nnodes));             % elem nodes
        vert = model.Nodes(nodes,1+(1:dim+1)).';  % corner vertices
        mater = model.Materials(elem(3), 3:4);  % materials (rho, c)
        % Assemble element matrix
        if dim == 1
            a = elem_a_1d_cyl(type, vert, w, N, dN, mater);
        end
        i = repmat(nodes.', 1, nnodes);         % row indices
        j = i.';                                % column indices
        ind = start(ielem)+(1:nnodes^2);        % index positions
        I(ind) = i(:);
        J(ind) = j(:);
        A(ind) = a(:);
        
        iElem = iElem+1;
        progbar(1, nElem, iElem);
    end % of cycle through elements of the given type
end % of cycle through elements

nNodes = size(model.Nodes,1);           % Number of nodes
% Assemble the sparse matrices
A = sparse(I, J, A, nNodes, nNodes);

% Make symmetric
A = 0.5 * (A + A.');

% DOF output as needed
if nargout > 1
    DOF = mode.Nodes(:, 1);
end

end % of function model_a

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

% Compute Gaussian quadrature weights and shape function samples
nnod = mod(type, 10);
switch type
    case 12
        [xi, w] = gaussquad1(3);
    case {23, 24}
        [xi, w] = gaussquad2(3, nnod);
end

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

end % of function elem_quad

function A = elem_a_1d_cyl(type, xnode, w, N, dNxi, material)

dim = floor(type/10);           %#ok<NASGU>

% Coordinate transform
J = dNxi * xnode.';             % Jacobian of the coordiate transform
% The jacobian determinant is the length of the tangential |J| = |dx/dxi|
j = sqrt(dot(J, J, 2));
r = N * xnode(1, :).';

% Numerical integration using Gaussian quadrature
q = w.*j.*r;
A = N.'*diag(q)*N * material(1)^2*material(2)^2;

end % of function elem_a_1d


function A = elem_a_2d(type, xnode, w, N, dNxi, material)   %#ok<INUSL>

% Coordinate transform
% Jacobian of the coordiate transform
% Each 2-row block of J contains the transposed jacobi
J = dNxi * xnode.';             
X = cross(J(1 : 2 : end, :), J(2 : 2: end, :), 2);
j = sqrt(dot(X, X, 2));

% Numerical integration using Gaussian quadrature
q = w.*j;
A = N.'*diag(q)*N * material(1)^2 * material(2)^2;

end