matrix_advection.m

function res = matrix_advection(x, model, df_info, i, j)
%function res = matrix_advection(x, model, df_info, i, j)

% IM 21.07.2016

gradient_hat_phi_i = evaluate_basis_function_derivative(df_info, x, i);
hat_phi_j = evaluate_basis_function(df_info, x, j);

JIT = df_info.grid.JIT;

v = cache_function(@(xarg, xmodel)xmodel.velocity(df_info.grid, ...
    1:df_info.grid.nelements, ...
    xarg, xmodel), x, model);

if ~iscell(v)
    % the following is vectorial of size nelements!
    res = - hat_phi_j * ( v(:,1) .* (JIT(:,1,1) * gradient_hat_phi_i(1) ...
        + JIT(:,1,2)* gradient_hat_phi_i(2)) ...
        + v(:,2) .* (JIT(:,2,1) * gradient_hat_phi_i(1) ...
        + JIT(:,2,2)* gradient_hat_phi_i(2)));
else
    res = cell(1,length(v));
    for q = 1:length(res)
        res{q} = - hat_phi_j * ( v{q}(:,1) .* (JIT(:,1,1) * gradient_hat_phi_i(1) ...
            + JIT(:,1,2)* gradient_hat_phi_i(2)) ...
            + v{q}(:,2) .* (JIT(:,2,1) * gradient_hat_phi_i(1) ...
            + JIT(:,2,2)* gradient_hat_phi_i(2)));
        
    end
end

end