operator.m

Go to the documentation of this file.

function varargout = operator(model, model_data)
%function varargout = operator(model, model_data)
% Complete Fem system assembly for lin_stat and stokes models.
%
% This function calls all necessary assembly routines for lin_stat
% and stokes models (possibly nonlinear, i.e. stationary
% Navier-Stokes) and returns a full Fe system with all boundary
% conditions incorporated. This function is used for
% example in lin_stat_model_default() for calling it in the
% detailed simulation, the reduced_data generation and also the
% reduced simulation. It also is called during a fixed-point
% algorithm for the nonlinear stokes model to evaluate the residual.
%
% Return values:
%  varargout: The first output represents the system matrix, the second
%  output (only if requested) the right-hand side vector.
%
% Required fields of model:
%  decomp_mode:  If 0, the usual Fem system is returned. If 1, the
%  parameter independent components are returned in cell arrays. If
%  2, the parametric coefficients are returned in column vectors.
%
% Optional fields of model:
%  uh: An instance of Fem.DiscFunc. When it is present, the
%  trilinear term of the Navier-Stokes equations is evaluated. (Can
%  be used more general for evaluating a nonlinear function).

% (note: strictly speaking (linear) Stokes problems should be
% contained in lin_stat ...)
% maybe check, if general nonlinear terms can be treated with this.

if model.decomp_mode < 2
    
    df_info = model_data.df_info;
    bc_info = model_data.bc_info;
        
    %% assemble system matrix
    if model.decomp_mode == 0
        mat_init = spalloc(df_info.ndofs, df_info.ndofs, 1);
    else
        mat_init = {};
    end
    
    A_vol = mat_init;
    A_bnd = mat_init;
    
    if model.has_volume_integral_matrix
        A_vol = Fem.Assembly.matrix_part(...
            @(varargin) Fem.Assembly.add_integral_kernels(...
            model.matrix_volume_int_kernel, varargin{:}), model, df_info);
    end
    
    if model.has_boundary_integral_matrix
        A_bnd = Fem.Assembly.matrix_part(model.matrix_boundary_int_kernel, ...
            model, df_info, bc_info);
    end
    
    % treat Dirichlet boundary condition
    if model.has_dirichlet_values
        dirichlet_gids = bc_info.dirichlet_gids;
        
        if model.decomp_mode == 0
            A_dir = sparse(dirichlet_gids, dirichlet_gids, ...
                ones(size(dirichlet_gids)), df_info.ndofs, df_info.ndofs);
            
            A_vol(dirichlet_gids, :) = 0;
            A_bnd(dirichlet_gids, :) = 0;
        else
            for q = 1:length(A_vol)
                A_vol{q}(dirichlet_gids, :) = 0;
            end
            for q = 1:length(A_bnd)
                A_bnd{q}(dirichlet_gids, :) = 0;
            end
        end
    end
    
    % add up
    if model.decomp_mode == 0
        A = A_vol + A_bnd;
        if model.has_dirichlet_values && ~isfield(model, 'uh') % only to linear term
            A = A + A_dir;% + eps*speye(df_info.ndofs); % prevents some condition problems
        end
    else
        
        A = [A_vol(:)', A_bnd(:)'];
    end
    
    %% assemble right hand side
    if nargout > 1
        
        if model.decomp_mode == 0
            vec_init = zeros(df_info.ndofs, 1);
        else
            vec_init = {};
        end
        
        r_vol = vec_init;
        r_bnd = vec_init;
        
        if model.has_volume_integral_rhs
            r_vol = Fem.Assembly.rhs_part(...
                @(varargin) Fem.Assembly.add_integral_kernels(...
                model.rhs_volume_int_kernel, varargin{:}), model, df_info);
        end
        
        if model.has_boundary_integral_rhs
            r_bnd = Fem.Assembly.rhs_part(model.rhs_boundary_int_kernel, ...
                model, df_info, bc_info);
        end
        
        % treat Dirichlet boundary condition
        r_dir = vec_init;
        
        if model.has_dirichlet_values
            
            % right hand side
            if model.decomp_mode == 0;
                
                if iscell(bc_info.dirichlet_dof_vector_components)
                    model.decomp_mode = 2;
                    r_dir_coeff = model.dirichlet_values([], [], [], [], model);
                    model.decomp_mode = 0;
                    
                    r_dir = lincomb_sequence(bc_info.dirichlet_dof_vector_components, ...
                        r_dir_coeff);
                else
                    r_dir = bc_info.dirichlet_dof_vector_components;
                end
                
                r_vol(dirichlet_gids) = 0;
                r_bnd(dirichlet_gids) = 0;
            else
                r_dir = bc_info.dirichlet_dof_vector_components;
                
                for q = 1:length(r_vol)
                    r_vol{q}(dirichlet_gids) = 0;
                    r_bnd{q}(dirichlet_gids) = 0;
                end
            end
        end
        
        % add up
        if model.decomp_mode == 0
            r = r_vol + r_bnd + r_dir;
        else
            
            if model.has_nonlinearity
                
                A_tmp = {};
                
                model.uh = Fem.DiscFunc([], model_data.df_info);
                for q = 1:length(r_dir)
                    model.uh.dofs = r_dir{q};
                    
                    A_nl = Fem.Assembly.operator(model, model_data);
                    A_tmp((q-1)*length(A_nl)+(1:length(A_nl))) = A_nl;
                end
                
                A = [A_tmp(:)', A(:)'];
            end
            
            r_tmp = cell(1, length(A)*length(r_dir));
            
            for q1 = 1:length(r_dir)
                for q2 = 1:length(A)
                    r_tmp{(q1-1)*length(A)+q2} = A{q2}*r_dir{q1};
                end
            end
            r = [r_vol(:)', r_bnd(:)', r_tmp(:)'];
        end
    end
    
else % model.decomp_mode = 2
    
    %% get coefficients, system matrix
    A_vol = [];
    A_bnd = [];
    
    if model.has_volume_integral_matrix
        A_vol = Fem.Assembly.eval_coefficients(model.matrix_volume_coeffs, model);
    end
    
    if model.has_boundary_integral_matrix
        A_bnd = Fem.Assembly.eval_coefficients(model.matrix_boundary_coeffs, model);
    end
    
    A = [A_vol(:); A_bnd(:)];
    
    if nargout > 1
        %% get coefficients, right hand side
        r_vol = [];
        r_bnd = [];
        
        if model.has_volume_integral_rhs
            r_vol = Fem.Assembly.eval_coefficients(model.rhs_volume_coeffs, model);
        end
        
        if model.has_boundary_integral_rhs
            r_bnd = Fem.Assembly.eval_coefficients(model.rhs_boundary_coeffs, model);
        end
        
        r_dir = [];
        
        if model.has_dirichlet_values
            r_dir = model.dirichlet_values([], [], [], [], model);
        end
        
        if model.has_nonlinearity
            A_nonlin = A_vol(model.matrix_nonlin_indices);
            A_tmp = A_nonlin(:)*r_dir(:)';
            A = [A_tmp(:); A(:)];
        end
        
        r_tmp = zeros(length(A)*length(r_dir), 1);
        for q = 1:length(r_dir)
            r_tmp((q-1)*length(A)+(1:length(A))) = -A(:)*r_dir(q);
        end
        r = [r_vol(:); r_bnd(:); r_tmp(:)];
    end
end % if model.decomp_mode < 2

varargout{1} = A;
if nargout > 1
    varargout{2} = r;
end

end