/[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 2059 by gross, Mon Nov 17 22:57:50 2008 UTC revision 2752 by caltinay, Wed Nov 18 06:01:18 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 23  Line 23 
23    
24  #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Finley_setError(IO_ERROR,"scan error while reading finley file"); return NULL;} }  #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Finley_setError(IO_ERROR,"scan error while reading finley file"); return NULL;} }
25    
 /**************************************************************/  
   
 /*  reads a mesh from a Finley file of name fname */  
   
 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)  
26    
27    Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
28  {  {
29        Paso_MPIInfo *mpi_info = NULL;
30    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );      dim_t numNodes, numDim, numEle, i0, i1;
31    dim_t numNodes, numDim, numEle, i0, i1;      Finley_Mesh *mesh_p=NULL;
32    index_t tag_key;      Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
33    Finley_Mesh *mesh_p=NULL;      char name[LenString_MAX],element_type[LenString_MAX],frm[20];
34    char name[LenString_MAX],element_type[LenString_MAX],frm[20];      char error_msg[LenErrorMsg_MAX];
35    char error_msg[LenErrorMsg_MAX];      FILE *fileHandle_p = NULL;
36    #ifdef Finley_TRACE      ElementTypeId typeID=NoRef;
37    double time0=Finley_timer();      int scan_ret, error_code;
38    #endif  
39    FILE *fileHandle_p = NULL;      Finley_resetError();
40    ElementTypeId typeID=NoType, faceTypeID=NoType, contactTypeID=NoType, pointTypeID=NoType;      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
41    int scan_ret;  
42          if (mpi_info->rank == 0) {
43    Finley_resetError();          /* get file handle */
44            fileHandle_p = fopen(fname, "r");
45    if (mpi_info->size > 1) {          if (fileHandle_p==NULL) {
46       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");              sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
47    } else {              Finley_setError(IO_ERROR,error_msg);
48       /* get file handle */              Paso_MPIInfo_free( mpi_info );
49       fileHandle_p = fopen(fname, "r");              return NULL;
      if (fileHandle_p==NULL) {  
        sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);  
        Finley_setError(IO_ERROR,error_msg);  
        Paso_MPIInfo_free( mpi_info );  
        return NULL;  
      }  
     
      /* read header */  
      sprintf(frm,"%%%d[^\n]",LenString_MAX-1);  
      scan_ret = fscanf(fileHandle_p, frm, name);  
      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
   
     
      /* get the nodes */  
     
      scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
      /* allocate mesh */  
      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);  
      if (Finley_noError()) {  
     
         /* read nodes */  
         Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
         if (Finley_noError()) {  
            if (1 == numDim) {  
                for (i0 = 0; i0 < numNodes; i0++) {  
                 scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],  
                        &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                        &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
            }  
            } else if (2 == numDim) {  
                     for (i0 = 0; i0 < numNodes; i0++) {  
                           scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],  
                                  &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);  
                   FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                 }  
            } else if (3 == numDim) {  
                     for (i0 = 0; i0 < numNodes; i0++) {  
                           scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],  
                                  &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);  
                   FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                 }  
            } /* if else else */  
         }  
         /* read elements */  
         if (Finley_noError()) {  
     
            scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
            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++) {  
                       scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);  
               FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                       mesh_p->Elements->Color[i0]=i0;  
                       mesh_p->Elements->Owner[i0]=0;  
                       for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {  
                            scan_ret = fscanf(fileHandle_p, " %d",  
                               &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                       } /* for i1 */  
                       scan_ret = fscanf(fileHandle_p, "\n");  
               FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                     } /* for i0 */  
                  }  
              }  
           }  
         }  
         /* get the face elements */  
         if (Finley_noError()) {  
              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
              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++) {  
                         scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                         mesh_p->FaceElements->Color[i0]=i0;  
                         mesh_p->FaceElements->Owner[i0]=0;  
                         for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {  
                              scan_ret = fscanf(fileHandle_p, " %d",  
                                 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);  
                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                         }   /* for i1 */  
                         scan_ret = fscanf(fileHandle_p, "\n");  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                       } /* for i0 */  
                    }  
                 }  
              }  
         }  
         /* get the Contact face element */  
         if (Finley_noError()) {  
              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
              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++) {  
                         scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                         mesh_p->ContactElements->Color[i0]=i0;  
                         mesh_p->ContactElements->Owner[i0]=0;  
                         for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {  
                             scan_ret = fscanf(fileHandle_p, " %d",  
                                &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);  
                 FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                         }   /* for i1 */  
                         scan_ret = fscanf(fileHandle_p, "\n");  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                       } /* for i0 */  
                   }  
                }  
              }  
         }    
         /* get the nodal element */  
         if (Finley_noError()) {  
              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
              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++) {  
                      scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                      mesh_p->Points->Color[i0]=i0;  
                      mesh_p->Points->Owner[i0]=0;  
                      for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
                          scan_ret = fscanf(fileHandle_p, " %d",  
                             &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                      }  /* for i1 */  
                      scan_ret = fscanf(fileHandle_p, "\n");  
              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                    } /* for i0 */  
                 }  
              }  
         }  
         /* get the name tags */  
         if (Finley_noError()) {  
            if (feof(fileHandle_p) == 0) {  
               scan_ret = fscanf(fileHandle_p, "%s\n", name);  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
               while (feof(fileHandle_p) == 0) {  
                    scan_ret = fscanf(fileHandle_p, "%s %d\n", name, &tag_key);  
            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
                    Finley_Mesh_addTagMap(mesh_p,name,tag_key);  
               }  
            }  
50          }          }
      }  
      /* close file */  
      fclose(fileHandle_p);  
     
      /*   resolve id's : */  
      /* rearrange elements: */  
     
      if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);  
      if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);  
     
      /* that's it */  
      #ifdef Finley_TRACE  
      printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);  
      #endif  
   }  
   
   /* that's it */  
   if (! Finley_noError()) {  
        Finley_Mesh_free(mesh_p);  
   }  
   Paso_MPIInfo_free( mpi_info );  
   return mesh_p;  
 }  
   
 Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)  
51    
52  {          /* read header */
53            sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
54    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );          scan_ret = fscanf(fileHandle_p, frm, name);
55    dim_t numNodes, numDim, numEle, i0, i1;          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
56    Finley_Mesh *mesh_p=NULL;  
57    char name[LenString_MAX],element_type[LenString_MAX],frm[20];          /* get the number of nodes */
58    char error_msg[LenErrorMsg_MAX];          scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
59    #ifdef Finley_TRACE          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
   double time0=Finley_timer();  
   #endif  
   FILE *fileHandle_p = NULL;  
   ElementTypeId typeID=NoType;  
   int scan_ret;  
   
   Finley_resetError();  
   
   if (mpi_info->rank == 0) {  
      /* get file handle */  
      fileHandle_p = fopen(fname, "r");  
      if (fileHandle_p==NULL) {  
        sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);  
        Finley_setError(IO_ERROR,error_msg);  
        Paso_MPIInfo_free( mpi_info );  
        return NULL;  
      }  
   
      /* read header */  
      sprintf(frm,"%%%d[^\n]",LenString_MAX-1);  
      scan_ret = fscanf(fileHandle_p, frm, name);  
      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
   
      /* get the number of nodes */  
      scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
   }  
   
 #ifdef PASO_MPI  
   /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/  
   if (mpi_info->size > 1) {  
     int temp1[3], error_code;  
     temp1[0] = numDim;  
     temp1[1] = numNodes;  
     temp1[2] = strlen(name) + 1;  
     error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);  
     if (error_code != MPI_SUCCESS) {  
       Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");  
       return NULL;  
     }  
     numDim = temp1[0];  
     numNodes = temp1[1];  
     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);  
     if (error_code != MPI_SUCCESS) {  
       Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");  
       return NULL;  
60      }      }
   }  
 #endif  
61    
62       /* allocate mesh */      #ifdef PASO_MPI
63       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);          /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
64       if (Finley_noError()) {          if (mpi_info->size > 1) {
65      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */              int temp1[3];
66      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;              if (mpi_info->rank == 0) {
67      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */                  temp1[0] = numDim;
68      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */                  temp1[1] = numNodes;
69                    temp1[2] = strlen(name) + 1;
70      /*              } else {
71        Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p                  temp1[0] = 0;
72        It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later                  temp1[1] = 0;
73        First chunk sent to CPU 1, second to CPU 2, ...                  temp1[2] = 1;
74        Last chunk stays on CPU 0 (the master)              }
75        The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message              MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
76      */              numDim = temp1[0];
77                numNodes = temp1[1];
78      if (mpi_info->rank == 0) {  /* Master */              error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
79        for (;;) {            /* Infinite loop */              if (error_code != MPI_SUCCESS) {
80          for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
         for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;  
         chunkNodes = 0;  
         for (i1=0; i1<chunkSize; i1++) {  
           if (totalNodes >= numNodes) break;    /* End of inner loop */  
               if (1 == numDim) {  
         scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",  
           &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],  
           &tempCoords[i1*numDim+0]);  
         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
               if (2 == numDim) {  
         scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",  
           &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],  
           &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);  
         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
               if (3 == numDim) {  
         scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",  
           &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],  
           &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);  
         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
           totalNodes++; /* When do we quit the infinite loop? */  
           chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */  
         }  
         if (chunkNodes > chunkSize) {  
               Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");  
               return NULL;  
         }  
 #ifdef PASO_MPI  
         /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
               int mpi_error;  
           tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */  
           mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);  
           if ( mpi_error != MPI_SUCCESS ) {  
                 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");  
                 return NULL;  
           }  
           mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);  
           if ( mpi_error != MPI_SUCCESS ) {  
                 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");  
81                  return NULL;                  return NULL;
82            }              }
         }  
 #endif  
         nextCPU++;  
         /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */  
         if (nextCPU > mpi_info->size) break; /* End infinite loop */  
       } /* Infinite loop */  
     }   /* End master */  
     else {  /* Worker */  
 #ifdef PASO_MPI  
       /* Each worker receives two messages */  
       MPI_Status status;  
           int mpi_error;  
       mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");  
             return NULL;  
       }  
       mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");  
             return NULL;  
       }  
       chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */  
 #endif  
     }   /* Worker */  
   
     /* Copy node data from tempMem to mesh_p */  
         Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);  
         if (Finley_noError()) {  
       for (i0=0; i0<chunkNodes; i0++) {  
         mesh_p->Nodes->Id[i0]               = tempInts[0+i0];  
         mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];  
         mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];  
         for (i1=0; i1<numDim; i1++) {  
           mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];  
         }  
           }  
83          }          }
84        #endif
85    
86      TMPMEMFREE(tempInts);      /* allocate mesh */
87      TMPMEMFREE(tempCoords);      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
88        
89        if (Finley_noError()) {
90            /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
91            int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0,  nextCPU=1;
92            int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
93            double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
94    
95            /*
96            Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
97            It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
98            First chunk sent to CPU 1, second to CPU 2, ...
99            Last chunk stays on CPU 0 (the master)
100            The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
101            */
102    
103            if (mpi_info->rank == 0) {  /* Master */
104                for (;;) {            /* Infinite loop */
105                    #pragma omp parallel for private (i0) schedule(static)
106                    for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
107                    
108                    #pragma omp parallel for private (i0) schedule(static)
109                    for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
110                    
111                    chunkNodes = 0;
112                    for (i1=0; i1<chunkSize; i1++) {
113                        if (totalNodes >= numNodes) break;    /* End of inner loop */
114                        if (1 == numDim) {
115                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
116                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
117                                                &tempCoords[i1*numDim+0]);
118                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
119                        }
120                        if (2 == numDim) {
121                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
122                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
123                                                &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
124                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
125                        }
126                        if (3 == numDim) {
127                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
128                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
129                                                &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
130                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
131                        }
132                        totalNodes++; /* When do we quit the infinite loop? */
133                        chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
134                    }
135                    if (chunkNodes > chunkSize) {
136                        Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
137                        return NULL;
138                    }
139                    #ifdef PASO_MPI
140                        /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
141                        if (nextCPU < mpi_info->size) {
142                            tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
143                            MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
144                            MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
145                        }
146                    #endif
147                    nextCPU++;
148                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
149                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
150                } /* Infinite loop */
151            }   /* End master */
152            else {  /* Worker */
153                #ifdef PASO_MPI
154                    /* Each worker receives two messages */
155                    MPI_Status status;
156                    MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
157                    MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
158                    chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
159                #endif
160            }   /* Worker */
161    
162          /* read elements */          /* Copy node data from tempMem to mesh_p */
163            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
164      /* Read the element typeID */          
165          if (Finley_noError()) {          if (Finley_noError()) {
166        if (mpi_info->rank == 0) {              #pragma omp parallel for private (i0, i1) schedule(static)
167              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              for (i0=0; i0<chunkNodes; i0++) {
168          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                  mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
169              typeID=Finley_RefElement_getTypeId(element_type);                  mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
170        }                  mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
171  #ifdef PASO_MPI                  for (i1=0; i1<numDim; i1++) {
172        if (mpi_info->size > 1) {                      mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
173          int temp1[2], mpi_error;                  }
174          temp1[0] = (int) typeID;              }
         temp1[1] = numEle;  
         mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);  
         if (mpi_error != MPI_SUCCESS) {  
           Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");  
           return NULL;  
         }  
         typeID = (ElementTypeId) temp1[0];  
         numEle = temp1[1];  
       }  
 #endif  
           if (typeID==NoType) {  
             sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);  
             Finley_setError(VALUE_ERROR, error_msg);  
           }  
     }  
   
       /* Allocate the ElementFile */  
       mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
       numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */  
   
       /* Read the element data */  
       if (Finley_noError()) {  
     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;  
     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */  
     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */  
     if (mpi_info->rank == 0) {  /* Master */  
       for (;;) {            /* Infinite loop */  
         for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
         chunkEle = 0;  
         for (i0=0; i0<chunkSize; i0++) {  
           if (totalEle >= numEle) break; /* End inner loop */  
           scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           for (i1 = 0; i1 < numNodes; i1++) {  
         scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
           scan_ret = fscanf(fileHandle_p, "\n");  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           totalEle++;  
           chunkEle++;  
         }  
 #ifdef PASO_MPI  
         /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
               int mpi_error;  
           tempInts[chunkSize*(2+numNodes)] = chunkEle;  
           mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);  
           if ( mpi_error != MPI_SUCCESS ) {  
                 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");  
                 return NULL;  
           }  
         }  
 #endif  
         nextCPU++;  
         /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */  
         if (nextCPU > mpi_info->size) break; /* End infinite loop */  
       } /* Infinite loop */  
     }   /* End master */  
     else {  /* Worker */  
 #ifdef PASO_MPI  
       /* Each worker receives one message */  
       MPI_Status status;  
       int mpi_error;  
       mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");  
             return NULL;  
       }  
       chunkEle = tempInts[chunkSize*(2+numNodes)];  
 #endif  
     }   /* Worker */  
   
     /* Copy Element data from tempInts to mesh_p */  
     Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);  
         mesh_p->Elements->minColor=0;  
         mesh_p->Elements->maxColor=chunkEle-1;  
         if (Finley_noError()) {  
           #pragma omp parallel for private (i0, i1)  
       for (i0=0; i0<chunkEle; i0++) {  
         mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];  
         mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];  
             mesh_p->Elements->Owner[i0]  =mpi_info->rank;  
             mesh_p->Elements->Color[i0] = i0;  
         for (i1 = 0; i1 < numNodes; i1++) {  
           mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];  
         }  
           }  
175          }          }
176            TMPMEMFREE(tempInts);
177            TMPMEMFREE(tempCoords);
178        }
179    
180      TMPMEMFREE(tempInts);      /* ***********************************  read elements ******************************************************************/
181        } /* end of Read the element data */      if (Finley_noError()) {
   
         /* read face elements */  
182    
183      /* Read the element typeID */          /* Read the element typeID */
184          if (Finley_noError()) {          if (mpi_info->rank == 0) {
185        if (mpi_info->rank == 0) {              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
186              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
187          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")              typeID=Finley_ReferenceElement_getTypeId(element_type);
188              typeID=Finley_RefElement_getTypeId(element_type);          }
189        }          #ifdef PASO_MPI
190  #ifdef PASO_MPI              if (mpi_info->size > 1) {
191        if (mpi_info->size > 1) {                  int temp1[2], mpi_error;
192          int temp1[2], mpi_error;                  temp1[0] = (int) typeID;
193          temp1[0] = (int) typeID;                  temp1[1] = numEle;
194          temp1[1] = numEle;                  mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
195          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);                  if (mpi_error != MPI_SUCCESS) {
196          if (mpi_error != MPI_SUCCESS) {                      Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
197            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");                      return NULL;
198            return NULL;                  }
199          }                  typeID = (ElementTypeId) temp1[0];
200          typeID = (ElementTypeId) temp1[0];                  numEle = temp1[1];
201          numEle = temp1[1];              }
202        }          #endif
203  #endif          if (typeID==NoRef) {
           if (typeID==NoType) {  
204              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
205              Finley_setError(VALUE_ERROR, error_msg);              Finley_setError(VALUE_ERROR, error_msg);
206            }            }
207      }      }
208    
209        /* Allocate the ElementFile */      /* Allocate the ElementFile */
210        mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);      if (Finley_noError()) {
211        numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */          refElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
212            mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
213        /* Read the face element data */          numNodes = mesh_p->Elements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
214        if (Finley_noError()) {      }
     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;  
     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */  
     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */  
     if (mpi_info->rank == 0) {  /* Master */  
       for (;;) {            /* Infinite loop */  
         for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
         chunkEle = 0;  
         for (i0=0; i0<chunkSize; i0++) {  
           if (totalEle >= numEle) break; /* End inner loop */  
           scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           for (i1 = 0; i1 < numNodes; i1++) {  
         scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
           scan_ret = fscanf(fileHandle_p, "\n");  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           totalEle++;  
           chunkEle++;  
         }  
 #ifdef PASO_MPI  
         /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
               int mpi_error;  
           tempInts[chunkSize*(2+numNodes)] = chunkEle;  
           mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);  
           if ( mpi_error != MPI_SUCCESS ) {  
                 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");  
                 return NULL;  
           }  
         }  
 #endif  
         nextCPU++;  
         /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */  
         if (nextCPU > mpi_info->size) break; /* End infinite loop */  
       } /* Infinite loop */  
     }   /* End master */  
     else {  /* Worker */  
 #ifdef PASO_MPI  
       /* Each worker receives one message */  
       MPI_Status status;  
       int mpi_error;  
       mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");  
             return NULL;  
       }  
       chunkEle = tempInts[chunkSize*(2+numNodes)];  
 #endif  
     }   /* Worker */  
   
     /* Copy Element data from tempInts to mesh_p */  
     Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);  
         mesh_p->FaceElements->minColor=0;  
         mesh_p->FaceElements->maxColor=chunkEle-1;  
         if (Finley_noError()) {  
           #pragma omp parallel for private (i0, i1)  
       for (i0=0; i0<chunkEle; i0++) {  
         mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];  
         mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];  
             mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;  
             mesh_p->FaceElements->Color[i0] = i0;  
         for (i1 = 0; i1 < numNodes; i1++) {  
           mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];  
         }  
           }  
         }  
215    
216      TMPMEMFREE(tempInts);      /* *************************** Read the element data *******************************************************************/
217        } /* end of Read the face element data */      if (Finley_noError()) {
218    
219          /* read contact elements */          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
220            int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
221            /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
222            if (mpi_info->rank == 0) {  /* Master */
223                for (;;) {            /* Infinite loop */
224                    #pragma omp parallel for private (i0) schedule(static)
225                    for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
226            
227                    chunkEle = 0;
228                    for (i0=0; i0<chunkSize; i0++) {
229                        if (totalEle >= numEle) break; /* End inner loop */
230                        scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
231                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
232                        for (i1 = 0; i1 < numNodes; i1++) {
233                            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
234                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
235                        }
236                        scan_ret = fscanf(fileHandle_p, "\n");
237                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
238                        totalEle++;
239                        chunkEle++;
240                    }
241                    #ifdef PASO_MPI
242                        /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
243                        if (nextCPU < mpi_info->size) {
244                            tempInts[chunkSize*(2+numNodes)] = chunkEle;
245                            MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
246                        }
247                    #endif
248                    nextCPU++;
249                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
250                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
251                } /* Infinite loop */
252            }   /* End master */
253            else {  /* Worker */
254                #ifdef PASO_MPI
255                    /* Each worker receives one message */
256                    MPI_Status status;
257                    MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
258                    chunkEle = tempInts[chunkSize*(2+numNodes)];
259                #endif
260            }   /* Worker */
261    
262            
263            Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
264    
265        
266            /* Copy Element data from tempInts to mesh_p */
267            if (Finley_noError()) {
268    
269                mesh_p->Elements->minColor=0;
270                mesh_p->Elements->maxColor=chunkEle-1;
271                #pragma omp parallel for private (i0, i1) schedule(static)
272                for (i0=0; i0<chunkEle; i0++) {
273                    mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
274                    mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
275                    mesh_p->Elements->Owner[i0]  =mpi_info->rank;
276                    mesh_p->Elements->Color[i0] = i0;
277                    for (i1 = 0; i1 < numNodes; i1++) {
278                        mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
279                    }
280                }
281            }
282    
283            TMPMEMFREE(tempInts);
284        }
285        /* ******************** end of Read the element data ******************************************************/
286    
287        /* ********************* read face elements ***************************************************************/
288        if (Finley_noError()) {
289            /* Read the element typeID */
290    
291      /* Read the element typeID */          if (mpi_info->rank == 0) {
292          if (Finley_noError()) {               scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
293        if (mpi_info->rank == 0) {               FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
294              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               typeID=Finley_ReferenceElement_getTypeId(element_type);
295          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")          }
296              typeID=Finley_RefElement_getTypeId(element_type);          #ifdef PASO_MPI
297        }              if (mpi_info->size > 1) {
298  #ifdef PASO_MPI                  int temp1[2];
299        if (mpi_info->size > 1) {                  temp1[0] = (int) typeID;
300          int temp1[2], mpi_error;                  temp1[1] = numEle;
301          temp1[0] = (int) typeID;                  MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
302          temp1[1] = numEle;                  typeID = (ElementTypeId) temp1[0];
303          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);                  numEle = temp1[1];
304          if (mpi_error != MPI_SUCCESS) {              }
305            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");          #endif
306            return NULL;          if (typeID==NoRef) {
         }  
         typeID = (ElementTypeId) temp1[0];  
         numEle = temp1[1];  
       }  
 #endif  
           if (typeID==NoType) {  
307              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
308              Finley_setError(VALUE_ERROR, error_msg);              Finley_setError(VALUE_ERROR, error_msg);
           }  
     }  
   
       /* Allocate the ElementFile */  
       mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
       numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */  
   
       /* Read the contact element data */  
       if (Finley_noError()) {  
     int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;  
     int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */  
     /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */  
     if (mpi_info->rank == 0) {  /* Master */  
       for (;;) {            /* Infinite loop */  
         for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
         chunkEle = 0;  
         for (i0=0; i0<chunkSize; i0++) {  
           if (totalEle >= numEle) break; /* End inner loop */  
           scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           for (i1 = 0; i1 < numNodes; i1++) {  
         scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);  
             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           }  
           scan_ret = fscanf(fileHandle_p, "\n");  
           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
           totalEle++;  
           chunkEle++;  
         }  
 #ifdef PASO_MPI  
         /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
               int mpi_error;  
           tempInts[chunkSize*(2+numNodes)] = chunkEle;  
           mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);  
           if ( mpi_error != MPI_SUCCESS ) {  
                 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");  
                 return NULL;  
           }  
         }  
 #endif  
         nextCPU++;  
         /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */  
         if (nextCPU > mpi_info->size) break; /* End infinite loop */  
       } /* Infinite loop */  
     }   /* End master */  
     else {  /* Worker */  
 #ifdef PASO_MPI  
       /* Each worker receives one message */  
       MPI_Status status;  
       int mpi_error;  
       mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");  
             return NULL;  
       }  
       chunkEle = tempInts[chunkSize*(2+numNodes)];  
 #endif  
     }   /* Worker */  
   
     /* Copy Element data from tempInts to mesh_p */  
     Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);  
         mesh_p->ContactElements->minColor=0;  
         mesh_p->ContactElements->maxColor=chunkEle-1;  
         if (Finley_noError()) {  
           #pragma omp parallel for private (i0, i1)  
       for (i0=0; i0<chunkEle; i0++) {  
         mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];  
         mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];  
             mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;  
             mesh_p->ContactElements->Color[i0] = i0;  
         for (i1 = 0; i1 < numNodes; i1++) {  
           mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];  
         }  
           }  
309          }          }
310            if (Finley_noError()) {
311                /* Allocate the ElementFile */
312                refFaceElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
313                mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
314                numNodes = mesh_p->FaceElements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
315            }
316    
317      TMPMEMFREE(tempInts);      }
       } /* end of Read the contact element data */  
318    
319          /* read nodal elements */      /* ********************** Read the face element data ********************************************************* */
320    
321      /* Read the element typeID */      if (Finley_noError()) {
322          if (Finley_noError()) {          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
323        if (mpi_info->rank == 0) {          int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
324              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
325          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")          if (mpi_info->rank == 0) {  /* Master */
326              typeID=Finley_RefElement_getTypeId(element_type);              for (;;) {            /* Infinite loop */
327        }                  #pragma omp parallel for private (i0) schedule(static)
328  #ifdef PASO_MPI                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
329        if (mpi_info->size > 1) {          
330          int temp1[2], mpi_error;                  chunkEle = 0;
331          temp1[0] = (int) typeID;                  for (i0=0; i0<chunkSize; i0++) {
332          temp1[1] = numEle;                      if (totalEle >= numEle) break; /* End inner loop */
333          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);                      scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
334          if (mpi_error != MPI_SUCCESS) {                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
335            Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");                      for (i1 = 0; i1 < numNodes; i1++) {
336            return NULL;                          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
337          }                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
338          typeID = (ElementTypeId) temp1[0];                      }
339          numEle = temp1[1];                      scan_ret = fscanf(fileHandle_p, "\n");
340        }                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
341  #endif                      totalEle++;
342            if (typeID==NoType) {                      chunkEle++;
343              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);                  }
344                    #ifdef PASO_MPI
345                        /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
346                        if (nextCPU < mpi_info->size) {
347                            tempInts[chunkSize*(2+numNodes)] = chunkEle;
348                            MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
349                        }
350                    #endif
351                    nextCPU++;
352                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
353                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
354                } /* Infinite loop */
355            }   /* End master */
356            else {  /* Worker */
357                #ifdef PASO_MPI
358                    /* Each worker receives one message */
359                    MPI_Status status;
360                    MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
361                    chunkEle = tempInts[chunkSize*(2+numNodes)];
362                    #endif
363            }   /* Worker */
364            
365            Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
366            
367            if (Finley_noError()) {
368                /* Copy Element data from tempInts to mesh_p */
369                
370                mesh_p->FaceElements->minColor=0;
371                mesh_p->FaceElements->maxColor=chunkEle-1;
372                #pragma omp parallel for private (i0, i1)
373                for (i0=0; i0<chunkEle; i0++) {
374                    mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
375                    mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
376                    mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
377                    mesh_p->FaceElements->Color[i0] = i0;
378                    for (i1 = 0; i1 < numNodes; i1++) {
379                        mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
380                    }
381                }
382            }
383    
384            TMPMEMFREE(tempInts);
385        }
386        /* ************************************* end of Read the face element data *************************************** */
387    
388    
389        /* ************************************* read contact elements ************************************************** */
390    
391        /* Read the element typeID */
392        if (Finley_noError()) {
393            if (mpi_info->rank == 0) {
394                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
395                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
396                typeID=Finley_ReferenceElement_getTypeId(element_type);
397            }
398            #ifdef PASO_MPI
399                if (mpi_info->size > 1) {
400                    int temp1[2];
401                    temp1[0] = (int) typeID;
402                    temp1[1] = numEle;
403                    MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
404                    typeID = (ElementTypeId) temp1[0];
405                    numEle = temp1[1];
406                }
407            #endif
408            if (typeID==NoRef) {
409                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
410              Finley_setError(VALUE_ERROR, error_msg);              Finley_setError(VALUE_ERROR, error_msg);
411            }           }
412      }      }
413    
414        /* Allocate the ElementFile */      if (Finley_noError()) {
415        mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          /* Allocate the ElementFile */
416        numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */          refContactElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
417            mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
418        /* Read the nodal element data */          numNodes = mesh_p->ContactElements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
419        if (Finley_noError()) {      }
420      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;      /* *************************** Read the contact element data ************************************************* */
421      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      if (Finley_noError()) {
422      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
423      if (mpi_info->rank == 0) {  /* Master */          int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
424        for (;;) {            /* Infinite loop */          /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
425          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;          if (mpi_info->rank == 0) {  /* Master */
426          chunkEle = 0;              for (;;) {            /* Infinite loop */
427          for (i0=0; i0<chunkSize; i0++) {                  #pragma omp parallel for private (i0) schedule(static)
428            if (totalEle >= numEle) break; /* End inner loop */                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
429            scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);        
430            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                  chunkEle = 0;
431            for (i1 = 0; i1 < numNodes; i1++) {                  for (i0=0; i0<chunkSize; i0++) {
432          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);                      if (totalEle >= numEle) break; /* End inner loop */
433              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                      scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
434            }                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
435            scan_ret = fscanf(fileHandle_p, "\n");                      for (i1 = 0; i1 < numNodes; i1++) {
436            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
437            totalEle++;                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
438            chunkEle++;                      }
439          }                      scan_ret = fscanf(fileHandle_p, "\n");
440  #ifdef PASO_MPI                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
441          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */                      totalEle++;
442          if (nextCPU < mpi_info->size) {                      chunkEle++;
443                int mpi_error;                  }
444            tempInts[chunkSize*(2+numNodes)] = chunkEle;                  #ifdef PASO_MPI
445            mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);                      /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
446            if ( mpi_error != MPI_SUCCESS ) {                      if (nextCPU < mpi_info->size) {
447                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");                          tempInts[chunkSize*(2+numNodes)] = chunkEle;
448                  return NULL;                          MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
449            }                      }
450          }                  #endif
451  #endif                  nextCPU++;
452          nextCPU++;                  /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
453          /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */                  if (nextCPU > mpi_info->size) break; /* End infinite loop */
454          if (nextCPU > mpi_info->size) break; /* End infinite loop */              } /* Infinite loop */
455        } /* Infinite loop */          }   /* End master */
456      }   /* End master */          else {  /* Worker */
457      else {  /* Worker */              #ifdef PASO_MPI
458  #ifdef PASO_MPI                  /* Each worker receives one message */
459        /* Each worker receives one message */                  MPI_Status status;
460        MPI_Status status;                  MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
461        int mpi_error;                  chunkEle = tempInts[chunkSize*(2+numNodes)] ;
462        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);              #endif
463        if ( mpi_error != MPI_SUCCESS ) {          }   /* Worker */
464              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");  
465              return NULL;          /* Copy Element data from tempInts to mesh_p */
466        }           Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
467        chunkEle = tempInts[chunkSize*(2+numNodes)];          
468  #endif           if (Finley_noError()) {        
469      }   /* Worker */               mesh_p->ContactElements->minColor=0;
470                 mesh_p->ContactElements->maxColor=chunkEle-1;
471      /* Copy Element data from tempInts to mesh_p */               #pragma omp parallel for private (i0, i1)
472      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);               for (i0=0; i0<chunkEle; i0++) {
473          mesh_p->Points->minColor=0;                   mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
474          mesh_p->Points->maxColor=chunkEle-1;                   mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
475          if (Finley_noError()) {                   mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
476            #pragma omp parallel for private (i0, i1)                   mesh_p->ContactElements->Color[i0] = i0;
477        for (i0=0; i0<chunkEle; i0++) {                   for (i1 = 0; i1 < numNodes; i1++) {
478          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];                       mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
479          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];                   }
480              mesh_p->Points->Owner[i0]  =mpi_info->rank;               }
481              mesh_p->Points->Color[i0] = i0;          }
482          for (i1 = 0; i1 < numNodes; i1++) {          TMPMEMFREE(tempInts);
483            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];      } /* end of Read the contact element data */
484          }  
485            }      /* ********************************* read nodal elements ****************************************************** */
486          }  
487        /* *******************************  Read the element typeID */
488    
489        if (Finley_noError()) {
490            if (mpi_info->rank == 0) {
491                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
492                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
493                typeID=Finley_ReferenceElement_getTypeId(element_type);
494            }
495            #ifdef PASO_MPI
496                if (mpi_info->size > 1) {
497                    int temp1[2];
498                    temp1[0] = (int) typeID;
499                    temp1[1] = numEle;
500                    MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
501                    typeID = (ElementTypeId) temp1[0];
502                    numEle = temp1[1];
503                }
504            #endif
505            if (typeID==NoRef) {
506                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
507                Finley_setError(VALUE_ERROR, error_msg);
508             }
509        }
510    
511      TMPMEMFREE(tempInts);      if (Finley_noError()) {
512        } /* end of Read the nodal element data */          /* Allocate the ElementFile */
513            refPoints= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
514            mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
515            numNodes = mesh_p->Points->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
516        }
517    
518        /* get the name tags */      /**********************************  Read the nodal element data **************************************************/
519        if (Finley_noError()) {      if (Finley_noError()) {
520            int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
521            int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
522            /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
523            if (mpi_info->rank == 0) {  /* Master */
524                for (;;) {            /* Infinite loop */
525                    #pragma omp parallel for private (i0) schedule(static)
526                    for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
527            
528                    chunkEle = 0;
529                    for (i0=0; i0<chunkSize; i0++) {
530                        if (totalEle >= numEle) break; /* End inner loop */
531                        scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
532                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
533                        for (i1 = 0; i1 < numNodes; i1++) {
534                            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
535                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
536                        }
537                        scan_ret = fscanf(fileHandle_p, "\n");
538                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
539                        totalEle++;
540                        chunkEle++;
541                    }
542                    #ifdef PASO_MPI
543                        /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
544                        if (nextCPU < mpi_info->size) {
545                            tempInts[chunkSize*(2+numNodes)] = chunkEle;
546                            MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
547                        }
548                    #endif
549                    nextCPU++;
550                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
551                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
552                } /* Infinite loop */
553            }   /* End master */
554            else {  /* Worker */
555                #ifdef PASO_MPI
556                    /* Each worker receives one message */
557                    MPI_Status status;
558                    MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
559                    chunkEle = tempInts[chunkSize*(2+numNodes)];
560                #endif
561            }   /* Worker */
562    
563            /* Copy Element data from tempInts to mesh_p */
564            Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
565            
566            if (Finley_noError()) {
567                mesh_p->Points->minColor=0;
568                mesh_p->Points->maxColor=chunkEle-1;
569                #pragma omp parallel for private (i0, i1) schedule(static)
570                for (i0=0; i0<chunkEle; i0++) {
571                    mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
572                    mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
573                    mesh_p->Points->Owner[i0]  =mpi_info->rank;
574                    mesh_p->Points->Color[i0] = i0;
575                    for (i1 = 0; i1 < numNodes; i1++) {
576                        mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
577                    }
578                }
579            }
580    
581            TMPMEMFREE(tempInts);
582        } /* ******************************** end of Read the nodal element data ************************************************************* */
583    
584        
585        /******************  get the name tags *****************************************/
586        if (Finley_noError()) {
587          char *remainder=0, *ptr;          char *remainder=0, *ptr;
588          size_t len;          size_t len=0;
589            #ifdef PASO_MPI
590                int len_i;
591            #endif
592          int tag_key;          int tag_key;
593            if (mpi_info->rank == 0) {  /* Master */
594                /* Read the word 'Tag' */
595                if (! feof(fileHandle_p)) {
596          if (mpi_info->rank == 0) {  /* Master */                  scan_ret = fscanf(fileHandle_p, "%s\n", name);
597        /* Read the word 'Tag' */                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
598        if (! feof(fileHandle_p)) {              }
599          scan_ret = fscanf(fileHandle_p, "%s\n", name);  
600          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")              #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
601        }                  remainder = NULL;
602                    len=0;
603  #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */                  while (1)
604                    {
605        remainder = NULL;                      size_t malloc_chunk = 1024;
606            len=0;                      size_t buff_size = 0;
607        while (1)                      int ch;
608            {  
609               size_t malloc_chunk = 1024;                      ch = fgetc(fileHandle_p);
610               size_t buff_size = 0;                      if( ch == '\r' )
611               int ch;                      {
612                            continue;
613               ch = fgetc(fileHandle_p);                      }
614               if( ch == '\r' )                      if( len+1 > buff_size )
615               {                      {
616                  continue;                          TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
617               }                      }
618               if( len+1 > buff_size )                      if( ch == EOF )
619               {                      {
620                  TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);                          /* hit EOF */
621               }                          remainder[len] = (char)0;
622               if( ch == EOF )                          break;
623               {                      }
624                  /* hit EOF */                      remainder[len] = (char)ch;
625                  remainder[len] = (char)0;                      len++;
626                  break;                  }
627               }              #else
628               remainder[len] = (char)ch;                  /* Read rest of file in one chunk, after using seek to find length */
629               len++;                  {
630            }                      long cur_pos, end_pos;
631  #else  
632        /* Read rest of file in one chunk, after using seek to find length */                      cur_pos = ftell(fileHandle_p);
633            {                      fseek(fileHandle_p, 0L, SEEK_END);
634               long cur_pos, end_pos;                      end_pos = ftell(fileHandle_p);
635                        fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
636               cur_pos = ftell(fileHandle_p);                      remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
637               fseek(fileHandle_p, 0L, SEEK_END);                      if (! feof(fileHandle_p))
638               end_pos = ftell(fileHandle_p);                      {
639               fseek(fileHandle_p, (long)cur_pos, SEEK_SET);                          scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
640               remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);                                               sizeof(char), fileHandle_p);
641               if (! feof(fileHandle_p))  
642               {                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
643                  scan_ret = fread(remainder, (size_t) end_pos-cur_pos,                          remainder[end_pos-cur_pos] = 0;
644                                   sizeof(char), fileHandle_p);                      }
645                    }
646                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")              #endif
647                  remainder[end_pos-cur_pos] = 0;              len = strlen(remainder);
648              }              while ((len>1) && isspace(remainder[--len])) {remainder[len]=0;}
649            }              len = strlen(remainder);
650  #endif              TMPMEMREALLOC(remainder,remainder,len+1,char);
651        len = strlen(remainder);          } /* Master */
652        while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;          #ifdef PASO_MPI
653        len = strlen(remainder);  
654            TMPMEMREALLOC(remainder,remainder,len+1,char);              len_i=(int) len;
655          }              MPI_Bcast (&len_i, 1, MPI_INT,  0, mpi_info->comm);
656  #ifdef PASO_MPI              len=(size_t) len_i;
657                if (mpi_info->rank != 0) {
658          if (MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm) != MPI_SUCCESS)                  remainder = TMPMEMALLOC(len+1, char);
659          {                  remainder[0] = 0;
660             Finley_setError(PASO_MPI_ERROR,              }
661                             "Finley_Mesh_read: broadcast of tag len failed");              if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
662             return NULL;                      MPI_SUCCESS)
663          }                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of remainder failed");
664      if (mpi_info->rank != 0) {          #endif
665        remainder = TMPMEMALLOC(len+1, char);  
666        remainder[0] = 0;          if (remainder[0]) {
667      }              ptr = remainder;
668                do {
669          if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=                  sscanf(ptr, "%s %d\n", name, &tag_key);
670              MPI_SUCCESS)                  if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
671          {                  ptr++;
672             Finley_setError(PASO_MPI_ERROR,              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
673                             "Finley_Mesh_read: broadcast of tags failed");          }
674             return NULL;          if (remainder)
675          }              TMPMEMFREE(remainder);
 #endif  
     if (remainder[0]) {  
           ptr = remainder;  
           do {  
             sscanf(ptr, "%s %d\n", name, &tag_key);  
             if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);  
             ptr++;  
           } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);  
           TMPMEMFREE(remainder);  
676      }      }
       }  
   
      }  
677    
678       /* close file */      /* close file */
679       if (mpi_info->rank == 0) fclose(fileHandle_p);      if (mpi_info->rank == 0) fclose(fileHandle_p);
680    
681       /*   resolve id's : */      /*   resolve id's : */
682       /* rearrange elements: */      /* rearrange elements: */
683        if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
684       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);      if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
685       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);  
686        /* that's it */
687       /* that's it */      if (! Finley_noError()) {
688       #ifdef Finley_TRACE          Finley_Mesh_free(mesh_p);
689       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);      }
690       #endif      /* free up memory */
691        Finley_ReferenceElementSet_dealloc(refPoints);
692    /* that's it */      Finley_ReferenceElementSet_dealloc(refContactElements);
693    if (! Finley_noError()) {      Finley_ReferenceElementSet_dealloc(refFaceElements);
694         Finley_Mesh_free(mesh_p);      Finley_ReferenceElementSet_dealloc(refElements);
695    }      Paso_MPIInfo_free( mpi_info );
696    Paso_MPIInfo_free( mpi_info );      return mesh_p;
   return mesh_p;  
697  }  }
698    

Legend:
Removed from v.2059  
changed lines
  Added in v.2752

  ViewVC Help
Powered by ViewVC 1.1.26