BoxWave

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='BoxWave';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=0.002;
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=0.5;%width, length, height, average radius
B.sampleL=0;%when L is zero, it is a 2-dimensional model
B.sampleH=0.25;
B.isSample=1;
B.platenStatus=[0;0;0;0;0;1];%topPlaten
B.BexpandRate=4;%boundary is 4-ball wider than sample
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;

%---------- gravity sedimentation
B.gravitySediment();%you may use B.gravitySediment(10); to increase sedimentation time (10)
B.compactSample(1);%input is compaction time
mfs.reduceGravity(d,3);%reduce the gravity of element
%------------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/BoxWave1.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
sampleFilter=mZ>B.sampleH*0.9|mZ<B.sampleH*0.1;
d.delElement(find(sampleFilter));%delete elements according to id
d.mo.setGPU('off');

%--------------assign new material
matTxt=load('Mats\WeakRock.txt');%load material file
Mats{1,1}=material('WeakRock',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;%assign new material
d.groupMat2Model({'sample'},1);%apply the new material
d.mo.mGZ(:)=0;%no gravitation
d.moveGroup('lefB',-B.sampleW*0.05,0,0);%move the lateral boundaries
d.moveGroup('rigB',B.sampleW*0.05,0,0);
d.removePrestress(0.1);%reduce the prestress
d.mo.setGPU('auto');
d.balance('Standard',2);%balance the model

d.show('mV');
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();

%in this code, we will learn how to generate seismic wave and make
%receivers to record data of seismic wave
clear;
load('TempModel/BoxWave2.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.isHeat=1;%calculate heat in the model
visRate=0.0001;
d.mo.mVis=d.mo.mVis*visRate;
d.mo.isCrack=1;

%---------------------define the source of the wave
%leftBlock is used to generate wave
mX=d.mo.aX(1:d.mNum);
leftFilter=mX<B.sampleW*0.05;
d.addGroup('leftBlock',find(leftFilter));
%generating sine wave on the leftLine group
totalBalanceNum=10000;%data number of the wave
totalT=totalBalanceNum*d.mo.dT;%total time of the wave

period=1e-4;%period of the wave
Ts=(0:totalBalanceNum)*d.mo.dT;
maxA=100;%maximum acceleration
Values=maxA*sin(Ts*(2*pi)/period);

waveProp='mAX';
d.addTimeProp('leftBlock',waveProp,Ts,Values);%assign the AZ to elements of LeftLine
d.addRecordProp('leftBlock',waveProp);%declare recording mAZ
%---------------------end define the source of the wave

%-------------------define the receiver
receiverNum=6;%receiver number
centerx=B.sampleW/receiverNum;%center position of the receiver
centery=0;
centerz=B.sampleH/2;
R=B.ballR*4;%radius of the receiver
gNames={};
prop1='mAX';%record the property 1
prop2='aX';%record the property 2
for i=1:receiverNum
    gName=['Receiver' num2str(i)];
    gNames=[gNames(:);gName];
    f.run('fun/defineSphereGroup.m',d,gName,centerx*i,centery,centerz,R);
    d.addRecordProp(gName,prop1);%declare recording mAZ
    d.addRecordProp(gName,prop2);%declare recording mAZ
end
figure;
subplot(2,1,1);
d.setGroupId();
d.showFilter('Group',gNames,'aR');
subplot(2,1,2);
plot(Ts,Values);xlabel('Ts [second]');ylabel('X acceleration of the leftBlock [m/s^2]');title('Wave on the leftBlock');
%-------------------end define the receiver

gpuStatus=d.mo.setGPU('auto');
totalCircle=30;
d.tic(totalCircle);
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
figure;
d.figureNumber=d.show('mV');
for i=1:totalCircle
    d.mo.setGPU(gpuStatus);
    d.balance('Standard',0.01);
    d.show('mV');
    d.clearData(1);
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.toc();%show the note of time
end
d.mo.setGPU('off');

%show the curves
curveFigure=figure;
for i=1:receiverNum
    subplot(6,2,i*2-1);
    d.status.show(['PROPReceiver' num2str(i) '_' prop1]);
    subplot(6,2,i*2);
    d.status.show(['PROPReceiver' num2str(i) '_' prop2]);
end
set(curveFigure, 'position', get(0,'ScreenSize'));

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();