-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathMeshGenerate.m
98 lines (79 loc) · 2.25 KB
/
MeshGenerate.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
function [p,t]=MeshGenerate(outline_src,hmax)
%use the program Mesh2d_v24 by Darren Engwirda to generate a triangle mesh
%and write the mesh into a text file
if outline_src==1
[nd,cnect]=outline_sms; %read from sms map file
elseif outline_src==2
[nd,cnect]=outline_cs; %read from cross-section file
elseif outline_src==3
[nd,cnect]=outline_obj; %read from a struct variable
else
disp('invalid input parameter');
return;
end
if nargin==1
hdata.hmax=100; %limit the cell size
else
hdata.hmax=hmax;
end
[p,t] = mesh2d(nd,cnect,hdata);
%------------------------------------------------------------
function [nd,cnect]=outline_sms
[nodes,arcs,polygons]=read_sms_map;
nd=nodes(:,2:3);
nplg=size(polygons,2);
cnect=[];
for i=1:1:nplg
PlgArcs=polygons{i};
narc_in_plg=size(PlgArcs,1);
c=zeros(narc_in_plg,2);
for j=1:1:narc_in_plg
ind=find(arcs(:,1)==PlgArcs(j)); %find the row index of the jth arc in a polygon
rnd=find(nodes(:,1)==arcs(ind,2));
c(j,1)=rnd;
rnd=find(nodes(:,1)==arcs(ind,3));
c(j,2)=rnd;
end
%check the connective relationship
for j=2:1:narc_in_plg
if (sum(c(j-1,:)==c(j,1))==0)&&(sum(c(j-1,:)==c(j,2))==0)
msgbox('invalid connective relationship');
else
if c(j-1,1)==c(j,1)
c(j-1,:)=fliplr(c(j-1,:));
elseif c(j-1,1)==c(j,2)
c(j-1:j,:)=fliplr(c(j-1:j,:));
elseif c(j-1,2)==c(j,2)
c(j,:)=fliplr(c(j,:));
end
end
end
cnect=[cnect,c];
end
%------------------------------------------------------------
function [nd,cnect]=outline_cs
cs_xyz=read_elevation;
nd=[];
ncs=size(cs_xyz,2);
for i=1:1:ncs
nd=[nd; cs_xyz(i).xyz(1,1:2)];
end
for i=ncs:-1:1
nd=[nd; cs_xyz(i).xyz(cs_xyz(i).npt, 1:2)];
end
nnd=size(nd,1);
cnect=[(1:nnd).',[(2:nnd).';1]];
%------------------------------------------------------------
function [nd,cnect]=outline_obj
[filename,path]=uigetfile('*.*');
load([path,filename],'CSold');
nd=[];
ncs=size(CSold,2);
for i=1:1:ncs
nd=[nd; CSold(i).xy(1,1:2)];
end
for i=ncs:-1:1
nd=[nd; CSold(i).xy(CSold(i).nodes, 1:2)];
end
nnd=size(nd,1);
cnect=[(1:nnd).',[(2:nnd).';1]];