comsol_thermal_block_operators.m

function [A, r] = comsol_thermal_block_operators(model,model_data)
%function [A, r] = comsol_thermal_block_operators(model,model_data)
%
% function returns the Operators of the comsol 2D thermal block model with
% affine decomposition!
%
% The Operators are extracted from COMSOL if model.use_comsol=1;
% Otherwise they are just forwarded from model_data where they are stored.
%
% Input:
%   - model
%   - model_data
%
% Output:
%   - A, r: components (model.decomp_mode=1) or coefficients (model.decomp_mode=2)
%           of the System A*u=r
%
% Oliver Zeeb, 25.08.2012



switch model.decomp_mode
    %%%%%%%%%%%%%%%%
    case 0 % full matrices!
        if model.use_comsol
            cmodel = model_data.comsol_model;
            cmodel.sol('sol1').feature('v1').set('scalemethod','none');
            cmodel.sol('sol1').feature('s1').active(false);
            str = mphmatrix(cmodel,'sol1','Out',{'K','L'});
            A = str.K;
            r = str.L;
            cmodel.sol('sol1').feature('v1').set('scalemethod','auto');
            cmodel.sol('sol1').feature('s1').active(true); 
        else %~model.use_comsol
            error('not yet implemented (since this case is most times unneeded...)')
        end
        
    %%%%%%%%%%%%%%%%    
    case 1 % components
        if model.use_comsol
            mu_old = get_mu(model);
            cmodel = model_data.comsol_model;
            cmodel.sol('sol1').feature('v1').set('scalemethod','none');
            cmodel.sol('sol1').feature('s1').active(false);
            A=cell(1,length(mu_old));
            %extract the matrices
            for k = 1:length(mu_old)
                mu = zeros(length(mu_old),1);
                mu(k) = 1;
                model = set_mu(model,mu);
                cmodel = comsol_set_mu_from_rbmatlab_model(model,cmodel);
                str = mphmatrix(cmodel,'sol1','Out',{'K','L'});
                A{k} = str.K;
                if k==1
                    r{k} = str.L;
                end
            end
            model=set_mu(model,mu_old);
            cmodel = comsol_set_mu_from_rbmatlab_model(model,cmodel);
            cmodel.sol('sol1').feature('v1').set('scalemethod','auto');
            cmodel.sol('sol1').feature('s1').active(true); 
            if model.comsol_get_eliminated_data %cut the matrices --> return eliminated matrices!
                for k=1:length(A)
                    A{k} = A{k}(model_data.ind_vectors.com_DOF_ind,model_data.ind_vectors.com_DOF_ind);
                end
                for k=1:length(r)
                    r{k} = r{k}(model_data.ind_vectors.com_DOF_ind);
                end
            end
        else %decomp_mode = 1, ~model.use_comsol
            if model.comsol_get_eliminated_data %return eliminated data
                A = model_data.operators.A_comp_eliminated;
                r = model_data.operators.f_comp_eliminated;
            else %return full data
                A = model_data.operators.A_comp_full;
                r = model_data.operators.f_comp_full;
            end
        end

    
    %%%%%%%%%%%%%%%%
    case 2 %coefficients
        A = get_mu(model);
        r = 1;
    otherwise
        error('unknown decomp_mode')
end