Skip to content

Commit

Permalink
Bits for the seismology students
Browse files Browse the repository at this point in the history
  • Loading branch information
@Timmmdavis
Timmmdavis committed Aug 30, 2019
1 parent dd54981 commit 4070265
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 106 deletions.
28 changes: 14 additions & 14 deletions 2dCode/MainFrame2D.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,13 @@
%Option A = Loading XY ascii data or manually creating fractures.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % %Single flat user defined fracture
% x=linspace(-0.5,0.5,15);
% y=zeros(1,numel(x));
% Pointsxy=[x;y]';
% mystruct.line1=(1:(length(Pointsxy(:,1))));
% Fnms=fieldnames(mystruct);
% nf=numel(Fnms);clear Fnms
% %Single flat user defined fracture
x=linspace(-0.5,0.5,15);
y=zeros(1,numel(x));
Pointsxy=[x;y]';
mystruct.line1=(1:(length(Pointsxy(:,1))));
Fnms=fieldnames(mystruct);
nf=numel(Fnms);clear Fnms



Expand All @@ -65,10 +65,10 @@
%Loads Points that represent the file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Pointsxy,mystruct]=m_shaperead('FracturesStoneHaven');
Fnms=fieldnames(mystruct);
nf=numel(Fnms);clear Fnms
Pointsxy=[Pointsxy(:,1),Pointsxy(:,3)];
% [Pointsxy,mystruct]=m_shaperead('FracturesStoneHaven');
% Fnms=fieldnames(mystruct);
% nf=numel(Fnms);clear Fnms
% Pointsxy=[Pointsxy(:,1),Pointsxy(:,3)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expand Down Expand Up @@ -280,10 +280,10 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Define the number of cells on the square grid
cells=16;
cells=25;
%How many extra cells you want to add away from the faults,
%reduce to 0 if having half space issues
padding=10;
padding=20;
% Create a vector that is just all the points in X Y.
PointsXY=[[P1(:,1);P2(:,1)],[P1(:,2);P2(:,2)]];
[maxgriX,mingriX,maxgriY,mingriY,sz]=MinMaxDataExtents2d(PointsXY,cells,padding);
Expand Down Expand Up @@ -343,7 +343,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Drawing and fixing Obs Point data just created.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dist=500;
Dist=0.05;
[X,Y]=NanTolDistLine2Pnt( X,Y,P1,P2,LineNormalVector,Dist );
%Drawing these and the location compared to the Obs Points
figure;
Expand Down
216 changes: 124 additions & 92 deletions 3dCode/AnalyticalTests/ComparisonScripts/Okada1985_RectangularDislocation_Test.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
%Option A = Loading Gocad ascii data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


string='OkadaFault_ON_2Faces.ts';
[ Points,Triangles ] = GoCadAsciiReader( string );

Expand All @@ -46,6 +47,8 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%


SecondSurface=0;
% Defining halfspace. 1 is on, 0 is off. The code will run much faster with
% this off as the calculations are simpler in a full space.
halfspace = 1;
Expand Down Expand Up @@ -111,6 +114,26 @@
rowcount = length(X(:,1));
colcount = length(X(1,:));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Option D = Load a secondary Fault Surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %Secondary Surface
% SecondSurface=1; %Flag set
% string='GoCadExportHector.ts';
% [ PointsObs,TrianglesObs ] = GoCadAsciiReader( string );
% PointsObs(:,2:4)=PointsObs(:,2:4)/0.5E4;
% PointsObs=[PointsObs(:,1),PointsObs(:,2)-max(PointsObs(:,2)),PointsObs(:,3)-max(PointsObs(:,3)),(PointsObs(:,4)-freesurface_height)];
% hold on; trisurf(TrianglesObs,PointsObs(:,2),PointsObs(:,3),PointsObs(:,4),'FaceAlpha',(.8),'facecolor', 'cyan');
% [MidPointObs,FaceNormalVectorObs] = MidPointCreate(PointsObs,TrianglesObs,0);
% X=MidPointObs(:,1);
% Y=MidPointObs(:,2);
% Z=MidPointObs(:,3);
% rowcount = length(X(:,1));
% colcount = length(X(1,:));


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%STEP 5: Calculate Stresses/Disps at defined observation points in XYZ.
Expand Down Expand Up @@ -139,99 +162,108 @@
[Ux,Uy,Uz] = CalculateDisplacementOnSurroundingPoints3d...
(Dss,Dds,Dn,nu,X,Y,Z,P1,P2,P3,halfspace);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%STEP 7: VISUALISATION AND ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if SecondSurface==0
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%STEP 7: VISUALISATION AND ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Calculating the depth of the centre of the fault using Pythagoras theorem
dip=70;
width=5;
TipDepth=1;
MidDeptha=sind(dip)*width;
MidDepth=MidDeptha/2+TipDepth;
%Running the Okada function, see: http://uk.mathworks.com/MATLABcentral/fileexchange/25982-okada--surface-
%deformation-due-to-a-finite-rectangular-source
%For details
[uE,uN,uZ,uNN,uNE,uEN,uEE] = Okada1985_RectangularDislocation(X,Y,MidDepth,60,dip,5,width,90,1,0,'plot') ;
hold on

uxy=-(uEN(:)+uNE(:))/2;

%Tri dislocations
[E1,E2,E1dir,E2dir]=EigCalc2d(exx,eyy,exy);
%Okada dislocations
[E1ok,E2ok,E1dirok,E2dirok]=EigCalc2d(-uEE(:),-uNN(:),-uxy);

[E1,E2,E1ok,E2ok ]=ReshapeData2d( rowcount,colcount,E1,E2,E1ok,E2ok );


%Plotting the displacement induced in the TDE code and its fault
figure;subplot(1,2,1),quiver3(X(:),Y(:),Z(:),Ux,Uy,Uz)
xlabel('x'); ylabel('y'); axis('equal'); title('DisplacementAndFaultTDECode')
hold on
trisurf(Triangles,Points(:,2),Points(:,3),Points(:,4),'FaceColor', 'cyan', 'faceAlpha',0.8);
hold off
%Drawing the displacement vectors induced by the Okada function
subplot(1,2,2),quiver3(X,Y,Z,uE,uN,uZ)
xlabel('x'); ylabel('y'); axis('equal'); title('Displacement OkadaCode')
hold off

%Reshaping vectors before calculating residual
Ux=reshape(Ux,(size(uE)));
Uy=reshape(Uy,(size(uE)));
Uz=reshape(Uz,(size(uE)));

%Calculating residuals
Resx=uE-Ux;
Resy=uN-Uy;
Resz=uZ-Uz;

figure;subplot(1,3,1),contourf(X,Y,Resx);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual ux');colorbar;
subplot(1,3,2),contourf(X,Y,Resy);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual uy');colorbar;
subplot(1,3,3),contourf(X,Y,Resz);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual uz');colorbar;


%Plotting Tri Dis displacements vs the Okada displacements
DrawDeformedGrid2d( X,Y,Ux,Uy,cmap2,(E1+E2),'Scale',5);
hax1=gca;no1=gcf;title('Triangular dislocations');
DrawDeformedGrid2d( X,Y,uE,uN,cmap2,(E1+E2),'Scale',5);
hax2=gca;no2=gcf;title('Rectangular dislocation');
hf2=figure;
s1=subplot(1,2,1);
s2=subplot(1,2,2);
pos1=get(s1,'Position');pos2=get(s2,'Position');
delete(s1);delete(s2);
hax_1=copyobj(hax1,hf2);
hax_2=copyobj(hax2,hf2);
set(hax_1, 'Position', pos1);colorbar;colormap(cmap2);
xlabel('x'); ylabel('y');gca;% title('Triangles');
set(hax_2, 'Position', pos2);colorbar;colormap(cmap2);
xlabel('x'); ylabel('y');gca;%title('gggg');
close(no1);close(no2);
titlesz=25;
fntsz=21;
ChangeFontSizes(fntsz,titlesz);
text(-25,13,'Ground displacements and dilatation, comparison of different analyical sources','Fontsize',25,'FontWeight','bold')

bad=isnan(Resx);Resx=Resx(~bad);%removing nans to get proper mean
fprintf('MaxResidualUx %i.\n',max(abs(Resx(:))))

bad=isnan(Resy);Resy=Resy(~bad);%removing nans to get proper mean
fprintf('MaxResidualUy %i.\n',max(abs(Resy(:))))

bad=isnan(Resz);Resz=Resz(~bad);%removing nans to get proper mean
fprintf('MaxResidualUz %i.\n',max(abs(Resz(:))))

P=max(abs(Resx(:)))+max(abs(Resy(:)))+max(abs(Resz(:)));
if any(P>1e-8) %Around the precision and errors introduced. Small values.
error('Calculated displacements are a long way from the analytical solutions')
else
disp('Everything looks good, tolerance checks max stress residuals and flags errors above 1e-8')
end

%Calculating residuals
Resex=exx--uEE(:);
Resey=exy-(-(uEN(:)+uNE(:))/2);

%Calculating the depth of the centre of the fault using Pythagoras theorem
dip=70;
width=5;
TipDepth=1;
MidDeptha=sind(dip)*width;
MidDepth=MidDeptha/2+TipDepth;
%Running the Okada function, see: http://uk.mathworks.com/MATLABcentral/fileexchange/25982-okada--surface-
%deformation-due-to-a-finite-rectangular-source
%For details
[uE,uN,uZ,uNN,uNE,uEN,uEE] = Okada1985_RectangularDislocation(X,Y,MidDepth,60,dip,5,width,90,1,0,'plot') ;
hold on

uxy=-(uEN(:)+uNE(:))/2;

%Tri dislocations
[E1,E2,E1dir,E2dir]=EigCalc2d(exx,eyy,exy);
%Okada dislocations
[E1ok,E2ok,E1dirok,E2dirok]=EigCalc2d(-uEE(:),-uNN(:),-uxy);

[E1,E2,E1ok,E2ok ]=ReshapeData2d( rowcount,colcount,E1,E2,E1ok,E2ok );


%Plotting the displacement induced in the TDE code and its fault
figure;subplot(1,2,1),quiver3(X(:),Y(:),Z(:),Ux,Uy,Uz)
xlabel('x'); ylabel('y'); axis('equal'); title('DisplacementAndFaultTDECode')
hold on
trisurf(Triangles,Points(:,2),Points(:,3),Points(:,4),'FaceColor', 'cyan', 'faceAlpha',0.8);
hold off
%Drawing the displacement vectors induced by the Okada function
subplot(1,2,2),quiver3(X,Y,Z,uE,uN,uZ)
xlabel('x'); ylabel('y'); axis('equal'); title('Displacement OkadaCode')
hold off

%Reshaping vectors before calculating residual
Ux=reshape(Ux,(size(uE)));
Uy=reshape(Uy,(size(uE)));
Uz=reshape(Uz,(size(uE)));

%Calculating residuals
Resx=uE-Ux;
Resy=uN-Uy;
Resz=uZ-Uz;

figure;subplot(1,3,1),contourf(X,Y,Resx);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual ux');colorbar;
subplot(1,3,2),contourf(X,Y,Resy);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual uy');colorbar;
subplot(1,3,3),contourf(X,Y,Resz);xlabel('x'); ylabel('y'); axis('equal'); title('Disp residual uz');colorbar;


%Plotting Tri Dis displacements vs the Okada displacements
DrawDeformedGrid2d( X,Y,Ux,Uy,cmap2,(E1+E2),'Scale',5);
hax1=gca;no1=gcf;title('Triangular dislocations');
DrawDeformedGrid2d( X,Y,uE,uN,cmap2,(E1+E2),'Scale',5);
hax2=gca;no2=gcf;title('Rectangular dislocation');
hf2=figure;
s1=subplot(1,2,1);
s2=subplot(1,2,2);
pos1=get(s1,'Position');pos2=get(s2,'Position');
delete(s1);delete(s2);
hax_1=copyobj(hax1,hf2);
hax_2=copyobj(hax2,hf2);
set(hax_1, 'Position', pos1);colorbar;colormap(cmap2);
xlabel('x'); ylabel('y');gca;% title('Triangles');
set(hax_2, 'Position', pos2);colorbar;colormap(cmap2);
xlabel('x'); ylabel('y');gca;%title('gggg');
close(no1);close(no2);
titlesz=25;
fntsz=21;
ChangeFontSizes(fntsz,titlesz);
text(-25,13,'Ground displacements and dilatation, comparison of different analyical sources','Fontsize',25,'FontWeight','bold')

bad=isnan(Resx);Resx=Resx(~bad);%removing nans to get proper mean
fprintf('MaxResidualUx %i.\n',max(abs(Resx(:))))

bad=isnan(Resy);Resy=Resy(~bad);%removing nans to get proper mean
fprintf('MaxResidualUy %i.\n',max(abs(Resy(:))))

bad=isnan(Resz);Resz=Resz(~bad);%removing nans to get proper mean
fprintf('MaxResidualUz %i.\n',max(abs(Resz(:))))

P=max(abs(Resx(:)))+max(abs(Resy(:)))+max(abs(Resz(:)));
if any(P>1e-8) %Around the precision and errors introduced. Small values.
error('Calculated displacements are a long way from the analytical solutions')
else
disp('Everything looks good, tolerance checks max stress residuals and flags errors above 1e-8')

Mu=ones(size(Sxx))*0.6; %Coeff Friction
Cohesion=zeros(size(Sxx)); %Coeff Friction
[CSS] = CalculateCoulombStressOnPlane(MidPointObs,FaceNormalVectorObs,...
Sxx,Syy,Szz,Sxy,Sxz,Syz,Mu,Cohesion,PointsObs,TrianglesObs,cmap2 );

end

%Calculating residuals
Resex=exx--uEE(:);
Resey=exy-(-(uEN(:)+uNE(:))/2);

0 comments on commit 4070265

Please sign in to comment.