Skip to content

Commit 6a04cce

Browse files
daisy20170101daisy20170101
committed
add scripts for mesh generating
1 parent 68410b1 commit 6a04cce

5 files changed

+439
-0
lines changed

‎Mesh/PlotVariable.m

+166
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
% read from accuacc.data to plot snapshot of nucleation
2+
% var = load('./lk22_h80_N25_T1/');
3+
% var = load('./lk22_h80_N25_T1/vardep-lk22_h80_N25_T1.dat');
4+
% var = load('./lk20_h90_N25_T2/vardep-lk20_h90_N25_T2.dat');
5+
% var = load('./lk 20_h90_N25_T2_slip/vardep-lk20_h90_N25_T2.dat');
6+
var = load('h90_N25_T2/vardep-h90_N25_T2.dat');
7+
% var = load('eq_h90_N25_T2/NG/vardep-lk20_h90_N25_T2.dat');
8+
9+
depth = var(:,1);
10+
eff = var(:,2);
11+
dc = var(:,3);
12+
pab = var(:,4);
13+
pa = var(:,5);
14+
15+
len=length(pa);
16+
17+
%% calculate accumulative slip on the fault
18+
% num = load('../Patch_num4.txt');
19+
data = load('400km_1km_smooth.gts'); % load meshing - read geometry
20+
nvex = data(1,1); ncell = data(1,3);
21+
vertex = data(2:nvex+1,:);
22+
cell = data(2+nvex:1+nvex+ncell,:);
23+
24+
%%
25+
load('/Users/duoli/Documents/Mexico/2014SSE/Contour/guerrero.mat');
26+
[x,y] = km2lonlat_mex(vertex(:,1)/1000,vertex(:,2)/1000,-55); % rotation back to longitude-latitude
27+
28+
figure;
29+
set(gcf,'position',[100 100 850 750]);
30+
subplot(2,2,1);
31+
hold on;
32+
box on;
33+
set(gca,'linewidth',1.2,'fontsize',12);
34+
% scatter(x(1:1:end,1),y(1:1:end,1),8,100*accslip(1:1:end),'filled');
35+
set(gca,'xlim',[-102.5, -99],'ylim',[16.5,19],'fontsize',12);
36+
% title('Mapview of a-b');
37+
xlabel('longitude');
38+
ylabel('latitude');
39+
colormap jet;
40+
41+
cl = colorbar ;
42+
% cl.Limits = [0 50] ;
43+
caxis([-0.0035 0.0035 ]) ;
44+
cl.Label.String = 'a-b' ;
45+
46+
tr = triangulation(cell(1:len,:),x,y,-vertex(:,3));
47+
trisurf(tr,pab,'edgecolor','none','facealpha',0.9);
48+
49+
% earthquake
50+
eq = load('/Users/duoli/Documents/Mexico//2014SSE/Contour/2014Eq_USGS.txt');
51+
plot(eq(1,1),eq(1,2),'pr','markerfacecolor','r','markersize',16); % 2014 Papanoa epicenter from USGS
52+
af = load('/Users/duoli/Documents/Mexico//2014SSE/Contour/2014Eq_aftershock.txt');
53+
% plot(af(:,1),af(:,2),'pk','markerfacecolor','y','markersize',16); % 2014 Papanoa epicenter from USGS
54+
% city
55+
city = load('/Users/duoli/Documents/Mexico//2014SSE/Contour/city2.txt');
56+
plot(city(:,1),city(:,2),'sk','markerfacecolor','k','markersize',6); % 2014 Papanoa epicenter from USGS
57+
58+
% set(gca,'Yscale','log','YLim',[1.91e20 1.98e20]);
59+
60+
% trench, coastline and depth contour
61+
trch = load('/Users/duoli/Documents/Mexico//2014SSE/Contour/MAT_trench.txt');
62+
[tr_xr,tr_yr] = lonlat2km_mex(trch(:,1),trch(:,2),0);
63+
plot(trch(:,1),trch(:,2),'-k','linewidth',2.0);
64+
% plot(trch(1:12:end,1),trch(1:12:end,2)+0.015,'<k','markerface','k');
65+
[ncst_xr,ncst_yr] = lonlat2km_mex(ncst(:,1),ncst(:,2),0); %% sphere plattern + rotate 70 degree.
66+
mapshow(ncst(:,1),ncst(:,2),'DisplayType','line','color','k','linewidth',1.2);
67+
68+
depcon = load('/Users/duoli/Documents/Mexico//2014SSE/Contour/contour_latlon2.txt');
69+
mapshow(depcon(:,1),depcon(:,2),'displaytype','line','color','w','linestyle','-.','linewidth',1.5);
70+
71+
72+
% legend('data','2014 Mw7.2','Aftershocks','GPS stations','trench','coastline','isodepth');
73+
74+
%
75+
subplot(2,2,2);
76+
hold on;
77+
box on;
78+
set(gca,'linewidth',1.2,'fontsize',12);
79+
% scatter(x(1:1:end,1),y(1:1:end,1),8,100*accslip(1:1:end),'filled');
80+
set(gca,'xlim',[-102.5, -99],'ylim',[16.5,19],'fontsize',12);
81+
% title('Mapview of effective normal stress');
82+
xlabel('longitude');
83+
ylabel('latitude');
84+
colormap jet;
85+
86+
cl = colorbar ;
87+
% cl.Limits = [0 50] ;
88+
caxis([0 50 ]) ;
89+
cl.Label.String = 'effective normal stress (MPa)' ;
90+
91+
trisurf(tr,eff/10,'edgecolor','none','facealpha',0.9);
92+
93+
% earthquake
94+
plot(eq(1,1),eq(1,2),'pr','markerfacecolor','r','markersize',16); % 2014 Papanoa epicenter from USGS
95+
% plot(af(:,1),af(:,2),'pk','markerfacecolor','y','markersize',16); % 2014 Papanoa epicenter from USGS
96+
% city
97+
plot(city(:,1),city(:,2),'sk','markerfacecolor','k','markersize',6); % 2014 Papanoa epicenter from USGS
98+
99+
% set(gca,'Yscale','log','YLim',[1.91e20 1.98e20]);
100+
101+
% trench, coastline and depth contour
102+
plot(trch(:,1),trch(:,2),'-k','linewidth',2.0);
103+
% plot(trch(1:12:end,1),trch(1:12:end,2)+0.015,'<k','markerface','k');
104+
mapshow(ncst(:,1),ncst(:,2),'DisplayType','line','color','k','linewidth',1.2);
105+
mapshow(depcon(:,1),depcon(:,2),'displaytype','line','color','w','linestyle','-.','linewidth',1.5);
106+
107+
% legend('data','2014 Mw7.2','Aftershocks','GPS stations','trench','coastline','isodepth');
108+
109+
subplot(2,2,3);
110+
hold on;
111+
box on;
112+
set(gca,'linewidth',1.2,'fontsize',12);
113+
% scatter(x(1:1:end,1),y(1:1:end,1),8,100*accslip(1:1:end),'filled');
114+
set(gca,'xlim',[-102.5, -99],'ylim',[16.5,19],'fontsize',12);
115+
xlabel('longitude');
116+
ylabel('latitude');
117+
colormap jet;
118+
119+
cl = colorbar ;
120+
% cl.Limits = [0 50] ;
121+
caxis([0 200 ]) ;
122+
cl.Label.String = 'dc (mm)' ;
123+
124+
trisurf(tr,dc,'edgecolor','none','facealpha',0.9);
125+
126+
% plot(af(:,1),af(:,2),'pk','markerfacecolor','y','markersize',16); % 2014 Papanoa epicenter from USGS
127+
% city
128+
plot(city(:,1),city(:,2),'sk','markerfacecolor','k','markersize',6); % 2014 Papanoa epicenter from USGS
129+
130+
% set(gca,'Yscale','log','YLim',[1.91e20 1.98e20]);
131+
132+
% trench, coastline and depth contour
133+
plot(trch(:,1),trch(:,2),'-k','linewidth',2.0);
134+
% plot(trch(1:12:end,1),trch(1:12:end,2)+0.015,'<k','markerface','k');
135+
mapshow(ncst(:,1),ncst(:,2),'DisplayType','line','color','k','linewidth',1.2);
136+
mapshow(depcon(:,1),depcon(:,2),'displaytype','line','color','w','linestyle','-.','linewidth',1.5);
137+
138+
% legend('data','2014 Mw7.2','Aftershocks','GPS stations','trench','coastline','isodepth');
139+
140+
subplot(2,2,4);
141+
hold on;
142+
box on;
143+
set(gca,'linewidth',1.2,'fontsize',12);
144+
% scatter(x(1:1:end,1),y(1:1:end,1),8,100*accslip(1:1:end),'filled');
145+
set(gca,'xlim',[-102.5, -99],'ylim',[16.5,19],'fontsize',12);
146+
xlabel('longitude');
147+
ylabel('latitude');
148+
colormap jet;
149+
150+
cl = colorbar ;
151+
% cl.Limits = [0 50] ;
152+
caxis([0.01 0.05]) ;
153+
cl.Label.String = 'a' ;
154+
155+
trisurf(tr,pa,'edgecolor','none','facealpha',0.9);
156+
% plot(af(:,1),af(:,2),'pk','markerfacecolor','y','markersize',16); % 2014 Papanoa epicenter from USGS
157+
plot(city(:,1),city(:,2),'sk','markerfacecolor','k','markersize',6); % 2014 Papanoa epicenter from USGS
158+
159+
% set(gca,'Yscale','log','YLim',[1.91e20 1.98e20]);
160+
% trench, coastline and depth contour
161+
plot(trch(:,1),trch(:,2),'-k','linewidth',2.0);
162+
% plot(trch(1:12:end,1),trch(1:12:end,2)+0.015,'<k','markerface','k');
163+
mapshow(ncst(:,1),ncst(:,2),'DisplayType','line','color','k','linewidth',1.2);
164+
mapshow(depcon(:,1),depcon(:,2),'displaytype','line','color','w','linestyle','-.','linewidth',1.5);
165+
166+
% legend('data','2014 Mw7.2','Aftershocks','GPS stations','trench','coastline','isodepth');

‎Mesh/ReadInpFile.m

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
%% read Inp file for mexico slab geometry. So Slow!! Use fortran ReadInp instead!!
2+
3+
% f1 = strcat('./mexicoslab_fine.inp');
4+
% n_node = 654; % Test for Coulomb stress.
5+
% n_ele = 1179;
6+
7+
% f1 = strcat('./400km_1km_smooth.inp');
8+
fnod = strcat('./400_1km_node');
9+
fele = strcat('./400_1km_ele');
10+
11+
n_node = 85884;
12+
n_ele = 170534;
13+
14+
[ind1,xx,yy,zz] = textread(fnod,'%d,%f,%f,%f\n',n_node);
15+
[ind2,n1,n2,n3] = textread(fele,'%d,%d,%d,%d\n',n_ele);
16+
17+
% figure;
18+
% hold on; box on;
19+
for i = 1:n_ele
20+
ele1 = find(ind1 == n1(i));
21+
ele2 = find(ind1 == n2(i));
22+
ele3 = find(ind1 == n3(i));
23+
xl(i,:) = [xx(ele1),xx(ele2),xx(ele3),xx(ele1)];
24+
yl(i,:) = [yy(ele1),yy(ele2),yy(ele3),yy(ele1)];
25+
zl(i,:) = [zz(ele1),zz(ele2),zz(ele3),zz(ele1)];
26+
nn(i,1:3) = [ele1,ele2,ele3];
27+
% plot3(xl(i,:),yl(i,:),zl(i,:),'-k');
28+
end
29+
30+
zz = -zz;
31+
32+
data1 = [xx,yy,zz];
33+
save 400km_1km_smooth.gts -ascii data1;
34+
save 400km_1km_smooth.gts -ascii -append nn;

‎Mesh/SelectPoint.m

+39
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
%% written by D. Li
2+
3+
fin = fopen('400km_1km_smooth.gts','r');
4+
% fin = fopen('triangular_mesh.gts','r'); % read in mesh file
5+
6+
nnum = textscan(fin,'%d %d %d\n',1);
7+
nvex = nnum{1}; nedge = nnum{2}; ncell = nnum{3};
8+
vex = textscan(fin,'%f %f %f\n',nvex);
9+
cell = textscan(fin,'%d %d %d\n',ncell);
10+
x = vex{1} ;
11+
y = vex{2};
12+
z = vex{3};
13+
x = x/1000;
14+
y = y/1000;
15+
z = z/1000;
16+
n1 = cell{1};
17+
n2 = cell{2};
18+
n3 = cell{3};
19+
20+
x1 = x(n1); y1=y(n1); z1=z(n1);
21+
x2 = x(n2); y2=y(n2); z2=z(n2);
22+
x3 = x(n3); y3=y(n3); z3=z(n3);
23+
24+
pp = [x1,y1,-z1];
25+
% find points at depths
26+
n1 = find(pp(:,3)> 9 & pp(:,3)< 10 & pp(:,2)>-132.0 & pp(:,2)<-130) ;
27+
n2 = find(pp(:,3)> 19 & pp(:,3)< 20 & pp(:,2)>-132.0 & pp(:,2)<-130) ;
28+
n3 = find(pp(:,3)> 29 & pp(:,3)< 30 & pp(:,2)>-132.0 & pp(:,2)<-130) ;
29+
n4 = find(pp(:,3)> 39 & pp(:,3)< 40 & pp(:,2)>-132.0 & pp(:,2)<-130) ;
30+
31+
n5 = find(pp(:,3)> 19 & pp(:,3)< 20 & pp(:,2)>30.0 & pp(:,2)<32) ;
32+
n6 = find(pp(:,3)> 19 & pp(:,3)< 20 & pp(:,2)>80.0 & pp(:,2)<82) ;
33+
n7 = find(pp(:,3)> 19 & pp(:,3)< 20 & pp(:,2)>-82.0 & pp(:,2)<-80) ;
34+
n8 = find(pp(:,3)> 19 & pp(:,3)< 20 & pp(:,2)>-32.0 & pp(:,2)<-30) ;
35+
36+
data1 = [n1(1),n2(1),n3(1),n4(1),n5(1),n6(1),n7(1),n8(1)];
37+
38+
save('ObvPoints.txt','-ascii','-int','data1');
39+
%%

‎Mesh/WriteGeoEqSSE.m

+95
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
%% write geometry, unit is meter.
2+
% generate triangular mesh of diameters approximating to 2000 m
3+
fout = fopen('geometry_400slab_4.jou','w+');
4+
fprintf(fout,'%s\n','${Units(si)}');
5+
fprintf(fout,'%s\n','reset');
6+
% fprintf(fout,'%s\n','undo off');
7+
8+
depth = [5 20 30 40 50 60];
9+
10+
for dd = 1:length(depth)
11+
fname = strcat('../../../2014SSE/CubitInput/csmooth2_',num2str(depth(dd)),'.txt');
12+
subgrd = load(fname);
13+
subgrd = subgrd(1:3:end,:);
14+
dnum = find(subgrd(:,2) < 170 & subgrd(:,2)> -230);
15+
subgrd = subgrd(dnum,:)*1e3;
16+
17+
% subgrd(:,3) = -subgrd(:,3)-60e3 ;
18+
dep = (depth(dd)-4.9)*1e3; % depth should in positive.
19+
20+
fprintf(fout,'%s%f%s%f%s%f%s\n','create vertex x {',subgrd(1,1),'} y {',subgrd(1,2),'} z {',dep,'}');
21+
fprintf(fout,'%s%d%s\n','${idPtTopeS',dd,'=Id("vertex")}');
22+
23+
for i = 2:length(subgrd(:,1))
24+
fprintf(fout,'%s%f%s%f%s%f%s\n','create vertex x {',subgrd(i,1),'} y {',subgrd(i,2),'} z {',dep,'}');
25+
end
26+
fprintf(fout,'%s%d%s\n','${idPtTopeN',dd,'=Id("vertex")}');
27+
% fprintf(fout,'%s\n','${idPtTopeN=Id("vertex")}');
28+
29+
fprintf(fout,'%s%d%s%d%s\n','create curve spline vertex {idPtTopeS',dd,'} to {idPtTopeN',dd,'}');
30+
fprintf(fout,'%s%d%s\n','curve {Id("curve")} name "curve',dd,'"');
31+
32+
% if dd >1
33+
% fprintf(fout,'%s%d%s%d%s\n','create curve spline vertex {idPtTopeS',dd,'} {idPtTopeS',dd-1,'}');
34+
% fprintf(fout,'%s%d%s\n','curve {Id("curve")} name "edgeS',dd,'"');
35+
% fprintf(fout,'%s%d%s%d%s\n','create curve spline vertex {idPtTopeN',dd,'} {idPtTopeN',dd-1,'}');
36+
% fprintf(fout,'%s%d%s\n','curve {Id("curve")} name "edgeN',dd,'"');
37+
% fprintf(fout,'%s%d%s%d%s%d%s%d\n','create surface curve curve',dd,' edgeN',dd,' curve',dd-1,' edgeS',dd);
38+
% fprintf(fout,'%s%d%s\n','surface {Id("surface")} name "surf',dd,'"');
39+
% end
40+
end
41+
% fprintf(fout,'%s%d %d %d %d %d %d %d %d %d %d %d\n','create surface curve ',1,2,3,4,5,6,7,8,9,10,11);
42+
% fprintf(fout,'%s\n','create curve spline vertex {idPtTopeN1} {idPtTopeN2} {idPtTopeN3} {idPtTopeN4} {idPtTopeN5} {idPtTopeN6} {idPtTopeN7} {idPtTopeN8} {idPtTopeN9} {idPtTopeN10}');
43+
% fprintf(fout,'%s\n','curve {Id("curve")} name "edgeN"');
44+
% fprintf(fout,'%s\n','create curve spline vertex {idPtTopeS1} {idPtTopeS2} {idPtTopeS3} {idPtTopeS4} {idPtTopeS5} {idPtTopeS6} {idPtTopeS7} {idPtTopeS8} {idPtTopeS9} {idPtTopeS10}');
45+
% fprintf(fout,'%s\n','curve {Id("curve")} name "edgeS"');
46+
% fprintf(fout,'%s%d%s\n','create curve spline vertex {idPtTopeN1} {idPtTopeN',dd,'}');
47+
% fprintf(fout,'%s\n','curve {Id("curve")} name "edgeN"');
48+
% fprintf(fout,'%s%d%s\n','create curve spline vertex {idPtTopeS',dd,'} {idPtTopeS1}');
49+
% fprintf(fout,'%s\n','curve {Id("curve")} name "edgeS"');
50+
% fprintf(fout,'%s%d%s\n','create surface curve curve1 curve',dd,' edgeN edgeS');
51+
% fprintf(fout,'%s\n','surface {Id("surface")} name "fault"');
52+
53+
% fprintf(fout,'%s\n','create surface vertex 100 120 130 on surface fault');
54+
% fprintf(fout,'%s\n','delete vertex all');
55+
% fprintf(fout,'%s\n','delete curve all');
56+
57+
%% create volume
58+
% fprintf(fout,'%s\n','create vertex x -300 y -800 z 0');
59+
% fprintf(fout,'%s\n','${box1=Id("vertex")}');
60+
% fprintf(fout,'%s\n','create vertex x 500 y -800 z 0');
61+
% fprintf(fout,'%s\n','${box2=Id("vertex")}');
62+
% fprintf(fout,'%s\n','create vertex x 500 y 300 z 0');
63+
% fprintf(fout,'%s\n','${box3=Id("vertex")}');
64+
% fprintf(fout,'%s\n','create vertex x -300 y 300 z 0');
65+
% fprintf(fout,'%s\n','${box4=Id("vertex")}');
66+
% fprintf(fout,'%s\n','create surface vertex {box1} {box2} {box3} {box4}');
67+
% fprintf(fout,'%s\n','surface {Id("surface")} name "freesurf"');
68+
%
69+
% fprintf(fout,'%s\n','sweep surface freesurf vector 0 0 -1 distance 200');
70+
% fprintf(fout,'%s\n','volume {Id("volume")} name "faultbody"');
71+
% fprintf(fout,'%s\n', 'Webcut volume faultbody sweep surface fault vector 0 0 1 Through_all');
72+
73+
% fprintf(fout,'%s\n','create vertex x -100 y -250 z 60');
74+
% fprintf(fout,'%s\n','create vertex x 300 y -250 z 60');
75+
% fprintf(fout,'%s\n','create vertex x -100 y 200 z 60');
76+
% fprintf(fout,'%s\n','create vertex x 300 y 200 z 60');
77+
% fprintf(fout,'%s\n','delete surface fault');
78+
79+
fprintf(fout,'%s\n','create surface skin curve all');
80+
fprintf(fout,'%s\n','delete curve all');
81+
fprintf(fout,'%s\n','delete vertex all');
82+
% fprintf(fout,'%s\n','export acis "surf_top.sat" overwrite');
83+
84+
% fprintf(fout,'%s\n','imprint all');
85+
% fprintf(fout,'%s\n','merge all');
86+
% fprintf(fout,'%s\n','rotate 45 about x');
87+
% fprintf(fout,'%s\n','rotate 30 about z');
88+
% fprintf(fout,'%s\n','save as "cascadia_mega.cub" overwrite');
89+
90+
91+
% fprintf(fout,'%s\n','${dx=0.5}');
92+
% fprintf(fout,'%s\n','surface fault scheme trimesh');
93+
% fprintf(fout,'%s\n','mesh surface fault');
94+
95+

0 commit comments

Comments
 (0)