function I = cell2node_interp(model)
%CELL2NODE_INTERP Interpolation from cell centers into node
%   I = CELL2NODE_INTERP(MODEL)
%       Creates the sparse interpolation matrix I that can interpolate
%       cell-based values into nodal values. The matrix I is N X E where N
%       is the number of nodes and E is the number of elements in MODEL.
%
%       The interpolation 

elem = drop_IDs(model);
nElem = size(elem,1);
nNodes = size(model.Nodes,1);

elem_types = [12; 23; 24; 34; 36; 38];
n_types = numel(elem_types);

cent = centnorm(model);

is = ismember(elem(:,2), elem_types);
n_values = sum(mod(elem(is,2), 10));

rows = zeros(n_values, 1);
cols = zeros(n_values, 1);
vals = zeros(n_values, 1);

start = 0;
for i_type = 1 : n_types
    type = elem_types(i_type);
    nnod = mod(type, 10);
    sel = (elem(:, 2) == type);
    n_type = sum(sel);
    if n_type == 0
        continue
    end
    % Node indices as a column vector
    node_inds = reshape(elem(sel, 4 + (1:nnod)).', [], 1);
    % Node coordinates
    node_crd = model.Nodes(node_inds, 2:4);
    % Vector from center to each node
    r = node_crd - kron(cent(sel, :), ones(nnod, 1));
    % Inverse distance from center to each node
    inv_d = 1 ./ sqrt(dot(r, r, 2));
    
    % Indices in the rows, cols, vals vectors
    inds = start + (1 : n_type*nnod);
    rows(inds) = node_inds;
    cols(inds) = kron(find(sel), ones(nnod, 1));
    vals(inds) = inv_d;
    
    start = start + n_type * nnod;
end

% Create a sparse matrix
I = sparse(rows, cols, vals, nNodes, nElem);
% Weighting matrix 
W = spdiags(1./sum(I,2), 0, nNodes, nNodes);

I = W * I;

end

