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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4257 - (hide annotations)
Wed Feb 27 03:42:40 2013 UTC (6 years, 1 month ago) by jfenwick
Original Path: branches/doubleplusgood/dudley/src/Mesh_readGmsh.c
File MIME type: text/plain
File size: 13230 byte(s)
Some simple experiments for c++ Finley

1 gross 935
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4154 * Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
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 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     * Development since 2012 by School of Earth Sciences
13     *
14     *****************************************************************************/
15 gross 935
16 jfenwick 3981 /************************************************************************************/
17 gross 935
18 jfenwick 3086 /* Dudley: read mesh */
19 gross 935
20 jfenwick 3981 /************************************************************************************/
21 gross 935
22     #include "Mesh.h"
23     #include <stdio.h>
24    
25 jfenwick 3086 #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 gross 935
27 jfenwick 3981 /************************************************************************************/
28 gross 935
29 jfenwick 3086 /* reads a mesh from a Dudley file of name fname */
30 gross 935
31 gross 2934 #define MAX_numNodes_gmsh 20
32 gross 935
33 jfenwick 3224 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 gross 935
37 jfenwick 3224 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 jfenwick 3086 #ifdef Dudley_TRACE
47 jfenwick 3224 double time0 = Dudley_timer();
48 jfenwick 1986 #endif
49 jfenwick 3224 FILE *fileHandle_p = NULL;
50 caltinay 3247 Dudley_ElementTypeId *element_type = NULL;
51 gross 935
52 jfenwick 3227 Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
53 jfenwick 3224 Dudley_resetError();
54     if (mpi_info->size > 1)
55     {
56     Dudley_setError(IO_ERROR, "reading GMSH with MPI is not supported yet.");
57 jfenwick 3227 Esys_MPIInfo_free(mpi_info);
58 jfenwick 3224 return NULL;
59     }
60     else
61     {
62 gross 2722
63 jfenwick 3224 /* allocate mesh */
64 gross 997
65 jfenwick 3224 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 jfenwick 3227 Esys_MPIInfo_free(mpi_info);
76 jfenwick 3224 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 caltinay 3247 Dudley_ElementTypeId final_element_type = Dudley_NoRef;
146     Dudley_ElementTypeId final_face_element_type = Dudley_NoRef;
147 jfenwick 3224 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 caltinay 3247 element_type = TMPMEMALLOC(totalNumElements, Dudley_ElementTypeId);
156 jfenwick 3224 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 caltinay 3247 element_type[e] = Dudley_Line2;
170 jfenwick 3224 element_dim = 1;
171     numNodesPerElement = 2;
172     break;
173 caltinay 3247 case 2: /* triangle order 1 */
174     element_type[e] = Dudley_Tri3;
175 jfenwick 3224 numNodesPerElement = 3;
176     element_dim = 2;
177     break;
178     case 4: /* tetrahedron order 1 */
179 caltinay 3247 element_type[e] = Dudley_Tet4;
180 jfenwick 3224 numNodesPerElement = 4;
181     element_dim = 3;
182     break;
183     case 15: /* point */
184 caltinay 3247 element_type[e] = Dudley_Point1;
185 jfenwick 3224 numNodesPerElement = 1;
186     element_dim = 0;
187     break;
188     default:
189 caltinay 3247 element_type[e] = Dudley_NoRef;
190 jfenwick 3224 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 caltinay 3247 if (final_element_type == Dudley_NoRef)
196 jfenwick 3224 {
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 caltinay 3247 if (final_face_element_type == Dudley_NoRef)
210 jfenwick 3224 {
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 caltinay 3247 if (final_element_type == Dudley_NoRef)
273 jfenwick 3224 {
274     if (numDim == 1)
275     {
276 caltinay 3247 final_element_type = Dudley_Line2;
277 jfenwick 3224 }
278     else if (numDim == 2)
279     {
280 caltinay 3247 final_element_type = Dudley_Tri3;
281 jfenwick 3224 }
282     else if (numDim == 3)
283     {
284 caltinay 3247 final_element_type = Dudley_Tet4;
285 jfenwick 3224 }
286     }
287 caltinay 3247 if (final_face_element_type == Dudley_NoRef)
288 jfenwick 3224 {
289     if (numDim == 1)
290     {
291 caltinay 3247 final_face_element_type = Dudley_Point1;
292 jfenwick 3224 }
293     else if (numDim == 2)
294     {
295 caltinay 3247 final_face_element_type = Dudley_Line2;
296 jfenwick 3224 }
297     else if (numDim == 3)
298     {
299 caltinay 3247 final_face_element_type = Dudley_Tri3;
300 jfenwick 3224 }
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 caltinay 3247 mesh_p->Points = Dudley_ElementFile_alloc(Dudley_Point1, mpi_info);
305 jfenwick 3224 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 jfenwick 3920 }
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 jfenwick 3921 if (! Dudley_noError()) break;
373 jfenwick 3920 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 jfenwick 3921 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 jfenwick 3920 name[strlen(name)-1]='\0';
380     Dudley_Mesh_addTagMap(mesh_p,&name[1],tag_key);
381     }
382     }
383 jfenwick 3224 /* 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 jfenwick 3227 Esys_MPIInfo_free(mpi_info);
418 jfenwick 3224 if (!Dudley_noError())
419     {
420     Dudley_Mesh_free(mesh_p);
421     return NULL;
422     }
423     else
424     {
425     return mesh_p;
426     }
427     }
428 gross 935 }

  ViewVC Help
Powered by ViewVC 1.1.26