/[escript]/branches/doubleplusgood/finley/src/Mesh_read.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/finley/src/Mesh_read.cpp

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

revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# 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        dim_t numNodes, numDim, numEle, i0, i1;
31        Finley_Mesh *mesh_p=NULL;
32        Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
33        char name[LenString_MAX],element_type[LenString_MAX],frm[20];
34        char error_msg[LenErrorMsg_MAX];
35        FILE *fileHandle_p = NULL;
36        ElementTypeId typeID=NoRef;
37        int scan_ret;
38    
39    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );      Finley_resetError();
40    dim_t numNodes, numDim, numEle, i0, i1;      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
   index_t tag_key;  
   Finley_Mesh *mesh_p=NULL;  
   char name[LenString_MAX],element_type[LenString_MAX],frm[20];  
   char error_msg[LenErrorMsg_MAX];  
   #ifdef Finley_TRACE  
   double time0=Finley_timer();  
   #endif  
   FILE *fileHandle_p = NULL;  
   ElementTypeId typeID=NoType, faceTypeID=NoType, contactTypeID=NoType, pointTypeID=NoType;  
   int scan_ret;  
     
   Finley_resetError();  
   
   if (mpi_info->size > 1) {  
      Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");  
   } else {  
      /* 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 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);  
               }  
            }  
         }  
      }  
      /* 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)  
   
 {  
41    
   Paso_MPIInfo *mpi_info = NULL;  
   dim_t numNodes, numDim, numEle, i0, i1;  
   Finley_Mesh *mesh_p=NULL;  
   char name[LenString_MAX],element_type[LenString_MAX],frm[20];  
   char error_msg[LenErrorMsg_MAX];  
   #ifdef Finley_TRACE  
   double time0=Finley_timer();  
   #endif  
   FILE *fileHandle_p = NULL;  
   ElementTypeId typeID=NoType;  
   int scan_ret;  
   
   Finley_resetError();  
   mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  
   
   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;  
42      if (mpi_info->rank == 0) {      if (mpi_info->rank == 0) {
43         temp1[0] = numDim;          /* get file handle */
44         temp1[1] = numNodes;          fileHandle_p = fopen(fname, "r");
45         temp1[2] = strlen(name) + 1;          if (fileHandle_p==NULL) {
46      } else {              sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
47         temp1[0] = 0;              Finley_setError(IO_ERROR,error_msg);
48         temp1[1] = 0;              Paso_MPIInfo_free( mpi_info );
        temp1[2] = 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;  
     }  
   }  
 #endif  
   
   
      /* allocate mesh */  
      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);  
      if (Finley_noError()) {  
     /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */  
     int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;  
     int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */  
     double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */  
   
     /*  
       Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p  
       It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later  
       First chunk sent to CPU 1, second to CPU 2, ...  
       Last chunk stays on CPU 0 (the master)  
       The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message  
     */  
   
     if (mpi_info->rank == 0) {  /* Master */  
       for (;;) {            /* Infinite loop */  
             #pragma omp parallel for private (i0) schedule(static)  
         for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;  
             #pragma omp parallel for private (i0) schedule(static)  
         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");  
                 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 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");  
49              return NULL;              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()) {  
           #pragma omp parallel for private (i0, i1) schedule(static)  
       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];  
         }  
           }  
50          }          }
51    
52      TMPMEMFREE(tempInts);          /* read header */
53      TMPMEMFREE(tempCoords);          sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
54            scan_ret = fscanf(fileHandle_p, frm, name);
55          /* read elements */          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
56    
57      /* Read the element typeID */          /* get the number of nodes */
58          if (Finley_noError()) {          scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
59        if (mpi_info->rank == 0) {          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
60              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);      }
         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")  
             typeID=Finley_RefElement_getTypeId(element_type);  
       }  
 #ifdef PASO_MPI  
       if (mpi_info->size > 1) {  
         int temp1[2], mpi_error;  
         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);  
           }  
     }  
61    
62        /* Allocate the ElementFile */      #ifdef PASO_MPI
63        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
64        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */          if (mpi_info->size > 1) {
65                int temp1[3];
66        /* Read the element data */              if (mpi_info->rank == 0) {
67        if (Finley_noError()) {                  temp1[0] = numDim;
68      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;                  temp1[1] = numNodes;
69      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */                  temp1[2] = strlen(name) + 1;
70      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */              } else {
71      if (mpi_info->rank == 0) {  /* Master */                  temp1[0] = 0;
72        for (;;) {            /* Infinite loop */                  temp1[1] = 0;
73              #pragma omp parallel for private (i0) schedule(static)                  temp1[2] = 1;
74          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;              }
75          chunkEle = 0;              MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
76          for (i0=0; i0<chunkSize; i0++) {              numDim = temp1[0];
77            if (totalEle >= numEle) break; /* End inner loop */              numNodes = temp1[1];
78            scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);              error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
79            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")              if (error_code != MPI_SUCCESS) {
80            for (i1 = 0; i1 < numNodes; i1++) {                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
         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");  
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 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) schedule(static)  
       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];  
         }  
           }  
83          }          }
84        #endif
85    
86        /* allocate mesh */
87        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      TMPMEMFREE(tempInts);          /* Copy node data from tempMem to mesh_p */
163        } /* end of Read the element data */          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
164            
165            if (Finley_noError()) {
166                #pragma omp parallel for private (i0, i1) schedule(static)
167                for (i0=0; i0<chunkNodes; i0++) {
168                    mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
169                    mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
170                    mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
171                    for (i1=0; i1<numDim; i1++) {
172                        mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
173                    }
174                }
175            }
176            TMPMEMFREE(tempInts);
177            TMPMEMFREE(tempCoords);
178        }
179    
180          /* read face elements */      /* ***********************************  read elements ******************************************************************/
181        if (Finley_noError()) {
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 */  
             #pragma omp parallel for private (i0) schedule(static)  
         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 */  
             #pragma omp parallel for private (i0) schedule(static)  
         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              #pragma omp parallel for private (i0) schedule(static)          if (mpi_info->rank == 0) {  /* Master */
426          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;              for (;;) {            /* Infinite loop */
427          chunkEle = 0;                  #pragma omp parallel for private (i0) schedule(static)
428          for (i0=0; i0<chunkSize; i0++) {                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
429            if (totalEle >= numEle) break; /* End inner loop */        
430            scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);                  chunkEle = 0;
431            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                  for (i0=0; i0<chunkSize; i0++) {
432            for (i1 = 0; i1 < numNodes; i1++) {                      if (totalEle >= numEle) break; /* End inner loop */
433          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);                      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")                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
435            }                      for (i1 = 0; i1 < numNodes; i1++) {
436            scan_ret = fscanf(fileHandle_p, "\n");                          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
437            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
438            totalEle++;                      }
439            chunkEle++;                      scan_ret = fscanf(fileHandle_p, "\n");
440          }                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
441  #ifdef PASO_MPI                      totalEle++;
442          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */                      chunkEle++;
443          if (nextCPU < mpi_info->size) {                  }
444                int mpi_error;                  #ifdef PASO_MPI
445            tempInts[chunkSize*(2+numNodes)] = chunkEle;                      /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
446            mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);                      if (nextCPU < mpi_info->size) {
447            if ( mpi_error != MPI_SUCCESS ) {                          tempInts[chunkSize*(2+numNodes)] = chunkEle;
448                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");                          MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
449                  return NULL;                      }
450            }                  #endif
451          }                  nextCPU++;
452  #endif                  /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
453          nextCPU++;                  if (nextCPU > mpi_info->size) break; /* End infinite loop */
454          /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */              } /* Infinite loop */
455          if (nextCPU > mpi_info->size) break; /* End infinite loop */          }   /* End master */
456        } /* Infinite loop */          else {  /* Worker */
457      }   /* End master */              #ifdef PASO_MPI
458      else {  /* Worker */                  /* Each worker receives one message */
459  #ifdef PASO_MPI                  MPI_Status status;
460        /* Each worker receives one message */                  MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
461        MPI_Status status;                  chunkEle = tempInts[chunkSize*(2+numNodes)] ;
462        int mpi_error;              #endif
463        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);          }   /* Worker */
464        if ( mpi_error != MPI_SUCCESS ) {  
465              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");          /* Copy Element data from tempInts to mesh_p */
466              return NULL;           Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
467        }          
468        chunkEle = tempInts[chunkSize*(2+numNodes)];           if (Finley_noError()) {        
469  #endif               mesh_p->ContactElements->minColor=0;
470      }   /* Worker */               mesh_p->ContactElements->maxColor=chunkEle-1;
471                 #pragma omp parallel for private (i0, i1)
472      /* Copy Element data from tempInts to mesh_p */               for (i0=0; i0<chunkEle; i0++) {
473      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);                   mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
474          mesh_p->Points->minColor=0;                   mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
475          mesh_p->Points->maxColor=chunkEle-1;                   mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
476          if (Finley_noError()) {                   mesh_p->ContactElements->Color[i0] = i0;
477            #pragma omp parallel for private (i0, i1) schedule(static)                   for (i1 = 0; i1 < numNodes; i1++) {
478        for (i0=0; i0<chunkEle; i0++) {                       mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
479          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];                   }
480          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];               }
481              mesh_p->Points->Owner[i0]  =mpi_info->rank;          }
482              mesh_p->Points->Color[i0] = i0;          TMPMEMFREE(tempInts);
483          for (i1 = 0; i1 < numNodes; i1++) {      } /* end of Read the contact element data */
484            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];  
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=0;          size_t len=0;
589  #ifdef PASO_MPI          #ifdef PASO_MPI
590          int len_i;              int len_i;
591  #endif          #endif
592          int tag_key;          int tag_key;
593          if (mpi_info->rank == 0) {  /* Master */          if (mpi_info->rank == 0) {  /* Master */
594        /* Read the word 'Tag' */              /* Read the word 'Tag' */
595        if (! feof(fileHandle_p)) {              if (! feof(fileHandle_p)) {
596          scan_ret = fscanf(fileHandle_p, "%s\n", name);                  scan_ret = fscanf(fileHandle_p, "%s\n", name);
597          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
598        }              }
599    
600  #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */              #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
601                    remainder = NULL;
602        remainder = NULL;                  len=0;
603            len=0;                  while (1)
604        while (1)                  {
605            {                      size_t malloc_chunk = 1024;
606               size_t malloc_chunk = 1024;                      size_t buff_size = 0;
607               size_t buff_size = 0;                      int ch;
608               int ch;  
609                        ch = fgetc(fileHandle_p);
610               ch = fgetc(fileHandle_p);                      if( ch == '\r' )
611               if( ch == '\r' )                      {
612               {                          continue;
613                  continue;                      }
614               }                      if( len+1 > buff_size )
615               if( len+1 > buff_size )                      {
616               {                          TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
617                  TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);                      }
618               }                      if( ch == EOF )
619               if( ch == EOF )                      {
620               {                          /* hit EOF */
621                  /* hit EOF */                          remainder[len] = (char)0;
622                  remainder[len] = (char)0;                          break;
623                  break;                      }
624               }                      remainder[len] = (char)ch;
625               remainder[len] = (char)ch;                      len++;
626               len++;                  }
627            }              #else
628  #else                  /* Read rest of file in one chunk, after using seek to find length */
629        /* Read rest of file in one chunk, after using seek to find length */                  {
630            {                      long cur_pos, end_pos;
631               long cur_pos, end_pos;  
632                        cur_pos = ftell(fileHandle_p);
633               cur_pos = ftell(fileHandle_p);                      fseek(fileHandle_p, 0L, SEEK_END);
634               fseek(fileHandle_p, 0L, SEEK_END);                      end_pos = ftell(fileHandle_p);
635               end_pos = ftell(fileHandle_p);                      fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
636               fseek(fileHandle_p, (long)cur_pos, SEEK_SET);                      remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
637               remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);                      if (! feof(fileHandle_p))
638               if (! feof(fileHandle_p))                      {
639               {                          scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
640                  scan_ret = fread(remainder, (size_t) end_pos-cur_pos,                                               sizeof(char), fileHandle_p);
641                                   sizeof(char), fileHandle_p);  
642                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
643                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")                          remainder[end_pos-cur_pos] = 0;
644                  remainder[end_pos-cur_pos] = 0;                      }
645              }                  }
646            }              #endif
647  #endif              len = strlen(remainder);
648        len = strlen(remainder);              while ((len>1) && isspace(remainder[--len])) {remainder[len]=0;}
649  /*    while (((--len)>0) && isspace((int)(remainder[len]))) {remainder[len]=0;}  */              len = strlen(remainder);
650        while ((len>1) && isspace(remainder[--len])) {remainder[len]=0;}              TMPMEMREALLOC(remainder,remainder,len+1,char);
651        len = strlen(remainder);          }
652            TMPMEMREALLOC(remainder,remainder,len+1,char);          #ifdef PASO_MPI
653          }  
654  #ifdef PASO_MPI              len_i=(int) len;
655                MPI_Bcast (&len_i, 1, MPI_INT,  0, mpi_info->comm);
656          len_i=(int) len;              len=(size_t) len_i;
657          if (MPI_Bcast (&len_i, 1, MPI_INT,  0, mpi_info->comm) != MPI_SUCCESS)              if (mpi_info->rank != 0) {
658          {                  remainder = TMPMEMALLOC(len+1, char);
659             Finley_setError(PASO_MPI_ERROR,                  remainder[0] = 0;
660                             "Finley_Mesh_read: broadcast of tag len failed");              }
661             return NULL;              MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
         }  
         len=(size_t) len_i;  
     if (mpi_info->rank != 0) {  
       remainder = TMPMEMALLOC(len+1, char);  
       remainder[0] = 0;  
     }  
         if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=  
662              MPI_SUCCESS)              MPI_SUCCESS)
663          {          #endif
            Finley_setError(PASO_MPI_ERROR,  
                            "Finley_Mesh_read: broadcast of tags failed");  
            return NULL;  
         }  
 #endif  
664    
665      if (remainder[0]) {          if (remainder[0]) {
666            ptr = remainder;              ptr = remainder;
667            do {              do {
668              sscanf(ptr, "%s %d\n", name, &tag_key);                  sscanf(ptr, "%s %d\n", name, &tag_key);
669              if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);                  if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
670              ptr++;                  ptr++;
671            } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
672            TMPMEMFREE(remainder);              TMPMEMFREE(remainder);
673            }
674      }      }
       }  
675    
676       }      /* close file */
677        if (mpi_info->rank == 0) fclose(fileHandle_p);
678    
679       /* close file */      /*   resolve id's : */
680       if (mpi_info->rank == 0) fclose(fileHandle_p);      /* rearrange elements: */
681        if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
682       /*   resolve id's : */      if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
683       /* rearrange elements: */  
684        /* that's it */
685       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);      if (! Finley_noError()) {
686       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);          Finley_Mesh_free(mesh_p);
687        }
688       /* that's it */      /* free up memory */
689       #ifdef Finley_TRACE      Finley_ReferenceElementSet_dealloc(refPoints);
690       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);      Finley_ReferenceElementSet_dealloc(refContactElements);
691       #endif      Finley_ReferenceElementSet_dealloc(refFaceElements);
692        Finley_ReferenceElementSet_dealloc(refElements);
693    /* that's it */      Paso_MPIInfo_free( mpi_info );
694    if (! Finley_noError()) {      return mesh_p;
        Finley_Mesh_free(mesh_p);  
   }  
   Paso_MPIInfo_free( mpi_info );  
   return mesh_p;  
695  }  }

Legend:
Removed from v.2548  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26