/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_readGmsh.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Mesh_readGmsh.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 14260 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "Mesh.h"
18
19 using escript::IOError;
20
21 namespace dudley {
22
23 #define MAX_numNodes_gmsh 20
24
25 /// reads a mesh from a gmsh file of name filename
26 Mesh* Mesh::readGmsh(escript::JMPI mpiInfo, const std::string& filename,
27 int numDim, bool optimize)
28 {
29 double version = 1.0;
30 int format = 0, size = sizeof(double), scan_ret;
31 dim_t numNodes, totalNumElements = 0, numTags = 0;
32 int numNodesPerElement = 0, numNodesPerElement2, element_dim = 0;
33 index_t gmsh_type, partition_id, itmp, elementary_id;
34 index_t numElements = 0, numFaceElements = 0, *id = NULL, *tag = NULL, *vertices = NULL;
35 char line[1024];
36 double rtmp0, rtmp1;
37
38 if (mpiInfo->size > 1)
39 throw DudleyException("reading GMSH with MPI is not supported yet.");
40
41 // allocate mesh
42 Mesh* mesh = new Mesh(filename, numDim, mpiInfo);
43
44 // open file
45 FILE* fileHandle = fopen(filename.c_str(), "r");
46 if (fileHandle == NULL) {
47 std::stringstream ss;
48 ss << "Opening gmsh file " << filename << " for reading failed.";
49 throw IOError(ss.str());
50 }
51
52 // start reading
53 while (1) {
54 // find line staring with $
55 do {
56 if (!fgets(line, sizeof(line), fileHandle))
57 break;
58 if (feof(fileHandle))
59 break;
60 } while (line[0] != '$');
61
62 if (feof(fileHandle))
63 break;
64
65 // format
66 if (!strncmp(&line[1], "MeshFormat", 10)) {
67 scan_ret = fscanf(fileHandle, "%lf %d %d\n", &version, &format, &size);
68 if (scan_ret == EOF)
69 throw IOError("readGmsh: early EOF while reading file");
70 }
71 // nodes
72 if (!strncmp(&line[1], "NOD", 3) || !strncmp(&line[1], "NOE", 3) || !strncmp(&line[1], "Nodes", 5)) {
73 scan_ret = fscanf(fileHandle, "%d", &numNodes);
74 if (scan_ret == EOF)
75 throw IOError("readGmsh: early EOF while reading file");
76 mesh->Nodes->allocTable(numNodes);
77 for (index_t i0 = 0; i0 < numNodes; i0++) {
78 if (1 == numDim) {
79 scan_ret = fscanf(fileHandle, "%d %le %le %le\n",
80 &mesh->Nodes->Id[i0],
81 &mesh->Nodes->Coordinates[INDEX2(0, i0, numDim)],
82 &rtmp0, &rtmp1);
83 if (scan_ret == EOF)
84 throw IOError("readGmsh: early EOF while reading file");
85 } else if (2 == numDim) {
86 scan_ret = fscanf(fileHandle, "%d %le %le %le\n",
87 &mesh->Nodes->Id[i0],
88 &mesh->Nodes->Coordinates[INDEX2(0, i0, numDim)],
89 &mesh->Nodes->Coordinates[INDEX2(1, i0, numDim)],
90 &rtmp0);
91 if (scan_ret == EOF)
92 throw IOError("readGmsh: early EOF while reading file");
93 } else if (3 == numDim) {
94 scan_ret = fscanf(fileHandle, "%d %le %le %le\n",
95 &mesh->Nodes->Id[i0],
96 &mesh->Nodes->Coordinates[INDEX2(0, i0, numDim)],
97 &mesh->Nodes->Coordinates[INDEX2(1, i0, numDim)],
98 &mesh->Nodes->Coordinates[INDEX2(2, i0, numDim)]);
99 if (scan_ret == EOF)
100 throw IOError("readGmsh: early EOF while reading file");
101 }
102 mesh->Nodes->globalDegreesOfFreedom[i0] = mesh->Nodes->Id[i0];
103 mesh->Nodes->Tag[i0] = 0;
104 }
105 }
106 // elements
107 else if (!strncmp(&line[1], "ELM", 3) || !strncmp(&line[1], "Elements", 8)) {
108 ElementTypeId final_element_type = Dudley_NoRef;
109 ElementTypeId final_face_element_type = Dudley_NoRef;
110 numElements = 0;
111 numFaceElements = 0;
112 scan_ret = fscanf(fileHandle, "%d", &totalNumElements);
113 if (scan_ret == EOF)
114 throw IOError("readGmsh: early EOF while reading file");
115
116 id = new index_t[totalNumElements];
117 tag = new index_t[totalNumElements];
118
119 ElementTypeId* element_type = new ElementTypeId[totalNumElements];
120 vertices = new index_t[totalNumElements * MAX_numNodes_gmsh];
121 // read all in
122 for (index_t e = 0; e < totalNumElements; e++) {
123 scan_ret = fscanf(fileHandle, "%d %d", &id[e], &gmsh_type);
124 if (scan_ret == EOF)
125 throw IOError("readGmsh: early EOF while reading file");
126 switch (gmsh_type) {
127 case 1: // line order 1
128 element_type[e] = Dudley_Line2;
129 element_dim = 1;
130 numNodesPerElement = 2;
131 break;
132 case 2: // triangle order 1
133 element_type[e] = Dudley_Tri3;
134 numNodesPerElement = 3;
135 element_dim = 2;
136 break;
137 case 4: // tetrahedron order 1
138 element_type[e] = Dudley_Tet4;
139 numNodesPerElement = 4;
140 element_dim = 3;
141 break;
142 case 15: // point
143 element_type[e] = Dudley_Point1;
144 numNodesPerElement = 1;
145 element_dim = 0;
146 break;
147 default: {
148 std::stringstream ss;
149 ss << "Unexpected gmsh element type " << gmsh_type
150 << " in mesh file " << filename;
151 throw IOError(ss.str());
152 }
153 }
154 if (element_dim == numDim) {
155 if (final_element_type == Dudley_NoRef) {
156 final_element_type = element_type[e];
157 } else if (final_element_type != element_type[e]) {
158 throw IOError("Dudley can handle a single type of "
159 "internal elements only.");
160 }
161 numElements++;
162 } else if (element_dim == numDim - 1) {
163 if (final_face_element_type == Dudley_NoRef) {
164 final_face_element_type = element_type[e];
165 } else if (final_face_element_type != element_type[e]) {
166 throw IOError("Dudley can handle a single type of "
167 "face elements only.");
168 }
169 numFaceElements++;
170 }
171
172 if (version <= 1.0) {
173 scan_ret = fscanf(fileHandle, "%d %d %d", &tag[e],
174 &elementary_id, &numNodesPerElement2);
175 if (scan_ret == EOF)
176 throw IOError("readGmsh: early EOF while reading file");
177 partition_id = 1;
178 if (numNodesPerElement2 != numNodesPerElement) {
179 std::stringstream ss;
180 ss << "Illegal number of nodes for element " << id[e]
181 << " in mesh file " << filename;
182 throw IOError(ss.str());
183 }
184 } else {
185 scan_ret = fscanf(fileHandle, "%d", &numTags);
186 if (scan_ret == EOF)
187 throw IOError("readGmsh: early EOF while reading file");
188 elementary_id = tag[e] = partition_id = 1;
189 numNodesPerElement2 = -1;
190 for (int j = 0; j < numTags; j++) {
191 scan_ret = fscanf(fileHandle, "%d", &itmp);
192 if (scan_ret == EOF)
193 throw IOError("readGmsh: early EOF while reading file");
194 if (j == 0) {
195 tag[e] = itmp;
196 } else if (j == 1) {
197 elementary_id = itmp;
198 } else if (j == 2) {
199 partition_id = itmp;
200 }
201 // ignore any other tags
202 }
203 }
204 for (int j = 0; j < numNodesPerElement; j++) {
205 scan_ret = fscanf(fileHandle, "%d",
206 &vertices[INDEX2(j, e, MAX_numNodes_gmsh)]);
207 if (scan_ret == EOF)
208 throw IOError("readGmsh: early EOF while reading file");
209 }
210 }
211 // all elements have been read, now we have to identify the
212 // dudley elements to define Elements and FaceElements
213 if (final_element_type == Dudley_NoRef) {
214 if (numDim == 1) {
215 final_element_type = Dudley_Line2;
216 } else if (numDim == 2) {
217 final_element_type = Dudley_Tri3;
218 } else if (numDim == 3) {
219 final_element_type = Dudley_Tet4;
220 }
221 }
222 if (final_face_element_type == Dudley_NoRef) {
223 if (numDim == 1) {
224 final_face_element_type = Dudley_Point1;
225 } else if (numDim == 2) {
226 final_face_element_type = Dudley_Line2;
227 } else if (numDim == 3) {
228 final_face_element_type = Dudley_Tri3;
229 }
230 }
231 mesh->Elements = new ElementFile(final_element_type, mpiInfo);
232 mesh->FaceElements = new ElementFile(final_face_element_type, mpiInfo);
233 mesh->Points = new ElementFile(Dudley_Point1, mpiInfo);
234 mesh->Elements->allocTable(numElements);
235 mesh->FaceElements->allocTable(numFaceElements);
236 mesh->Points->allocTable(0);
237 mesh->Elements->minColor = 0;
238 mesh->Elements->maxColor = numElements - 1;
239 mesh->FaceElements->minColor = 0;
240 mesh->FaceElements->maxColor = numFaceElements - 1;
241 mesh->Points->minColor = 0;
242 mesh->Points->maxColor = 0;
243 numElements = 0;
244 numFaceElements = 0;
245 for (index_t e = 0; e < totalNumElements; e++) {
246 if (element_type[e] == final_element_type) {
247 mesh->Elements->Id[numElements] = id[e];
248 mesh->Elements->Tag[numElements] = tag[e];
249 mesh->Elements->Color[numElements] = numElements;
250 mesh->Elements->Owner[numElements] = 0;
251 for (int j = 0; j < mesh->Elements->numNodes; ++j) {
252 mesh->Elements->Nodes[INDEX2(j, numElements,
253 mesh->Elements->numNodes)] =
254 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
255 }
256 numElements++;
257 } else if (element_type[e] == final_face_element_type) {
258 mesh->FaceElements->Id[numFaceElements] = id[e];
259 mesh->FaceElements->Tag[numFaceElements] = tag[e];
260 mesh->FaceElements->Color[numFaceElements] = numFaceElements;
261 mesh->FaceElements->Owner[numFaceElements] = 0;
262 for (int j = 0; j < mesh->FaceElements->numNodes; ++j)
263 {
264 mesh->FaceElements->Nodes[INDEX2(j, numFaceElements,
265 mesh->FaceElements->numNodes)] =
266 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
267 }
268 numFaceElements++;
269 }
270 }
271 // and clean up
272 delete[] id;
273 delete[] tag;
274 delete[] element_type;
275 delete[] vertices;
276 } // end reading elements
277 // name tags (thanks to Antoine Lefebvre, antoine.lefebvre2@mail.mcgill.ca )
278 else if (!strncmp(&line[1], "PhysicalNames", 13)) {
279 char name[1024];
280 index_t tag_key;
281 scan_ret = fscanf(fileHandle, "%d", &numTags);
282 if (scan_ret == EOF)
283 throw IOError("readGmsh: early EOF while reading file");
284 for (int i0 = 0; i0 < numTags; i0++) {
285 scan_ret = fscanf(fileHandle, "%d %d %s\n", &itmp, &tag_key, name);
286 if (scan_ret == EOF)
287 throw IOError("readGmsh: early EOF while reading file");
288 if (itmp != 2)
289 throw IOError("readGmsh: expecting two entries per physical name.");
290 if (strlen(name) < 3)
291 throw IOError("readGmsh: illegal tagname (\" missing?)");
292 name[strlen(name)-1]='\0';
293 mesh->addTagMap(&name[1], tag_key);
294 }
295 }
296 // search for end of data block
297 do {
298 if (!fgets(line, sizeof(line), fileHandle)) {
299 std::stringstream ss;
300 ss << "Unexpected end of file in " << filename;
301 throw IOError(ss.str());
302 }
303 if (feof(fileHandle)) {
304 std::stringstream ss;
305 ss << "Unexpected end of file in " << filename;
306 throw IOError(ss.str());
307 }
308 } while (line[0] != '$');
309 }
310
311 fclose(fileHandle);
312 mesh->resolveNodeIds();
313 mesh->prepare(optimize);
314 return mesh;
315 }
316
317 } // namespace dudley
318

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_readGmsh.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_readGmsh.cpp:2682-2741 /branches/pasowrap/dudley/src/Mesh_readGmsh.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_readGmsh.cpp:3871-3891 /branches/restext/finley/src/Mesh_readGmsh.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Mesh_readGmsh.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_readGmsh.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Mesh_readGmsh.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_readGmsh.cpp:3517-3974 /release/3.0/finley/src/Mesh_readGmsh.cpp:2591-2601 /release/4.0/dudley/src/Mesh_readGmsh.cpp:5380-5406 /trunk/dudley/src/Mesh_readGmsh.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Mesh_readGmsh.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26