MonteCarlo

Step 1Step 2Step 3Step 4Step 5

%build the geometrical model
clear;
fs.randSeed(2);%build random model
B=obj_Box;%build a box object
B.name='BoxMonteCarlo';
B.GPUstatus='auto';
B.ballR=1;
B.isClump=0;
B.distriRate=0.2; 
B.sampleW=60;
B.sampleL=0;
B.sampleH=60;
B.type='topPlaten';
B.setType();
B.buildInitialModel();%B.show();
B.setUIoutput();
d=B.d;

%--------------end initial model------------
B.gravitySediment(0.1);
B.compactSample(1);%input is compaction time
%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model

d.clearData(1);%clear dependent data
d.recordCalHour('BoxModel1Finish');
save(['TempModel/' B.name '1.mat'],'B','d');
save(['TempModel/' B.name '1R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
d.show('aR');

%set the material of the model
clear
load('TempModel/BoxMonteCarlo1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo

%---------cut the model to make slope
C=Tool_Cut(d);%cut the model
lSurf=load('slope/layer surface.txt');%load the surface data
C.addSurf(lSurf);%add the surfaces to the cut
C.setLayer({'sample'},[1,2,3,4]);%set layers according geometrical data
gNames={'lefPlaten';'rigPlaten';'botPlaten';'layer1';'layer2';'layer3'};
d.makeModelByGroups(gNames);
%---------end cut the model to make slope

%----------set material of model
matTxt=load('Mats\Soil1.txt');matTxt(3:4)=matTxt(3:4)*10;
Mats{1,1}=material('Soil1',matTxt,B.ballR);
Mats{1,1}.Id=1;
matTxt2=load('Mats\Soil2.txt');matTxt2(3:4)=matTxt2(3:4)*10;
Mats{2,1}=material('Soil2',matTxt2,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;
%----------end set material of model

%---------assign material to layers and balance the model
d.setGroupMat('layer2','Soil2');
d.groupMat2Model({'sample','layer2'});
d.balanceBondedModel();
d.show('aR');
%---------end assign material to layers and balance the model

%---------save the data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxStep2Finish');
save(['TempModel/' B.name '2.mat'],'B','d');
save(['TempModel/' B.name '2R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

clear;
load('TempModel/BoxMonteCarlo2.mat');
d2=B.d;
load('TempModel/BoxMonteCarlo2.mat');
B.setUIoutput();

modelNum=5;%defines the number of slope in one model
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%d.setModel();%reset the initial status of the model
d.status=modelStatus(d);%initialize model status, which records running information

d.is2D=0;
disY=max(d.aR)*4;
d.GROUP.sample1=d.GROUP.sample;
d.SET.modelNum=modelNum;
%the slope model is added one by one here @@@@
%we can make a new object including all slopes, and add it to d
for i=2:modelNum
    newSampleId=mfs.addModel(d,d2,0,disY*(i-1),0);
    d.GROUP.(['sample' num2str(i)])=newSampleId;
end
%---------cut the model to make slope
C=Tool_Cut(d);%cut the model
lSurf=load('slope/layer surface.txt');%load the surface data
C.addSurf(lSurf);%add the surfaces to the cut
C.setLayer({'sample'},[1,2,3,4]);%set layers according geometrical data
%---------end cut the model to make slope

%---------assign material to layers and balance the model
d.setGroupMat('layer2','Soil2');
d.groupMat2Model({'sample','layer2'});
d.balanceBondedModel();
d.show('aR');
%---------end assign material to layers and balance the model

d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxStep3Finish');
save(['TempModel/' B.name '2_3.mat'],'B','d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

%clear;
fs.randSeed(1);%
totalCircle=2;%define the number of random simulation
load('TempModel/BoxMonteCarlo2_3.mat');
%----------define the parameters of the slopes
aBFRates=(2+rand(d.SET.modelNum,1)*8)/20;%rates of d.mo.aBF, tensile strength
aMUpRates=(2+rand(totalCircle,1)*8)/20;%rate of d.mo.aMUp, friction coefficient

fName=['data/step/' B.name  num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
for i=1:totalCircle
    load('TempModel/BoxMonteCarlo2_3.mat');
    B.setUIoutput();
    d=B.d;
    d.calculateData();
    d.mo.setGPU('off');
    d.getModel();%d.setModel();%reset the initial status of the model
    d.resetStatus();%initialize model status, which records running information

    d.mo.isHeat=1;%calculate heat in the model
    visRate=0.001;
    d.mo.mVis=d.mo.mVis*visRate;
    d.setStandarddT();
    sId=d.GROUP.sample;
    d.mo.aMUp(sId)=d.mo.aMUp(sId)*aMUpRates(i);%assign a different friction coefficient@@@

    d.SET.aMUpRate=aMUpRates(i);%record the setup of this model
    d.SET.aBFRates=aBFRates;
    for j=1:d.SET.modelNum
        sName=['sample' num2str(j)];
        sId=d.GROUP.(sName);
        d.mo.aBF(sId)=d.mo.aBF(sId)*aBFRates(j);%assign a different break force (tensile strength)@@@
    end
    d.mo.setGPU('auto');
    d.balance('Standard',0.5);%we may reduce the value to 0.2 to increase the speed

    d.clearData(1);%save the data
    save([fName num2str(i) '.mat']);
    d.calculateData();
end

%show the displacement of each slope, and mark them
d.show('Displacement');
title(['Displacement, aMUpRate: ' num2str(d.SET.aMUpRate)]);
for i=1:d.SET.modelNum
    maxSZ=max(d.mo.aZ(d.GROUP.sample));
    sY=d.aY(d.GROUP.(['sample' num2str(i)]));
    text(0.05*maxSZ,sY(1),0.1*maxSZ+maxSZ,num2str(i));
end

%the file is used to analyse the data of MonteCarlo3
clear
MCdatas=[];
for dataI=1:100
    load(['data/step/BoxMonteCarlo0.2-0.2loopNum' num2str(dataI) '.mat']);%load the saved file
    d.setData();
    movedRates=zeros(d.SET.modelNum,1);
    for i=1:d.SET.modelNum
        sampleId=d.GROUP.(['sample' num2str(i)]);
        sDisplacement=d.data.Displacement(sampleId);
        sFilter=sDisplacement>B.ballR*2;
        movedRates(i)=sum(sFilter)/length(sFilter);
    end
    movedRateLimit=0.01;
    movedFilter=movedRates>movedRateLimit;

    MCdata=[aBFRates,d.SET.aMUpRate*ones(size(aBFRates)),movedFilter];
    MCdatas=[MCdatas;MCdata];
end
save('slope/MCdatas3.mat','MCdatas');