rhs_part.m

Go to the documentation of this file.

function vec = rhs_part(int_kernel, model, df_info, bc_info)
%function vec = rhs_part(int_kernel, model, df_info, bc_info)
% Fem right-hand side assembly.
%
% Assembly of volume or boundary Fem right hand side. This function
% is used in Fem.Assembly.operator().
%
% boundary integrals are computed if bc_info is provided. Then int_kernel
% is a cell array of length bc_info.nbint_parts.

% IM 16.07.2013


%% initialize
is_boundary = (nargin > 3);

nint_parts = 1;
if is_boundary
    nint_parts = bc_info.nbint_parts;
end

vec = zeros(df_info.ndofs, 1);
vec_cell = {};


%% assemble each integral part
for j = 1:nint_parts
    
    % if available
    if ~is_boundary || ~isempty(int_kernel{j})
        
        if is_boundary
            lin_ind = sub2ind(size(df_info.grid.EL), bc_info.elinds{j}, bc_info.edgeinds{j});
            c = df_info.grid.EL(lin_ind);
        else
            c = df_info.detDF;
            
        end
        
        for i = 1:df_info.ndofs_per_element
            
            if is_boundary
                func = @(x) int_kernel{j}(x, model, df_info, i, bc_info.elinds{j}, bc_info.edgeinds{j});
                res = intervalquadrature(model.qdeg, func);
                
                gids = df_info.get_global_dof_index(bc_info.elinds{j}, i);
            else
                func = @(x) int_kernel(x, model, df_info, i);
                res = triaquadrature(model.qdeg, func);
                
                gids = df_info.get_global_dof_index(1:df_info.grid.nelements, i);
            end
            
            if ~isempty(res)
                
                % depending on decomp_mode the output is a single vector or a cell
                % array of vectors
                if ~iscell(res)
                    
                    res = c .* res;
                    
                    vec_tmp = sparse(gids, ones(length(gids), 1), res, df_info.ndofs, 1);
                    vec = vec + vec_tmp;
                else
                    
                    if ~iscell(vec)
                        vec = cell(1, length(res));
                        for q = 1:length(res)
                            vec{q} = zeros(df_info.ndofs, 1);
                        end
                    end
                    
                    for q = 1:length(res)
                        
                        res{q} = c .* res{q};
                        
                        vec_tmp = sparse(gids, ones(length(gids), 1), res{q}, df_info.ndofs, 1);
                        vec{q} = vec{q} + vec_tmp;
                    end
                end % iscell
            end % isempty
        end
        
        if is_boundary
            % concatenate cell arrays
            if iscell(vec)
                
                vec_cell = [vec_cell, vec(:)'];
                vec = [];
            end
        end
        
    end
    
end

if ~isempty(vec_cell)
    vec = vec_cell;
end

end