matrix_part.m

Go to the documentation of this file.

function Mat = matrix_part(int_kernel, model, df_info1, varargin)
%function Mat = matrix_part(int_kernel, model, df_info1, varargin)
% Fem Matrix assembly.
%
% Assembly of volume or boundary FE matrix, possibility of
% providing two different function spaces on same grid or on grids with
% common boundary in the case of boundary integrals. This function
% is used in Fem.Assembly.operator().
%
% varargin possibilities:
% - {}
% - df_info2
% - bc_info1
% - bc_info1, df_info2, bc_info2
%
% boundary integrals are computed if bc_info(s) are provided.

% IM 16.07.2013


%% initialize
is_boundary = 0;
spaces_identical = 1;

if nargin > 3
    
    if nargin > 4
        
        is_boundary = 1;
        spaces_identical = 0;
        
        bc_info1 = varargin{1};
        df_info2 = varargin{2};
        bc_info2 = varargin{3};
    else
        
        if isa(varargin{1}, 'Fem.IFemInfo')
            
            spaces_identical = 0;
            
            df_info2 = varargin{1};
        else
            
            is_boundary = 1;
            
            bc_info1 = varargin{1};
            df_info2 = df_info1;
            bc_info2 = bc_info1;
        end
    end
else
    
    df_info2 = df_info1;
end

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

Mat = spalloc(df_info1.ndofs, df_info2.ndofs, 1);
Mat_cell = {};


%% assemble each integral part
for k = 1:nint_parts
    
    % if available
    if ~is_boundary || ~isempty(int_kernel{k})
        
        if is_boundary
            lin_ind = sub2ind(size(df_info1.grid.EL), bc_info1.elinds{k}, bc_info1.edgeinds{k});
            c = df_info1.grid.EL(lin_ind);
        else
            c = df_info1.detDF;
        end
        
        for i = 1:df_info1.ndofs_per_element
            
            for j = 1:df_info2.ndofs_per_element
                
                if is_boundary
                    if spaces_identical
                        func = @(x) int_kernel{k}(x, model, df_info1, i, j, bc_info1.elinds{k}, bc_info1.edgeinds{k});
                    else
                        func = @(x) int_kernel{k}(x, model, df_info1, df_info2, i, j, ...
                            bc_info1.elinds{k}, bc_info1.edgeinds{k}, bc_info2.elinds{k}, bc_info2.edgeinds{k});
                    end
                    res = intervalquadrature(model.qdeg, func);
                    
                    gids_i = df_info1.get_global_dof_index(bc_info1.elinds{k}, i);
                    gids_j = df_info2.get_global_dof_index(bc_info2.elinds{k}, j);
                else
                    if spaces_identical
                        func = @(x) int_kernel(x, model, df_info1, i, j);
                    else
                        func = @(x) int_kernel(x, model, df_info1, df_info2, i, j);
                    end
                    res = triaquadrature(model.qdeg, func);
                    
                    gids_i = df_info1.get_global_dof_index(':', i);
                    gids_j = df_info2.get_global_dof_index(':', j);
                end
                
                if ~isempty(res)
                    
                    % depending on decomp_mode the output is a single matrix or a cell
                    % array of matrices
                    if ~iscell(res)
                        
                        res = c .* res;
                        
                        Mat_tmp = sparse(gids_i, gids_j, res, df_info1.ndofs, df_info2.ndofs);
                        Mat = Mat + Mat_tmp;
                        
                    else
                        
                        if ~iscell(Mat)
                            Mat = cell(1, length(res));
                            for q = 1:length(res)
                                Mat{q} = spalloc(df_info1.ndofs, df_info2.ndofs, 1);
                            end
                        end
                        
                        for q = 1:length(res)
                            
                            res{q} = c .* res{q};
                            
                            Mat_tmp = sparse(gids_i, gids_j, res{q}, df_info1.ndofs, df_info2.ndofs);
                            Mat{q} = Mat{q} + Mat_tmp;
                        end
                        
                    end % iscell
                end % isempty
            end            
        end
        
        if is_boundary
            % concatenate cell arrays
            if iscell(Mat)
                
                Mat_cell = [Mat_cell, Mat(:)'];
                Mat = [];
            end
        end
        
    end
    
end

if ~isempty(Mat_cell)
    Mat = Mat_cell;
end

cache('clear');

end