/[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 1771 - (hide annotations)
Mon Sep 8 22:47:55 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 36015 byte(s)
Removed failing test until we fix a bug.
Fixed several compiler warnings.

1 jgs 150
2 ksteube 1312 /* $Id$ */
3 jgs 82
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 jgs 82
16     /**************************************************************/
17    
18 ksteube 1312 /* Finley: read mesh */
19 jgs 82
20     /**************************************************************/
21    
22 ksteube 1771 #include <ctype.h>
23 jgs 82 #include "Mesh.h"
24    
25     /**************************************************************/
26    
27     /* reads a mesh from a Finley file of name fname */
28    
29 ksteube 1312 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize)
30 jgs 82
31 ksteube 1312 {
32    
33     Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
34 jgs 123 dim_t numNodes, numDim, numEle, i0, i1;
35 gross 1044 index_t tag_key;
36 bcumming 730 Finley_Mesh *mesh_p=NULL;
37 jgs 150 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
38     char error_msg[LenErrorMsg_MAX];
39 jgs 82 double time0=Finley_timer();
40 gross 1028 FILE *fileHandle_p = NULL;
41     ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
42 btully 1170
43 jgs 150 Finley_resetError();
44 jgs 82
45 ksteube 1312 if (mpi_info->size > 1) {
46     Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");
47     } else {
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     Paso_MPIInfo_free( mpi_info );
54     return NULL;
55 gross 1044 }
56 ksteube 1312
57     /* read header */
58     sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
59     fscanf(fileHandle_p, frm, name);
60    
61     /* get the nodes */
62    
63     fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
64     /* allocate mesh */
65     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
66     if (Finley_noError()) {
67    
68     /* read nodes */
69     Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
70     if (Finley_noError()) {
71     if (1 == numDim) {
72     for (i0 = 0; i0 < numNodes; i0++)
73     fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
74     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
75     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
76     } else if (2 == numDim) {
77     for (i0 = 0; i0 < numNodes; i0++)
78     fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
79     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
80     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
81     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
82     } else if (3 == numDim) {
83     for (i0 = 0; i0 < numNodes; i0++)
84     fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
85     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
86     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
87     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
88     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
89     } /* if else else */
90     }
91     /* read elements */
92     if (Finley_noError()) {
93    
94     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
95     typeID=Finley_RefElement_getTypeId(element_type);
96     if (typeID==NoType) {
97     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
98     Finley_setError(VALUE_ERROR,error_msg);
99     } else {
100     /* read the elements */
101     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
102     if (Finley_noError()) {
103     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
104     mesh_p->Elements->minColor=0;
105     mesh_p->Elements->maxColor=numEle-1;
106     if (Finley_noError()) {
107     for (i0 = 0; i0 < numEle; i0++) {
108     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
109     mesh_p->Elements->Color[i0]=i0;
110 gross 1739 mesh_p->Elements->Owner[i0]=0;
111 ksteube 1312 for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
112     fscanf(fileHandle_p, " %d",
113     &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
114     } /* for i1 */
115     fscanf(fileHandle_p, "\n");
116     } /* for i0 */
117     }
118     }
119     }
120     }
121     /* get the face elements */
122     if (Finley_noError()) {
123     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
124     faceTypeID=Finley_RefElement_getTypeId(element_type);
125     if (faceTypeID==NoType) {
126     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
127     Finley_setError(VALUE_ERROR,error_msg);
128     } else {
129     mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
130     if (Finley_noError()) {
131     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
132     if (Finley_noError()) {
133     mesh_p->FaceElements->minColor=0;
134     mesh_p->FaceElements->maxColor=numEle-1;
135     for (i0 = 0; i0 < numEle; i0++) {
136     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
137     mesh_p->FaceElements->Color[i0]=i0;
138 gross 1739 mesh_p->FaceElements->Owner[i0]=0;
139 ksteube 1312 for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
140     fscanf(fileHandle_p, " %d",
141     &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
142     } /* for i1 */
143     fscanf(fileHandle_p, "\n");
144     } /* for i0 */
145     }
146     }
147     }
148     }
149     /* get the Contact face element */
150     if (Finley_noError()) {
151     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
152     contactTypeID=Finley_RefElement_getTypeId(element_type);
153     if (contactTypeID==NoType) {
154     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
155     Finley_setError(VALUE_ERROR,error_msg);
156     } else {
157     mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
158     if (Finley_noError()) {
159     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
160     if (Finley_noError()) {
161     mesh_p->ContactElements->minColor=0;
162     mesh_p->ContactElements->maxColor=numEle-1;
163     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 gross 1739 mesh_p->ContactElements->Owner[i0]=0;
167 ksteube 1312 for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
168     fscanf(fileHandle_p, " %d",
169     &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
170     } /* for i1 */
171     fscanf(fileHandle_p, "\n");
172     } /* for i0 */
173     }
174     }
175     }
176     }
177     /* get the nodal element */
178     if (Finley_noError()) {
179     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
180     pointTypeID=Finley_RefElement_getTypeId(element_type);
181     if (pointTypeID==NoType) {
182     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
183     Finley_setError(VALUE_ERROR,error_msg);
184     }
185     mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
186     if (Finley_noError()) {
187     Finley_ElementFile_allocTable(mesh_p->Points, numEle);
188     if (Finley_noError()) {
189     mesh_p->Points->minColor=0;
190     mesh_p->Points->maxColor=numEle-1;
191     for (i0 = 0; i0 < numEle; i0++) {
192     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
193     mesh_p->Points->Color[i0]=i0;
194 gross 1739 mesh_p->Points->Owner[i0]=0;
195 ksteube 1312 for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
196     fscanf(fileHandle_p, " %d",
197     &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
198     } /* for i1 */
199     fscanf(fileHandle_p, "\n");
200     } /* for i0 */
201     }
202     }
203     }
204     /* get the name tags */
205     if (Finley_noError()) {
206     if (feof(fileHandle_p) == 0) {
207     fscanf(fileHandle_p, "%s\n", name);
208     while (feof(fileHandle_p) == 0) {
209     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
210     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
211     }
212     }
213     }
214     }
215     /* close file */
216     fclose(fileHandle_p);
217    
218     /* resolve id's : */
219     /* rearrange elements: */
220    
221     if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
222     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
223    
224     /* that's it */
225     #ifdef Finley_TRACE
226     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
227     #endif
228 gross 1044 }
229 jgs 82
230     /* that's it */
231 ksteube 1312 if (! Finley_noError()) {
232     Finley_Mesh_free(mesh_p);
233 gross 964 }
234 ksteube 1312 Paso_MPIInfo_free( mpi_info );
235 jgs 82 return mesh_p;
236     }
237 ksteube 1360
238 ksteube 1402 Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order, bool_t optimize)
239 ksteube 1360
240     {
241    
242     Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
243     dim_t numNodes, numDim, numEle, i0, i1;
244     Finley_Mesh *mesh_p=NULL;
245     char name[LenString_MAX],element_type[LenString_MAX],frm[20];
246     char error_msg[LenErrorMsg_MAX];
247     double time0=Finley_timer();
248     FILE *fileHandle_p = NULL;
249 ksteube 1725 ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
250 ksteube 1744 Finley_TagMap* tag_map;
251 phornby 1646 index_t tag_key;
252    
253 ksteube 1360 Finley_resetError();
254    
255     if (mpi_info->rank == 0) {
256     /* get file handle */
257     fileHandle_p = fopen(fname, "r");
258     if (fileHandle_p==NULL) {
259     sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
260     Finley_setError(IO_ERROR,error_msg);
261     Paso_MPIInfo_free( mpi_info );
262     return NULL;
263     }
264 ksteube 1402
265 ksteube 1360 /* read header */
266     sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
267     fscanf(fileHandle_p, frm, name);
268 ksteube 1402
269     /* get the number of nodes */
270 ksteube 1360 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
271     }
272    
273     #ifdef PASO_MPI
274 ksteube 1713 /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
275     if (mpi_info->size > 1) {
276 ksteube 1360 int temp1[3], error_code;
277     temp1[0] = numDim;
278     temp1[1] = numNodes;
279     temp1[2] = strlen(name) + 1;
280     error_code = MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
281     if (error_code != MPI_SUCCESS) {
282     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
283     return NULL;
284     }
285     numDim = temp1[0];
286     numNodes = temp1[1];
287     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
288     if (error_code != MPI_SUCCESS) {
289     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
290     return NULL;
291     }
292     }
293     #endif
294    
295     /* allocate mesh */
296     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
297     if (Finley_noError()) {
298 ksteube 1725 /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
299 phornby 1646 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
300 ksteube 1713 int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t); /* Stores the integer message data */
301     double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
302 ksteube 1360
303     /*
304 ksteube 1713 Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
305 ksteube 1402 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
306 ksteube 1360 First chunk sent to CPU 1, second to CPU 2, ...
307     Last chunk stays on CPU 0 (the master)
308     The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
309     */
310    
311     if (mpi_info->rank == 0) { /* Master */
312     for (;;) { /* Infinite loop */
313 ksteube 1713 for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
314     for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
315 ksteube 1360 chunkNodes = 0;
316     for (i1=0; i1<chunkSize; i1++) {
317 ksteube 1732 if (totalNodes >= numNodes) break; /* End of inner loop */
318 ksteube 1360 if (1 == numDim)
319     fscanf(fileHandle_p, "%d %d %d %le\n",
320 ksteube 1713 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
321 ksteube 1360 &tempCoords[i1*numDim+0]);
322     if (2 == numDim)
323     fscanf(fileHandle_p, "%d %d %d %le %le\n",
324 ksteube 1713 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
325 ksteube 1360 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
326     if (3 == numDim)
327     fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
328 ksteube 1713 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
329 ksteube 1360 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
330 ksteube 1713 totalNodes++; /* When do we quit the infinite loop? */
331     chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
332 ksteube 1360 }
333 ksteube 1713 if (chunkNodes > chunkSize) {
334     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
335     return NULL;
336     }
337     #ifdef PASO_MPI
338     /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
339 ksteube 1360 if (nextCPU < mpi_info->size) {
340 phornby 1646 int mpi_error;
341 ksteube 1713 tempInts[chunkSize*3] = chunkNodes; /* The message has one more int to send chunkNodes */
342     mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
343 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
344     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
345     return NULL;
346     }
347 ksteube 1713 mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
348 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
349     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
350     return NULL;
351     }
352     }
353 ksteube 1713 #endif
354 ksteube 1732 nextCPU++;
355     /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
356     if (nextCPU > mpi_info->size) break; /* End infinite loop */
357 ksteube 1360 } /* Infinite loop */
358     } /* End master */
359     else { /* Worker */
360     #ifdef PASO_MPI
361     /* Each worker receives two messages */
362     MPI_Status status;
363 ksteube 1660 int mpi_error;
364 ksteube 1713 mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
365 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
366     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
367     return NULL;
368     }
369 ksteube 1713 mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
370 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
371     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
372     return NULL;
373     }
374 ksteube 1713 chunkNodes = tempInts[chunkSize*3]; /* How many nodes are in this workers chunk? */
375 ksteube 1360 #endif
376     } /* Worker */
377 ksteube 1402
378     /* Copy node data from tempMem to mesh_p */
379     Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
380     if (Finley_noError()) {
381     for (i0=0; i0<chunkNodes; i0++) {
382     mesh_p->Nodes->Id[i0] = tempInts[0+i0];
383 ksteube 1713 mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize+i0];
384     mesh_p->Nodes->Tag[i0] = tempInts[chunkSize*2+i0];
385 ksteube 1402 for (i1=0; i1<numDim; i1++) {
386     mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
387     }
388     }
389 ksteube 1360 }
390    
391     TMPMEMFREE(tempInts);
392     TMPMEMFREE(tempCoords);
393    
394 ksteube 1402 /* read elements */
395 ksteube 1360
396 ksteube 1402 /* Read the element typeID */
397     if (Finley_noError()) {
398     if (mpi_info->rank == 0) {
399     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
400     typeID=Finley_RefElement_getTypeId(element_type);
401     }
402     #ifdef PASO_MPI
403 ksteube 1713 if (mpi_info->size > 1) {
404 ksteube 1725 int temp1[2], mpi_error;
405 ksteube 1660 temp1[0] = (int) typeID;
406 ksteube 1402 temp1[1] = numEle;
407     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
408     if (mpi_error != MPI_SUCCESS) {
409     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
410     return NULL;
411     }
412 ksteube 1660 typeID = (ElementTypeId) temp1[0];
413 ksteube 1402 numEle = temp1[1];
414     }
415     #endif
416     if (typeID==NoType) {
417     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
418     Finley_setError(VALUE_ERROR, error_msg);
419     }
420     }
421    
422 ksteube 1732 /* Allocate the ElementFile */
423 ksteube 1402 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
424     numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
425    
426 ksteube 1732 /* Read the element data */
427 ksteube 1402 if (Finley_noError()) {
428 ksteube 1725 int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
429     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
430     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
431 ksteube 1402 if (mpi_info->rank == 0) { /* Master */
432     for (;;) { /* Infinite loop */
433 ksteube 1725 for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
434 ksteube 1402 chunkEle = 0;
435     for (i0=0; i0<chunkSize; i0++) {
436 ksteube 1732 if (totalEle >= numEle) break; /* End inner loop */
437 ksteube 1402 fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
438     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
439     fscanf(fileHandle_p, "\n");
440     totalEle++;
441     chunkEle++;
442     }
443 ksteube 1725 #ifdef PASO_MPI
444     /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
445 ksteube 1402 if (nextCPU < mpi_info->size) {
446 phornby 1646 int mpi_error;
447 ksteube 1725 tempInts[chunkSize*(2+numNodes)] = chunkEle;
448     mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
449 ksteube 1402 if ( mpi_error != MPI_SUCCESS ) {
450     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
451     return NULL;
452     }
453     }
454 ksteube 1725 #endif
455 ksteube 1732 nextCPU++;
456     /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
457     if (nextCPU > mpi_info->size) break; /* End infinite loop */
458 ksteube 1402 } /* Infinite loop */
459     } /* End master */
460     else { /* Worker */
461     #ifdef PASO_MPI
462 ksteube 1725 /* Each worker receives one message */
463 ksteube 1402 MPI_Status status;
464 ksteube 1660 int mpi_error;
465 ksteube 1725 mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
466 ksteube 1402 if ( mpi_error != MPI_SUCCESS ) {
467     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
468     return NULL;
469     }
470 ksteube 1725 chunkEle = tempInts[chunkSize*(2+numNodes)];
471 ksteube 1402 #endif
472     } /* Worker */
473 ksteube 1725
474 ksteube 1402 /* Copy Element data from tempInts to mesh_p */
475     Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
476     mesh_p->Elements->minColor=0;
477     mesh_p->Elements->maxColor=chunkEle-1;
478     if (Finley_noError()) {
479     #pragma omp parallel for private (i0, i1)
480     for (i0=0; i0<chunkEle; i0++) {
481     mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
482     mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
483 gross 1739 mesh_p->Elements->Owner[i0] =mpi_info->rank;
484 ksteube 1725 mesh_p->Elements->Color[i0] = i0;
485 ksteube 1402 for (i1 = 0; i1 < numNodes; i1++) {
486     mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
487     }
488     }
489     }
490    
491     TMPMEMFREE(tempInts);
492 ksteube 1732 } /* end of Read the element data */
493    
494     /* read face elements */
495    
496     /* Read the element typeID */
497     if (Finley_noError()) {
498     if (mpi_info->rank == 0) {
499     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
500     typeID=Finley_RefElement_getTypeId(element_type);
501     }
502     #ifdef PASO_MPI
503     if (mpi_info->size > 1) {
504     int temp1[2], mpi_error;
505     temp1[0] = (int) typeID;
506     temp1[1] = numEle;
507     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
508     if (mpi_error != MPI_SUCCESS) {
509     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
510     return NULL;
511     }
512     typeID = (ElementTypeId) temp1[0];
513     numEle = temp1[1];
514     }
515     #endif
516     if (typeID==NoType) {
517     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
518     Finley_setError(VALUE_ERROR, error_msg);
519     }
520     }
521    
522     /* Allocate the ElementFile */
523     mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524     numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
525    
526     /* Read the face element data */
527     if (Finley_noError()) {
528     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
529     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
530     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
531     if (mpi_info->rank == 0) { /* Master */
532     for (;;) { /* Infinite loop */
533     for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
534     chunkEle = 0;
535     for (i0=0; i0<chunkSize; i0++) {
536     if (totalEle >= numEle) break; /* End inner loop */
537     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
538     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
539     fscanf(fileHandle_p, "\n");
540     totalEle++;
541     chunkEle++;
542     }
543     #ifdef PASO_MPI
544     /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
545     if (nextCPU < mpi_info->size) {
546     int mpi_error;
547     tempInts[chunkSize*(2+numNodes)] = chunkEle;
548     mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
549     if ( mpi_error != MPI_SUCCESS ) {
550     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
551     return NULL;
552     }
553     }
554     #endif
555     nextCPU++;
556     /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
557     if (nextCPU > mpi_info->size) break; /* End infinite loop */
558     } /* Infinite loop */
559     } /* End master */
560     else { /* Worker */
561     #ifdef PASO_MPI
562     /* Each worker receives one message */
563     MPI_Status status;
564     int mpi_error;
565     mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
566     if ( mpi_error != MPI_SUCCESS ) {
567     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
568     return NULL;
569     }
570     chunkEle = tempInts[chunkSize*(2+numNodes)];
571     #endif
572     } /* Worker */
573    
574     /* Copy Element data from tempInts to mesh_p */
575     Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
576     mesh_p->FaceElements->minColor=0;
577     mesh_p->FaceElements->maxColor=chunkEle-1;
578     if (Finley_noError()) {
579     #pragma omp parallel for private (i0, i1)
580     for (i0=0; i0<chunkEle; i0++) {
581     mesh_p->FaceElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
582     mesh_p->FaceElements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
583 gross 1739 mesh_p->FaceElements->Owner[i0] =mpi_info->rank;
584 ksteube 1732 mesh_p->FaceElements->Color[i0] = i0;
585     for (i1 = 0; i1 < numNodes; i1++) {
586     mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
587     }
588     }
589     }
590    
591     TMPMEMFREE(tempInts);
592     } /* end of Read the face element data */
593    
594     /* read contact elements */
595    
596     /* Read the element typeID */
597     if (Finley_noError()) {
598     if (mpi_info->rank == 0) {
599     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
600     typeID=Finley_RefElement_getTypeId(element_type);
601     }
602     #ifdef PASO_MPI
603     if (mpi_info->size > 1) {
604     int temp1[2], mpi_error;
605     temp1[0] = (int) typeID;
606     temp1[1] = numEle;
607     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
608     if (mpi_error != MPI_SUCCESS) {
609     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
610     return NULL;
611     }
612     typeID = (ElementTypeId) temp1[0];
613     numEle = temp1[1];
614     }
615     #endif
616     if (typeID==NoType) {
617     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
618     Finley_setError(VALUE_ERROR, error_msg);
619     }
620     }
621    
622     /* Allocate the ElementFile */
623     mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
624     numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
625    
626     /* Read the contact element data */
627     if (Finley_noError()) {
628     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
629     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
630     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
631     if (mpi_info->rank == 0) { /* Master */
632     for (;;) { /* Infinite loop */
633     for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
634     chunkEle = 0;
635     for (i0=0; i0<chunkSize; i0++) {
636     if (totalEle >= numEle) break; /* End inner loop */
637     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
638     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
639     fscanf(fileHandle_p, "\n");
640     totalEle++;
641     chunkEle++;
642     }
643     #ifdef PASO_MPI
644     /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
645     if (nextCPU < mpi_info->size) {
646     int mpi_error;
647     tempInts[chunkSize*(2+numNodes)] = chunkEle;
648     mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
649     if ( mpi_error != MPI_SUCCESS ) {
650     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
651     return NULL;
652     }
653     }
654     #endif
655     nextCPU++;
656     /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
657     if (nextCPU > mpi_info->size) break; /* End infinite loop */
658     } /* Infinite loop */
659     } /* End master */
660     else { /* Worker */
661     #ifdef PASO_MPI
662     /* Each worker receives one message */
663     MPI_Status status;
664     int mpi_error;
665     mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
666     if ( mpi_error != MPI_SUCCESS ) {
667     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
668     return NULL;
669     }
670     chunkEle = tempInts[chunkSize*(2+numNodes)];
671     #endif
672     } /* Worker */
673    
674     /* Copy Element data from tempInts to mesh_p */
675     Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
676     mesh_p->ContactElements->minColor=0;
677     mesh_p->ContactElements->maxColor=chunkEle-1;
678     if (Finley_noError()) {
679     #pragma omp parallel for private (i0, i1)
680     for (i0=0; i0<chunkEle; i0++) {
681     mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
682     mesh_p->ContactElements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
683 gross 1739 mesh_p->ContactElements->Owner[i0] =mpi_info->rank;
684 ksteube 1732 mesh_p->ContactElements->Color[i0] = i0;
685     for (i1 = 0; i1 < numNodes; i1++) {
686     mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
687     }
688     }
689     }
690    
691     TMPMEMFREE(tempInts);
692     } /* end of Read the contact element data */
693    
694     /* read nodal elements */
695    
696     /* Read the element typeID */
697     if (Finley_noError()) {
698     if (mpi_info->rank == 0) {
699     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
700     typeID=Finley_RefElement_getTypeId(element_type);
701     }
702     #ifdef PASO_MPI
703     if (mpi_info->size > 1) {
704     int temp1[2], mpi_error;
705     temp1[0] = (int) typeID;
706     temp1[1] = numEle;
707     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
708     if (mpi_error != MPI_SUCCESS) {
709     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
710     return NULL;
711     }
712     typeID = (ElementTypeId) temp1[0];
713     numEle = temp1[1];
714     }
715     #endif
716     if (typeID==NoType) {
717     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
718     Finley_setError(VALUE_ERROR, error_msg);
719     }
720     }
721    
722     /* Allocate the ElementFile */
723     mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
724     numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
725    
726     /* Read the nodal element data */
727     if (Finley_noError()) {
728     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
729     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
730     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
731     if (mpi_info->rank == 0) { /* Master */
732     for (;;) { /* Infinite loop */
733     for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
734     chunkEle = 0;
735     for (i0=0; i0<chunkSize; i0++) {
736     if (totalEle >= numEle) break; /* End inner loop */
737     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
738     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
739     fscanf(fileHandle_p, "\n");
740     totalEle++;
741     chunkEle++;
742     }
743     #ifdef PASO_MPI
744     /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
745     if (nextCPU < mpi_info->size) {
746     int mpi_error;
747     tempInts[chunkSize*(2+numNodes)] = chunkEle;
748     mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
749     if ( mpi_error != MPI_SUCCESS ) {
750     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
751     return NULL;
752     }
753     }
754     #endif
755     nextCPU++;
756     /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
757     if (nextCPU > mpi_info->size) break; /* End infinite loop */
758     } /* Infinite loop */
759     } /* End master */
760     else { /* Worker */
761     #ifdef PASO_MPI
762     /* Each worker receives one message */
763     MPI_Status status;
764     int mpi_error;
765     mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
766     if ( mpi_error != MPI_SUCCESS ) {
767     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
768     return NULL;
769     }
770     chunkEle = tempInts[chunkSize*(2+numNodes)];
771     #endif
772     } /* Worker */
773    
774     /* Copy Element data from tempInts to mesh_p */
775     Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
776     mesh_p->Points->minColor=0;
777     mesh_p->Points->maxColor=chunkEle-1;
778     if (Finley_noError()) {
779     #pragma omp parallel for private (i0, i1)
780     for (i0=0; i0<chunkEle; i0++) {
781     mesh_p->Points->Id[i0] = tempInts[i0*(2+numNodes)+0];
782     mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
783 gross 1739 mesh_p->Points->Owner[i0] =mpi_info->rank;
784 ksteube 1732 mesh_p->Points->Color[i0] = i0;
785     for (i1 = 0; i1 < numNodes; i1++) {
786     mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
787     }
788     }
789     }
790    
791     TMPMEMFREE(tempInts);
792     } /* end of Read the nodal element data */
793    
794 ksteube 1744 /* get the name tags */
795 ksteube 1754 if (Finley_noError()) {
796     char *remainder, *ptr;
797 ksteube 1745 int tag_key, len, error_code;
798 ksteube 1754 long cur_pos, end_pos;
799 ksteube 1744 if (mpi_info->rank == 0) { /* Master */
800     /* Read the word 'Tag' */
801 ksteube 1754 if (! feof(fileHandle_p)) fscanf(fileHandle_p, "%s\n", name);
802     /* Read rest of file in one chunk, after using seek to find length */
803     cur_pos = ftell(fileHandle_p);
804     fseek(fileHandle_p, 0L, SEEK_END);
805     end_pos = ftell(fileHandle_p);
806     fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
807     remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
808     if (! feof(fileHandle_p)) fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
809     remainder[end_pos-cur_pos] = 0;
810     len = strlen(remainder);
811     while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
812     len = strlen(remainder);
813 ksteube 1744 }
814     #ifdef PASO_MPI
815     error_code = MPI_Bcast (&len, 1, MPI_INT, 0, mpi_info->comm);
816     if (error_code != MPI_SUCCESS) {
817     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
818     return NULL;
819     }
820 ksteube 1754 if (mpi_info->rank != 0) {
821     remainder = TMPMEMALLOC(len+1, char);
822     remainder[0] = 0;
823     }
824 ksteube 1744 error_code = MPI_Bcast (remainder, len+1, MPI_CHAR, 0, mpi_info->comm);
825     if (error_code != MPI_SUCCESS) {
826     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tags failed");
827     return NULL;
828     }
829     #endif
830 ksteube 1754 if (remainder[0]) {
831     ptr = remainder;
832     do {
833     sscanf(ptr, "%s %d\n", name, &tag_key);
834     if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
835     ptr++;
836     } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
837     TMPMEMFREE(remainder);
838     }
839 ksteube 1744 }
840 ksteube 1732
841 ksteube 1725 }
842 ksteube 1402
843 ksteube 1360 /* close file */
844 ksteube 1402 if (mpi_info->rank == 0) fclose(fileHandle_p);
845    
846 ksteube 1360 /* resolve id's : */
847     /* rearrange elements: */
848 ksteube 1402
849 ksteube 1360 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
850     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
851 ksteube 1402
852 ksteube 1360 /* that's it */
853     #ifdef Finley_TRACE
854     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
855     #endif
856    
857     /* that's it */
858     if (! Finley_noError()) {
859     Finley_Mesh_free(mesh_p);
860     }
861     Paso_MPIInfo_free( mpi_info );
862     return mesh_p;
863     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26