/[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 1725 by ksteube, Tue Aug 26 03:31:49 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 242  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 247  Finley_Mesh* Finley_Mesh_read_MPI(char*
247    double time0=Finley_timer();    double time0=Finley_timer();
248    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
249    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
250      Finley_TagMap* tag_map;
251    index_t tag_key;    index_t tag_key;
252    
253    Finley_resetError();    Finley_resetError();
# Line 308  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 343  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 368  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 375  Finley_Mesh* Finley_Mesh_read_MPI(char*
375  #endif  #endif
376      }   /* Worker */      }   /* Worker */
377    
 #if 0  
         /* 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");  
         }  
             printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);  
 #endif  
   
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 427  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 chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
429      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
# Line 440  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 433  Finley_Mesh* Finley_Mesh_read_MPI(char*
433          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
434          chunkEle = 0;          chunkEle = 0;
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");
# Line 457  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 450  Finley_Mesh* Finley_Mesh_read_MPI(char*
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            }            }
           nextCPU++;  
453          }          }
454  #endif  #endif
455          if (totalEle >= numEle) break; /* End infinite loop */          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 */
# Line 477  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 471  Finley_Mesh* Finley_Mesh_read_MPI(char*
471  #endif  #endif
472      }   /* Worker */      }   /* Worker */
473    
 #if 0  
     /* Display the temp mem for debugging */  
     for (i0=0; i0<chunkSize*(numNodes+2); i0++) {  
       printf(" %2d", tempInts[i0]);  
       if (i0%(numNodes+2)==numNodes+2-1) printf("\n");  
     }  
 #endif  
   
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);
476          mesh_p->Elements->minColor=0;          mesh_p->Elements->minColor=0;
# Line 494  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 480  Finley_Mesh* Finley_Mesh_read_MPI(char*
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;              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];
# Line 502  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 489  Finley_Mesh* Finley_Mesh_read_MPI(char*
489          }          }
490    
491      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
492          } /* end of Read the element data */
493    
494            /* read face elements */
495    
496        /* Read the element typeID */
497            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          /* Allocate the ElementFile */
523          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
525    
526          /* Read the face element data */
527          if (Finley_noError()) {
528        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
529        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
530        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
531        if (mpi_info->rank == 0) {  /* Master */
532          for (;;) {            /* Infinite loop */
533            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
534            chunkEle = 0;
535            for (i0=0; i0<chunkSize; i0++) {
536              if (totalEle >= numEle) break; /* End inner loop */
537              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
538              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
539              fscanf(fileHandle_p, "\n");
540              totalEle++;
541              chunkEle++;
542            }
543    #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
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          /* 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        /* 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        TMPMEMFREE(tempInts);
692          } /* end of Read the contact element data */
693    
694            /* read nodal elements */
695    
696        /* Read the element typeID */
697            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          /* Allocate the ElementFile */
723          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
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        /* Copy Element data from tempInts to mesh_p */
775        Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
776            mesh_p->Points->minColor=0;
777            mesh_p->Points->maxColor=chunkEle-1;
778            if (Finley_noError()) {
779              #pragma omp parallel for private (i0, i1)
780          for (i0=0; i0<chunkEle; i0++) {
781            mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
782            mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
783                mesh_p->Points->Owner[i0]  =mpi_info->rank;
784                mesh_p->Points->Color[i0] = i0;
785            for (i1 = 0; i1 < numNodes; i1++) {
786              mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
787            }
788              }
789            }
790    
791        TMPMEMFREE(tempInts);
792          } /* end of Read the nodal element data */
793    
794          /* get the name tags */
795          if (Finley_noError()) {
796            char *remainder, *ptr;
797            int tag_key, len, error_code;
798            long cur_pos, end_pos;
799            if (mpi_info->rank == 0) {  /* Master */
800          /* Read the word 'Tag' */
801          if (! feof(fileHandle_p)) fscanf(fileHandle_p, "%s\n", name);
802          /* Read rest of file in one chunk, after using seek to find length */
803              cur_pos = ftell(fileHandle_p);
804              fseek(fileHandle_p, 0L, SEEK_END);
805              end_pos = ftell(fileHandle_p);
806              fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
807          remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
808          if (! feof(fileHandle_p)) fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
809          remainder[end_pos-cur_pos] = 0;
810          len = strlen(remainder);    
811          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
812          len = strlen(remainder);
813            }
814    #ifdef PASO_MPI
815            error_code = MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm);
816            if (error_code != MPI_SUCCESS) {
817              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
818              return NULL;
819            }
820        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
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 */
# Line 511  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 846  Finley_Mesh* Finley_Mesh_read_MPI(char*
846       /*   resolve id's : */       /*   resolve id's : */
847       /* rearrange elements: */       /* rearrange elements: */
848    
      return mesh_p; /* ksteube temp return for debugging */  
   
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);
851    

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

  ViewVC Help
Powered by ViewVC 1.1.26