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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1170 - (show 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 /*
2 ************************************************************
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 */
12
13 /**************************************************************/
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 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize_labeling) {
31
32 dim_t numNodes, numDim, numEle, i0, i1;
33 index_t tag_key;
34 Finley_Mesh *mesh_p=NULL;
35 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
36 char error_msg[LenErrorMsg_MAX];
37 double time0=Finley_timer();
38 FILE *fileHandle_p = NULL;
39 ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
40
41 Finley_resetError();
42 #ifdef PASO_MPI
43 /* TODO */
44 Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");`
45 return NULL;
46 #endif
47
48 /* get file handle */
49 fileHandle_p = fopen(fname, "r");
50 if (fileHandle_p==NULL) {
51 sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
52 Finley_setError(IO_ERROR,error_msg);
53 return NULL;
54 }
55
56 /* read header */
57 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58 fscanf(fileHandle_p, frm, name);
59
60 /* get the nodes */
61
62 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63 /* allocate mesh */
64 mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order);
65 if (! Finley_noError()) return NULL;
66
67 Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
68 if (! Finley_noError()) return NULL;
69
70 if (1 == numDim) {
71 for (i0 = 0; i0 < numNodes; i0++)
72 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 } else if (2 == numDim) {
76 for (i0 = 0; i0 < numNodes; i0++)
77 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 } else if (3 == numDim) {
82 for (i0 = 0; i0 < numNodes; i0++)
83 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 } /* if else else */
89
90 /* get the element type */
91
92 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
93 typeID=Finley_RefElement_getTypeId(element_type);
94 if (typeID==NoType) {
95 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
96 Finley_setError(VALUE_ERROR,error_msg);
97 return NULL;
98 }
99 /* read the elements */
100 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
101 Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
102 mesh_p->Elements->minColor=0;
103 mesh_p->Elements->maxColor=numEle-1;
104 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 fscanf(fileHandle_p, " %d",
109 &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 faceTypeID=Finley_RefElement_getTypeId(element_type);
117 faceTypeID=Finley_RefElement_getTypeId(element_type);
118 if (faceTypeID==NoType) {
119 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
120 Finley_setError(VALUE_ERROR,error_msg);
121 return NULL;
122 }
123 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order);
124 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
125 mesh_p->FaceElements->minColor=0;
126 mesh_p->FaceElements->maxColor=numEle-1;
127 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 fscanf(fileHandle_p, " %d",
132 &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 contactTypeID=Finley_RefElement_getTypeId(element_type);
140 if (contactTypeID==NoType) {
141 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
142 Finley_setError(VALUE_ERROR,error_msg);
143 return NULL;
144 }
145 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order);
146 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
147 mesh_p->ContactElements->minColor=0;
148 mesh_p->ContactElements->maxColor=numEle-1;
149 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 fscanf(fileHandle_p, " %d",
154 &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 pointTypeID=Finley_RefElement_getTypeId(element_type);
162 if (pointTypeID==NoType) {
163 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
164 Finley_setError(VALUE_ERROR,error_msg);
165 return NULL;
166 }
167 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order);
168 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
169 mesh_p->Points->minColor=0;
170 mesh_p->Points->maxColor=numEle-1;
171 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 fscanf(fileHandle_p, " %d",
176 &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 /* 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 /* close file */
189
190 fclose(fileHandle_p);
191
192 /* resolve id's : */
193
194 if (Finley_noError()) {
195 Finley_Mesh_resolveNodeIds(mesh_p);
196 }
197
198 /* rearrange elements: */
199
200 if (Finley_noError()) {
201 Finley_Mesh_prepare(mesh_p);
202 }
203
204 /* optimize node labeling*/
205
206 if (Finley_noError()) {
207 if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
208 }
209
210 /* that's it */
211 #ifdef Finley_TRACE
212 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
213 #endif
214 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 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