comsol_ThermalBlock3D_calculate_inner_product_matrices.m

function inner_product_matrices = comsol_ThermalBlock3D_calculate_inner_product_matrices(model,model_data)
%function inner_product_matrices = comsol_ThermalBlock3D_calculate_inner_product_matrices(model,model_data)
%
% This function computes the inner product matrices (L2 and H1) for the 
% comsol 3D Thermal Block model! It can be used for all thermal block models,
% no matter how many Parameters are specified.
% For the L2-Matrix, the function sets the reaction coefficient a of the
% model = 1 and extracts the stiffness matrix from comsol. For the
% H1-Matrix, the reaction coefficient a and the diffusion coefficient c
% must be = 1:  a=c=1;
% After extracting the Matrices this function sets the parameters in the 
% Comsol model to the values they had before this function changed them.
%
% Input: model, model_data
% Output: struct containing the inner product matrices
%
% Oliver Zeeb, 08/2013

comsol_model = model_data.comsol_model;

%get all tags of the coefficient forms (so that this file works with all
%thermal block models no matter which dimension they are
all_feature_tags = cellstr(char(comsol_model.physics('c').feature.tags));
nr_cfeq=0;
for k = 1:length(all_feature_tags)
    feature_type = char(comsol_model.physics('c').feature(all_feature_tags{k}).getType());
    if strcmp(feature_type, 'CoefficientFormPDE')
        nr_cfeq=nr_cfeq+1;
        cfeq_tag{nr_cfeq} = all_feature_tags{k};
    end
end

if nr_cfeq==0
    warning('No CoefficientFormPDE found in the model!')
    keyboard
end

%get all original parameters
for k=1:nr_cfeq
    c_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('c')));
    a_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('a')));
    f_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('f')));
    ea_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('ea')));
    da_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('da')));
    al_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('al')));
    be_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('be')));
    ga_coeff{k} = cellstr(char(comsol_model.physics('c').feature(cfeq_tag{k}).getStringMatrix('ga')));
end
    
%reset the coefficients to a=1, all others=0 for all coefficient forms!--> L2-inner-product-matrix
for k=1:nr_cfeq    
    comsol_model.physics('c').feature(cfeq_tag{k}).set('c', '0');
    comsol_model.physics('c').feature(cfeq_tag{k}).set('a', '1');
    comsol_model.physics('c').feature(cfeq_tag{k}).set('f', '0');
    comsol_model.physics('c').feature(cfeq_tag{k}).set('ea', '0');
    comsol_model.physics('c').feature(cfeq_tag{k}).set('da', '0');
    comsol_model.physics('c').feature(cfeq_tag{k}).set('al', {'0' '0' '0'});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('be', {'0' '0' '0'});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('ga', {'0' '0' '0'});
end
    
%get L2 matrix:
comsol_model.sol('sol1').feature('v1').set('scalemethod','none');
comsol_model.sol('sol1').feature('s1').active(false);
str = mphmatrix(comsol_model,'sol1','Out',{'K','Kc'},'initmethod','init', 'initsol', 'zero');
inner_product_matrices.L2 = str.K;
inner_product_matrices.L2_eliminated = str.Kc;
    
    %set c=a=1, all others 0 --> H1-inner-product-matrix
for k=1:nr_cfeq
    comsol_model.physics('c').feature(cfeq_tag{k}).set('c', '1');
end
    
%get H1 matrix:
comsol_model.sol('sol1').feature('v1').set('scalemethod','none');
comsol_model.sol('sol1').feature('s1').active(false);
str = mphmatrix(comsol_model,'sol1','Out',{'K','Kc'},'initmethod','init', 'initsol', 'zero');
inner_product_matrices.H1 = str.K;
inner_product_matrices.H1_eliminated = str.Kc;
    

%change back to original values:
for k=1:nr_cfeq
    comsol_model.physics('c').feature(cfeq_tag{k}).set('c', c_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('a', a_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('f', f_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('ea', ea_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('da', da_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('al', al_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('be', be_coeff{k});
    comsol_model.physics('c').feature(cfeq_tag{k}).set('ga', ga_coeff{k});
end