-
Notifications
You must be signed in to change notification settings - Fork 5
/
load_gmsh.m
198 lines (191 loc) · 6.94 KB
/
load_gmsh.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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
function mesh = load_gmsh(filename)
% Reads a mesh in msh format, version 1 or 2
% Copyright (C) 10/2007 R Lorphèvre (r.lorphevre@ulg.ac.be)
% Based on load_gmsh supplied with gmsh-2.0 and load_gmsh2 from JP
% Moitinho de Almeida (moitinho@civil.ist.utl.pt)
% number of nodes in function of the element type
msh.NODES_PER_TYPE_OF_ELEMENT = [
2 3 4 4 8 6 5 3 6 9 10 27 18 14 1 8 20 15 13 ];
% The format 2 don't sort the elements by reg phys but by
% reg-elm. If this classification is important for your program,
% use this (after calling this function):
%
% [OldRowNumber, NewRowNumber] = sort(OldMatrix(:,SortColumn));
% NewMatrix = OldMatrix(NewRowNumber,:);
%
% Change the name of OldMatrix and NewMatrix with the name of yours
% SortColumn by the number of the last column
mesh = [];
mesh.MIN = zeros(3, 1);
mesh.MAX = zeros(3, 1);
fid = fopen(filename, 'r');
disp (' ')
while 1
endoffile = 0;
while 1
tline = fgetl(fid);
if feof(fid), endoffile=1, break, end
if (tline(1) == '$' )
if tline(2) == 'N' && tline(3) == 'O'
fileformat = 1 ;
break
end
if tline(2) == 'M' && tline(3) == 'e'
fileformat = 2;
tline = fgetl(fid);
disp('Mesh Type')
disp (tline)
tline = fgetl(fid);
if (tline(1) == '$' && tline(2) == 'E'&& tline(3) == 'n')
tline = fgetl(fid);
break
else
disp (' This program can only read ASCII mesh file')
disp (' of format 1 or 2 from GMSH, try again?')
endoffile=1
break
end
end
if tline(2) == 'E' && (tline(3) == 'L' || tline(3) == 'l' )
break
end
end
end
if endoffile == 1, break, end
if tline(2) == 'N' && ( tline(3) == 'O' || tline(3) == 'o' )
disp('reading nodes')
mesh.nbNod = fscanf(fid, '%d', 1);
mesh.POS = zeros(mesh.nbNod, 3);
for(I=1:mesh.nbNod)
iNod = fscanf(fid, '%d', 1);
X = fscanf(fid, '%g', 3);
IDS(iNod) = I;
if (I == 1)
mesh.MIN = X;
mesh.MAX = X;
else
if(mesh.MAX(1) < X(1)) mesh.MAX(1) = X(1); end
if(mesh.MAX(2) < X(2)) mesh.MAX(2) = X(2); end
if(mesh.MAX(3) < X(3)) mesh.MAX(3) = X(3); end
if(mesh.MIN(1) > X(1)) mesh.MIN(1) = X(1); end
if(mesh.MIN(2) > X(2)) mesh.MIN(2) = X(2); end
if(mesh.MIN(3) > X(3)) mesh.MIN(3) = X(3); end
end
mesh.POS(I, 1) = X(1);
mesh.POS(I, 2) = X(2);
mesh.POS(I, 3) = X(3);
end
tline = fgetl(fid);
disp('nodes have been read')
elseif tline(2) == 'E' && ( tline(3) == 'L' || tline(3) == 'l')
disp('reading elements')
mesh.nbElm = fscanf(fid, '%d', 1);
if (fileformat == 1)
nbinfo = 5;
tags = 3;
end
if (fileformat == 2)
nbinfo = 4;
tags = 4;
end
mesh.ELE_INFOS = zeros(mesh.nbElm, nbinfo);
mesh.nbPoints = 0;
mesh.nbLines = 0;
mesh.nbTriangles = 0;
mesh.nbQuads = 0;
mesh.nbTets = 0;
mesh.nbHexas = 0;
mesh.nbPrisms = 0;
mesh.nbPyramids = 0;
% comment next 8 lines to get "tight" arrays (will slow down reading)
mesh.POINTS = zeros(mesh.nbElm, 2);
mesh.LINES = zeros(mesh.nbElm, 3);
mesh.TRIANGLES = zeros(mesh.nbElm, 4);
mesh.QUADS = zeros(mesh.nbElm, 5);
mesh.TETS = zeros(mesh.nbElm, 5);
mesh.HEXAS = zeros(mesh.nbElm, 9);
mesh.PRISMS = zeros(mesh.nbElm, 7);
mesh.PYRAMIDS = zeros(mesh.nbElm, 6);
for(I = 1:mesh.nbElm)
mesh.ELE_INFOS(I, :) = fscanf(fid, '%d', nbinfo);
if (fileformat == 1)
% take the mesh.ELE_INFOS(I, 5) nodes of the element
NODES_ELEM = fscanf(fid, '%d', mesh.ELE_INFOS(I, 5));
end
if (fileformat == 2)
mesh.ELE_TAGS(I,:) = fscanf(fid, '%d', (mesh.ELE_INFOS(I,3)-1));
% take the msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) nodes of the element
NODES_ELEM = fscanf(fid, '%d', msh.NODES_PER_TYPE_OF_ELEMENT (mesh.ELE_INFOS(I, 2)) );
end
if(mesh.ELE_INFOS(I, 2) == 15) %% point
mesh.nbPoints = mesh.nbPoints + 1;
mesh.POINTS(mesh.nbPoints, 1) = IDS(NODES_ELEM(1));
mesh.POINTS(mesh.nbPoints, 2) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 1) %% line
mesh.nbLines = mesh.nbLines + 1;
mesh.LINES(mesh.nbLines, 1) = IDS(NODES_ELEM(1));
mesh.LINES(mesh.nbLines, 2) = IDS(NODES_ELEM(2));
mesh.LINES(mesh.nbLines, 3) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 2) %% triangle
mesh.nbTriangles = mesh.nbTriangles + 1;
mesh.TRIANGLES(mesh.nbTriangles, 1) = IDS(NODES_ELEM(1));
mesh.TRIANGLES(mesh.nbTriangles, 2) = IDS(NODES_ELEM(2));
mesh.TRIANGLES(mesh.nbTriangles, 3) = IDS(NODES_ELEM(3));
mesh.TRIANGLES(mesh.nbTriangles, 4) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 3) %% quadrangle
mesh.nbQuads = mesh.nbQuads + 1;
mesh.QUADS(mesh.nbQuads, 1) = IDS(NODES_ELEM(1));
mesh.QUADS(mesh.nbQuads, 2) = IDS(NODES_ELEM(2));
mesh.QUADS(mesh.nbQuads, 3) = IDS(NODES_ELEM(3));
mesh.QUADS(mesh.nbQuads, 4) = IDS(NODES_ELEM(4));
mesh.QUADS(mesh.nbQuads, 5) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 4) %% tetrahedron
mesh.nbTets = mesh.nbTets + 1;
mesh.TETS(mesh.nbTets, 1) = IDS(NODES_ELEM(1));
mesh.TETS(mesh.nbTets, 2) = IDS(NODES_ELEM(2));
mesh.TETS(mesh.nbTets, 3) = IDS(NODES_ELEM(3));
mesh.TETS(mesh.nbTets, 4) = IDS(NODES_ELEM(4));
mesh.TETS(mesh.nbTets, 5) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 5) %% hexahedron
mesh.nbHexas = mesh.nbHexas + 1;
mesh.HEXAS(mesh.nbHexas, 1) = IDS(NODES_ELEM(1));
mesh.HEXAS(mesh.nbHexas, 2) = IDS(NODES_ELEM(2));
mesh.HEXAS(mesh.nbHexas, 3) = IDS(NODES_ELEM(3));
mesh.HEXAS(mesh.nbHexas, 4) = IDS(NODES_ELEM(4));
mesh.HEXAS(mesh.nbHexas, 5) = IDS(NODES_ELEM(5));
mesh.HEXAS(mesh.nbHexas, 6) = IDS(NODES_ELEM(6));
mesh.HEXAS(mesh.nbHexas, 7) = IDS(NODES_ELEM(7));
mesh.HEXAS(mesh.nbHexas, 8) = IDS(NODES_ELEM(8));
mesh.HEXAS(mesh.nbHexas, 9) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 6) %% prism
mesh.nbPrisms = mesh.nbPrisms + 1;
mesh.PRISMS(mesh.nbPrisms, 1) = IDS(NODES_ELEM(1));
mesh.PRISMS(mesh.nbPrisms, 2) = IDS(NODES_ELEM(2));
mesh.PRISMS(mesh.nbPrisms, 3) = IDS(NODES_ELEM(3));
mesh.PRISMS(mesh.nbPrisms, 4) = IDS(NODES_ELEM(4));
mesh.PRISMS(mesh.nbPrisms, 5) = IDS(NODES_ELEM(5));
mesh.PRISMS(mesh.nbPrisms, 6) = IDS(NODES_ELEM(6));
mesh.PRISMS(mesh.nbPrisms, 7) = mesh.ELE_INFOS(I, tags);
end
if(mesh.ELE_INFOS(I, 2) == 7) %% pyramid
mesh.nbPyramids = mesh.nbPyramids + 1;
mesh.PYRAMIDS(mesh.nbPyramids, 1) = IDS(NODES_ELEM(1));
mesh.PYRAMIDS(mesh.nbPyramids, 2) = IDS(NODES_ELEM(2));
mesh.PYRAMIDS(mesh.nbPyramids, 3) = IDS(NODES_ELEM(3));
mesh.PYRAMIDS(mesh.nbPyramids, 4) = IDS(NODES_ELEM(4));
mesh.PYRAMIDS(mesh.nbPyramids, 5) = IDS(NODES_ELEM(5));
mesh.PYRAMIDS(mesh.nbPyramids, 6) = IDS(NODES_ELEM(6));
mesh.PYRAMIDS(mesh.nbPyramids, 7) = mesh.ELE_INFOS(I, tags);
end
end
tline = fgetl(fid);
disp('elements have been read')
end
end
fclose(fid);