/[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 1713 - (hide annotations)
Wed Aug 20 07:03:45 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 29540 byte(s)
MPI read of mesh.fly: reduced memory requirements

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     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
110     fscanf(fileHandle_p, " %d",
111     &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
112     } /* for i1 */
113     fscanf(fileHandle_p, "\n");
114     } /* for i0 */
115     }
116     }
117     }
118     }
119     /* get the face elements */
120     if (Finley_noError()) {
121     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
122     faceTypeID=Finley_RefElement_getTypeId(element_type);
123     if (faceTypeID==NoType) {
124     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
125     Finley_setError(VALUE_ERROR,error_msg);
126     } else {
127     mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
128     if (Finley_noError()) {
129     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
130     if (Finley_noError()) {
131     mesh_p->FaceElements->minColor=0;
132     mesh_p->FaceElements->maxColor=numEle-1;
133     for (i0 = 0; i0 < numEle; i0++) {
134     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
135     mesh_p->FaceElements->Color[i0]=i0;
136     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
137     fscanf(fileHandle_p, " %d",
138     &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
139     } /* for i1 */
140     fscanf(fileHandle_p, "\n");
141     } /* for i0 */
142     }
143     }
144     }
145     }
146     /* get the Contact face element */
147     if (Finley_noError()) {
148     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
149     contactTypeID=Finley_RefElement_getTypeId(element_type);
150     if (contactTypeID==NoType) {
151     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
152     Finley_setError(VALUE_ERROR,error_msg);
153     } else {
154     mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
155     if (Finley_noError()) {
156     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
157     if (Finley_noError()) {
158     mesh_p->ContactElements->minColor=0;
159     mesh_p->ContactElements->maxColor=numEle-1;
160     for (i0 = 0; i0 < numEle; i0++) {
161     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
162     mesh_p->ContactElements->Color[i0]=i0;
163     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
164     fscanf(fileHandle_p, " %d",
165     &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
166     } /* for i1 */
167     fscanf(fileHandle_p, "\n");
168     } /* for i0 */
169     }
170     }
171     }
172     }
173     /* get the nodal element */
174     if (Finley_noError()) {
175     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
176     pointTypeID=Finley_RefElement_getTypeId(element_type);
177     if (pointTypeID==NoType) {
178     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
179     Finley_setError(VALUE_ERROR,error_msg);
180     }
181     mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
182     if (Finley_noError()) {
183     Finley_ElementFile_allocTable(mesh_p->Points, numEle);
184     if (Finley_noError()) {
185     mesh_p->Points->minColor=0;
186     mesh_p->Points->maxColor=numEle-1;
187     for (i0 = 0; i0 < numEle; i0++) {
188     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
189     mesh_p->Points->Color[i0]=i0;
190     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
191     fscanf(fileHandle_p, " %d",
192     &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
193     } /* for i1 */
194     fscanf(fileHandle_p, "\n");
195     } /* for i0 */
196     }
197     }
198     }
199     /* get the name tags */
200     if (Finley_noError()) {
201     if (feof(fileHandle_p) == 0) {
202     fscanf(fileHandle_p, "%s\n", name);
203     while (feof(fileHandle_p) == 0) {
204     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
205     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
206     }
207     }
208     }
209     }
210     /* close file */
211     fclose(fileHandle_p);
212    
213     /* resolve id's : */
214     /* rearrange elements: */
215    
216     if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
217     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
218    
219     /* that's it */
220     #ifdef Finley_TRACE
221     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
222     #endif
223 gross 1044 }
224 jgs 82
225     /* that's it */
226 ksteube 1312 if (! Finley_noError()) {
227     Finley_Mesh_free(mesh_p);
228 gross 964 }
229 ksteube 1312 Paso_MPIInfo_free( mpi_info );
230 jgs 82 return mesh_p;
231     }
232 ksteube 1360
233 ksteube 1402 Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order, bool_t optimize)
234 ksteube 1360
235     {
236    
237     Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
238     dim_t numNodes, numDim, numEle, i0, i1;
239     Finley_Mesh *mesh_p=NULL;
240     char name[LenString_MAX],element_type[LenString_MAX],frm[20];
241     char error_msg[LenErrorMsg_MAX];
242     double time0=Finley_timer();
243     FILE *fileHandle_p = NULL;
244 phornby 1646 ElementTypeId typeID;
245     ElementTypeId faceTypeID, contactTypeID, pointTypeID;
246     index_t tag_key;
247    
248 ksteube 1360 Finley_resetError();
249    
250     if (mpi_info->rank == 0) {
251     /* get file handle */
252     fileHandle_p = fopen(fname, "r");
253     if (fileHandle_p==NULL) {
254     sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
255     Finley_setError(IO_ERROR,error_msg);
256     Paso_MPIInfo_free( mpi_info );
257     return NULL;
258     }
259 ksteube 1402
260 ksteube 1360 /* read header */
261     sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
262     fscanf(fileHandle_p, frm, name);
263 ksteube 1402
264     /* get the number of nodes */
265 ksteube 1360 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
266     }
267    
268     #ifdef PASO_MPI
269 ksteube 1713 /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
270     if (mpi_info->size > 1) {
271 ksteube 1360 int temp1[3], error_code;
272     temp1[0] = numDim;
273     temp1[1] = numNodes;
274     temp1[2] = strlen(name) + 1;
275     error_code = MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
276     if (error_code != MPI_SUCCESS) {
277     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
278     return NULL;
279     }
280     numDim = temp1[0];
281     numNodes = temp1[1];
282     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
283     if (error_code != MPI_SUCCESS) {
284     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
285     return NULL;
286     }
287     }
288     #endif
289    
290     /* allocate mesh */
291     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
292     if (Finley_noError()) {
293 ksteube 1713 /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large for that */
294 phornby 1646 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
295 ksteube 1713 int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t); /* Stores the integer message data */
296     double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
297 ksteube 1360
298     /*
299 ksteube 1713 Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
300 ksteube 1402 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
301 ksteube 1360 First chunk sent to CPU 1, second to CPU 2, ...
302     Last chunk stays on CPU 0 (the master)
303     The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
304     */
305    
306     if (mpi_info->rank == 0) { /* Master */
307     for (;;) { /* Infinite loop */
308 ksteube 1713 for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
309     for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
310 ksteube 1360 chunkNodes = 0;
311     for (i1=0; i1<chunkSize; i1++) {
312 ksteube 1713 if (totalNodes >= numNodes) break; /* Maybe end the infinite loop */
313 ksteube 1360 if (1 == numDim)
314     fscanf(fileHandle_p, "%d %d %d %le\n",
315 ksteube 1713 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
316 ksteube 1360 &tempCoords[i1*numDim+0]);
317     if (2 == numDim)
318     fscanf(fileHandle_p, "%d %d %d %le %le\n",
319 ksteube 1713 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
320 ksteube 1360 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
321     if (3 == numDim)
322     fscanf(fileHandle_p, "%d %d %d %le %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], &tempCoords[i1*numDim+2]);
325 ksteube 1713 totalNodes++; /* When do we quit the infinite loop? */
326     chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
327 ksteube 1360 }
328 ksteube 1713 if (chunkNodes > chunkSize) {
329     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
330     return NULL;
331     }
332     #ifdef PASO_MPI
333     /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
334 ksteube 1360 if (nextCPU < mpi_info->size) {
335 phornby 1646 int mpi_error;
336 ksteube 1713 tempInts[chunkSize*3] = chunkNodes; /* The message has one more int to send chunkNodes */
337     mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
338 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
339     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
340     return NULL;
341     }
342 ksteube 1713 mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
343 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
344     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
345     return NULL;
346     }
347     nextCPU++;
348     }
349 ksteube 1713 #endif
350     if (totalNodes >= numNodes) break; /* Maybe end the infinite loop */
351 ksteube 1360 } /* Infinite loop */
352     } /* End master */
353     else { /* Worker */
354     #ifdef PASO_MPI
355     /* Each worker receives two messages */
356     MPI_Status status;
357 ksteube 1660 int mpi_error;
358 ksteube 1713 mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
359 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
360     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
361     return NULL;
362     }
363 ksteube 1713 mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
364 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
365     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
366     return NULL;
367     }
368 ksteube 1713 chunkNodes = tempInts[chunkSize*3]; /* How many nodes are in this workers chunk? */
369 ksteube 1360 #endif
370     } /* Worker */
371 ksteube 1402
372 ksteube 1713 #if 1
373 ksteube 1402 /* Display the temp mem for debugging */
374 ksteube 1713 printf("ksteube Nodes tempInts\n");
375     for (i0=0; i0<chunkSize*3; i0++) {
376 ksteube 1360 printf(" %2d", tempInts[i0]);
377 ksteube 1713 if (i0%chunkSize==chunkSize-1) printf("\n");
378 ksteube 1360 }
379     printf("ksteube tempCoords:\n");
380 ksteube 1713 for (i0=0; i0<chunkSize*numDim; i0++) {
381 ksteube 1360 printf(" %20.15e", tempCoords[i0]);
382     if (i0%numDim==numDim-1) printf("\n");
383     }
384     #endif
385    
386 ksteube 1713 printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);
387    
388 ksteube 1402 /* Copy node data from tempMem to mesh_p */
389     Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
390     if (Finley_noError()) {
391     for (i0=0; i0<chunkNodes; i0++) {
392     mesh_p->Nodes->Id[i0] = tempInts[0+i0];
393 ksteube 1713 mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize+i0];
394     mesh_p->Nodes->Tag[i0] = tempInts[chunkSize*2+i0];
395 ksteube 1402 for (i1=0; i1<numDim; i1++) {
396     mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
397     }
398     }
399 ksteube 1360 }
400    
401     TMPMEMFREE(tempInts);
402     TMPMEMFREE(tempCoords);
403    
404 ksteube 1402 /* read elements */
405 ksteube 1360
406 ksteube 1402 /* Read the element typeID */
407     if (Finley_noError()) {
408     if (mpi_info->rank == 0) {
409     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
410     typeID=Finley_RefElement_getTypeId(element_type);
411     }
412     #ifdef PASO_MPI
413 ksteube 1713 if (mpi_info->size > 1) {
414 ksteube 1660 int temp1[3], mpi_error;
415     temp1[0] = (int) typeID;
416 ksteube 1402 temp1[1] = numEle;
417     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
418     if (mpi_error != MPI_SUCCESS) {
419     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
420     return NULL;
421     }
422 ksteube 1660 typeID = (ElementTypeId) temp1[0];
423 ksteube 1402 numEle = temp1[1];
424     }
425     #endif
426     if (typeID==NoType) {
427     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
428     Finley_setError(VALUE_ERROR, error_msg);
429     }
430     }
431    
432     /* Read the element data */
433     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
434     numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
435    
436     if (Finley_noError()) {
437     int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
438 phornby 1646 int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439 ksteube 1402 if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
440     if (mpi_info->rank == 0) { /* Master */
441     for (;;) { /* Infinite loop */
442     chunkEle = 0;
443     for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
444     for (i0=0; i0<chunkSize; i0++) {
445     if (totalEle >= numEle) break; /* End infinite loop */
446     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
447     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
448     fscanf(fileHandle_p, "\n");
449     totalEle++;
450     chunkEle++;
451     }
452     /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
453     if (nextCPU < mpi_info->size) {
454     #ifdef PASO_MPI
455 phornby 1646 int mpi_error;
456    
457 ksteube 1402 tempInts[numEle*(2+numNodes)] = chunkEle;
458     printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
459     mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
460     if ( mpi_error != MPI_SUCCESS ) {
461     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
462     return NULL;
463     }
464     #endif
465     nextCPU++;
466     }
467     if (totalEle >= numEle) break; /* End infinite loop */
468     } /* Infinite loop */
469     } /* End master */
470     else { /* Worker */
471     #ifdef PASO_MPI
472     /* Each worker receives two messages */
473     MPI_Status status;
474 ksteube 1660 int mpi_error;
475 ksteube 1402 printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
476     mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
477     if ( mpi_error != MPI_SUCCESS ) {
478     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
479     return NULL;
480     }
481     chunkEle = tempInts[numEle*(2+numNodes)];
482     #endif
483     } /* Worker */
484 ksteube 1713 #if 0
485 ksteube 1402 /* Display the temp mem for debugging */
486 ksteube 1713 printf("ksteube Elements tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
487 ksteube 1402 for (i0=0; i0<numEle*(numNodes+2); i0++) {
488     printf(" %2d", tempInts[i0]);
489     if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
490     }
491     #endif
492    
493     /* Copy Element data from tempInts to mesh_p */
494     Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
495     mesh_p->Elements->minColor=0;
496     mesh_p->Elements->maxColor=chunkEle-1;
497     if (Finley_noError()) {
498     #pragma omp parallel for private (i0, i1)
499     for (i0=0; i0<chunkEle; i0++) {
500     mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
501     mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
502     for (i1 = 0; i1 < numNodes; i1++) {
503     mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
504     }
505     }
506     }
507    
508     TMPMEMFREE(tempInts);
509     }
510    
511    
512     #if 0 /* this is the original code for reading elements */
513 ksteube 1360 /* read elements */
514     if (Finley_noError()) {
515 ksteube 1402
516 ksteube 1360 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
517     typeID=Finley_RefElement_getTypeId(element_type);
518     if (typeID==NoType) {
519     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
520     Finley_setError(VALUE_ERROR,error_msg);
521     } else {
522     /* read the elements */
523     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524     if (Finley_noError()) {
525     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
526     mesh_p->Elements->minColor=0;
527     mesh_p->Elements->maxColor=numEle-1;
528     if (Finley_noError()) {
529     for (i0 = 0; i0 < numEle; i0++) {
530     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
531     mesh_p->Elements->Color[i0]=i0;
532     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
533     fscanf(fileHandle_p, " %d",
534     &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
535     } /* for i1 */
536     fscanf(fileHandle_p, "\n");
537     } /* for i0 */
538     }
539     }
540     }
541     }
542 ksteube 1402 #endif
543    
544     printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
545    
546     #if 1
547    
548     /* Define other structures to keep mesh_write from crashing */
549 ksteube 1637 /* Change the typeid from NoType later */
550 ksteube 1402
551 ksteube 1637 mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
552 ksteube 1402 Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
553    
554 ksteube 1637 mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
555 ksteube 1402 Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
556    
557 ksteube 1637 mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
558 ksteube 1402 Finley_ElementFile_allocTable(mesh_p->Points, 0);
559    
560     #endif
561    
562     #if 0 /* comment out the rest of the un-implemented crap for now */
563 ksteube 1360 /* get the face elements */
564     if (Finley_noError()) {
565     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
566     faceTypeID=Finley_RefElement_getTypeId(element_type);
567     if (faceTypeID==NoType) {
568     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
569     Finley_setError(VALUE_ERROR,error_msg);
570     } else {
571     mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
572     if (Finley_noError()) {
573     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
574     if (Finley_noError()) {
575     mesh_p->FaceElements->minColor=0;
576     mesh_p->FaceElements->maxColor=numEle-1;
577     for (i0 = 0; i0 < numEle; i0++) {
578     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
579     mesh_p->FaceElements->Color[i0]=i0;
580     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
581     fscanf(fileHandle_p, " %d",
582     &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
583     } /* for i1 */
584     fscanf(fileHandle_p, "\n");
585     } /* for i0 */
586     }
587     }
588     }
589     }
590     /* get the Contact face element */
591     if (Finley_noError()) {
592     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
593     contactTypeID=Finley_RefElement_getTypeId(element_type);
594     if (contactTypeID==NoType) {
595     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
596     Finley_setError(VALUE_ERROR,error_msg);
597     } else {
598     mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
599     if (Finley_noError()) {
600     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
601     if (Finley_noError()) {
602     mesh_p->ContactElements->minColor=0;
603     mesh_p->ContactElements->maxColor=numEle-1;
604     for (i0 = 0; i0 < numEle; i0++) {
605     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
606     mesh_p->ContactElements->Color[i0]=i0;
607     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
608     fscanf(fileHandle_p, " %d",
609     &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
610     } /* for i1 */
611     fscanf(fileHandle_p, "\n");
612     } /* for i0 */
613     }
614     }
615     }
616 ksteube 1402 }
617 ksteube 1360 /* get the nodal element */
618     if (Finley_noError()) {
619     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
620     pointTypeID=Finley_RefElement_getTypeId(element_type);
621     if (pointTypeID==NoType) {
622     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
623     Finley_setError(VALUE_ERROR,error_msg);
624     }
625     mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
626     if (Finley_noError()) {
627     Finley_ElementFile_allocTable(mesh_p->Points, numEle);
628     if (Finley_noError()) {
629     mesh_p->Points->minColor=0;
630     mesh_p->Points->maxColor=numEle-1;
631     for (i0 = 0; i0 < numEle; i0++) {
632     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
633     mesh_p->Points->Color[i0]=i0;
634     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
635     fscanf(fileHandle_p, " %d",
636     &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
637     } /* for i1 */
638     fscanf(fileHandle_p, "\n");
639     } /* for i0 */
640     }
641     }
642     }
643     /* get the name tags */
644     if (Finley_noError()) {
645     if (feof(fileHandle_p) == 0) {
646     fscanf(fileHandle_p, "%s\n", name);
647     while (feof(fileHandle_p) == 0) {
648     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
649     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
650     }
651     }
652     }
653 ksteube 1402 #endif /* comment out the rest of the un-implemented crap for now */
654 ksteube 1360 }
655     /* close file */
656 ksteube 1402 if (mpi_info->rank == 0) fclose(fileHandle_p);
657    
658 ksteube 1360 /* resolve id's : */
659     /* rearrange elements: */
660 ksteube 1402
661 ksteube 1360 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
662     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
663 ksteube 1402 return mesh_p; /* ksteube temp return for debugging */
664    
665 ksteube 1360 /* that's it */
666     #ifdef Finley_TRACE
667     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
668     #endif
669    
670     /* that's it */
671     if (! Finley_noError()) {
672     Finley_Mesh_free(mesh_p);
673     }
674     Paso_MPIInfo_free( mpi_info );
675     return mesh_p;
676     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26