/[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 1771 by ksteube, Mon Sep 8 22:47:55 2008 UTC
# Line 19  Line 19 
19    
20  /**************************************************************/  /**************************************************************/
21    
22    #include <ctype.h>
23  #include "Mesh.h"  #include "Mesh.h"
24    
25  /**************************************************************/  /**************************************************************/
# Line 106  Finley_Mesh* Finley_Mesh_read(char* fnam Line 107  Finley_Mesh* Finley_Mesh_read(char* fnam
107                      for (i0 = 0; i0 < numEle; i0++) {                      for (i0 = 0; i0 < numEle; i0++) {
108                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
109                        mesh_p->Elements->Color[i0]=i0;                        mesh_p->Elements->Color[i0]=i0;
110                          mesh_p->Elements->Owner[i0]=0;
111                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
112                             fscanf(fileHandle_p, " %d",                             fscanf(fileHandle_p, " %d",
113                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
# Line 133  Finley_Mesh* Finley_Mesh_read(char* fnam Line 135  Finley_Mesh* Finley_Mesh_read(char* fnam
135                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
136                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
137                          mesh_p->FaceElements->Color[i0]=i0;                          mesh_p->FaceElements->Color[i0]=i0;
138                            mesh_p->FaceElements->Owner[i0]=0;
139                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
140                               fscanf(fileHandle_p, " %d",                               fscanf(fileHandle_p, " %d",
141                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
# Line 160  Finley_Mesh* Finley_Mesh_read(char* fnam Line 163  Finley_Mesh* Finley_Mesh_read(char* fnam
163                        for (i0 = 0; i0 < numEle; i0++) {                        for (i0 = 0; i0 < numEle; i0++) {
164                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
165                          mesh_p->ContactElements->Color[i0]=i0;                          mesh_p->ContactElements->Color[i0]=i0;
166                            mesh_p->ContactElements->Owner[i0]=0;
167                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
168                              fscanf(fileHandle_p, " %d",                              fscanf(fileHandle_p, " %d",
169                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
# Line 187  Finley_Mesh* Finley_Mesh_read(char* fnam Line 191  Finley_Mesh* Finley_Mesh_read(char* fnam
191                     for (i0 = 0; i0 < numEle; i0++) {                     for (i0 = 0; i0 < numEle; i0++) {
192                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
193                       mesh_p->Points->Color[i0]=i0;                       mesh_p->Points->Color[i0]=i0;
194                         mesh_p->Points->Owner[i0]=0;
195                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
196                           fscanf(fileHandle_p, " %d",                           fscanf(fileHandle_p, " %d",
197                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
# Line 241  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 246  Finley_Mesh* Finley_Mesh_read_MPI(char*
246    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
247    double time0=Finley_timer();    double time0=Finley_timer();
248    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
249    ElementTypeId typeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
250    ElementTypeId faceTypeID, contactTypeID, pointTypeID;    Finley_TagMap* tag_map;
251    index_t tag_key;    index_t tag_key;
252    
253    Finley_resetError();    Finley_resetError();
# Line 290  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 295  Finley_Mesh* Finley_Mesh_read_MPI(char*
295       /* allocate mesh */       /* allocate mesh */
296       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
297       if (Finley_noError()) {       if (Finley_noError()) {
298      /* 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 */
299      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;
300      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 */
301      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
# Line 309  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 314  Finley_Mesh* Finley_Mesh_read_MPI(char*
314          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
315          chunkNodes = 0;          chunkNodes = 0;
316          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
317            if (totalNodes >= numNodes) break;    /* Maybe end the infinite loop */            if (totalNodes >= numNodes) break;    /* End of inner loop */
318                if (1 == numDim)                if (1 == numDim)
319          fscanf(fileHandle_p, "%d %d %d %le\n",          fscanf(fileHandle_p, "%d %d %d %le\n",
320            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
# Line 344  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 349  Finley_Mesh* Finley_Mesh_read_MPI(char*
349                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
350                  return NULL;                  return NULL;
351            }            }
           nextCPU++;  
352          }          }
353  #endif  #endif
354          if (totalNodes >= numNodes) break;  /* Maybe end the infinite loop */          nextCPU++;
355            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
356            if (nextCPU > mpi_info->size) break; /* End infinite loop */
357        } /* Infinite loop */        } /* Infinite loop */
358      }   /* End master */      }   /* End master */
359      else {  /* Worker */      else {  /* Worker */
# Line 369  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 375  Finley_Mesh* Finley_Mesh_read_MPI(char*
375  #endif  #endif
376      }   /* Worker */      }   /* Worker */
377    
 #if 1  
         /* Display the temp mem for debugging */  
         printf("ksteube Nodes tempInts\n");  
         for (i0=0; i0<chunkSize*3; i0++) {  
           printf(" %2d", tempInts[i0]);  
           if (i0%chunkSize==chunkSize-1) printf("\n");  
         }  
         printf("ksteube tempCoords:\n");  
         for (i0=0; i0<chunkSize*numDim; i0++) {  
           printf(" %20.15e", tempCoords[i0]);  
           if (i0%numDim==numDim-1) printf("\n");  
         }  
 #endif  
   
         printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);  
   
378      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
379          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
380          if (Finley_noError()) {          if (Finley_noError()) {
# Line 411  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 401  Finley_Mesh* Finley_Mesh_read_MPI(char*
401        }        }
402  #ifdef PASO_MPI  #ifdef PASO_MPI
403        if (mpi_info->size > 1) {        if (mpi_info->size > 1) {
404          int temp1[3], mpi_error;          int temp1[2], mpi_error;
405          temp1[0] = (int) typeID;          temp1[0] = (int) typeID;
406          temp1[1] = numEle;          temp1[1] = numEle;
407          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 429  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 419  Finley_Mesh* Finley_Mesh_read_MPI(char*
419            }            }
420      }      }
421    
422        /* Read the element data */        /* Allocate the ElementFile */
423        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
424        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 */
425    
426          /* Read the element data */
427        if (Finley_noError()) {        if (Finley_noError()) {
428      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;
429      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) */
430      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 */
431      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
432        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
433            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
434          chunkEle = 0;          chunkEle = 0;
         for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
435          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
436            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End inner loop */
437            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]);
438            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
439            fscanf(fileHandle_p, "\n");            fscanf(fileHandle_p, "\n");
440            totalEle++;            totalEle++;
441            chunkEle++;            chunkEle++;
442          }          }
         /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
443  #ifdef PASO_MPI  #ifdef PASO_MPI
444            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
445            if (nextCPU < mpi_info->size) {
446                int mpi_error;                int mpi_error;
447              tempInts[chunkSize*(2+numNodes)] = chunkEle;
448            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);  
449            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
450                  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");
451                  return NULL;                  return NULL;
452            }            }
 #endif  
           nextCPU++;  
453          }          }
454          if (totalEle >= numEle) break; /* End infinite loop */  #endif
455            nextCPU++;
456            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
457            if (nextCPU > mpi_info->size) break; /* End infinite loop */
458        } /* Infinite loop */        } /* Infinite loop */
459      }   /* End master */      }   /* End master */
460      else {  /* Worker */      else {  /* Worker */
461  #ifdef PASO_MPI  #ifdef PASO_MPI
462        /* Each worker receives two messages */        /* Each worker receives one message */
463        MPI_Status status;        MPI_Status status;
464        int mpi_error;        int mpi_error;
465  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);  
466        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
467              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");
468              return NULL;              return NULL;
469        }        }
470        chunkEle = tempInts[numEle*(2+numNodes)];        chunkEle = tempInts[chunkSize*(2+numNodes)];
471  #endif  #endif
472      }   /* Worker */      }   /* Worker */
 #if 0  
     /* Display the temp mem for debugging */  
     printf("ksteube Elements tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);  
     for (i0=0; i0<numEle*(numNodes+2); i0++) {  
       printf(" %2d", tempInts[i0]);  
       if (i0%(numNodes+2)==numNodes+2-1) printf("\n");  
     }  
 #endif  
473    
474      /* Copy Element data from tempInts to mesh_p */      /* Copy Element data from tempInts to mesh_p */
475      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
# Line 499  printf("ksteube CPU=%d/%d recv on %d\n", Line 480  printf("ksteube CPU=%d/%d recv on %d\n",
480        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
481          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
482          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
483                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
484                mesh_p->Elements->Color[i0] = i0;
485          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
486            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];
487          }          }
# Line 506  printf("ksteube CPU=%d/%d recv on %d\n", Line 489  printf("ksteube CPU=%d/%d recv on %d\n",
489          }          }
490    
491      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
492        }        } /* end of Read the element data */
493    
494            /* read face elements */
495    
496  #if 0 /* this is the original code for reading elements */      /* Read the element typeID */
         /* read elements */  
497          if (Finley_noError()) {          if (Finley_noError()) {
498          if (mpi_info->rank == 0) {
499                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
500                typeID=Finley_RefElement_getTypeId(element_type);
501          }
502    #ifdef PASO_MPI
503          if (mpi_info->size > 1) {
504            int temp1[2], mpi_error;
505            temp1[0] = (int) typeID;
506            temp1[1] = numEle;
507            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
508            if (mpi_error != MPI_SUCCESS) {
509              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
510              return NULL;
511            }
512            typeID = (ElementTypeId) temp1[0];
513            numEle = temp1[1];
514          }
515    #endif
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              }
520        }
521    
522             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        /* Allocate the ElementFile */
523             typeID=Finley_RefElement_getTypeId(element_type);        mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524             if (typeID==NoType) {        numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
525               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
526               Finley_setError(VALUE_ERROR,error_msg);        /* Read the face element data */
527             } else {        if (Finley_noError()) {
528               /* read the elements */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
529               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
530               if (Finley_noError()) {      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
531                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);      if (mpi_info->rank == 0) {  /* Master */
532                   mesh_p->Elements->minColor=0;        for (;;) {            /* Infinite loop */
533                   mesh_p->Elements->maxColor=numEle-1;          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
534                   if (Finley_noError()) {          chunkEle = 0;
535                      for (i0 = 0; i0 < numEle; i0++) {          for (i0=0; i0<chunkSize; i0++) {
536                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);            if (totalEle >= numEle) break; /* End inner loop */
537                        mesh_p->Elements->Color[i0]=i0;            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
538                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
539                             fscanf(fileHandle_p, " %d",            fscanf(fileHandle_p, "\n");
540                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);            totalEle++;
541                        } /* for i1 */            chunkEle++;
542                        fscanf(fileHandle_p, "\n");          }
543                      } /* for i0 */  #ifdef PASO_MPI
544                   }          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
545               }          if (nextCPU < mpi_info->size) {
546                  int mpi_error;
547              tempInts[chunkSize*(2+numNodes)] = chunkEle;
548              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
549              if ( mpi_error != MPI_SUCCESS ) {
550                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
551                    return NULL;
552              }
553            }
554    #endif
555            nextCPU++;
556            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
557            if (nextCPU > mpi_info->size) break; /* End infinite loop */
558          } /* Infinite loop */
559        }   /* End master */
560        else {  /* Worker */
561    #ifdef PASO_MPI
562          /* Each worker receives one message */
563          MPI_Status status;
564          int mpi_error;
565          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
566          if ( mpi_error != MPI_SUCCESS ) {
567                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
568                return NULL;
569          }
570          chunkEle = tempInts[chunkSize*(2+numNodes)];
571    #endif
572        }   /* Worker */
573    
574        /* Copy Element data from tempInts to mesh_p */
575        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
576            mesh_p->FaceElements->minColor=0;
577            mesh_p->FaceElements->maxColor=chunkEle-1;
578            if (Finley_noError()) {
579              #pragma omp parallel for private (i0, i1)
580          for (i0=0; i0<chunkEle; i0++) {
581            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
582            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
583                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
584                mesh_p->FaceElements->Color[i0] = i0;
585            for (i1 = 0; i1 < numNodes; i1++) {
586              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
587            }
588            }            }
589          }          }
590    
591        TMPMEMFREE(tempInts);
592          } /* end of Read the face element data */
593    
594            /* read contact elements */
595    
596        /* Read the element typeID */
597            if (Finley_noError()) {
598          if (mpi_info->rank == 0) {
599                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
600                typeID=Finley_RefElement_getTypeId(element_type);
601          }
602    #ifdef PASO_MPI
603          if (mpi_info->size > 1) {
604            int temp1[2], mpi_error;
605            temp1[0] = (int) typeID;
606            temp1[1] = numEle;
607            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
608            if (mpi_error != MPI_SUCCESS) {
609              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
610              return NULL;
611            }
612            typeID = (ElementTypeId) temp1[0];
613            numEle = temp1[1];
614          }
615  #endif  #endif
616              if (typeID==NoType) {
617                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
618                Finley_setError(VALUE_ERROR, error_msg);
619              }
620        }
621    
622          /* Allocate the ElementFile */
623          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
624          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
625    
626  printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);        /* Read the contact element data */
627          if (Finley_noError()) {
628        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
629        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
630        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
631        if (mpi_info->rank == 0) {  /* Master */
632          for (;;) {            /* Infinite loop */
633            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
634            chunkEle = 0;
635            for (i0=0; i0<chunkSize; i0++) {
636              if (totalEle >= numEle) break; /* End inner loop */
637              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
638              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
639              fscanf(fileHandle_p, "\n");
640              totalEle++;
641              chunkEle++;
642            }
643    #ifdef PASO_MPI
644            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
645            if (nextCPU < mpi_info->size) {
646                  int mpi_error;
647              tempInts[chunkSize*(2+numNodes)] = chunkEle;
648              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
649              if ( mpi_error != MPI_SUCCESS ) {
650                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
651                    return NULL;
652              }
653            }
654    #endif
655            nextCPU++;
656            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
657            if (nextCPU > mpi_info->size) break; /* End infinite loop */
658          } /* Infinite loop */
659        }   /* End master */
660        else {  /* Worker */
661    #ifdef PASO_MPI
662          /* Each worker receives one message */
663          MPI_Status status;
664          int mpi_error;
665          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
666          if ( mpi_error != MPI_SUCCESS ) {
667                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
668                return NULL;
669          }
670          chunkEle = tempInts[chunkSize*(2+numNodes)];
671    #endif
672        }   /* Worker */
673    
674  #if 1      /* Copy Element data from tempInts to mesh_p */
675        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
676            mesh_p->ContactElements->minColor=0;
677            mesh_p->ContactElements->maxColor=chunkEle-1;
678            if (Finley_noError()) {
679              #pragma omp parallel for private (i0, i1)
680          for (i0=0; i0<chunkEle; i0++) {
681            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
682            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
683                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
684                mesh_p->ContactElements->Color[i0] = i0;
685            for (i1 = 0; i1 < numNodes; i1++) {
686              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
687            }
688              }
689            }
690    
691    /* Define other structures to keep mesh_write from crashing */      TMPMEMFREE(tempInts);
692    /* Change the typeid from NoType later */        } /* end of Read the contact element data */
693    
694    mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);          /* read nodal elements */
   Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);  
695    
696    mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);      /* Read the element typeID */
697    Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);          if (Finley_noError()) {
698          if (mpi_info->rank == 0) {
699                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
700                typeID=Finley_RefElement_getTypeId(element_type);
701          }
702    #ifdef PASO_MPI
703          if (mpi_info->size > 1) {
704            int temp1[2], mpi_error;
705            temp1[0] = (int) typeID;
706            temp1[1] = numEle;
707            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
708            if (mpi_error != MPI_SUCCESS) {
709              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
710              return NULL;
711            }
712            typeID = (ElementTypeId) temp1[0];
713            numEle = temp1[1];
714          }
715    #endif
716              if (typeID==NoType) {
717                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
718                Finley_setError(VALUE_ERROR, error_msg);
719              }
720        }
721    
722    mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);        /* Allocate the ElementFile */
723    Finley_ElementFile_allocTable(mesh_p->Points, 0);        mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
724          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
725    
726          /* Read the nodal element data */
727          if (Finley_noError()) {
728        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
729        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
730        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
731        if (mpi_info->rank == 0) {  /* Master */
732          for (;;) {            /* Infinite loop */
733            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
734            chunkEle = 0;
735            for (i0=0; i0<chunkSize; i0++) {
736              if (totalEle >= numEle) break; /* End inner loop */
737              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
738              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
739              fscanf(fileHandle_p, "\n");
740              totalEle++;
741              chunkEle++;
742            }
743    #ifdef PASO_MPI
744            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
745            if (nextCPU < mpi_info->size) {
746                  int mpi_error;
747              tempInts[chunkSize*(2+numNodes)] = chunkEle;
748              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
749              if ( mpi_error != MPI_SUCCESS ) {
750                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
751                    return NULL;
752              }
753            }
754  #endif  #endif
755            nextCPU++;
756            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
757            if (nextCPU > mpi_info->size) break; /* End infinite loop */
758          } /* Infinite loop */
759        }   /* End master */
760        else {  /* Worker */
761    #ifdef PASO_MPI
762          /* Each worker receives one message */
763          MPI_Status status;
764          int mpi_error;
765          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
766          if ( mpi_error != MPI_SUCCESS ) {
767                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
768                return NULL;
769          }
770          chunkEle = tempInts[chunkSize*(2+numNodes)];
771    #endif
772        }   /* Worker */
773    
774  #if 0 /* comment out the rest of the un-implemented crap for now */      /* Copy Element data from tempInts to mesh_p */
775          /* get the face elements */      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
776          if (Finley_noError()) {          mesh_p->Points->minColor=0;
777               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          mesh_p->Points->maxColor=chunkEle-1;
              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 */  
778          if (Finley_noError()) {          if (Finley_noError()) {
779               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);            #pragma omp parallel for private (i0, i1)
780               contactTypeID=Finley_RefElement_getTypeId(element_type);        for (i0=0; i0<chunkEle; i0++) {
781               if (contactTypeID==NoType) {          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
782                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
783                 Finley_setError(VALUE_ERROR,error_msg);              mesh_p->Points->Owner[i0]  =mpi_info->rank;
784               } else {              mesh_p->Points->Color[i0] = i0;
785                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          for (i1 = 0; i1 < numNodes; i1++) {
786                 if (Finley_noError()) {            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
787                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);          }
788                     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 */  
                   }  
                }  
              }  
789          }          }
790          /* get the nodal element */  
791          if (Finley_noError()) {      TMPMEMFREE(tempInts);
792               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        } /* end of Read the nodal element data */
793               pointTypeID=Finley_RefElement_getTypeId(element_type);  
794               if (pointTypeID==NoType) {        /* get the name tags */
795                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);        if (Finley_noError()) {
796                 Finley_setError(VALUE_ERROR,error_msg);          char *remainder, *ptr;
797               }          int tag_key, len, error_code;
798               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          long cur_pos, end_pos;
799               if (Finley_noError()) {          if (mpi_info->rank == 0) {  /* Master */
800                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);        /* Read the word 'Tag' */
801                  if (Finley_noError()) {        if (! feof(fileHandle_p)) fscanf(fileHandle_p, "%s\n", name);
802                     mesh_p->Points->minColor=0;        /* Read rest of file in one chunk, after using seek to find length */
803                     mesh_p->Points->maxColor=numEle-1;            cur_pos = ftell(fileHandle_p);
804                     for (i0 = 0; i0 < numEle; i0++) {            fseek(fileHandle_p, 0L, SEEK_END);
805                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);            end_pos = ftell(fileHandle_p);
806                       mesh_p->Points->Color[i0]=i0;            fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
807                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {        remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
808                           fscanf(fileHandle_p, " %d",        if (! feof(fileHandle_p)) fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
809                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);        remainder[end_pos-cur_pos] = 0;
810                       }  /* for i1 */        len = strlen(remainder);    
811                       fscanf(fileHandle_p, "\n");        while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
812                     } /* for i0 */        len = strlen(remainder);
                 }  
              }  
813          }          }
814          /* get the name tags */  #ifdef PASO_MPI
815          if (Finley_noError()) {          error_code = MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm);
816             if (feof(fileHandle_p) == 0) {          if (error_code != MPI_SUCCESS) {
817                fscanf(fileHandle_p, "%s\n", name);            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
818                while (feof(fileHandle_p) == 0) {            return NULL;
819                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);          }
820                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);      if (mpi_info->rank != 0) {
821                }        remainder = TMPMEMALLOC(len+1, char);
822             }        remainder[0] = 0;
823        }
824            error_code = MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm);
825            if (error_code != MPI_SUCCESS) {
826              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tags failed");
827              return NULL;
828          }          }
829  #endif /* comment out the rest of the un-implemented crap for now */  #endif
830        if (remainder[0]) {
831              ptr = remainder;
832              do {
833                sscanf(ptr, "%s %d\n", name, &tag_key);
834                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
835                ptr++;
836              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
837              TMPMEMFREE(remainder);
838        }
839          }
840    
841       }       }
842    
843       /* close file */       /* close file */
844       if (mpi_info->rank == 0) fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
845    
# Line 660  printf("ksteube CPU=%d/%d Element typeID Line 848  printf("ksteube CPU=%d/%d Element typeID
848    
849       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
850       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 */  
851    
852       /* that's it */       /* that's it */
853       #ifdef Finley_TRACE       #ifdef Finley_TRACE

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

  ViewVC Help
Powered by ViewVC 1.1.26