function llmatlb

clc;
clear;
close all

model=mphload('domain_activation_llmatlab.mph');

domInd=[2 3 4 5];

ht=model.physics('ht');

for i=1:8
    k=mod(i,4);
    
    if k==0
        k=4;
    end
    
    ht.selection.set([1 domInd(k)]);
    ht.feature('init2').selection.set(domInd(k));
    
    model.study('std1').run;
    
    if i==1
        cpt1=model.result.dataset.create('cpt1','CutPoint3D');
        cpt1.set('pointx','0 L/2 L');
        cpt1.set('pointy','0 L/2 L');
        cpt1.set('pointz','L/10');
        pg1=model.result.create('pg1','PlotGroup1D');
        pg1.set('data','cpt1');
        ptgr1=model.result('pg1').feature.create('ptgr1','PointGraph');
        ptgr1.set('lenged','on');
        
        pg2=model.result.create('pg2','PlotGroup3D');
        surf1=pg2.feature.create('surf1','Surface');
        surf1.set('rangecoloractive','on');
        surf1.set('rangecolormax','336');
        surf1.set('rangecolormin','293.15');
        
        ht.feature('init1').set('T',1,'T');
        
        v1=model.sol('sol1').feature('v1');
        v1.set('initsol','sol1');    
    end 

    figure(1)
    mphplot(model,'pg1','rangenum',1);
    hold on

    figure(2)
    subplot(4,2,i);
    mphplot(model,'pg2');
    
    time=mphglobal(model,'t','solnum','end');
    model.param.set('t0',time);
    
    disp(sprintf('End of iteration No. %d', i));    
end