BoxCrash

Step 1Step 2Step 3

%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box;%declare a box object
B.name='BoxCrash';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=5*0.1;
B.isShear=0;
B.isClump=0;%if isClump=1, particles are composed of several balls
B.distriRate=0.2;%define distribution of ball radius, 
B.sampleW=500;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=300;
B.BexpandRate=4;%boundary is 4-ball wider than sample
B.PexpandRate=0;
B.type='topPlaten';%add a top platen to compact model
B.isSample=1;
%B.type='TriaxialCompression';
B.setType();
B.buildInitialModel();%B.show();
B.setUIoutput();

d=B.d;%d.breakGroup('sample');d.breakGroup('lefPlaten');
%you may change the size distribution of elements here, e.g. d.mo.aR=d.aR*0.95;

d.mo.setGPU('auto');
%1. The simpleContact is the default contact model
%d.mo.balanceCommand='ContactModel.simpleContact(obj);';
%2. The normalContact only consider the normal force of element
%d.mo.balanceCommand='ContactModel.normalContact(obj);';
%3. The normalContact model is defined in a function file
d.mo.setBalanceFunction('fun/normalContact.m');%user-defined normal model
%--------------end initial model------------

%---------- gravity sedimentation
B.gravitySediment();%you may use B.gravitySediment(10); to increase sedimentation time (10)
%d.show('mV');return;
%B.compactSample(1);%input is compaction time
%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model

d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('Step1Finish');
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();

%set the material of the model
clear
load('TempModel/BoxCrash1.mat');
B.setUIoutput();%set the output
d=B.d;
d.calculateData();%calculate data
d.mo.setGPU('off');%close the GPU calculation
d.getModel();%get xyz from d.mo
%---------------delele elements on the top
mZ=d.mo.aZ(1:d.mNum);%get the Z of elements
topLayerFilter=mZ>max(mZ)*0.5;
d.delElement(find(topLayerFilter));%delete elements according to id

%--------------assign new material
matTxt=load('Mats\WeakRock.txt');%load material file
Mats{1,1}=material('WeakRock',matTxt,B.ballR);
Mats{1,1}.Id=1;
matTxt2=load('Mats\StrongRock.txt');
Mats{2,1}=material('StrongRock',matTxt2,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;%assign new material
d.groupMat2Model({'sample'},1);%apply the new material

%----------make disc sample------------
sampleId=d.GROUP.sample;
sX=d.aX(sampleId);sZ=d.aZ(sampleId);sR=d.aR(sampleId);
discCX=mean(sX);discCZ=mean(sZ);
discR=20;
discFilter=(d.aX-discCX).^2+(d.aZ-discCZ).^2<discR^2;
%d.setData();d.data.showFilter=discFilter;d.show('aR');
d.addGroup('Disc0',find(discFilter));%add a new group
discObj=d.group2Obj('Disc0');

discId=d.addElement('StrongRock',discObj);%mat Id, obj
d.addGroup('Disc',discId);%add a new group
disZ=max(sZ+sR)-min(discObj.Z-discObj.R);
d.moveGroup('Disc',0,0,disZ);%move the group

d.balanceBondedModel0();%balance the bonded model without friction
d.show('aMatId');
d.clearData(1);%clear dependent data
d.recordCalHour('Step2Finish');
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();

ufs.setTitle('MatDEM陨石撞击地面实时模拟');
clear;
load('TempModel/BoxCrash2.mat');
d.calculateData();
d.mo.setGPU('off');
B.setUIoutput();
d=B.d;
d.getModel();%d.setModel();%reset the initial status of the model
d.resetStatus();%initialize model status, which records running information
d.mo.isShear=0;
d.mo.isHeat=1;%calculate heat in the model
visRate=0.0001;
d.mo.mVis=d.mo.mVis*visRate;
discId=d.GROUP.Disc;
d.setStandarddT();
d.mo.dT=d.mo.dT*4;%increase step time to increase computing speed

d.mo.mVZ(discId)=-2000;
d.mo.mVX(discId)=-1000;

d.showB=0;
d.status.legendLocation='West';

%----------set the drawing of result during iterations
setappdata(0,'simpleFigure',1);
setappdata(0,'ballRate',0.01);
showType='StressZZ';
%----------end set the drawing of result during iterations

gpuStatus=d.mo.setGPU('auto');
totalCircle=50;
d.tic(totalCircle);
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
for i=1:totalCircle
    d.balance('Standard',0.0005);
    d.figureNumber=d.show(showType);%result will be shown in one figure
    pause(0.05);
end
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxCrush3Finish');
save(['TempModel/' B.name '3.mat'],'d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();