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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4261 - (show annotations)
Wed Feb 27 06:09:33 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 13230 byte(s)
Initial all c++ build.
But ... there are now reinterpret_cast<>'s
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by 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 since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 /************************************************************************************/
17
18 /* Dudley: read mesh */
19
20 /************************************************************************************/
21
22 #include "Mesh.h"
23 #include <stdio.h>
24
25 #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Dudley_setError(IO_ERROR,"scan error while reading dudley file"); return NULL;} }
26
27 /************************************************************************************/
28
29 /* reads a mesh from a Dudley file of name fname */
30
31 #define MAX_numNodes_gmsh 20
32
33 Dudley_Mesh *Dudley_Mesh_readGmsh(char *fname, index_t numDim, index_t order, index_t reduced_order, bool_t optimize,
34 bool_t useMacroElements)
35 {
36
37 double version = 1.0;
38 int format = 0, size = sizeof(double), scan_ret;
39 dim_t numNodes, totalNumElements = 0, numTags = 0, numNodesPerElement = 0, numNodesPerElement2, element_dim = 0;
40 index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;
41 index_t numElements = 0, numFaceElements = 0, *id = NULL, *tag = NULL, *vertices = NULL;
42 Dudley_Mesh *mesh_p = NULL;
43 char line[LenString_MAX + 1];
44 char error_msg[LenErrorMsg_MAX];
45 double rtmp0, rtmp1;
46 #ifdef Dudley_TRACE
47 double time0 = Dudley_timer();
48 #endif
49 FILE *fileHandle_p = NULL;
50 Dudley_ElementTypeId *element_type = NULL;
51
52 Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
53 Dudley_resetError();
54 if (mpi_info->size > 1)
55 {
56 Dudley_setError(IO_ERROR, "reading GMSH with MPI is not supported yet.");
57 Esys_MPIInfo_free(mpi_info);
58 return NULL;
59 }
60 else
61 {
62
63 /* allocate mesh */
64
65 mesh_p = Dudley_Mesh_alloc(fname, numDim, mpi_info);
66 if (!Dudley_noError())
67 return NULL;
68
69 /* get file handle */
70 fileHandle_p = fopen(fname, "r");
71 if (fileHandle_p == NULL)
72 {
73 sprintf(error_msg, "Opening Gmsh file %s for reading failed.", fname);
74 Dudley_setError(IO_ERROR, error_msg);
75 Esys_MPIInfo_free(mpi_info);
76 return NULL;
77 }
78
79 /* start reading */
80 while (1)
81 {
82 if (!Dudley_noError())
83 break;
84 /* find line staring with $ */
85 do
86 {
87 if (!fgets(line, sizeof(line), fileHandle_p))
88 break;
89 if (feof(fileHandle_p))
90 break;
91 }
92 while (line[0] != '$');
93
94 if (feof(fileHandle_p))
95 break;
96
97 /* format */
98 if (!strncmp(&line[1], "MeshFormat", 10))
99 {
100 scan_ret = fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
101 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
102 }
103 /* nodes are read */
104 if (!strncmp(&line[1], "NOD", 3) || !strncmp(&line[1], "NOE", 3) || !strncmp(&line[1], "Nodes", 5))
105 {
106
107 scan_ret = fscanf(fileHandle_p, "%d", &numNodes);
108 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
109 if (!Dudley_noError())
110 break;
111 Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
112 if (!Dudley_noError())
113 break;
114 for (i0 = 0; i0 < numNodes; i0++)
115 {
116 if (1 == numDim)
117 {
118 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
119 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)], &rtmp0, &rtmp1);
120 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
121 }
122 else if (2 == numDim)
123 {
124 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
125 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
126 &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)], &rtmp0);
127 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
128 }
129 else if (3 == numDim)
130 {
131 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
132 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
133 &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)],
134 &mesh_p->Nodes->Coordinates[INDEX2(2, i0, numDim)]);
135 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
136 }
137 mesh_p->Nodes->globalDegreesOfFreedom[i0] = mesh_p->Nodes->Id[i0];
138 mesh_p->Nodes->Tag[i0] = 0;
139 }
140 }
141 /* elements */
142 else if (!strncmp(&line[1], "ELM", 3) || !strncmp(&line[1], "Elements", 8))
143 {
144
145 Dudley_ElementTypeId final_element_type = Dudley_NoRef;
146 Dudley_ElementTypeId final_face_element_type = Dudley_NoRef;
147 numElements = 0;
148 numFaceElements = 0;
149 scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
150 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
151
152 id = TMPMEMALLOC(totalNumElements, index_t);
153 tag = TMPMEMALLOC(totalNumElements, index_t);
154
155 element_type = TMPMEMALLOC(totalNumElements, Dudley_ElementTypeId);
156 vertices = TMPMEMALLOC(totalNumElements * MAX_numNodes_gmsh, index_t);
157 if (!
158 (Dudley_checkPtr(id) || Dudley_checkPtr(tag) || Dudley_checkPtr(element_type)
159 || Dudley_checkPtr(vertices)))
160 {
161 /* read all in */
162 for (e = 0; e < totalNumElements; e++)
163 {
164 scan_ret = fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
165 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
166 switch (gmsh_type)
167 {
168 case 1: /* line order 1 */
169 element_type[e] = Dudley_Line2;
170 element_dim = 1;
171 numNodesPerElement = 2;
172 break;
173 case 2: /* triangle order 1 */
174 element_type[e] = Dudley_Tri3;
175 numNodesPerElement = 3;
176 element_dim = 2;
177 break;
178 case 4: /* tetrahedron order 1 */
179 element_type[e] = Dudley_Tet4;
180 numNodesPerElement = 4;
181 element_dim = 3;
182 break;
183 case 15: /* point */
184 element_type[e] = Dudley_Point1;
185 numNodesPerElement = 1;
186 element_dim = 0;
187 break;
188 default:
189 element_type[e] = Dudley_NoRef;
190 sprintf(error_msg, "Unexected gmsh element type %d in mesh file %s.", gmsh_type, fname);
191 Dudley_setError(IO_ERROR, error_msg);
192 }
193 if (element_dim == numDim)
194 {
195 if (final_element_type == Dudley_NoRef)
196 {
197 final_element_type = element_type[e];
198 }
199 else if (final_element_type != element_type[e])
200 {
201 sprintf(error_msg, "Dudley can handle a single type of internal elements only.");
202 Dudley_setError(IO_ERROR, error_msg);
203 break;
204 }
205 numElements++;
206 }
207 else if (element_dim == numDim - 1)
208 {
209 if (final_face_element_type == Dudley_NoRef)
210 {
211 final_face_element_type = element_type[e];
212 }
213 else if (final_face_element_type != element_type[e])
214 {
215 sprintf(error_msg, "Dudley can handle a single type of face elements only.");
216 Dudley_setError(IO_ERROR, error_msg);
217 break;
218 }
219 numFaceElements++;
220 }
221
222 if (version <= 1.0)
223 {
224 scan_ret = fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
225 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
226 partition_id = 1;
227 if (numNodesPerElement2 != numNodesPerElement)
228 {
229 sprintf(error_msg, "Illegal number of nodes for element %d in mesh file %s.", id[e],
230 fname);
231 Dudley_setError(IO_ERROR, error_msg);
232 }
233 }
234 else
235 {
236 scan_ret = fscanf(fileHandle_p, "%d", &numTags);
237 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
238 elementary_id = tag[e] = partition_id = 1;
239 numNodesPerElement2 = -1;
240 for (j = 0; j < numTags; j++)
241 {
242 scan_ret = fscanf(fileHandle_p, "%d", &itmp);
243 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
244 if (j == 0)
245 {
246 tag[e] = itmp;
247 }
248 else if (j == 1)
249 {
250 elementary_id = itmp;
251 }
252 else if (j == 2)
253 {
254 partition_id = itmp;
255 }
256 /* ignore any other tags */
257 }
258 }
259 if (!Dudley_noError())
260 break;
261 for (j = 0; j < numNodesPerElement; j++)
262 {
263 scan_ret = fscanf(fileHandle_p, "%d", &vertices[INDEX2(j, e, MAX_numNodes_gmsh)]);
264 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
265 }
266 }
267 /* all elements have been read, now we have to identify the elements for dudley */
268
269 if (Dudley_noError())
270 {
271 /* first we have to identify the elements to define Elementis and FaceElements */
272 if (final_element_type == Dudley_NoRef)
273 {
274 if (numDim == 1)
275 {
276 final_element_type = Dudley_Line2;
277 }
278 else if (numDim == 2)
279 {
280 final_element_type = Dudley_Tri3;
281 }
282 else if (numDim == 3)
283 {
284 final_element_type = Dudley_Tet4;
285 }
286 }
287 if (final_face_element_type == Dudley_NoRef)
288 {
289 if (numDim == 1)
290 {
291 final_face_element_type = Dudley_Point1;
292 }
293 else if (numDim == 2)
294 {
295 final_face_element_type = Dudley_Line2;
296 }
297 else if (numDim == 3)
298 {
299 final_face_element_type = Dudley_Tri3;
300 }
301 }
302 mesh_p->Elements = Dudley_ElementFile_alloc(final_element_type, mpi_info);
303 mesh_p->FaceElements = Dudley_ElementFile_alloc(final_face_element_type, mpi_info);
304 mesh_p->Points = Dudley_ElementFile_alloc(Dudley_Point1, mpi_info);
305 if (Dudley_noError())
306 {
307 Dudley_ElementFile_allocTable(mesh_p->Elements, numElements);
308 Dudley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
309 Dudley_ElementFile_allocTable(mesh_p->Points, 0);
310 if (Dudley_noError())
311 {
312 mesh_p->Elements->minColor = 0;
313 mesh_p->Elements->maxColor = numElements - 1;
314 mesh_p->FaceElements->minColor = 0;
315 mesh_p->FaceElements->maxColor = numFaceElements - 1;
316 mesh_p->Points->minColor = 0;
317 mesh_p->Points->maxColor = 0;
318 numElements = 0;
319 numFaceElements = 0;
320 for (e = 0; e < totalNumElements; e++)
321 {
322 if (element_type[e] == final_element_type)
323 {
324 mesh_p->Elements->Id[numElements] = id[e];
325 mesh_p->Elements->Tag[numElements] = tag[e];
326 mesh_p->Elements->Color[numElements] = numElements;
327 mesh_p->Elements->Owner[numElements] = 0;
328 for (j = 0; j < mesh_p->Elements-> /*referenceElementSet-> */ numNodes; ++j)
329 {
330 mesh_p->Elements->Nodes[INDEX2
331 (j, numElements,
332 mesh_p->
333 Elements-> /*referenceElementSet-> */ numNodes)] =
334 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
335 }
336 numElements++;
337 }
338 else if (element_type[e] == final_face_element_type)
339 {
340 mesh_p->FaceElements->Id[numFaceElements] = id[e];
341 mesh_p->FaceElements->Tag[numFaceElements] = tag[e];
342 mesh_p->FaceElements->Color[numFaceElements] = numFaceElements;
343 mesh_p->FaceElements->Owner[numFaceElements] = 0;
344 for (j = 0; j < mesh_p->FaceElements-> /*referenceElementSet-> */ numNodes; ++j)
345 {
346 mesh_p->FaceElements->Nodes[INDEX2
347 (j, numFaceElements,
348 mesh_p->
349 FaceElements-> /*referenceElementSet-> */
350 numNodes)] =
351 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
352 }
353 numFaceElements++;
354 }
355 }
356 }
357 }
358 }
359 }
360 /* and clean up */
361 TMPMEMFREE(id);
362 TMPMEMFREE(tag);
363 TMPMEMFREE(element_type);
364 TMPMEMFREE(vertices);
365 }
366 /* name tags (thanks to Antoine Lefebvre, antoine.lefebvre2@mail.mcgill.ca ) */
367 else if (!strncmp(&line[1], "PhysicalNames", 13)) {
368 char name[LenString_MAX+1];
369 index_t tag_key;
370 scan_ret = fscanf(fileHandle_p, "%d", &numTags);
371 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
372 if (! Dudley_noError()) break;
373 for (i0 = 0; i0 < numTags; i0++) {
374 scan_ret = fscanf(fileHandle_p, "%d %d %s\n", &itmp, &tag_key, name);
375 FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
376 if (! (itmp == 2)) Dudley_setError(IO_ERROR,"Dudley_Mesh_readGmsh: expecting two entries per physical name.");
377 if ( strlen(name) < 3 ) Dudley_setError(IO_ERROR,"Dudley_Mesh_readGmsh: illegal tagname (\" missing?)");
378 if (! Dudley_noError()) break;
379 name[strlen(name)-1]='\0';
380 Dudley_Mesh_addTagMap(mesh_p,&name[1],tag_key);
381 }
382 }
383 /* serach for end of data block */
384 do
385 {
386 if (!fgets(line, sizeof(line), fileHandle_p))
387 {
388 sprintf(error_msg, "Unexected end of file in %s", fname);
389 Dudley_setError(IO_ERROR, error_msg);
390 }
391 if (feof(fileHandle_p))
392 {
393 sprintf(error_msg, "Unexected end of file in %s", fname);
394 Dudley_setError(IO_ERROR, error_msg);
395 }
396 if (!Dudley_noError())
397 break;
398 }
399 while (line[0] != '$');
400 }
401
402 /* close file */
403 fclose(fileHandle_p);
404 /* clean up */
405 if (!Dudley_noError())
406 {
407 Dudley_Mesh_free(mesh_p);
408 return NULL;
409 }
410 /* resolve id's : */
411 if (Dudley_noError())
412 Dudley_Mesh_resolveNodeIds(mesh_p);
413 /* rearrange elements: */
414 if (Dudley_noError())
415 Dudley_Mesh_prepare(mesh_p, optimize);
416 /* free up memory */
417 Esys_MPIInfo_free(mpi_info);
418 if (!Dudley_noError())
419 {
420 Dudley_Mesh_free(mesh_p);
421 return NULL;
422 }
423 else
424 {
425 return mesh_p;
426 }
427 }
428 }

  ViewVC Help
Powered by ViewVC 1.1.26