/[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 1745 - (hide annotations)
Wed Sep 3 00:41:24 2008 UTC (11 years ago) by ksteube
File MIME type: text/plain
File size: 35482 byte(s)
MPI mesh_read was crashing if tags not present

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26