/[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 1028 - (hide annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/plain
File size: 7930 byte(s)
modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
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 jgs 123 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {
31 jgs 82
32 jgs 123 dim_t numNodes, numDim, numEle, i0, i1;
33 bcumming 730 Finley_Mesh *mesh_p=NULL;
34 jgs 150 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
35     char error_msg[LenErrorMsg_MAX];
36 jgs 82 double time0=Finley_timer();
37 gross 1028 FILE *fileHandle_p = NULL;
38     ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
39 bcumming 730
40 jgs 150 Finley_resetError();
41 jgs 82
42     /* get file handle */
43 gross 1028 fileHandle_p = fopen(fname, "r");
44 jgs 82 if (fileHandle_p==NULL) {
45 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
46 jgs 150 Finley_setError(IO_ERROR,error_msg);
47 jgs 82 return NULL;
48     }
49    
50     /* read header */
51 jgs 150 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
52     fscanf(fileHandle_p, frm, name);
53 gross 934
54 gross 939 /* get the nodes */
55    
56     fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
57 gross 934 /* allocate mesh */
58     #ifndef PASO_MPI
59     mesh_p = Finley_Mesh_alloc(name,numDim,order);
60 gross 939 if (! Finley_noError()) return NULL;
61 gross 934 #else
62     /* TODO */
63     #endif
64 jgs 82
65 bcumming 730 #ifndef PASO_MPI
66 gross 934 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
67 jgs 150 if (! Finley_noError()) return NULL;
68 bcumming 730 #else
69     /* TODO */
70     #endif
71 jgs 82
72     if (1 == numDim) {
73     for (i0 = 0; i0 < numNodes; i0++)
74 gross 934 fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
75     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
76     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
77 jgs 82 } else if (2 == numDim) {
78     for (i0 = 0; i0 < numNodes; i0++)
79 gross 934 fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
80     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
81     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
82     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
83 jgs 82 } else if (3 == numDim) {
84     for (i0 = 0; i0 < numNodes; i0++)
85 gross 934 fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
86     &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
87     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
88     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
89     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
90 jgs 82 } /* if else else */
91    
92     /* get the element type */
93    
94     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
95 gross 1028 typeID=Finley_RefElement_getTypeId(element_type);
96 jgs 82 if (typeID==NoType) {
97 gross 1028 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
98 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
99 jgs 82 return NULL;
100     }
101     /* read the elements */
102 bcumming 730 #ifndef PASO_MPI
103 jgs 82 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);
104 bcumming 730 #else
105     /* TODO */
106     #endif
107 jgs 82 Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
108 jgs 123 mesh_p->Elements->minColor=0;
109     mesh_p->Elements->maxColor=numEle-1;
110 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
111     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
112     mesh_p->Elements->Color[i0]=i0;
113     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
114 jgs 126 fscanf(fileHandle_p, " %d",
115 jgs 82 &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
116     } /* for i1 */
117     fscanf(fileHandle_p, "\n");
118     } /* for i0 */
119    
120     /* get the face elements */
121     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
122 gross 1028 faceTypeID=Finley_RefElement_getTypeId(element_type);
123     faceTypeID=Finley_RefElement_getTypeId(element_type);
124 jgs 82 if (faceTypeID==NoType) {
125 gross 1028 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
126 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
127 jgs 82 return NULL;
128     }
129 bcumming 730 #ifndef PASO_MPI
130 jgs 82 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);
131 bcumming 730 #else
132     /* TODO */
133     #endif
134 jgs 82 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
135 jgs 123 mesh_p->FaceElements->minColor=0;
136     mesh_p->FaceElements->maxColor=numEle-1;
137 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
138     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
139     mesh_p->FaceElements->Color[i0]=i0;
140     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
141 jgs 126 fscanf(fileHandle_p, " %d",
142 jgs 82 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
143     } /* for i1 */
144     fscanf(fileHandle_p, "\n");
145     } /* for i0 */
146    
147     /* get the Contact face element */
148     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
149 gross 1028 contactTypeID=Finley_RefElement_getTypeId(element_type);
150 jgs 82 if (contactTypeID==NoType) {
151 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
152 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
153 jgs 82 return NULL;
154     }
155 bcumming 730 #ifndef PASO_MPI
156 jgs 82 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);
157 bcumming 730 #else
158     /* TODO */
159     #endif
160 jgs 82 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
161 jgs 123 mesh_p->ContactElements->minColor=0;
162     mesh_p->ContactElements->maxColor=numEle-1;
163 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
164     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
165     mesh_p->ContactElements->Color[i0]=i0;
166     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
167 jgs 126 fscanf(fileHandle_p, " %d",
168 jgs 82 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
169     } /* for i1 */
170     fscanf(fileHandle_p, "\n");
171     } /* for i0 */
172    
173     /* get the nodal element */
174     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
175 gross 1028 pointTypeID=Finley_RefElement_getTypeId(element_type);
176 jgs 82 if (pointTypeID==NoType) {
177 gross 1028 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
178 jgs 150 Finley_setError(VALUE_ERROR,error_msg);
179 jgs 82 return NULL;
180     }
181 bcumming 730 #ifndef PASO_MPI
182 jgs 82 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);
183 bcumming 730 #else
184     /* TODO */
185     #endif
186 jgs 82 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
187 jgs 123 mesh_p->Points->minColor=0;
188     mesh_p->Points->maxColor=numEle-1;
189 jgs 82 for (i0 = 0; i0 < numEle; i0++) {
190     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
191     mesh_p->Points->Color[i0]=i0;
192     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
193 jgs 126 fscanf(fileHandle_p, " %d",
194 jgs 82 &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
195     } /* for i1 */
196     fscanf(fileHandle_p, "\n");
197     } /* for i0 */
198    
199    
200     /* close file */
201    
202     fclose(fileHandle_p);
203    
204     /* resolve id's : */
205 jgs 126
206 jgs 82 Finley_Mesh_resolveNodeIds(mesh_p);
207 jgs 126
208 jgs 82 /* rearrange elements: */
209 jgs 126
210 jgs 82 Finley_Mesh_prepare(mesh_p);
211 jgs 126
212 jgs 82 /* that's it */
213 jgs 149 #ifdef Finley_TRACE
214 jgs 82 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
215 jgs 149 #endif
216 gross 964 if (Finley_noError()) {
217     if ( ! Finley_Mesh_isPrepared(mesh_p)) {
218     Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
219     }
220     } else {
221     Finley_Mesh_dealloc(mesh_p);
222     }
223 jgs 82 return mesh_p;
224     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26