ThermalBlock_3D_5x5x5_LinLagr_Einfluss_unten_8Params.m

function out = model
%
% ThermalBlock_3D_5x5x5_LinLagr_Einfluss_unten_8Params.m
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Properties of the model here!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%mu_init = [1,1,1]; % [1,3,8,2] --> sets parameters: mu1=1, mu2=3, mu3=8, mu4=2
mu_init = ones(8,1);

nr_blocks = 125;

%zufallszahlen fuer das modell generieren! verwende seed, damit immer die
%gleichen zahlen beim model generieren verwendet werden
seed = 1;
rng(seed);
block_mu_index  = randi([1 length(mu_init)],nr_blocks,1); %random


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generation of the Comsol Model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');

model.modelPath('/Users/oliverzeeb/Documents/MATLAB/rbmatlab_COMSOL/COMSOL_rbmatlab/models/ThermalBlock/3D_stat');

model.name('ThermalBlock_3D_5x5x5_LinLagr_Einfluss_unten_8Params.mph');

model.baseSystem('none');

for k=1:length(mu_init)
    model.param.set(['mu',num2str(k)], num2str(mu_init(k)));
end

model.modelNode.create('mod1');

model.geom.create('geom1', 3);
%Bloecke generieren
for k=1:nr_blocks
    eval_string = ['model.geom(''geom1'').feature.create(''blk' num2str(k),  ''', ''Block'');'];
    eval(eval_string);
end
% model.geom('geom1').feature.create('blk1', 'Block');
% model.geom('geom1').feature.create('blk2', 'Block');
% ...

%groesse der bloecke setzen
counter = zeros(3,1);
for k=1:nr_blocks
    %position festlegen
    pos_string = ['model.geom(''geom1'').feature(''blk',  num2str(k)  , ''').set(''pos'', {'''...
        num2str(counter(1)) '/5'' ''' num2str(counter(2)), '/5'' ''', num2str(counter(3)) '/5''});'];
    eval(pos_string);

    size_string = ['model.geom(''geom1'').feature(''blk', num2str(k), ''').set(''size'', {''1/5'' ''1/5'' ''1/5''});'];
    eval(size_string);
    
    if k<nr_blocks
        counter = increase_counter(5,1,counter);
    end
end
% model.geom('geom1').feature('blk1').set('pos', {'0' '0' '0'});
% model.geom('geom1').feature('blk1').set('size', {'1/4' '1/4' '1/4'});
% ...

model.geom('geom1').run;

model.physics.create('c', 'CoefficientFormPDE', 'geom1');
for k=2:nr_blocks
    create_string = ['model.physics(''c'').feature.create(''cfeq' num2str(k) ''', ''CoefficientFormPDE'', 3);'];
    eval(create_string);
    selection_string = ['model.physics(''c'').feature(''cfeq' num2str(k) ''').selection.set([' num2str(k) ']);'];
    eval(selection_string);
end
% model.physics('c').feature.create('cfeq2', 'CoefficientFormPDE', 3);
% model.physics('c').feature('cfeq2').selection.set([2]);
% ...

model.physics('c').feature.create('dir1', 'DirichletBoundary', 2);
model.physics('c').feature('dir1').selection.set([16, 32, 48, 64, 80, 101, 117, 133, 149, 165,...
    186, 202, 218, 234, 250, 271, 287, 303, 319, 335, 356, 372, 388, 404, 420]);
model.physics('c').feature.create('flux1', 'FluxBoundary', 2);
model.physics('c').feature('flux1').selection.set([3, 19, 35, 51, 67, 88, 104, 120, 136, 152,...
    173, 189, 205, 221, 237, 258, 274, 290, 306, 322, 343, 359, 375, 391, 407]);

model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').feature.create('ftet1', 'FreeTet');

model.view('view1').set('renderwireframe', true);

model.physics('c').prop('ShapeProperty').set('order', '1');

% koeffizienten der PDE setzen
for k=1:nr_blocks
    set_f_string = ['model.physics(''c'').feature(''cfeq' num2str(k) ''').set(''f'', ''0'');'];
    eval(set_f_string)
    set_da_string = ['model.physics(''c'').feature(''cfeq' num2str(k) ''').set(''da'', ''0'');'];
    eval(set_da_string)
end
model.physics('c').feature('flux1').set('g', '50');
for k=1:length(block_mu_index)
    mu_string = ['mu',num2str(block_mu_index(k))];
    model.physics('c').feature(['cfeq',num2str(k)]).set('c', {mu_string '0' '0' '0' mu_string '0' '0' '0' mu_string});
end
% model.physics('c').feature('cfeq1').set('c', {'mu1' '0' '0' '0' 'mu1' '0' '0' '0' 'mu1'});
% model.physics('c').feature('cfeq1').set('f', '0');
% model.physics('c').feature('cfeq1').set('da', '0');
% ...

%model.mesh('mesh1').run;

model.frame('material1').sorder(1);

model.study.create('std1');
model.study('std1').feature.create('stat', 'Stationary');

model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').attach('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('s1', 'Stationary');
model.sol('sol1').feature('s1').feature.create('fc1', 'FullyCoupled');
model.sol('sol1').feature('s1').feature.remove('fcDef');

% Ausgabefuntional
model.result.numerical.create('int1', 'IntSurface');
model.result.numerical('int1').selection.set([3, 19, 35, 51, 67, 88, 104, 120, 136, 152,...
    173, 189, 205, 221, 237, 258, 274, 290, 306, 322, 343, 359, 375, 391, 407]);
model.result.numerical('int1').set('probetag', 'none');

model.result.create('pg1', 'PlotGroup3D');
model.result('pg1').feature.create('slc1', 'Slice');

model.sol('sol1').feature('st1').name('Compile Equations: Stationary {stat}');
model.sol('sol1').feature('st1').set('studystep', 'stat');
model.sol('sol1').feature('v1').set('control', 'stat');
model.sol('sol1').feature('s1').set('control', 'stat');
%model.sol('sol1').runAll;

model.result('pg1').feature('slc1').set('unit', '1');

out = model;

end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function counter = increase_counter(basis, pos,counter)
    if mod(counter(pos)+1,basis) == 0 %aktuelle pos auf null setzen, nächste pos um 1 erhöhen!
        counter(pos) = 0;
        counter = increase_counter(basis, pos+1,counter);
    else 
        counter(pos) = counter(pos)+1;
    end
end