/[escript]/branches/domexper/dudley/src/Mesh_readGmsh.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_readGmsh.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1814 - (hide annotations)
Tue Sep 30 00:21:21 2008 UTC (10 years, 6 months ago) by gross
Original Path: trunk/finley/src/Mesh_readGmsh.c
File MIME type: text/plain
File size: 15386 byte(s)
Contact elements have the right tyoe now.
1 gross 935
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 gross 935
14 ksteube 1811
15 gross 935 /**************************************************************/
16    
17 ksteube 1312 /* Finley: read mesh */
18 gross 935
19     /**************************************************************/
20    
21     #include "Mesh.h"
22     #include <stdio.h>
23    
24    
25     /**************************************************************/
26    
27     /* reads a mesh from a Finley file of name fname */
28    
29     #define MAX_numNodes_gmsh 10
30    
31 ksteube 1312 Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize) {
32 gross 935
33     double version = 1.0;
34     int format = 0, size = sizeof(double);
35     dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;
36 phornby 1628 index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;
37 gross 1028 index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;
38 gross 935 Finley_Mesh *mesh_p=NULL;
39 gross 1028 char line[LenString_MAX+1];
40 gross 935 char error_msg[LenErrorMsg_MAX];
41     double rtmp0, rtmp1;
42     double time0=Finley_timer();
43 gross 1028 FILE * fileHandle_p = NULL;
44     ElementTypeId* element_type=NULL;
45 gross 935
46 ksteube 1312 Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
47 gross 935 Finley_resetError();
48 ksteube 1312 if (mpi_info->size>1) {
49     Finley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");
50     Paso_MPIInfo_free( mpi_info );
51 gross 935 return NULL;
52 ksteube 1312 } else {
53 gross 997
54 ksteube 1312 /* allocate mesh */
55 gross 935
56 ksteube 1312 mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order, mpi_info);
57     if (! Finley_noError()) return NULL;
58    
59     /* get file handle */
60     fileHandle_p = fopen(fname, "r");
61     if (fileHandle_p==NULL) {
62     sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
63     Finley_setError(IO_ERROR,error_msg);
64     Paso_MPIInfo_free( mpi_info );
65     return NULL;
66     }
67    
68     /* start reading */
69     while(1) {
70     if (! Finley_noError()) break;
71     /* find line staring with $ */
72     do {
73     if( ! fgets(line, sizeof(line), fileHandle_p) ) break;
74     if(feof(fileHandle_p)) break;
75     } while(line[0] != '$');
76    
77     if (feof(fileHandle_p)) break;
78    
79    
80     /* format */
81     if (!strncmp(&line[1], "MeshFormat", 10)) {
82     fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
83     }
84     /* nodes are read */
85     if ( !strncmp(&line[1], "NOD", 3) ||
86     !strncmp(&line[1], "NOE", 3) ||
87     !strncmp(&line[1], "Nodes", 5) ) {
88    
89     fscanf(fileHandle_p, "%d", &numNodes);
90     if (! Finley_noError()) break;
91     Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
92     if (! Finley_noError()) break;
93     for (i0 = 0; i0 < numNodes; i0++) {
94     if (1 == numDim) {
95     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
96     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
97     &rtmp0,
98     &rtmp1);
99     } else if (2 == numDim) {
100     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
101     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
102     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
103     &rtmp0);
104     } else if (3 == numDim) {
105     fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
106     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
107     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
108     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
109    
110    
111     }
112     mesh_p->Nodes->globalDegreesOfFreedom[i0]=mesh_p->Nodes->Id[i0];
113     mesh_p->Nodes->Tag[i0]=0;
114 gross 935 }
115 ksteube 1312 }
116     /* elements */
117     else if(!strncmp(&line[1], "ELM", 3) ||
118     !strncmp(&line[1], "Elements", 8)) {
119    
120     ElementTypeId final_element_type = NoType;
121     ElementTypeId final_face_element_type = NoType;
122 gross 1814 ElementTypeId contact_element_type = NoType;
123 ksteube 1312 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 gross 1814 if (final_face_element_type == Line2) {
268     contact_element_type=Line2_Contact;
269     } else if (final_face_element_type == Line3) {
270     contact_element_type=Line3_Contact;
271     } else if (final_face_element_type == Tri3) {
272     contact_element_type=Tri3_Contact;
273     } else if (final_face_element_type == Tri6) {
274     contact_element_type=Tri6_Contact;
275     } else {
276     contact_element_type=Point1_Contact;
277     }
278    
279 ksteube 1312 mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);
280     mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);
281 gross 1814 mesh_p->ContactElements=Finley_ElementFile_alloc(contact_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);
282 ksteube 1312 mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order, mpi_info);
283     if (Finley_noError()) {
284     Finley_ElementFile_allocTable(mesh_p->Elements, numElements);
285     Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
286     Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
287     Finley_ElementFile_allocTable(mesh_p->Points, 0);
288     if (Finley_noError()) {
289     mesh_p->Elements->minColor=0;
290     mesh_p->Elements->maxColor=numElements-1;
291     mesh_p->FaceElements->minColor=0;
292     mesh_p->FaceElements->maxColor=numFaceElements-1;
293     mesh_p->ContactElements->minColor=0;
294     mesh_p->ContactElements->maxColor=0;
295     mesh_p->Points->minColor=0;
296     mesh_p->Points->maxColor=0;
297     numElements=0;
298     numFaceElements=0;
299     for(e = 0; e < totalNumElements; e++) {
300     if (element_type[e] == final_element_type) {
301     mesh_p->Elements->Id[numElements]=id[e];
302     mesh_p->Elements->Tag[numElements]=tag[e];
303     mesh_p->Elements->Color[numElements]=numElements;
304 gross 1740 mesh_p->Elements->Owner[numElements]=0;
305 ksteube 1312 for (j = 0; j< mesh_p->Elements->ReferenceElement->Type->numNodes; ++j) {
306     mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
307     }
308     numElements++;
309     } else if (element_type[e] == final_face_element_type) {
310     mesh_p->FaceElements->Id[numFaceElements]=id[e];
311     mesh_p->FaceElements->Tag[numFaceElements]=tag[e];
312     mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;
313 gross 1740 mesh_p->FaceElements->Owner[numElements]=0;
314 ksteube 1312 for (j = 0; j< mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j) {
315     mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
316     }
317     numFaceElements++;
318     }
319 gross 935 }
320 ksteube 1312 }
321 gross 935 }
322 ksteube 1312 }
323 gross 935 }
324 ksteube 1312 /* and clean up */
325     TMPMEMFREE(id);
326     TMPMEMFREE(tag);
327     TMPMEMFREE(element_type);
328     TMPMEMFREE(vertices);
329 gross 935 }
330 ksteube 1312 /* serach for end of data block */
331     do {
332     if (!fgets(line, sizeof(line), fileHandle_p)) {
333     sprintf(error_msg,"Unexected end of file in %s",fname);
334     Finley_setError(IO_ERROR,error_msg);
335     }
336     if (feof(fileHandle_p)) {
337     sprintf(error_msg,"Unexected end of file in %s",fname);
338     Finley_setError(IO_ERROR,error_msg);
339     }
340     if (! Finley_noError()) break;
341     } while(line[0] != '$');
342     }
343    
344     /* close file */
345     fclose(fileHandle_p);
346     /* clean up */
347     if (! Finley_noError()) {
348     Finley_Mesh_free(mesh_p);
349     return NULL;
350     }
351     /* resolve id's : */
352     Finley_Mesh_resolveNodeIds(mesh_p);
353     /* rearrange elements: */
354     Finley_Mesh_prepare(mesh_p, optimize);
355     /* that's it */
356     #ifdef Finley_TRACE
357     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
358     #endif
359     Paso_MPIInfo_free( mpi_info );
360     return mesh_p;
361 gross 935 }
362     }

  ViewVC Help
Powered by ViewVC 1.1.26