/[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 1402 - (hide annotations)
Thu Jan 31 00:17:49 2008 UTC (11 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 28671 byte(s)
now implemented to read/fan out elements

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     index_t tag_key;
240     Finley_Mesh *mesh_p=NULL;
241     char name[LenString_MAX],element_type[LenString_MAX],frm[20];
242     char error_msg[LenErrorMsg_MAX];
243     double time0=Finley_timer();
244     FILE *fileHandle_p = NULL;
245     ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
246 ksteube 1402
247 ksteube 1360 Finley_resetError();
248    
249     if (mpi_info->rank == 0) {
250     /* get file handle */
251     fileHandle_p = fopen(fname, "r");
252     if (fileHandle_p==NULL) {
253     sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
254     Finley_setError(IO_ERROR,error_msg);
255     Paso_MPIInfo_free( mpi_info );
256     return NULL;
257     }
258 ksteube 1402
259 ksteube 1360 /* read header */
260     sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
261     fscanf(fileHandle_p, frm, name);
262 ksteube 1402
263     /* get the number of nodes */
264 ksteube 1360 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
265     }
266    
267     #ifdef PASO_MPI
268     /* MPI Broadcast numDim, numNodes, name */
269     if (mpi_info->size > 0) {
270     int temp1[3], error_code;
271     temp1[0] = numDim;
272     temp1[1] = numNodes;
273     temp1[2] = strlen(name) + 1;
274     error_code = MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
275     if (error_code != MPI_SUCCESS) {
276     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
277     return NULL;
278     }
279     numDim = temp1[0];
280     numNodes = temp1[1];
281     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
282     if (error_code != MPI_SUCCESS) {
283     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
284     return NULL;
285     }
286     }
287     #endif
288    
289     /* allocate mesh */
290     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
291     if (Finley_noError()) {
292 ksteube 1402 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1, mpi_error;
293 ksteube 1360 int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
294     double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
295    
296     /*
297 ksteube 1402 Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
298     It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
299 ksteube 1360 First chunk sent to CPU 1, second to CPU 2, ...
300     Last chunk stays on CPU 0 (the master)
301     The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
302     */
303    
304     if (mpi_info->rank == 0) { /* Master */
305     for (;;) { /* Infinite loop */
306     chunkNodes = 0;
307 ksteube 1402 for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
308 ksteube 1360 for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
309     for (i1=0; i1<chunkSize; i1++) {
310     if (totalNodes >= numNodes) break;
311     if (1 == numDim)
312     fscanf(fileHandle_p, "%d %d %d %le\n",
313     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
314     &tempCoords[i1*numDim+0]);
315     if (2 == numDim)
316     fscanf(fileHandle_p, "%d %d %d %le %le\n",
317     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
318     &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
319     if (3 == numDim)
320     fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
321     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
322     &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
323     totalNodes++;
324     chunkNodes++;
325     }
326     /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
327     if (nextCPU < mpi_info->size) {
328     #ifdef PASO_MPI
329     tempInts[numNodes*3] = chunkNodes;
330 ksteube 1402 /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
331     mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
332 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
333     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
334     return NULL;
335     }
336 ksteube 1402 mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
337 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
338     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
339     return NULL;
340     }
341     #endif
342     nextCPU++;
343     }
344     if (totalNodes >= numNodes) break;
345     } /* Infinite loop */
346     } /* End master */
347     else { /* Worker */
348     #ifdef PASO_MPI
349     /* Each worker receives two messages */
350     MPI_Status status;
351 ksteube 1402 mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
352 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
353     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
354     return NULL;
355     }
356 ksteube 1402 mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
357 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
358     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
359     return NULL;
360     }
361     chunkNodes = tempInts[numNodes*3];
362     #endif
363     } /* Worker */
364 ksteube 1402
365 ksteube 1360 #if 0
366 ksteube 1402 /* Display the temp mem for debugging */
367 ksteube 1360 printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
368     for (i0=0; i0<numNodes*3; i0++) {
369     printf(" %2d", tempInts[i0]);
370     if (i0%numNodes==numNodes-1) printf("\n");
371     }
372     printf("ksteube tempCoords:\n");
373     for (i0=0; i0<chunkNodes*numDim; i0++) {
374     printf(" %20.15e", tempCoords[i0]);
375     if (i0%numDim==numDim-1) printf("\n");
376     }
377     #endif
378    
379 ksteube 1402 printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
380     /* Copy node data from tempMem to mesh_p */
381     Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
382     if (Finley_noError()) {
383     for (i0=0; i0<chunkNodes; i0++) {
384     mesh_p->Nodes->Id[i0] = tempInts[0+i0];
385     mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[numNodes+i0];
386     mesh_p->Nodes->Tag[i0] = tempInts[numNodes*2+i0];
387     for (i1=0; i1<numDim; i1++) {
388     mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
389     }
390     }
391 ksteube 1360 }
392    
393     TMPMEMFREE(tempInts);
394     TMPMEMFREE(tempCoords);
395    
396 ksteube 1402 /* read elements */
397 ksteube 1360
398 ksteube 1402 /* Read the element typeID */
399     if (Finley_noError()) {
400     if (mpi_info->rank == 0) {
401     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
402     typeID=Finley_RefElement_getTypeId(element_type);
403     }
404     #ifdef PASO_MPI
405     if (mpi_info->size > 0) {
406     int temp1[3];
407     temp1[0] = typeID;
408     temp1[1] = numEle;
409     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
410     if (mpi_error != MPI_SUCCESS) {
411     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
412     return NULL;
413     }
414     typeID = temp1[0];
415     numEle = temp1[1];
416     }
417     #endif
418     if (typeID==NoType) {
419     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
420     Finley_setError(VALUE_ERROR, error_msg);
421     }
422     }
423    
424     /* Read the element data */
425     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
426     numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
427    
428     if (Finley_noError()) {
429     int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
430     int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1, mpi_error;
431     if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
432     if (mpi_info->rank == 0) { /* Master */
433     for (;;) { /* Infinite loop */
434     chunkEle = 0;
435     for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
436     for (i0=0; i0<chunkSize; i0++) {
437     if (totalEle >= numEle) break; /* End infinite loop */
438     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
439     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
440     fscanf(fileHandle_p, "\n");
441     totalEle++;
442     chunkEle++;
443     }
444     /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
445     if (nextCPU < mpi_info->size) {
446     #ifdef PASO_MPI
447     tempInts[numEle*(2+numNodes)] = chunkEle;
448     printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
449     mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
450     if ( mpi_error != MPI_SUCCESS ) {
451     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
452     return NULL;
453     }
454     #endif
455     nextCPU++;
456     }
457     if (totalEle >= numEle) break; /* End infinite loop */
458     } /* Infinite loop */
459     } /* End master */
460     else { /* Worker */
461     #ifdef PASO_MPI
462     /* Each worker receives two messages */
463     MPI_Status status;
464     printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
465     mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
466     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     chunkEle = tempInts[numEle*(2+numNodes)];
471     #endif
472     } /* Worker */
473     #if 1
474     /* Display the temp mem for debugging */
475     printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
476     for (i0=0; i0<numEle*(numNodes+2); i0++) {
477     printf(" %2d", tempInts[i0]);
478     if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
479     }
480     #endif
481    
482     /* Copy Element data from tempInts to mesh_p */
483     Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
484     mesh_p->Elements->minColor=0;
485     mesh_p->Elements->maxColor=chunkEle-1;
486     if (Finley_noError()) {
487     #pragma omp parallel for private (i0, i1)
488     for (i0=0; i0<chunkEle; i0++) {
489     mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
490     mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
491     for (i1 = 0; i1 < numNodes; i1++) {
492     mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
493     }
494     }
495     }
496    
497     TMPMEMFREE(tempInts);
498     }
499    
500    
501     #if 0 /* this is the original code for reading elements */
502 ksteube 1360 /* read elements */
503     if (Finley_noError()) {
504 ksteube 1402
505 ksteube 1360 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
506     typeID=Finley_RefElement_getTypeId(element_type);
507     if (typeID==NoType) {
508     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
509     Finley_setError(VALUE_ERROR,error_msg);
510     } else {
511     /* read the elements */
512     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
513     if (Finley_noError()) {
514     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
515     mesh_p->Elements->minColor=0;
516     mesh_p->Elements->maxColor=numEle-1;
517     if (Finley_noError()) {
518     for (i0 = 0; i0 < numEle; i0++) {
519     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
520     mesh_p->Elements->Color[i0]=i0;
521     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
522     fscanf(fileHandle_p, " %d",
523     &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
524     } /* for i1 */
525     fscanf(fileHandle_p, "\n");
526     } /* for i0 */
527     }
528     }
529     }
530     }
531 ksteube 1402 #endif
532    
533     printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
534    
535     #if 1
536    
537     /* Define other structures to keep mesh_write from crashing */
538    
539     mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
540     Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
541    
542     mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
543     Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
544    
545     mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
546     Finley_ElementFile_allocTable(mesh_p->Points, 0);
547    
548     #endif
549    
550     #if 0 /* comment out the rest of the un-implemented crap for now */
551 ksteube 1360 /* get the face elements */
552     if (Finley_noError()) {
553     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
554     faceTypeID=Finley_RefElement_getTypeId(element_type);
555     if (faceTypeID==NoType) {
556     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
557     Finley_setError(VALUE_ERROR,error_msg);
558     } else {
559     mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
560     if (Finley_noError()) {
561     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
562     if (Finley_noError()) {
563     mesh_p->FaceElements->minColor=0;
564     mesh_p->FaceElements->maxColor=numEle-1;
565     for (i0 = 0; i0 < numEle; i0++) {
566     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
567     mesh_p->FaceElements->Color[i0]=i0;
568     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
569     fscanf(fileHandle_p, " %d",
570     &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
571     } /* for i1 */
572     fscanf(fileHandle_p, "\n");
573     } /* for i0 */
574     }
575     }
576     }
577     }
578     /* get the Contact face element */
579     if (Finley_noError()) {
580     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
581     contactTypeID=Finley_RefElement_getTypeId(element_type);
582     if (contactTypeID==NoType) {
583     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
584     Finley_setError(VALUE_ERROR,error_msg);
585     } else {
586     mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
587     if (Finley_noError()) {
588     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
589     if (Finley_noError()) {
590     mesh_p->ContactElements->minColor=0;
591     mesh_p->ContactElements->maxColor=numEle-1;
592     for (i0 = 0; i0 < numEle; i0++) {
593     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
594     mesh_p->ContactElements->Color[i0]=i0;
595     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
596     fscanf(fileHandle_p, " %d",
597     &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
598     } /* for i1 */
599     fscanf(fileHandle_p, "\n");
600     } /* for i0 */
601     }
602     }
603     }
604 ksteube 1402 }
605 ksteube 1360 /* get the nodal element */
606     if (Finley_noError()) {
607     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
608     pointTypeID=Finley_RefElement_getTypeId(element_type);
609     if (pointTypeID==NoType) {
610     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
611     Finley_setError(VALUE_ERROR,error_msg);
612     }
613     mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
614     if (Finley_noError()) {
615     Finley_ElementFile_allocTable(mesh_p->Points, numEle);
616     if (Finley_noError()) {
617     mesh_p->Points->minColor=0;
618     mesh_p->Points->maxColor=numEle-1;
619     for (i0 = 0; i0 < numEle; i0++) {
620     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
621     mesh_p->Points->Color[i0]=i0;
622     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
623     fscanf(fileHandle_p, " %d",
624     &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
625     } /* for i1 */
626     fscanf(fileHandle_p, "\n");
627     } /* for i0 */
628     }
629     }
630     }
631     /* get the name tags */
632     if (Finley_noError()) {
633     if (feof(fileHandle_p) == 0) {
634     fscanf(fileHandle_p, "%s\n", name);
635     while (feof(fileHandle_p) == 0) {
636     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
637     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
638     }
639     }
640     }
641 ksteube 1402 #endif /* comment out the rest of the un-implemented crap for now */
642 ksteube 1360 }
643     /* close file */
644 ksteube 1402 if (mpi_info->rank == 0) fclose(fileHandle_p);
645    
646 ksteube 1360 /* resolve id's : */
647     /* rearrange elements: */
648 ksteube 1402
649 ksteube 1360 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
650     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
651 ksteube 1402 return mesh_p; /* ksteube temp return for debugging */
652    
653 ksteube 1360 /* that's it */
654     #ifdef Finley_TRACE
655     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
656     #endif
657    
658     /* that's it */
659     if (! Finley_noError()) {
660     Finley_Mesh_free(mesh_p);
661     }
662     Paso_MPIInfo_free( mpi_info );
663     return mesh_p;
664     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26