/[escript]/trunk/finley/src/Mesh_read.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC
# Line 230  Finley_Mesh* Finley_Mesh_read(char* fnam Line 230  Finley_Mesh* Finley_Mesh_read(char* fnam
230    return mesh_p;    return mesh_p;
231  }  }
232    
233  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
234    
235  {  {
236    
237    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
238    dim_t numNodes, numDim, numEle, i0, i1;    dim_t numNodes, numDim, numEle, i0, i1;
   index_t tag_key;  
239    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
240    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
241    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
242    double time0=Finley_timer();    double time0=Finley_timer();
243    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
244    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID;
245      
246      /* these are in unimplemented code below */
247    #if 0
248      index_t tag_key;
249      ElementTypeId faceTypeID, contactTypeID, pointTypeID;
250    #endif
251    
252    Finley_resetError();    Finley_resetError();
253    
254    if (mpi_info->rank == 0) {    if (mpi_info->rank == 0) {
# Line 255  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 260  Finley_Mesh* Finley_Mesh_read_MPI(char*
260         Paso_MPIInfo_free( mpi_info );         Paso_MPIInfo_free( mpi_info );
261         return NULL;         return NULL;
262       }       }
263      
264       /* read header */       /* read header */
265       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
266       fscanf(fileHandle_p, frm, name);       fscanf(fileHandle_p, frm, name);
267      
268       /* get the nodes */       /* get the number of nodes */
     
269       fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);       fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
270    }    }
271    
# Line 286  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 290  Finley_Mesh* Finley_Mesh_read_MPI(char*
290      }      }
291    }    }
292  #endif  #endif
   printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes);  
293    
294       /* allocate mesh */       /* allocate mesh */
295       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
296       if (Finley_noError()) {       if (Finley_noError()) {
297      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, nextCPU=1, mpi_error;      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
298    #ifdef PASO_MPI
299        int mpi_error;
300    #endif
301      int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);      int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
302      double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);      double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
303    
304      /*      /*
305        Read a chunk of nodes, send to worker CPU if necessary, copy chunk into local mesh_p        Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
306        It doesn't matter that a CPU has the wrong nodes, this is sorted out later        It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
307        First chunk sent to CPU 1, second to CPU 2, ...        First chunk sent to CPU 1, second to CPU 2, ...
308        Last chunk stays on CPU 0 (the master)        Last chunk stays on CPU 0 (the master)
309        The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message        The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
# Line 306  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 312  Finley_Mesh* Finley_Mesh_read_MPI(char*
312      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
313        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
314          chunkNodes = 0;          chunkNodes = 0;
315          for (i0=0; i0<numNodes*3; i0++) tempInts[i0] = -1;          for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
316          for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;          for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
317          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
318            if (totalNodes >= numNodes) break;            if (totalNodes >= numNodes) break;
# Line 329  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 335  Finley_Mesh* Finley_Mesh_read_MPI(char*
335          if (nextCPU < mpi_info->size) {          if (nextCPU < mpi_info->size) {
336  #ifdef PASO_MPI  #ifdef PASO_MPI
337            tempInts[numNodes*3] = chunkNodes;            tempInts[numNodes*3] = chunkNodes;
338            mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 8727, mpi_info->comm);            /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
339              mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
340            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
341                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
342                  return NULL;                  return NULL;
343            }            }
344            mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 8728, mpi_info->comm);            mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
345            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
346                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
347                  return NULL;                  return NULL;
# Line 349  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 356  Finley_Mesh* Finley_Mesh_read_MPI(char*
356  #ifdef PASO_MPI  #ifdef PASO_MPI
357        /* Each worker receives two messages */        /* Each worker receives two messages */
358        MPI_Status status;        MPI_Status status;
359        mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 8727, mpi_info->comm, &status);        mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
360        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
361              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
362              return NULL;              return NULL;
363        }        }
364        mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 8728, mpi_info->comm, &status);        mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
365        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
366              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
367              return NULL;              return NULL;
# Line 362  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 369  Finley_Mesh* Finley_Mesh_read_MPI(char*
369        chunkNodes = tempInts[numNodes*3];        chunkNodes = tempInts[numNodes*3];
370  #endif  #endif
371      }   /* Worker */      }   /* Worker */
372  printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d chunkNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes, chunkNodes);  
373  #if 0  #if 0
374            /* Display the temp mem for debugging */
375          printf("ksteube tempInts totalNodes=%d:\n", totalNodes);          printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
376          for (i0=0; i0<numNodes*3; i0++) {          for (i0=0; i0<numNodes*3; i0++) {
377            printf(" %2d", tempInts[i0]);            printf(" %2d", tempInts[i0]);
# Line 376  printf("ksteube CPU=%d/%d name='%s' numD Line 384  printf("ksteube CPU=%d/%d name='%s' numD
384          }          }
385  #endif  #endif
386    
387      /* Copy tempMem into mesh_p */      printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
388          Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);      /* Copy node data from tempMem to mesh_p */
389      for (i0=0; i0<numNodes; i0++) {          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
390        mesh_p->Nodes->Id[i0]                 = tempInts[0+i0];          if (Finley_noError()) {
391        mesh_p->Nodes->globalDegreesOfFreedom[i0]     = tempInts[numNodes+i0];        for (i0=0; i0<chunkNodes; i0++) {
392        mesh_p->Nodes->Tag[i0]                = tempInts[numNodes*2+i0];          mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
393        for (i1=0; i1<numDim; i1++) {          mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[numNodes+i0];
394          mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]    = tempCoords[i0*numDim+i1];          mesh_p->Nodes->Tag[i0]              = tempInts[numNodes*2+i0];
395        }          for (i1=0; i1<numDim; i1++) {
396              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
397            }
398              }
399          }          }
400    
401      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
402      TMPMEMFREE(tempCoords);      TMPMEMFREE(tempCoords);
403    
404          return mesh_p; /* ksteube temp for debugging */          /* read elements */
405    
406        /* 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          if (mpi_info->size > 0) {
414            int temp1[3];
415            temp1[0] = typeID;
416            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            typeID = temp1[0];
423            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        int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439    #ifdef PASO_MPI
440        int mpi_error;
441    #endif
442        if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
443        if (mpi_info->rank == 0) {  /* Master */
444          for (;;) {            /* Infinite loop */
445            chunkEle = 0;
446            for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
447            for (i0=0; i0<chunkSize; i0++) {
448              if (totalEle >= numEle) break; /* End infinite loop */
449              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
450              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
451              fscanf(fileHandle_p, "\n");
452              totalEle++;
453              chunkEle++;
454            }
455            /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
456            if (nextCPU < mpi_info->size) {
457    #ifdef PASO_MPI
458              tempInts[numEle*(2+numNodes)] = chunkEle;
459    printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
460              mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
461              if ( mpi_error != MPI_SUCCESS ) {
462                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
463                    return NULL;
464              }
465    #endif
466              nextCPU++;
467            }
468            if (totalEle >= numEle) break; /* End infinite loop */
469          } /* Infinite loop */
470        }   /* End master */
471        else {  /* Worker */
472    #ifdef PASO_MPI
473          /* Each worker receives two messages */
474          MPI_Status status;
475    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    #if 1
485        /* Display the temp mem for debugging */
486        printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
487        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          /* read elements */          /* read elements */
514          if (Finley_noError()) {          if (Finley_noError()) {
515      
516             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
517             typeID=Finley_RefElement_getTypeId(element_type);             typeID=Finley_RefElement_getTypeId(element_type);
518             if (typeID==NoType) {             if (typeID==NoType) {
# Line 421  printf("ksteube CPU=%d/%d name='%s' numD Line 539  printf("ksteube CPU=%d/%d name='%s' numD
539               }               }
540            }            }
541          }          }
542    #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    
550      mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
551      Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
552    
553      mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
554      Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
555    
556      mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
557      Finley_ElementFile_allocTable(mesh_p->Points, 0);
558    
559    #endif
560    
561    #if 0 /* comment out the rest of the un-implemented crap for now */
562          /* get the face elements */          /* get the face elements */
563          if (Finley_noError()) {          if (Finley_noError()) {
564               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
# Line 474  printf("ksteube CPU=%d/%d name='%s' numD Line 612  printf("ksteube CPU=%d/%d name='%s' numD
612                    }                    }
613                 }                 }
614               }               }
615          }            }
616          /* get the nodal element */          /* get the nodal element */
617          if (Finley_noError()) {          if (Finley_noError()) {
618               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
# Line 511  printf("ksteube CPU=%d/%d name='%s' numD Line 649  printf("ksteube CPU=%d/%d name='%s' numD
649                }                }
650             }             }
651          }          }
652    #endif /* comment out the rest of the un-implemented crap for now */
653       }       }
654       /* close file */       /* close file */
655       fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
656      
657       /*   resolve id's : */       /*   resolve id's : */
658       /* rearrange elements: */       /* rearrange elements: */
659      
660       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
661       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
662           return mesh_p; /* ksteube temp return for debugging */
663    
664       /* that's it */       /* that's it */
665       #ifdef Finley_TRACE       #ifdef Finley_TRACE
666       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);

Legend:
Removed from v.1388  
changed lines
  Added in v.1628

  ViewVC Help
Powered by ViewVC 1.1.26