comsol_calculate_inner_product_matrices_default.m

function inner_product_matrices = comsol_calculate_inner_product_matrices_default(model,model_data)
%function [L2, L2_full, H1, H1_full] = comsol_get_inner_product_matrices(comsol_model)
%
% function returns the L2 and H1 inner product matrices using a comsol
% coefficient form pde model. This is done by setting the parameters of the
% coefficient form to 
% a=1, all others=0 --> L2-matrix
% a=c=1, all others=0 --> H1-matrix
%
% input:
% comsol_model: a comsol coefficient PDE form model
%
% output:
% L2, L2_full, H1, H1_full: L2 and H1 inner product matrices (full and
% eliminated versions)
%
% How it works::
% 1.) get all the original parameter and store them.
% 2.) set the reactin (L2) or reaction and diffusion (H1) parameters to 1
%     all others to zero
% 3.) extract the stiffness matrix str.K (and str.Kc)
% 4.) set back all the parameters to their original values
%     so that the comsol model is the same as it was when this function is
%     called.
%
% Oliver Zeeb
% 2013/06/05



comsol_model=model_data.comsol_model;
ptag = model.comsol_tags.physics;

% get number of cfeq's used in the comsol model
feature_tags = char(comsol_model.physics(ptag).feature.tags);
nr_cfeq = 0;
for k=1:length(feature_tags)
    if strfind(feature_tags(k,:), 'cfeq')
        nr_cfeq = nr_cfeq + 1;
    end
end



%get all original parameters
for k=1:nr_cfeq
    c_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('c')));
    a_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('a')));
    f_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('f')));
    ea_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('ea')));
    da_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('da')));
    al_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('al')));
    be_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('be')));
    ga_coeff{k} = cellstr(char(comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).getStringMatrix('ga')));
    
    
    %reset the coefficients to a=1, all others=0 --> L2-inner-product-matrix
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('c', mat2cell(num2str(zeros(size(c_coeff{k}))), ones(size(c_coeff{k},1),1)));
    %comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('c', {'0' '0' '0' '0'});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('a', '1');
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('f', '0');
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('ea', '0');
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('da', '0');
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('al', mat2cell(num2str(zeros(size(al_coeff{k}))), ones(size(al_coeff{k},1),1)));
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('be', mat2cell(num2str(zeros(size(be_coeff{k}))), ones(size(be_coeff{k},1),1)));
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('ga', mat2cell(num2str(zeros(size(ga_coeff{k}))), ones(size(ga_coeff{k},1),1)));
end


%get L2 matrix:
str = mphmatrix(comsol_model,'sol1','Out',{'K','Kc'});
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(ptag).feature(['cfeq' num2str(k)]).set('c', {'1'});
    %comsol_model.physics(ptag).feature('cfeq1').set('c', {'1' '0' '0' '1'});
end
str = mphmatrix(comsol_model,'sol1','Out',{'K','Kc'});
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(ptag).feature(['cfeq' num2str(k)]).set('c', c_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('a', a_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('f', f_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('ea', ea_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('da', da_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('al', al_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('be', be_coeff{k});
    comsol_model.physics(ptag).feature(['cfeq' num2str(k)]).set('ga', ga_coeff{k});
end