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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26