comsol_ThermalBlock3D_operators.m

function [A, r] = comsol_ThermalBlock3D_operators(model,model_data)
%function [A, r] = comsol_ThermalBlock3D_operators(model,model_data)
%
% function returns the operators of the comsol thermal block model with
% affine decomposition!
%
% Oliver Zeeb, 25.08.2012

if model.use_comsol
    switch model.decomp_mode
        %%%%%%%%%%%%%%%%
        case 0 % full matrices!
            warning('decomp_mode 0 NOCH TESTEN!!!');
            keyboard;

            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); 
            if model.comsol_get_eliminated_data
                A = A(model_data.ind_vectors.com_DOF_ind,model_data.ind_vectors.com_DOF_ind);
                r = r(model_data.ind_vectors.com_DOF_ind);
            end

        %%%%%%%%%%%%%%%%    
        case 1 % components
            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));
            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);
                cmodel.sol('sol1').feature('s1').active(false);
                str = mphmatrix(cmodel,'sol1','Out',{'K','L'},'initmethod','init', 'initsol', 'zero');
                A{k} = str.K;
                if k==1
                    r{k} = str.L;
                end
            end
            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

            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); 

        %%%%%%%%%%%%%%%%
        case 2 %coefficients
            A = get_mu(model);
            r = 1;
        otherwise
            error('unknown decomp_mode')
    end
    
else % use only matlab, i.e. saved matrices!
    switch model.decomp_mode
        case 0
            warning('decomp_mode 0 not yet implemented (because unneeded!!!)')
            keyboard
        case 1
            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
        case 2
            A = get_mu(model);
            r = 1;
        otherwise
            error('unknown decomp_mode')
    end
end