function F = model_f(model)

%% 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
nDim = floor(Elements(Elements(:,2) < 100,2)/10);         % Number of dimensions

% number of matrix entries per elem
% We have nDim^2 DOF for each node in the trial space,
% And 1 DOF for each node in the test space
nVal = nNode.^2.*nDim.^2;          

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

F  = zeros(N,1);    % preallocating forcing vector values

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

% Go through elements by element type
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
    [w, N, dN] = elem_quad(type);       %#ok<ASGLU> % 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
        f = elem_f(type, vert, w, dN);
        i = repmat(nodes.', 1, nnodes*dim^2);         % row indices
        j = reshape(repmat(dim^2*(nodes-1),nnodes*dim^2,1),nnodes,nnodes*dim^2);  % column indices
        sh = repmat(1:dim^2,nnodes,nnodes);
        j = j + sh;
        
        ind = start(ielem)+(1:(nnodes^2*dim^2));        % index positions
        I(ind) = i(:);
        J(ind) = j(:);
        F(ind) = f(:);
        iElem = iElem+1;
        progbar(1, nElem, iElem);
    end
end
    
%% Generate the sparse matrices from the arrays
nNodes = size(model.Nodes,1);
%m = ismember(1:nNodes,[I;iI]);
%vNodes = sum(m);
%index = zeros(vNodes,1);
%index(m) = 1:vNodes;
%I = index(I);
%J = index(J);
F = sparse(I, J, F, nNodes, dim^2*nNodes);


end


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, nnod);
    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

% Assemble force interpolation vector for element
function F = elem_f(type, xnode, w, dNxi)

dim = floor(type/10);       % Element (local) dimensions
ngauss = length(w);         % Number of gaussian points
nnode = mod(type,10);       % Number of element nodes

% Preallocate for shape function derivatives
dNxe = zeros(dim*ngauss, dim^2*nnode);

% 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,:);
    % Temporary
    %dNxe(ind(1),1:4:end) = dNx(ind(1),:);
    %dNxe(ind(2),2:4:end) = dNx(ind(2),:);
    %dNxe(ind(1),3:4:end) = dNx(ind(1),:);
    %dNxe(ind(2),4:4:end) = dNx(ind(2),:);
    dNxe(ind(1),1:4:end) = dNx(ind(1),:);
    dNxe(ind(2),2:4:end) = dNx(ind(1),:);
    dNxe(ind(1),3:4:end) = dNx(ind(2),:);
    dNxe(ind(2),4:4:end) = dNx(ind(2),:);
    j(n) = abs(det(Jac));
end

% Numerical integration with Gaussian quadrature
q = reshape(repmat((w.*j).',dim, 1),dim*ngauss,1);

% The extended version of the shape function derivatives

% Fill the extended matrix
% for i = 1:dim
%     for j = 1:dim
%         dNxe(i:dim:end,(i+(j-1)*dim):dim^2:end) = dNx(j:dim:end,:);
%     end
% end

F = -dNx.' * diag(q) * dNxe;
end