matrix_robin.m

function res = matrix_robin(x, model, df_info, i, j, elind, edgeind)
%function res = matrix_robin(x, model, df_info, i, j, elind, edgeind)
%
% auxiliary function for integral kernel for A_robin
% integral kernel on all robin-edges simultaneously
% here x is a scalar value on unit interval
% f(x) = (alpha/beta + b^T n) hatphi_i hatphi_j
% multiplication with |edge| is realized in caller after
% quadrature.
% This function can handle cell-array valued data

% B. Haasdonk 22.2.2011
% IM 21.07.2016

res = zeros(0,1);
grid = df_info.grid;
for local_edge_id = 1:3
    lcoord = llocal2local(grid,local_edge_id,x); % local coord for x on edge 1
    hat_phi_i = evaluate_basis_function(df_info,lcoord,i);
    hat_phi_j = evaluate_basis_function(df_info,lcoord,j);
    elinds = elind(edgeind==local_edge_id);
    
    if model.decomp_mode == 0
        res = [res; (model.robin_beta(grid,elinds,local_edge_id,x, ...
            model).^(-1)) .* ...
            model.robin_alpha(grid,elinds,local_edge_id,x, ...
            model) * hat_phi_i * hat_phi_j ];
    else
        % number of components = 1
        if ~iscell(res)
            res = {zeros(0, 1)};
        end
        
        res{1} = [res{1}; ...
            ( (model.robin_beta(grid,elinds,local_edge_id,x, ...
            model).^(-1)) .* ...
            model.robin_alpha(grid,elinds,local_edge_id,x, ...
            model) ...
            ) * hat_phi_i * hat_phi_j ]; % local model!
    end
    
    if model.has_advection
        b = cache_function(@(xelind, xarg, xmodel)...
                           xmodel.velocity(grid, xelind, ...
                                           xarg, xmodel), ...
                           elinds, lcoord, model);
        
        if ~iscell(b)
            res(end-numel(elinds)+1:end) ...
                = res(end-numel(elinds)+1:end) ...
                      + ( b(:,1) .* grid.NX(elinds,local_edge_id) ...
                          + b(:,2) .* grid.NY(elinds,local_edge_id) ...
                          ) * hat_phi_i * hat_phi_j; % local model!
        else % iscell => components mode
            
            % number of components = 1+ Q_velocity y
            if numel(res) < length(b)+1
                res = [res, cell(1,length(b))];
                for q = 1:length(b)
                    res{q+1} = zeros(0,1);
                end
            end
            
            for q = 1:length(b)
                res{q+1} = [res{q+1};...
                    ( b{q}(:,1).*grid.NX(elinds,local_edge_id) ...
                    + b{q}(:,2).*grid.NY(elinds,local_edge_id) ...
                    ) * hat_phi_i * hat_phi_j; ...
                    ];
            end
            %    keyboard;
        end
        
    end %model.has_advection
end %local_edge_id

end