3DSlope

Step 1Step 2Step 3Step 4Step 5

%converting scattered XYZ data to digital elevation model, which can be
%used to define a surface and to cut discrete element model
clear;
d1=load('slope/XYZData.txt');%scattered XYZ data

X1=d1(:,1);Y1=d1(:,2);Z1=d1(:,3);
[X1,Y1]=mfs.rotateIJ(X1,Y1,90);%rotate the surface data
x1=min(X1);x2=max(X1);%find the span of the data
y1=min(Y1);y2=max(Y1);

%make function for scattered data, see Matlab scatteredInterpolant
F1=scatteredInterpolant(X1,Y1,Z1,'natural','nearest');
gSide=20;%side of grid
[gX,gY]=meshgrid(x1:gSide:x2,y1:gSide:y2);%X, Y mesh coordinates
gZ=F1(gX,gY);%calculate Z by using function of scattered data

figure
surface(gX,gY,gZ,gZ);
fs.general3Dset();
view(120,30);

%-------------This file is used to create surface of model layers
clear;
r=10;%radius of element, in this example, it is between 5-20 m
surfPackNum=4;%number of surface element along depth direction
surfPackNum2=1;%additional surface element

botPackNum=1.5;topPackNum=1.5;%
load('slope/slopeSurface.mat');%you may use user_L3lpsXYZ2Surf.m to generate similar data
%Digital Elevation Model of model is saved in S, see Matlab commands
%"scatteredInterpolant" "meshgrid" "surf"
gSide=S.Y(2)-S.Y(1);
%--------------cut slope data;
cutX1=floor(250/gSide);cutX2=floor(150/gSide);%cut the digital elevation  data
cutY1=floor(400/gSide);cutY2=floor(300/gSide);
S.X=S.X(1+cutY1:end-cutY2,1+cutX1:end-cutX2);
S.Y=S.Y(1+cutY1:end-cutY2,1+cutX1:end-cutX2);
S.Z1=S.Z1(1+cutY1:end-cutY2,1+cutX1:end-cutX2);%Z1 is original height
S.Z2=S.Z2(1+cutY1:end-cutY2,1+cutX1:end-cutX2);%Z2 is the height after landslide
S.X=S.X-min(S.X(:));
S.Y=S.Y-min(S.Y(:));
S.dZ=S.Z2-S.Z1;
%--------------end cut slope data;

[imH,imW]=size(S.X);
S.Z=S.Z1;
S0min.X=S.X;S0min.Y=S.Y;
S0min.Z=min(S.Z1,S.Z2);
S0.X=S.X;S0.Y=S.Y;
S0.Z=S.Z1;

S2=S0;
oldSlopeT=30/3;
slopeOldFilter=mfs.image2RegionFilter('slope/slopeOld.png',imH,imW);
se1=strel('disk',40);%
slopeOldFilter2=imerode(slopeOldFilter,se1);
slopeOldFilter3=imerode(slopeOldFilter2,se1);
%imshow(slopeOldFilter3);return;
S2.Z(slopeOldFilter)=S2.Z(slopeOldFilter)-oldSlopeT;
S2.Z(slopeOldFilter2)=S2.Z(slopeOldFilter2)-oldSlopeT;
S2.Z(slopeOldFilter3)=S2.Z(slopeOldFilter3)-oldSlopeT;

S_bot=mfs.moveMeshGrid(S0min,-r*2*(surfPackNum+botPackNum));
%thickAreaFilter=mfs.loadImageRegion('slope/thickarea.png',imH,imW);
%S_bot.Z(~thickAreaFilter)=S0.Z(~thickAreaFilter)-2*r*2;
S1=mfs.moveMeshGrid(S_bot,r*2*botPackNum);

S_top=mfs.moveMeshGrid(S0,r*2*(topPackNum+surfPackNum2));
topRate=0.3;%add elements on top
dZ=max(S_top.Z(:))-min(S_top.Z(:));
topFilter=S_top.Z>max(S_top.Z(:))-dZ*topRate;
topAddH=100;
topZ=S_top.Z(topFilter);
dTopZ=topAddH*(topZ-min(topZ))/(dZ*topRate);

S_top.Z(topFilter)=S_top.Z(topFilter)+dTopZ;
S_top0=mfs.moveMeshGrid(S_top,-r*2*topPackNum);

S_source=S0;
S_source.Z=S.Z2;
sourceFilter=mfs.image2RegionFilter('slope/slopesource.png',imH,imW);
%imshow(sourceFilter);return;
S_source.Z(~sourceFilter)=max(S_source.Z(:))+100;
S_source.Z(sourceFilter)=S_source.Z(sourceFilter)-10;

V=S0;%original surface
surface(V.X,V.Y,V.Z,-1000*ones(size(V.X)));
hold all;
V=S_source;%slope source
%surface(V.X,V.Y,V.Z,mx.Z2-mx.Z1);
%surface(V.X,V.Y,S.Z2,S.Z2-S.Z1);
V=S1;%surface above bottom surface
%surface(V.X,V.Y,V.Z,50*ones(size(V.X)));
V=S2;%surface of bottom side of soil,soil2
%surface(V.X,V.Y,V.Z,1000*ones(size(V.X)));
V=S_bot;%bottom surface
%surface(V.X,V.Y,V.Z,0*ones(size(V.X)));
V=S_top;%top surface
%surface(V.X,V.Y,V.Z,10*ones(size(V.X)));
V=S_top0;%bottom side of top surface
%surface(V.X,V.Y,V.Z,30*ones(size(V.X)));

save('data/slopeSurface2.mat','S_bot','S_top','S_top0','S0','S1','S2','S_source','r');
%colorbar
fs.general3Dset();

%-------------------user_mxSlope1.m;
%build the geometrical model;
clear;
load('data/slopeSurface2.mat');
fs.randSeed(1);%build random model
B=obj_Box;%build a box object
B.name='3DSlope';
%--------------initial model------------
B.GPUstatus='auto';
boxWidth=2800;boxLength=1500;boxHeight=1800;ballR=r;
distriRate=0.2;%define distribution of ball radius, 
isClump=0;
%--------------initial model------------;
B.isUI=1;%when run the code in UI_command, isUI=1
B.ballR=ballR;
B.isClump=isClump;
B.distriRate=distriRate;
B.sampleW=boxWidth;
B.sampleL=boxLength;
B.sampleH=boxHeight;
B.platenStatus(:)=0;%no platen in the model
B.buildModel();
B.createSample();%create balls in the box
B.sample.R=B.sample.R*2^(1/12);%R deviation between close packing and cube packing is 2^(1/6)
S_Bbot=S_bot;
S_Bbot.Z=S_Bbot.Z-ballR*4;
S_Btop=S_top;%top limit of boundary
S_Btop.Z=S_Btop.Z+1000;%boundary limit greater
B.addSurf(S_bot);%add the bottom surface
B.addSurf(S_top);%add the top surface

B.addSurf(S_Bbot);%add the bot surface of boundary
B.addSurf(S_Btop);%add the top surface of boundary
B.cutGroup({'sample','botB','topB'},1,2);%cut sample and remove top boundary
B.cutGroup({'lefB','rigB','froB','bacB'},3,4);

B.finishModel();%built the geometric model
B.setSoftMat();%set soft balls to increase the speed of computing

%-------------end initial parameters and build model--------------;
%------------generate random-packing balls-------------;
B.d=B.exportModel();%tranform model data to build object
B.d.mo.isShear=0;
%B.R=B.R*0.5;B.show();

d=B.d;d.showB=1;%d.breakGroup('sample');d.breakGroup('lefPlaten');
%d.show('aR');return

C=Tool_Cut(d);%use the tool to cut sample and get layers
C.addSurf(S_bot);
C.addSurf(S1);
C.addSurf(S2);
C.addSurf(S0);
C.addSurf(S_top0);
C.addSurf(S_top);
C.addSurf(S_source);
C.setLayer({'sample'},[1,2,3,4,5,6]);%set layers according geometrical data
%layer1 is between surface 1 and 2, etc

gNames={'layer1';'layer2';'layer3';'layer4';'layer5'};
d.makeModelByGroups(gNames);%make the model by using the layers.
d.defineWallElement('layer1');
d.mo.aR(d.GROUP.layer1)=B.ballR*1.3;

mo=d.mo;
mo.isFix=1;%fix coordinates;
gId_top=d.getGroupId('layer5');%get element Id of group
mo.FixXId=gId_top;%fix X coordinate
mo.FixYId=gId_top;
mo.FixZId=gId_top;

nBall=d.mo.nBall;
bcFilter=sum(nBall>d.mNum&nBall~=d.aNum,2)>0;%boundary connected ball filter
gFilter=zeros(size(bcFilter));
gFilter(gId_top)=1;
mo.aR(gId_top)=B.ballR;
mo.aR(gFilter&(~bcFilter))=B.ballR*1.3;%not a boundary connected ball
d.setClump('layer5');

%d.show('aR');return;

B.uniformGRate=1;%the parameter used in B.gravitySediment
B.gravitySediment(0.5);
d.mo.FixZId=[];%unfix all Z coordinates

d.mo.dT=d.mo.dT*4;
d.balance('Standard');
d.mo.dT=d.mo.dT/4;

%--------------end initial model------------;
%B.compactSample(6);%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('Box3DSlope1Finish');
save(['TempModel/' B.name '1.mat'],'B','d','C');
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/3DSlope1.mat');
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo

%------------set the material of the model;
load('Mats\Mat_mxRock2.mat');%load a trained material
%you may use user_BoxMatTraining to train a material
Mats{1,1}=Mat_mxRock2;
Mats{1,1}.Id=1;
matTxt2=load('Mats\mxSoil.txt');%load a un-trained material file
Mats{2,1}=material('mxSoil',matTxt2,B.ballR);
Mats{2,1}.Id=2;
d.Mats=Mats;

%--------------redefine layers and apply new material
%because layer elements moved, the layers should be redefined.
d.addGroup('slopeBottom',d.GROUP.layer1);
d.protectGroup('slopeBottom','on');%protect the group, of which the element will not be removed
d.delGroup({'layer1';'layer2';'layer3';'layer4';'layer5'});
C.layerNum=0;
%original layer1 is changed to slopeBottom, and be protected, so new layer1
%is between surfaces 1 and 3
C.setLayer({'sample'},[2,3,4,5]);%set layers according to geometrical data
gNames={'layer1';'layer2';'layer3'};
d.makeModelByGroups(gNames);
d.setGroupMat('layer2','mxSoil');
d.groupMat2Model(gNames,1);

%d.showB=3;d.show('aR');return;
%-----------------define wall area to increase speed of code
sX=d.mo.aX(1:d.mNum);sY=d.mo.aY(1:d.mNum);

imH=302;imW=561;%height and width of image
%read the image and change the size,image is in black and white color
regionFilter=mfs.image2RegionFilter('slope/slopepack.png',imH,imW);%white is true
sFilter=mfs.applyRegionFilter(regionFilter,sX,sY);
sFilter=~sFilter;
sId=find(sFilter);
sId(sId>d.mNum)=[];
d.addGroup('slopeWall',sId);

d.defineWallElement('slopeWall');
d.protectGroup('slopeWall','on');
%-----------------end define wall area to increase speed of code
%d.showB=3;d.show('aR');return;
%---------------balance the model
d.balanceBondedModel0(0.5);
d.mo.mVis=d.mo.mVis*5;
d.balanceBondedModel(0.5);%bond all elements and balance the model
d.mo.setGPU('off');

%--------------balance model 2, cut the model to get final shape
d.delGroup({'layer1';'layer2';'layer3'});
C.layerNum=0;
C.setLayer({'sample'},[1,3,4]);%set layers according to geometrical data
C.setLayer({'sample'},[7,4]);%set layers according to geometrical data
gNames={'layer1';'layer2';'layer3'};
d.makeModelByGroups(gNames);

d.setGroupMat('layer2','mxSoil');
d.groupMat2Model(gNames,1);
d.mo.zeroBalance();
%d.showB=3;d.show('aR');return;
%----------------reduce friction and balance model

aKS=d.mo.aKS;
d.mo.aKS(:)=1;d.mo.setKNKS();%reduce the stiffness of element
d.balanceBondedModel0();
d.mo.aKS=aKS;d.mo.setKNKS();%restore the stiffness of element
d.mo.mVis=d.mo.mVis*5;
d.balanceBondedModel();

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

%-------------------user_mxSlope3.m;
clear;
load('TempModel/3DSlope2.mat');
d.calculateData();
B.setUIoutput();
d.mo.setGPU('off');
%d.showB=3;d.show('StressXX');return;

%setup of calculation
d.getModel();%d.setModel();%reset the initial status of the model
d.mo.aHeat(:)=0;
d.status=modelStatus(d);%initialize model status, which records running information
d.breakGroupOuter({'layer3'});%break the outer connection of the group
d.breakGroup({'layer3'});%break the outer connection of the group

d.mo.isHeat=1;%calculate heat in the model
visRate=0.00001;%ok
d.mo.mVis=d.mo.mVis*visRate;
d.setStandarddT();

fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
gpuStatus=d.mo.setGPU('auto');
totalCircle=10;
d.tic(totalCircle);
for i=1:totalCircle
    d.balance('Standard');
    d.clearData(1);
    d.mo.setGPU('off');
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.mo.setGPU(gpuStatus);
    d.toc();%show the note of time;
end
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Box3DSlope3Finish');
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();
d.show('aR');