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/lps1.txt');%scattered XYZ data before sliding
d2=load('slope/lps2.txt');%scattered XYZ data after sliding
gSide=5;%side of grid
baseZ=1150;%the base level of the data

%defined the coordinates of the digital elevation data
X1=d1(:,1);Y1=d1(:,2);Z1=d1(:,3)-baseZ;
X2=d2(:,1);Y2=d2(:,2);Z2=d2(:,3)-baseZ;
angle=-70;%rotate the data
[X1,Y1]=mfs.rotateIJ(X1,Y1,angle);
[X2,Y2]=mfs.rotateIJ(X2,Y2,angle);

%move the data
minX=min(X1);minY=min(Y1);
X1=X1-minX;Y1=Y1-minY;
X2=X2-minX;Y2=Y2-minY;
minX=min(X1);minY=min(Y1);
maxX=max(X1);maxY=max(Y1);

%interpret 3D surface according to the data
F1=scatteredInterpolant(X1,Y1,Z1,'natural','nearest');
F2=scatteredInterpolant(X2,Y2,Z2,'natural','nearest');

%get the mesh of the surface
[gX,gY]=meshgrid(minX:gSide:maxX,minY:gSide:maxY);
gZ1=F1(gX,gY);
gZ2=F2(gX,gY);
dgZ=gZ2-gZ1;
S.X=gX;S.Y=gY;
S.Z1=gZ1;S.Z2=gZ2;
save('TempModel/lps_Slope.mat');

figure
surface(gX,gY,gZ1,dgZ);
fs.general3Dset();
colorbar

%-------------This file is used to create surface of model layers
%-----basic setup of the layer surface
clear;
ballR=15;%radius of element, it is between 5-20 m. default value is 5, use 15 to increase speed
additionalDepth=15;%determine additional elements
surfPackNum=ceil(additionalDepth/(ballR*2));%number of surface element along depth direction

surfPackNum2=1;%additional surface element
shellTNum=6;%thickness of wall elements
botPackNum=2;
topPackNum=1.5;%elements of bottom and top shell
%-----end basic setup of the layer surface

%------------record the basic setup parameters
C=Tool_Cut();%use the tool to cut sample and get layers
C.SET.additionalDepth=additionalDepth;%additional elements
C.SET.surfPackNum=surfPackNum;%number of surface element along depth direction

C.SET.surfPackNum2=surfPackNum2;%additional surface element
C.SET.shellTNum=shellTNum;%thickness of wall elements
C.SET.botPackNum=botPackNum;
C.SET.topPackNum=topPackNum;%elements of bottom and top shell
C.SET.ballR=ballR;
%------------end record the basic setup parameters

load('TempModel/lps_Slope.mat');%run user_L3lpsXYZ2Surf.m to get the 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(0/gSide);cutX2=floor(120/gSide);%cut the digital elevation  data
cutY1=floor(160/gSide);cutY2=floor(200/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;
S0.X=S.X;S0.Y=S.Y;S0.Z=S.Z1;
S0.name='S0';
S0min=S;S0min.name='S0min';
S0min.Z=min(S.Z1,S.Z2);

S_bot=mfs.moveMeshGrid(S0min,-ballR*2*(botPackNum+surfPackNum));
S_bot.name='Sbot';
S_bot0=mfs.moveMeshGrid(S_bot,ballR*2*botPackNum);%i.e. S1 of previous code
S_bot0.name='Sbot0';

S_top=mfs.moveMeshGrid(S0,ballR*2*(topPackNum+surfPackNum2));
S_top.name='Stop';
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,-ballR*2*topPackNum);
S_top0.name='Stop0';

S_source=S0;
S_source.Z=S.Z2;
[imH,imW]=size(S.X);
sourceFilter=mfs.image2RegionFilter('slope/lps_slopesource.png',imH,imW);
%imshow(sourceFilter);return;
S_source.Z(~sourceFilter)=max(S_source.Z(:))+100;
S_source.Z(sourceFilter)=S_source.Z(sourceFilter)-ballR;
S_source.name='Ssource';

C.addSurf(S_bot,'S_bot');
C.addSurf(S_bot0,'S_bot0');
C.addSurf(S0,'S0');
C.addSurf(S_top0,'S_top0');
C.addSurf(S_top,'S_top');
C.addSurf(S_source,'S_source');
C.addSurf(S0min,'S0min');

C.SET.sampleThickness=max(abs(S.dZ(:)))+additionalDepth;
C.SET.S=S;

V=S_source;surface(V.X,V.Y,S.Z2,S.Z2-S.Z1);fs.general3Dset();colorbar

save('TempModel/lps_slopeSurface2.mat','C');

%-------------------user_mxSlope1.m;
%build the geometrical model;
clear;
load('TempModel/lps_slopeSurface2.mat');

fs.randSeed(1);%build random model
B=obj_Box;%build a box object
B.name='lps_3DSlope';
%--------------initial model------------
B.GPUstatus='auto';
S_top=C.SurfData.S_top;
boxWidth=max(S_top.X(:));
boxLength=max(S_top.Y(:));
boxHeight=max(S_top.Z(:))*1.1;%increase the height of the model
ballR=C.SET.ballR;
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.SET.shellHRate=0.1;%record the ratio of shell height to the height of 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=C.SurfData.S_bot;
S_Bbot.Z=S_Bbot.Z-ballR*4;
S_Btop=C.SurfData.S_top;%top limit of boundary
S_Btop.Z=S_Btop.Z+(max(S_top.Z)-min(S_top.Z))*0.8;%boudary limit greater

B.addSurf(C.SurfData.S_bot);%add the bottom surface
B.addSurf(C.SurfData.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 platen
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.d=d;
gNames={'botShell','slopeBody','addedBody','topShell'};
I=C.SurfId;
C.setLayer({'sample'},[I.S_bot,I.S_bot0,I.S0,I.S_top0,I.S_top],gNames);%set layers according geometrical data
%layer1 is between surface 1 and 2, etc
d.makeModelByGroups(gNames);%make the model by using the layers.
%defined bottom shell
d.defineWallElement('botShell');
d.mo.aR(d.GROUP.botShell)=B.ballR*1.4;

%--------------defined top shell
mo=d.mo;
mo.isFix=1;%fix coordinates;
gId_top=d.getGroupId('topShell');%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=false(size(bcFilter));
gFilter(gId_top)=true;
mo.aR(gId_top)=B.ballR;
mo.aR(gFilter&(~bcFilter))=B.ballR*1.3;%not a boundary connected ball
d.setClump('topShell');
%--------------end defined top shell

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

B.uniformGRate=1;%the parameter used in B.gravitySediment
sampleRate=C.SET.sampleThickness/B.sampleH;
B.gravitySediment(sampleRate);
d.mo.FixZId=[];%unfix all Z coordinates
%d.show();return;

d.mo.dT=d.mo.dT*4;%increase the step time to increase the computing speed
d.balance('Standard',sampleRate*2);
d.mo.dT=d.mo.dT/4;
%------------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('Step1Finished');
save(['TempModel/' B.name '1.mat'],'B','d','C','-v7.3');
save(['TempModel/' B.name '1R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat'],'-v7.3');
d.calculateData();

%set the material of the model
clear;
load('TempModel/lps_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;
matTxt2=load('Mats\lpsSoil.txt');%load a un-trained material file
Mats{1,1}=material('lpsSoil',matTxt2,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;

%--------------redefine layers and apply new material
%because layer elements moved, the layers should be redefined.
d.protectGroup('botShell','on');%protect the group, of which the element will not be removed in makeModelByGroups
gNames={'slopeBody','addedBody','topShell'};
d.delGroup(gNames);
%original layer1 is changed to slopeBottom, and be protected, so new layer1
%is between surfaces 1 and 3
gNames={'slopeBody','addedBody'};
ID=C.SurfId;
C.setLayer({'sample'},[ID.S_bot0,ID.S0,ID.S_top0],gNames);%set layers according to geometrical data
d.makeModelByGroups(gNames);

%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=256;imW=596;%height and width of image
%read the image and change the size,image is in black and white color
regionFilter=mfs.image2RegionFilter('slope/lps_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');
%-----------------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
C.layerNum=0;
ID=C.SurfId;
C.setLayer({'sample'},[ID.S_bot,ID.S0],{'slopeBody'});%set layers according to geometrical data
C.setLayer({'sample'},[ID.S_source,ID.S0],{'slopeSource'});%set layers according to geometrical data
C.setLayer({'slopeWall'},[ID.S_bot,ID.S0],{'slopeWall'});%set layers according to geometrical data
gNames={'slopeBody';'slopeSource';'slopeWall'};

d.makeModelByGroups(gNames);

%---------------added in MatDEM_v1.42, reduce wall to shell
clearWall=1;%if the code result in error, change it to 0
if clearWall==1
    %--------------------------added in MatDEM_v1.42
    S0_shell=mfs.moveMeshGrid(C.SurfData.S0,-C.SET.ballR*2*C.SET.shellTNum);
    %--------------------------end added in MatDEM_v1.42
    C.addSurf(S0_shell,'S0_shell');

    delFilter=false(d.aNum,1);
    delFilter(d.GROUP.slopeWall)=true;
    delFilter(d.GROUP.botShell)=true;

    dilatedTNum=C.SET.shellTNum;
    se1=strel('disk',ceil(B.ballR*2*dilatedTNum*imW/B.sampleW));%dilated pixel
    protectedRegionFilter=imdilate(regionFilter,se1);
    protectedFilter=mfs.applyRegionFilter(protectedRegionFilter,d.mo.aX,d.mo.aY);

    ID=C.SurfId;
    C.setLayer({'slopeWall','botShell'},[ID.S0_shell,ID.S0],{'surfaceShell'});%set layers according to geometrical data
    protectedFilter(d.GROUP.surfaceShell)=true;

    delFilter(protectedFilter)=false;
    d.delElement(find(delFilter));
end
%---------------end added in MatDEM_v1.42
d.setGroupMat('slopeBody','lpsSoil');
d.groupMat2Model(gNames,1);
d.mo.zeroBalance();
%d.showB=3;d.show('Displacement');return;
%----------------reduce friction and balance model

visRate=B.sampleH/C.SET.sampleThickness;
mVis=d.mo.mVis;
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*visRate;
d.balanceBondedModel();

d.mo.mVis=d.mo.mVis*C.SET.shellTNum;
for i=1:5
d.balance(400);
d.mo.mVX(:)=0;d.mo.mVY(:)=0;d.mo.mVZ(:)=0;
d.mo.bFilter(:)=true;
d.mo.zeroBalance();
end
d.mo.mVis=mVis;

%-----------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','-v7.3');
save(['TempModel/' B.name '2R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat'],'-v7.3');
d.calculateData();

%-------------------user_mxSlope3.m;
clear;
load('TempModel/lps_3DSlope2.mat');
d.calculateData();
B.setUIoutput();
d.mo.setGPU('off');

%-----------setup of the 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('slopeSource');%break the outer connection of the group
d.breakGroup('slopeSource');%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;%use lower viscosity
d.setStandarddT();%reset the step time
d.mo.dT=d.mo.dT*4;%increase step time to increase the computing speed
%-----------end setup of the calculation

%---------setup of numerical simulation
fName=['data/step/' B.name num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);
gpuStatus=d.mo.setGPU('auto');
totalCircle=50;
d.tic(totalCircle);
%---------start the numerical simulation
for i=1:totalCircle
    d.balance('Standard',0.5);
    d.clearData(1);
    d.mo.setGPU('off');
    save([fName num2str(i) '.mat'],'-v7.3');
    d.calculateData();
    d.mo.setGPU(gpuStatus);
    d.toc();%show the note of time;
end

%----------save the data
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('Box3DSlope3Finish');
save(['TempModel/' B.name '3.mat'],'d','-v7.3');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat'],'-v7.3');
d.calculateData();
d.show('aR');