BoxCrash

Step 1Step 2Step 3Step 4

%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;
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

%---------------delete 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;
load('Mats/Mat_mxRock2.mat');%load the trained material
Mats{2,1}=Mat_mxRock2;
Mat_mxRock2.name='StrongRock';%set the name of the material
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();

%Excavation process example. Please run the user_BoxCrash1.m first
%this code shows how to kill elements
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
B.name='BoxCrashHole';

%---------------delele elements on the top
d.delElement(d.GROUP.topPlaten);%delete elements according to id
d.balanceBondedModel0(0.1);%balance the bonded model without friction

d.mo.bFilter(:)=true;
d.mo.nBondRate=d.mo.nBondRate*0.1;
d.mo.zeroBalance();
d.resetStatus();
%----------make disc sample------------
gpuStatus=d.mo.setGPU('auto');
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
discR=20;
for showCircle=1:30
    mX=d.mo.aX(1:d.mNum);
    mZ=d.mo.aZ(1:d.mNum);
    discCX=showCircle*10;
    discCZ=mean(mZ);
    discFilter=(mX-discCX).^2+(mZ-discCZ).^2<discR^2;
    d.killElement(find(discFilter));%add a new group

    d.balance('Standard',0.1);
    save([fName num2str(showCircle) '.mat']);
    d.show('Displacement');
end

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;

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

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

%----------set the drawing of result during iterations
%ufs.simpleFigure('on');%simplify the figure drawing
showType='StressZZ';
figureNumber=d.show(showType);
d.figureNumber=figureNumber;
%----------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.show(showType);
    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();