function [Mw, Kwt, DOF] = model_wall_mkt_cyl(model)
% MODEL_WALL_MKT_CYL Cylindrical mass, stiffness matrices of a viscous wall
%   [MW, KWT, DOF] = MODEL_WALL_MKT_CYL(MODEL) Assembles wall mass and 
%       stiffness matrices MW and KWT for viscothermal acoustic simulations 
%       of the finite element model MODEL.
%       The input model should contain the wall boundary of a cylindrically
%       symmetric 2D geometry. (r variable: x, symmetry axis: y)
%   Inputs:
%       MODEL   - Finite element model of the wall
%   Outputs:
%       MW      - Wall mass matrix
%       KWT     - Wall stiffness matrix containing the tangential gradients
%       DOF     - Degrees of freedom vector
%
%   References:
%       [1] M. Berggren, A. Bernland, D. Noreland. "Acoustic boundary 
%           layers as boundary conditions." Journal of Computational 
%           Physics 371 (2018), pp. 633-650
%           DOI: 10.1016/j.jcp.2018.06.005

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

% 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

Mw  = zeros(N,1);           % preallocating mass values
Kwt  = zeros(N,1);          % preallocating stiffness 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
            [m, k] = elem_mkt_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(:);
        Mw(ind) = m(:);
        Kwt(ind) = k(:);
        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
% sel = ismember(1:nNodes, I);            % Logical array of selected nodes
% vNodes = sum(sel);                      % Number of value nodes
% index = zeros(vNodes,1);                % Indices
% index(sel) = 1:vNodes;
% 
% I = index(I);
% J = index(J);
% 
% % The DOF vector contains the node IDs corresponding to matrix rows
DOF = model.Nodes(:,1);

% Assemble the sparse matrices
Mw = sparse(I, J, Mw, nNodes, nNodes);
Kwt = sparse(I, J, Kwt, nNodes, nNodes);

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

end % of function model_wall_mkt

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 [M, K] = elem_mkt_1d_cyl(type, xnode, w, N, dNxi, material)

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

% Coordinate transform
J = dNxi * xnode.';             % Jacobian of the coordiate transform
dNtx = zeros(size(dNxi));       % Tangential gradient in the xyz space
j = zeros(ngauss,1);            % Jacobi determinant
r = N * xnode(1, :).';

for n = 1 : ngauss
    % The jacobian determinant is the length of the tangential:
    %   |J| = |dx/dxi|
    j(n) = norm(J(n, :));
    % Assemble the tangential transformation matrix
    T = J(n, :) / norm(J(n, :));
    D = J(n, :) * J(n, :).';
    % Tangential derivatives
    dNtx(n, :) = T * J(n, :).' * (D.' \ dNxi(n, :));
end

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

end % of function elem_mkt_1d_cyl


