function [M, K, G] = model_mkg(model)
%MODEL_MKG Compute mass and stiffness matrices of the FE model
%   [M K G DOF] = MODEL_MK(MODEL) Computes the mass matrix M, stiffness
%   matrix K and damping matrix C of a FE/IE model MODEL. The output vector
%   DOF relates the elements of the output matrices to the degrees of
%   freedom of the model.
%   [M K DOF] = MODEL_MK(MODEL)
%   [I J M K C DOF] = MODEL_MK(MODEL, 'ind')
%   [I J M K DOF] = MODE_MK(MODEL, 'ind')
% See also: MODEL_MK_CYL, MODEL_A, MODEL_A_CYL

% Peter Fiala and Peter Rucz
% 2009 - 2011

% 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);
% PML workaround
Elements = Elements(Elements(:,2) < 200,:);

%% Preallocating space for the acoustic 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
nNode = sum(Elements(Elements(:,2) < 100,5:end)~=0, 2);   % number of nodes per element
nVal = nNode.^2;                        % number of matrix entries per elem
dims = floor(Elements(:,2)/10);

udim = unique(dims);
if numel(udim) ~= 1
    error('model_mkg:not_unique_dimension', ...
        'Unique dimension needed!');
end
Start = [0; cumsum(nVal)];              % starting indices of elements
StartG = [0; cumsum(nVal.*dims)];

N = Start(end); % total number of matrix entries (non overlapped)
NG = StartG(end);
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

IG = zeros(NG, 1);
JG = zeros(NG, 1);
G  = zeros(NG, 1);


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

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);              % 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
    startg = StartG(typeind);
    [w, N, dN] = elem_quad(type);       % Gaussian quadarture for this type
    for ielem = 1 : nelem
        elem = elements(ielem,:);
        nodes = elem(4+(1:nnodes));             % elem nodes
        vert = model.Nodes(nodes,1+(1:dim)).';  % corner vertices
        % Assemble element matrix
        [m, k, g] = elem_mkg(type, vert, w, N, dN, model.Materials(elem(3),3:4));
        i = repmat(nodes.', 1, nnodes);         % row indices
        j = i.';                                % column indices
        
        ig = repmat(nodes.', dim, nnodes);
        jg = bsxfun(@plus, dim*(nodes.'-1), (1 : dim));
        jg = kron(reshape(jg.', [], 1), ones(nnodes, 1));
        
        ind = start(ielem)+(1:nnodes^2);        % index positions
        I(ind) = i(:);
        J(ind) = j(:);
        M(ind) = m(:);
        K(ind) = k(:);
        
        indg = startg(ielem)+(1:nnodes^2*dim);
        IG(indg) = ig(:);
        JG(indg) = jg(:);
        G(indg) = g(:);
        
        iElem = iElem+1;
        progbar(1, nElem, iElem);
    end
end

%% Generate the sparse matrices from the arrays
nNodes = size(model.Nodes,1);               % Number of nodes

M = sparse(I, J, M, nNodes, nNodes);
K = sparse(I, J, K, nNodes, nNodes);
G = sparse(IG, JG, G, nNodes, nNodes*udim);

end % of function model_mkg


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(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
[N, dNxi] = shapefun(xi, type);
dNxi = reshape(shiftdim(dNxi,2),numel(dNxi)/nnod,nnod);
end

function [M, K, G] = elem_mkg(type, xnode, w, N, dNxi, material)

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

%% Coordinate transform
J = dNxi * xnode.';      % Jacobian of the coordiate transform
dNx = zeros(size(dNxi)); % Gradient in xyz space
j = zeros(ngauss,1);
for n = 1 : ngauss
    ind = dim*(n-1)+(1:dim);
    Jac = J(ind,:);
    dNx(ind,:) = (Jac) \ dNxi(ind,:);
    j(n) = abs(det(Jac));
end

%% Numerical integration with Gaussian quadrature
q = w.*j;
M = N.'*diag(q)*N * material(1);
G = zeros(size(N,2), dim*size(N,2));
for i = 1 : dim
    G(:, i : dim : end) = N.' * diag(q) * dNx(i:dim:end,:) * material(1);
end
q = reshape(repmat(q.',dim, 1),dim*ngauss,1);
K = dNx.' * diag(q) * dNx * (material(1)*material(2)^2);
end