OperatorsDefault.m

classdef OperatorsDefault < handle
    %classdef OperatorsDefault < handle
    % Fem operators can be obtained from stored components.
    %
    % This class helps storing Fem system components for accelerating the
    % operator evaluation, if it is done multiple times. Since the
    % compontents are accessible, it is also useful for the reduced data
    % generation. Note that this is compatible with the function
    % Fem.Assembly.operator.
    %
    % The components can be assembled with the method
    % assemble_components(). If "operators" is a Fem.OperatorsDefault
    % object, the components are stored in operators.A_comp and
    % operators.f_comp. To obtain the Fem system simply call [A, f] =
    % operators(model, model_data).
    %
    % In stokes problems, it can only be used for linear terms.
    
    properties (SetAccess=private)
        
        A_comp = {};
        f_comp = {};
        componentsAssembled = 0;
    end
    
    methods
        
        function this = OperatorsDefault()
            %function this = OperatorsDefault()
            % Only initializing new object.
        end
        
        function varargout = evaluate(this, model, model_data)
            %function varargout = evaluate(this, model, model_data)
            % Operator evaluation. Returns a full Fem system.
            %
            % If the components are stored, the system is computed via
            % linear combinations of the components, else the assembly
            % routine of model is called.
            %
            % Return values:
            %  varargout: The first output represents the system matrix,
            %  the second output (only if requested) the right-hand side
            %  vector. The field model.decomp_mode has no influence.
            
            if this.componentsAssembled
                model.decomp_mode = 2;
                [A_coeff, f_coeff] = model.operators(model, model_data);
                % the following is needed because dirichlet dofs are
                % ignored in assembly of components for reduced data
                if model.has_nonlinearity
                    Qr_dir = length(model_data.bc_info.dirichlet_dof_vector_components);
                    QA_nonlin = length(model.matrix_nonlin_indices);                    
                    A_coeff(1:Qr_dir*QA_nonlin) = 0;
                end
                varargout{1} = lincomb_sequence(this.A_comp, A_coeff);
                dirichlet_gids = model_data.bc_info.dirichlet_gids;
                varargout{1} = varargout{1} + sparse(dirichlet_gids, dirichlet_gids, ...
                    ones(size(dirichlet_gids)), model_data.df_info.ndofs, ...
                    model_data.df_info.ndofs) + eps*speye(model_data.df_info.ndofs);
                if nargout > 1
                    r_dir_coeff = model.dirichlet_values([], [], [], [], model);
                    r_dir = lincomb_sequence(...
                        model_data.bc_info.dirichlet_dof_vector_components, ...
                        r_dir_coeff);                    
                    f_coeff(end+1-(1:length(r_dir_coeff)*length(A_coeff))) = 0;
                    varargout{2} = lincomb_sequence(this.f_comp, f_coeff) + r_dir;
                end
            else
                model.decomp_mode = 0;
                [varargout{1:nargout}] = model.operators(model, model_data);
            end                    
        end
        
        function this = assemble_components(this, model, model_data)
            %function this = assemble_components(this, model, model_data)
            % Calls the assembly routine of model and stores all
            % components.
            model.decomp_mode = 1;
            [this.A_comp, this.f_comp] = model.operators(model, model_data);
            this.componentsAssembled = 1;
        end
        
        function this = clear_components(this)
            %function this = clear_components(this)
            % Deletes all components.
            this.A_comp = {};
            this.f_comp = {};
            this.componentsAssembled = 0;
        end
        
        function varargout = subsref(this, S)
            %function varargout = subsref(this, S)
            % Overloads subsref.
            %
            % This enables an operator evaluation via operators(model,
            % model_data), if "operators" is a Fem.OperatorsDefault object.
            if strcmp(S(1).type, '()')
                [varargout{1:nargout}] = evaluate(this, S(1).subs{:});
            else
                [varargout{1:nargout}] = builtin('subsref', this, S);
            end
        end
        
    end %methods
    
end %classdef