/[escript]/trunk/finley/src/Mesh_read.c
ViewVC logotype

Annotation of /trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1170 - (hide annotations)
Fri May 25 05:24:57 2007 UTC (12 years, 5 months ago) by btully
File MIME type: text/plain
File size: 8434 byte(s)
Added new methods for integration of Tri's and Tet's For Order 2 and 3.
NB: There is a conflict in the testing of the getNormal() function. On consultation with Lutz, we have concluded this is an issue with the test rather than the new methods.

1 jgs 150 /*
2 elspeth 616 ************************************************************
3     * Copyright 2006 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 jgs 150 */
12    
13 jgs 82 /**************************************************************/
14    
15     /* Finley: read mesh */
16    
17     /**************************************************************/
18    
19     /* Author: gross@access.edu.au */
20     /* Version: $Id$ */
21    
22     /**************************************************************/
23    
24     #include "Mesh.h"
25    
26     /**************************************************************/
27    
28     /* reads a mesh from a Finley file of name fname */
29    
30 gross 1062 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize_labeling) {
31 jgs 82
32 jgs 123 dim_t numNodes, numDim, numEle, i0, i1;
33 gross 1044 index_t tag_key;
34 bcumming 730 Finley_Mesh *mesh_p=NULL;
35 jgs 150 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
36     char error_msg[LenErrorMsg_MAX];
37 jgs 82 double time0=Finley_timer();
38 gross 1028 FILE *fileHandle_p = NULL;
39     ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
40 btully 1170
41 jgs 150 Finley_resetError();
42 gross 1044 #ifdef PASO_MPI
43     /* TODO */
44     Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");`
45 gross 1156 return NULL;
46 gross 1044 #endif
47 jgs 82
48     /* get file handle */
49 gross 1028 fileHandle_p = fopen(fname, "r");
50 jgs 82 if (fileHandle_p==NULL) {
51 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
52 jgs 150 Finley_setError(IO_ERROR,error_msg);
53 jgs 82 return NULL;
54     }
55    
56     /* read header */
57 jgs 150 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58     fscanf(fileHandle_p, frm, name);
59 gross 934
60 gross 939 /* get the nodes */
61    
62     fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63 gross 934 /* allocate mesh */
64 gross 1062 mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order);
65 gross 939 if (! Finley_noError()) return NULL;
66 jgs 82
67 gross 934 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
68 jgs 150 if (! Finley_noError()) return NULL;
69 jgs 82
70     if (1 == numDim) {
71     for (i0 = 0; i0 < numNodes; i0++)
72 gross 934 fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
73     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
74     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
75 jgs 82 } else if (2 == numDim) {
76     for (i0 = 0; i0 < numNodes; i0++)
77 gross 934 fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
78     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
79     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
80     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
81 jgs 82 } else if (3 == numDim) {
82     for (i0 = 0; i0 < numNodes; i0++)
83 gross 934 fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
84     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
85     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
86     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
87     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
88 jgs 82 } /* if else else */
89    
90     /* get the element type */
91    
92     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
93 gross 1028 typeID=Finley_RefElement_getTypeId(element_type);
94 jgs 82 if (typeID==NoType) {
95 gross 1028 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
96 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
97 jgs 82 return NULL;
98     }
99     /* read the elements */
100 gross 1062 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
101 jgs 82 Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
102 jgs 123 mesh_p->Elements->minColor=0;
103     mesh_p->Elements->maxColor=numEle-1;
104 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
105     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
106     mesh_p->Elements->Color[i0]=i0;
107     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
108 jgs 126 fscanf(fileHandle_p, " %d",
109 jgs 82 &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
110     } /* for i1 */
111     fscanf(fileHandle_p, "\n");
112     } /* for i0 */
113    
114     /* get the face elements */
115     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
116 gross 1028 faceTypeID=Finley_RefElement_getTypeId(element_type);
117     faceTypeID=Finley_RefElement_getTypeId(element_type);
118 jgs 82 if (faceTypeID==NoType) {
119 gross 1028 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
120 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
121 jgs 82 return NULL;
122     }
123 gross 1062 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order);
124 jgs 82 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
125 jgs 123 mesh_p->FaceElements->minColor=0;
126     mesh_p->FaceElements->maxColor=numEle-1;
127 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
128     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
129     mesh_p->FaceElements->Color[i0]=i0;
130     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
131 jgs 126 fscanf(fileHandle_p, " %d",
132 jgs 82 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
133     } /* for i1 */
134     fscanf(fileHandle_p, "\n");
135     } /* for i0 */
136    
137     /* get the Contact face element */
138     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
139 gross 1028 contactTypeID=Finley_RefElement_getTypeId(element_type);
140 jgs 82 if (contactTypeID==NoType) {
141 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
142 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
143 jgs 82 return NULL;
144     }
145 gross 1062 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order);
146 jgs 82 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
147 jgs 123 mesh_p->ContactElements->minColor=0;
148     mesh_p->ContactElements->maxColor=numEle-1;
149 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
150     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
151     mesh_p->ContactElements->Color[i0]=i0;
152     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
153 jgs 126 fscanf(fileHandle_p, " %d",
154 jgs 82 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
155     } /* for i1 */
156     fscanf(fileHandle_p, "\n");
157     } /* for i0 */
158    
159     /* get the nodal element */
160     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
161 gross 1028 pointTypeID=Finley_RefElement_getTypeId(element_type);
162 jgs 82 if (pointTypeID==NoType) {
163 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
164 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
165 jgs 82 return NULL;
166     }
167 gross 1062 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order);
168 jgs 82 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
169 jgs 123 mesh_p->Points->minColor=0;
170     mesh_p->Points->maxColor=numEle-1;
171 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
172     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
173     mesh_p->Points->Color[i0]=i0;
174     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
175 jgs 126 fscanf(fileHandle_p, " %d",
176 jgs 82 &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
177     } /* for i1 */
178     fscanf(fileHandle_p, "\n");
179     } /* for i0 */
180 gross 1044 /* get the name tags */
181     if (feof(fileHandle_p) == 0) {
182     fscanf(fileHandle_p, "%s\n", name);
183     while (feof(fileHandle_p) == 0) {
184     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
185     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
186     }
187     }
188 jgs 82 /* close file */
189    
190     fclose(fileHandle_p);
191    
192     /* resolve id's : */
193 jgs 126
194 gross 1156 if (Finley_noError()) {
195     Finley_Mesh_resolveNodeIds(mesh_p);
196     }
197 jgs 126
198 jgs 82 /* rearrange elements: */
199 jgs 126
200 gross 1156 if (Finley_noError()) {
201     Finley_Mesh_prepare(mesh_p);
202     }
203 jgs 126
204 gross 1156 /* optimize node labeling*/
205    
206     if (Finley_noError()) {
207     if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
208     }
209    
210 jgs 82 /* that's it */
211 jgs 149 #ifdef Finley_TRACE
212 jgs 82 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
213 jgs 149 #endif
214 gross 964 if (Finley_noError()) {
215     if ( ! Finley_Mesh_isPrepared(mesh_p)) {
216     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
217     }
218     } else {
219     Finley_Mesh_dealloc(mesh_p);
220     }
221 jgs 82 return mesh_p;
222     }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26