stokes_matrix_fd_nl.m

Go to the documentation of this file.

function res = stokes_matrix_fd_nl(x, model, df_info, i, j)
%function res = stokes_matrix_fd_nl(x, model, df_info, i, j)
% additional trilinear term in Frechet derivative of global operator

% IM 02.12.2104


vel_ndofs_per_element = df_info.values{1}.values{1}.ndofs_per_element;

% detect velocity indices
if all([i j] <= vel_ndofs_per_element*2) && sum([i j] <= vel_ndofs_per_element) ~= 1
    
    if isfield(model, 'uh')
        
        hat_phi_i = evaluate_basis_function(df_info, x, i);
        hat_phi_j = evaluate_basis_function(df_info, x, j);
        
        grad_uh_vals = cache_function(@ ...
                                      (xarg, xmodel)xmodel.uh ...
                                      .evaluate_derivative(':', ...
                                                          xarg), x, model);
        d = cache_function(@(xarg, xmodel)xmodel.divergence_tensor(df_info.grid, ...
                                                          1: ...
                                                          df_info ...
                                                          .grid.nelements, ...
                                                          xarg, ...
                                                          xmodel), x, model);
        
        grad_uh1_phi_i = ...
            grad_uh_vals(:, 1, 1) * hat_phi_i(1) + ...
            grad_uh_vals(:, 1, 2) * hat_phi_i(2);
        grad_uh2_phi_i = ...
            grad_uh_vals(:, 2, 1) * hat_phi_i(1) + ...
            grad_uh_vals(:, 2, 2) * hat_phi_i(2);
        
        if ~iscell(d)
            
            res = hat_phi_j(1) * ...
                (grad_uh1_phi_i .* d(:, 1) + grad_uh2_phi_i .* d(:, 2)) + ...
                hat_phi_j(2) * ...
                (grad_uh1_phi_i .* d(:, 3) + grad_uh2_phi_i .* d(:, 4));
        else
            
            res = cell(1,length(d));
            for q = 1:length(d)
                
                res{q} = hat_phi_j(1) * ...
                    (grad_uh1_phi_i .* d{q}(:, 1) + grad_uh2_phi_i .* d{q}(:, 2)) + ...
                    hat_phi_j(2) * ...
                    (grad_uh1_phi_i .* d{q}(:, 3) + grad_uh2_phi_i .* d{q}(:, 4));
            end
            
        end
    end
    
else    
    res = [];
end

end