/[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 1387 - (hide annotations)
Fri Jan 11 07:45:26 2008 UTC (11 years, 3 months ago) by trankine
Original Path: temp/finley/src/Mesh_readGmsh.c
File MIME type: text/plain
File size: 14694 byte(s)
Restore the trunk that existed before the windows changes were committed to the (now moved to branches) old trunk.
1 gross 935
2 ksteube 1312 /* $Id$ */
3 gross 935
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 gross 935
16     /**************************************************************/
17    
18 ksteube 1312 /* Finley: read mesh */
19 gross 935
20     /**************************************************************/
21    
22     #include "Mesh.h"
23     #include <stdio.h>
24    
25    
26     /**************************************************************/
27    
28     /* reads a mesh from a Finley file of name fname */
29    
30     #define MAX_numNodes_gmsh 10
31    
32 ksteube 1312 Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize) {
33 gross 935
34     double version = 1.0;
35     int format = 0, size = sizeof(double);
36     dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;
37     index_t e, i0, j, gmsh_type, partition_id, itmp, final_element_type, elementary_id;
38 gross 1028 index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;
39 gross 935 Finley_Mesh *mesh_p=NULL;
40 gross 1028 char line[LenString_MAX+1];
41 gross 935 char error_msg[LenErrorMsg_MAX];
42     double rtmp0, rtmp1;
43     double time0=Finley_timer();
44 gross 1028 FILE * fileHandle_p = NULL;
45     ElementTypeId* element_type=NULL;
46 gross 935
47 ksteube 1312 Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
48 gross 935 Finley_resetError();
49 ksteube 1312 if (mpi_info->size>1) {
50     Finley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");
51     Paso_MPIInfo_free( mpi_info );
52 gross 935 return NULL;
53 ksteube 1312 } else {
54 gross 997
55 ksteube 1312 /* allocate mesh */
56 gross 935
57 ksteube 1312 mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order, mpi_info);
58     if (! Finley_noError()) return NULL;
59    
60     /* get file handle */
61     fileHandle_p = fopen(fname, "r");
62     if (fileHandle_p==NULL) {
63     sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
64     Finley_setError(IO_ERROR,error_msg);
65     Paso_MPIInfo_free( mpi_info );
66     return NULL;
67     }
68    
69     /* start reading */
70     while(1) {
71     if (! Finley_noError()) break;
72     /* find line staring with $ */
73     do {
74     if( ! fgets(line, sizeof(line), fileHandle_p) ) break;
75     if(feof(fileHandle_p)) break;
76     } while(line[0] != '$');
77    
78     if (feof(fileHandle_p)) break;
79    
80    
81     /* format */
82     if (!strncmp(&line[1], "MeshFormat", 10)) {
83     fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
84     }
85     /* nodes are read */
86     if ( !strncmp(&line[1], "NOD", 3) ||
87     !strncmp(&line[1], "NOE", 3) ||
88     !strncmp(&line[1], "Nodes", 5) ) {
89    
90     fscanf(fileHandle_p, "%d", &numNodes);
91     if (! Finley_noError()) break;
92     Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
93     if (! Finley_noError()) break;
94     for (i0 = 0; i0 < numNodes; i0++) {
95     if (1 == numDim) {
96     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
97     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
98     &rtmp0,
99     &rtmp1);
100     } else if (2 == numDim) {
101     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
102     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
103     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
104     &rtmp0);
105     } else if (3 == numDim) {
106     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
107     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
108     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
109     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
110    
111    
112     }
113     mesh_p->Nodes->globalDegreesOfFreedom[i0]=mesh_p->Nodes->Id[i0];
114     mesh_p->Nodes->Tag[i0]=0;
115 gross 935 }
116 ksteube 1312 }
117     /* elements */
118     else if(!strncmp(&line[1], "ELM", 3) ||
119     !strncmp(&line[1], "Elements", 8)) {
120    
121     ElementTypeId final_element_type = NoType;
122     ElementTypeId final_face_element_type = NoType;
123     numElements=0;
124     numFaceElements=0;
125     fscanf(fileHandle_p, "%d", &totalNumElements);
126    
127     id=TMPMEMALLOC(totalNumElements,index_t);
128     tag=TMPMEMALLOC(totalNumElements,index_t);
129    
130    
131     element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);
132     vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);
133     if (! (Finley_checkPtr(id) || Finley_checkPtr(tag) || Finley_checkPtr(element_type) || Finley_checkPtr(vertices) ) ) {
134     /* read all in */
135     for(e = 0; e < totalNumElements; e++) {
136     fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
137     switch (gmsh_type) {
138     case 1: /* line order 1 */
139     element_type[e]=Line2;
140     element_dim=1;
141     numNodesPerElement=2;
142     break;
143     case 2: /* traingle order 1 */
144     element_type[e]=Tri3;
145     numNodesPerElement= 3;
146     element_dim=2;
147     break;
148     case 3: /* quadrilateral order 1 */
149     element_type[e]=Rec4;
150     numNodesPerElement= 4;
151     element_dim=2;
152     break;
153     case 4: /* tetrahedron order 1 */
154     element_type[e]=Tet4;
155     numNodesPerElement= 4;
156     element_dim=3;
157     break;
158     case 5: /* hexahedron order 1 */
159     element_type[e]=Hex8;
160     numNodesPerElement= 8;
161     element_dim=3;
162     break;
163     case 8: /* line order 2 */
164     element_type[e]=Line3;
165     numNodesPerElement= 3;
166     element_dim=1;
167     break;
168     case 9: /* traingle order 2 */
169     element_type[e]=Tri6;
170     numNodesPerElement= 6;
171     element_dim=2;
172     break;
173     case 10: /* quadrilateral order 2 */
174     element_type[e]=Rec9;
175     numNodesPerElement= 9;
176     element_dim=2;
177     break;
178     case 11: /* tetrahedron order 2 */
179     element_type[e]=Tet10;
180     numNodesPerElement= 10;
181     element_dim=3;
182     break;
183     case 15 : /* point */
184     element_type[e]=Point1;
185     numNodesPerElement= 1;
186     element_dim=0;
187     break;
188     default:
189     element_type[e]=NoType;
190     sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);
191     Finley_setError(IO_ERROR,error_msg);
192     }
193     if (element_dim == numDim) {
194     if (final_element_type == NoType) {
195     final_element_type = element_type[e];
196     } else if (final_element_type != element_type[e]) {
197     sprintf(error_msg,"Finley can handle a single type of internal elements only.");
198     Finley_setError(IO_ERROR,error_msg);
199     break;
200     }
201     numElements++;
202     } else if (element_dim == numDim-1) {
203     if (final_face_element_type == NoType) {
204     final_face_element_type = element_type[e];
205     } else if (final_face_element_type != element_type[e]) {
206     sprintf(error_msg,"Finley can handle a single type of face elements only.");
207     Finley_setError(IO_ERROR,error_msg);
208     break;
209     }
210     numFaceElements++;
211     }
212    
213     if(version <= 1.0){
214     fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
215     partition_id = 1;
216     if (numNodesPerElement2 != numNodesPerElement) {
217     sprintf(error_msg,"Illegal number of nodes for element %d in mesh file %s.",id[e],fname);
218     Finley_setError(IO_ERROR,error_msg);
219     }
220     } else {
221     fscanf(fileHandle_p, "%d", &numTags);
222     elementary_id = tag[e] = partition_id = 1;
223     numNodesPerElement2=-1;
224     for(j = 0; j < numTags; j++){
225     fscanf(fileHandle_p, "%d", &itmp);
226     if (j == 0) {
227     tag[e] = itmp;
228     } else if (j == 1) {
229     elementary_id = itmp;
230     } else if (j == 2) {
231     partition_id = itmp;
232     }
233     /* ignore any other tags */
234     }
235     }
236     if (! Finley_noError()) break;
237     for(j = 0; j < numNodesPerElement; j++) fscanf(fileHandle_p, "%d", &vertices[INDEX2(j,e,MAX_numNodes_gmsh)]);
238     /* for tet10 the last two nodes need to be swapped */
239     if (element_type[e]==Tet10) {
240     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];
241     vertices[INDEX2(9,e,MAX_numNodes_gmsh)]=vertices[INDEX2(8,e,MAX_numNodes_gmsh)];
242     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;
243     }
244     }
245     /* all elements have been read, now we have to identify the elements for finley */
246    
247     if (Finley_noError()) {
248     /* first we have to identify the elements to define Elementis and FaceElements */
249 gross 935 if (final_element_type == NoType) {
250 ksteube 1312 if (numDim==1) {
251     final_element_type=Line2;
252     } else if (numDim==2) {
253     final_element_type=Tri3;
254     } else if (numDim==3) {
255     final_element_type=Tet4;
256     }
257 gross 935 }
258     if (final_face_element_type == NoType) {
259 ksteube 1312 if (numDim==1) {
260     final_face_element_type=Point1;
261     } else if (numDim==2) {
262     final_face_element_type=Line2;
263     } else if (numDim==3) {
264     final_face_element_type=Tri3;
265     }
266 gross 935 }
267 ksteube 1312 mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);
268     mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);
269     mesh_p->ContactElements=Finley_ElementFile_alloc(Point1_Contact,mesh_p->order, mesh_p->reduced_order, mpi_info);
270     mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order, mpi_info);
271     if (Finley_noError()) {
272     Finley_ElementFile_allocTable(mesh_p->Elements, numElements);
273     Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
274     Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
275     Finley_ElementFile_allocTable(mesh_p->Points, 0);
276     if (Finley_noError()) {
277     mesh_p->Elements->minColor=0;
278     mesh_p->Elements->maxColor=numElements-1;
279     mesh_p->FaceElements->minColor=0;
280     mesh_p->FaceElements->maxColor=numFaceElements-1;
281     mesh_p->ContactElements->minColor=0;
282     mesh_p->ContactElements->maxColor=0;
283     mesh_p->Points->minColor=0;
284     mesh_p->Points->maxColor=0;
285     numElements=0;
286     numFaceElements=0;
287     for(e = 0; e < totalNumElements; e++) {
288     if (element_type[e] == final_element_type) {
289     mesh_p->Elements->Id[numElements]=id[e];
290     mesh_p->Elements->Tag[numElements]=tag[e];
291     mesh_p->Elements->Color[numElements]=numElements;
292     for (j = 0; j< mesh_p->Elements->ReferenceElement->Type->numNodes; ++j) {
293     mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
294     }
295     numElements++;
296     } else if (element_type[e] == final_face_element_type) {
297     mesh_p->FaceElements->Id[numFaceElements]=id[e];
298     mesh_p->FaceElements->Tag[numFaceElements]=tag[e];
299     mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;
300     for (j = 0; j< mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j) {
301     mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
302     }
303     numFaceElements++;
304     }
305 gross 935 }
306 ksteube 1312 }
307 gross 935 }
308 ksteube 1312 }
309 gross 935 }
310 ksteube 1312 /* and clean up */
311     TMPMEMFREE(id);
312     TMPMEMFREE(tag);
313     TMPMEMFREE(element_type);
314     TMPMEMFREE(vertices);
315 gross 935 }
316 ksteube 1312 /* serach for end of data block */
317     do {
318     if (!fgets(line, sizeof(line), fileHandle_p)) {
319     sprintf(error_msg,"Unexected end of file in %s",fname);
320     Finley_setError(IO_ERROR,error_msg);
321     }
322     if (feof(fileHandle_p)) {
323     sprintf(error_msg,"Unexected end of file in %s",fname);
324     Finley_setError(IO_ERROR,error_msg);
325     }
326     if (! Finley_noError()) break;
327     } while(line[0] != '$');
328     }
329    
330     /* close file */
331     fclose(fileHandle_p);
332     /* clean up */
333     if (! Finley_noError()) {
334     Finley_Mesh_free(mesh_p);
335     return NULL;
336     }
337     /* resolve id's : */
338     Finley_Mesh_resolveNodeIds(mesh_p);
339     /* rearrange elements: */
340     Finley_Mesh_prepare(mesh_p, optimize);
341     /* that's it */
342     #ifdef Finley_TRACE
343     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
344     #endif
345     Paso_MPIInfo_free( mpi_info );
346     return mesh_p;
347 gross 935 }
348     }

  ViewVC Help
Powered by ViewVC 1.1.26