/[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 1062 - (hide annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years, 1 month ago) by gross
Original Path: trunk/finley/src/Mesh_readGmsh.c
File MIME type: text/plain
File size: 13913 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
1 gross 935 /*
2     ************************************************************
3     * Copyright 2006, 2007 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 3.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11     */
12    
13     /**************************************************************/
14    
15     /* Finley: read mesh */
16    
17     /**************************************************************/
18    
19     /* Author: gross@access.edu.au */
20     /* Version: $Id: Mesh_read.c 730 2006-05-15 04:03:49Z bcumming $ */
21    
22     /**************************************************************/
23    
24     #include "Mesh.h"
25     #include <stdio.h>
26    
27    
28     /**************************************************************/
29    
30     /* reads a mesh from a Finley file of name fname */
31    
32     #define MAX_numNodes_gmsh 10
33    
34     Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize_labeling) {
35    
36     double version = 1.0;
37     int format = 0, size = sizeof(double);
38     dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;
39     index_t e, i0, j, gmsh_type, partition_id, itmp, final_element_type, elementary_id;
40 gross 1028 index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;
41 gross 935 Finley_Mesh *mesh_p=NULL;
42 gross 1028 char line[LenString_MAX+1];
43 gross 935 char error_msg[LenErrorMsg_MAX];
44     double rtmp0, rtmp1;
45     double time0=Finley_timer();
46 gross 1028 FILE * fileHandle_p = NULL;
47     ElementTypeId* element_type=NULL;
48 gross 935
49     Finley_resetError();
50     #if PASO_MPI
51     sprintf(error_msg,"reading GMSH with MPI is not supported yet.");
52     Finley_setError(IO_ERROR,error_msg);
53     return NULL;
54    
55     #else
56     /* allocate mesh */
57 gross 997
58 gross 1062 mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order);
59 gross 935 if (! Finley_noError()) return NULL;
60    
61     /* get file handle */
62 gross 1028 fileHandle_p = fopen(fname, "r");
63 gross 935 if (fileHandle_p==NULL) {
64     sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
65     Finley_setError(IO_ERROR,error_msg);
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->degreeOfFreedom[i0]=mesh_p->Nodes->Id[i0];
114     mesh_p->Nodes->Tag[i0]=0;
115     }
116     }
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 gross 1028 numElements=0;
124     numFaceElements=0;
125 gross 935 fscanf(fileHandle_p, "%d", &totalNumElements);
126    
127 gross 1028 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 gross 935 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 gross 1028 /* ignore any other tags */
234 gross 935 }
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     }
239     /* all elements have been read, now we have to identify the elements for finley */
240 gross 997
241 gross 935 if (Finley_noError()) {
242 gross 1003 /* first we have to identify the elements to define Elementis and FaceElements */
243     if (final_element_type == NoType) {
244     if (numDim==1) {
245     final_element_type=Line2;
246     } else if (numDim==2) {
247     final_element_type=Tri3;
248     } else if (numDim==2) {
249     final_element_type=Tet4;
250     }
251     }
252     if (final_face_element_type == NoType) {
253     if (numDim==1) {
254     final_face_element_type=Point1;
255     } else if (numDim==2) {
256     final_face_element_type=Line2;
257     } else if (numDim==2) {
258     final_face_element_type=Tri3;
259     }
260     }
261 gross 1062 mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order);
262     mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order);
263     mesh_p->ContactElements=Finley_ElementFile_alloc(Point1_Contact,mesh_p->order, mesh_p->reduced_order);
264     mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order);
265 gross 935 if (Finley_noError()) {
266     Finley_ElementFile_allocTable(mesh_p->Elements, numElements);
267     Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
268     Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
269     Finley_ElementFile_allocTable(mesh_p->Points, 0);
270     if (Finley_noError()) {
271     mesh_p->Elements->minColor=0;
272     mesh_p->Elements->maxColor=numElements-1;
273     mesh_p->FaceElements->minColor=0;
274     mesh_p->FaceElements->maxColor=numFaceElements-1;
275     mesh_p->ContactElements->minColor=0;
276     mesh_p->ContactElements->maxColor=0;
277     mesh_p->Points->minColor=0;
278     mesh_p->Points->maxColor=0;
279     numElements=0;
280     numFaceElements=0;
281     for(e = 0; e < totalNumElements; e++) {
282     if (element_type[e] == final_element_type) {
283     mesh_p->Elements->Id[numElements]=id[e];
284     mesh_p->Elements->Tag[numElements]=tag[e];
285     mesh_p->Elements->Color[numElements]=numElements;
286     for (j = 0; j< mesh_p->Elements->ReferenceElement->Type->numNodes; ++j)
287     mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=
288     vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
289     numElements++;
290     } else if (element_type[e] == final_face_element_type) {
291     mesh_p->FaceElements->Id[numFaceElements]=id[e];
292     mesh_p->FaceElements->Tag[numFaceElements]=tag[e];
293     mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;
294     for (j = 0; j< mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j)
295     mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=
296     vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
297     numFaceElements++;
298     }
299     }
300     }
301     }
302     }
303     }
304     /* and clean up */
305     TMPMEMFREE(id);
306     TMPMEMFREE(tag);
307     TMPMEMFREE(element_type);
308     TMPMEMFREE(vertices);
309     }
310     /* serach for end of data block */
311     do {
312     if (!fgets(line, sizeof(line), fileHandle_p)) {
313     sprintf(error_msg,"Unexected end of file in %s",fname);
314     Finley_setError(IO_ERROR,error_msg);
315     }
316     if (feof(fileHandle_p)) {
317 gross 940 sprintf(error_msg,"Unexected end of file in %s",fname);
318 gross 935 Finley_setError(IO_ERROR,error_msg);
319     }
320     if (! Finley_noError()) break;
321     } while(line[0] != '$');
322     }
323    
324     /* close file */
325     fclose(fileHandle_p);
326     /* clean up */
327     if (! Finley_noError()) {
328     Finley_Mesh_dealloc(mesh_p);
329     return NULL;
330     }
331     /* resolve id's : */
332     Finley_Mesh_resolveNodeIds(mesh_p);
333     /* rearrange elements: */
334     Finley_Mesh_prepare(mesh_p);
335     /* that's it */
336     #ifdef Finley_TRACE
337     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
338     #endif
339     #endif
340 gross 964 if (Finley_noError()) {
341     if ( ! Finley_Mesh_isPrepared(mesh_p)) {
342     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
343     }
344     }
345 gross 935 return mesh_p;
346     }

  ViewVC Help
Powered by ViewVC 1.1.26