/[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 1713 by ksteube, Wed Aug 20 07:03:45 2008 UTC revision 1725 by ksteube, Tue Aug 26 03:31:49 2008 UTC
# Line 241  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 241  Finley_Mesh* Finley_Mesh_read_MPI(char*
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;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
   ElementTypeId faceTypeID, contactTypeID, pointTypeID;  
245    index_t tag_key;    index_t tag_key;
246    
247    Finley_resetError();    Finley_resetError();
# Line 290  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 289  Finley_Mesh* Finley_Mesh_read_MPI(char*
289       /* allocate mesh */       /* allocate mesh */
290       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
291       if (Finley_noError()) {       if (Finley_noError()) {
292      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large for that */      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
293      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
294      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
295      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
# Line 369  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 368  Finley_Mesh* Finley_Mesh_read_MPI(char*
368  #endif  #endif
369      }   /* Worker */      }   /* Worker */
370    
371  #if 1  #if 0
372          /* Display the temp mem for debugging */          /* Display the temp mem for debugging */
373          printf("ksteube Nodes tempInts\n");          printf("ksteube Nodes tempInts\n");
374          for (i0=0; i0<chunkSize*3; i0++) {          for (i0=0; i0<chunkSize*3; i0++) {
# Line 381  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 380  Finley_Mesh* Finley_Mesh_read_MPI(char*
380            printf(" %20.15e", tempCoords[i0]);            printf(" %20.15e", tempCoords[i0]);
381            if (i0%numDim==numDim-1) printf("\n");            if (i0%numDim==numDim-1) printf("\n");
382          }          }
383                printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);
384  #endif  #endif
385    
         printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);  
   
386      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
387          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
388          if (Finley_noError()) {          if (Finley_noError()) {
# Line 411  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 409  Finley_Mesh* Finley_Mesh_read_MPI(char*
409        }        }
410  #ifdef PASO_MPI  #ifdef PASO_MPI
411        if (mpi_info->size > 1) {        if (mpi_info->size > 1) {
412          int temp1[3], mpi_error;          int temp1[2], mpi_error;
413          temp1[0] = (int) typeID;          temp1[0] = (int) typeID;
414          temp1[1] = numEle;          temp1[1] = numEle;
415          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
# Line 434  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 432  Finley_Mesh* Finley_Mesh_read_MPI(char*
432        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
433    
434        if (Finley_noError()) {        if (Finley_noError()) {
435      int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
436      int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
437      if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
438      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
439        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
440            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
441          chunkEle = 0;          chunkEle = 0;
         for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
442          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
443            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End infinite loop */
444            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
# Line 449  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 447  Finley_Mesh* Finley_Mesh_read_MPI(char*
447            totalEle++;            totalEle++;
448            chunkEle++;            chunkEle++;
449          }          }
         /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
450  #ifdef PASO_MPI  #ifdef PASO_MPI
451            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
452            if (nextCPU < mpi_info->size) {
453                int mpi_error;                int mpi_error;
454              tempInts[chunkSize*(2+numNodes)] = chunkEle;
455            tempInts[numEle*(2+numNodes)] = chunkEle;            mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);  
           mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);  
456            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
457                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
458                  return NULL;                  return NULL;
459            }            }
 #endif  
460            nextCPU++;            nextCPU++;
461          }          }
462    #endif
463          if (totalEle >= numEle) break; /* End infinite loop */          if (totalEle >= numEle) break; /* End infinite loop */
464        } /* Infinite loop */        } /* Infinite loop */
465      }   /* End master */      }   /* End master */
466      else {  /* Worker */      else {  /* Worker */
467  #ifdef PASO_MPI  #ifdef PASO_MPI
468        /* Each worker receives two messages */        /* Each worker receives one message */
469        MPI_Status status;        MPI_Status status;
470        int mpi_error;        int mpi_error;
471  printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
       mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);  
472        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
473              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
474              return NULL;              return NULL;
475        }        }
476        chunkEle = tempInts[numEle*(2+numNodes)];        chunkEle = tempInts[chunkSize*(2+numNodes)];
477  #endif  #endif
478      }   /* Worker */      }   /* Worker */
479    
480  #if 0  #if 0
481      /* Display the temp mem for debugging */      /* Display the temp mem for debugging */
482      printf("ksteube Elements tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);      for (i0=0; i0<chunkSize*(numNodes+2); i0++) {
     for (i0=0; i0<numEle*(numNodes+2); i0++) {  
483        printf(" %2d", tempInts[i0]);        printf(" %2d", tempInts[i0]);
484        if (i0%(numNodes+2)==numNodes+2-1) printf("\n");        if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
485      }      }
# Line 499  printf("ksteube CPU=%d/%d recv on %d\n", Line 494  printf("ksteube CPU=%d/%d recv on %d\n",
494        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
495          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
496          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
497                mesh_p->Elements->Color[i0] = i0;
498          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
499            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
500          }          }
# Line 507  printf("ksteube CPU=%d/%d recv on %d\n", Line 503  printf("ksteube CPU=%d/%d recv on %d\n",
503    
504      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
505        }        }
   
   
 #if 0 /* this is the original code for reading elements */  
         /* read elements */  
         if (Finley_noError()) {  
   
            fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
            typeID=Finley_RefElement_getTypeId(element_type);  
            if (typeID==NoType) {  
              sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
              Finley_setError(VALUE_ERROR,error_msg);  
            } else {  
              /* read the elements */  
              mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
              if (Finley_noError()) {  
                  Finley_ElementFile_allocTable(mesh_p->Elements, numEle);  
                  mesh_p->Elements->minColor=0;  
                  mesh_p->Elements->maxColor=numEle-1;  
                  if (Finley_noError()) {  
                     for (i0 = 0; i0 < numEle; i0++) {  
                       fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);  
                       mesh_p->Elements->Color[i0]=i0;  
                       for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {  
                            fscanf(fileHandle_p, " %d",  
                               &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
                       } /* for i1 */  
                       fscanf(fileHandle_p, "\n");  
                     } /* for i0 */  
                  }  
              }  
           }  
         }  
 #endif  
   
 printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);  
   
 #if 1  
   
   /* Define other structures to keep mesh_write from crashing */  
   /* Change the typeid from NoType later */  
   
   mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);  
   Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);  
   
   mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);  
   Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);  
   
   mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);  
   Finley_ElementFile_allocTable(mesh_p->Points, 0);  
   
 #endif  
   
 #if 0 /* comment out the rest of the un-implemented crap for now */  
         /* get the face elements */  
         if (Finley_noError()) {  
              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
              faceTypeID=Finley_RefElement_getTypeId(element_type);  
              if (faceTypeID==NoType) {  
                sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);  
                Finley_setError(VALUE_ERROR,error_msg);  
              } else {  
                 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
                 if (Finley_noError()) {  
                    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
                    if (Finley_noError()) {  
                       mesh_p->FaceElements->minColor=0;  
                       mesh_p->FaceElements->maxColor=numEle-1;  
                       for (i0 = 0; i0 < numEle; i0++) {  
                         fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);  
                         mesh_p->FaceElements->Color[i0]=i0;  
                         for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {  
                              fscanf(fileHandle_p, " %d",  
                                 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);  
                         }   /* for i1 */  
                         fscanf(fileHandle_p, "\n");  
                       } /* for i0 */  
                    }  
                 }  
              }  
         }  
         /* get the Contact face element */  
         if (Finley_noError()) {  
              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
              contactTypeID=Finley_RefElement_getTypeId(element_type);  
              if (contactTypeID==NoType) {  
                sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);  
                Finley_setError(VALUE_ERROR,error_msg);  
              } else {  
                mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
                if (Finley_noError()) {  
                    Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);  
                    if (Finley_noError()) {  
                       mesh_p->ContactElements->minColor=0;  
                       mesh_p->ContactElements->maxColor=numEle-1;  
                       for (i0 = 0; i0 < numEle; i0++) {  
                         fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);  
                         mesh_p->ContactElements->Color[i0]=i0;  
                         for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {  
                             fscanf(fileHandle_p, " %d",  
                                &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);  
                         }   /* for i1 */  
                         fscanf(fileHandle_p, "\n");  
                       } /* for i0 */  
                   }  
                }  
              }  
         }  
         /* get the nodal element */  
         if (Finley_noError()) {  
              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
              pointTypeID=Finley_RefElement_getTypeId(element_type);  
              if (pointTypeID==NoType) {  
                sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);  
                Finley_setError(VALUE_ERROR,error_msg);  
              }  
              mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
              if (Finley_noError()) {  
                 Finley_ElementFile_allocTable(mesh_p->Points, numEle);  
                 if (Finley_noError()) {  
                    mesh_p->Points->minColor=0;  
                    mesh_p->Points->maxColor=numEle-1;  
                    for (i0 = 0; i0 < numEle; i0++) {  
                      fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
                      mesh_p->Points->Color[i0]=i0;  
                      for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
                          fscanf(fileHandle_p, " %d",  
                             &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
                      }  /* for i1 */  
                      fscanf(fileHandle_p, "\n");  
                    } /* for i0 */  
                 }  
              }  
         }  
         /* get the name tags */  
         if (Finley_noError()) {  
            if (feof(fileHandle_p) == 0) {  
               fscanf(fileHandle_p, "%s\n", name);  
               while (feof(fileHandle_p) == 0) {  
                    fscanf(fileHandle_p, "%s %d\n", name, &tag_key);  
                    Finley_Mesh_addTagMap(mesh_p,name,tag_key);  
               }  
            }  
         }  
 #endif /* comment out the rest of the un-implemented crap for now */  
506       }       }
507    
508       /* close file */       /* close file */
509       if (mpi_info->rank == 0) fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
510    
511       /*   resolve id's : */       /*   resolve id's : */
512       /* rearrange elements: */       /* rearrange elements: */
513    
514         return mesh_p; /* ksteube temp return for debugging */
515    
516       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
517       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
      return mesh_p; /* ksteube temp return for debugging */  
518    
519       /* that's it */       /* that's it */
520       #ifdef Finley_TRACE       #ifdef Finley_TRACE

Legend:
Removed from v.1713  
changed lines
  Added in v.1725

  ViewVC Help
Powered by ViewVC 1.1.26