/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_read.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Mesh_read.cpp

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

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 16  Line 16 
16    
17  #include "Mesh.h"  #include "Mesh.h"
18    
19  namespace dudley {  using escript::IOError;
20    
21  Dudley_Mesh *Dudley_Mesh_read(char *fname, index_t order, index_t reduced_order, bool optimize)  namespace {
 {  
     dim_t numNodes, numDim=0, numEle, i0, i1;  
     Dudley_Mesh *mesh_p = NULL;  
     char name[1024], element_type[1024], frm[20];  
     char error_msg[1024];  
     FILE *fileHandle_p = NULL;  
     Dudley_ElementTypeId typeID = Dudley_NoRef;  
     int scan_ret;  
22    
23      /* No! Bad! take a parameter for this */  using namespace dudley;
     escript::JMPI mpi_info = escript::makeInfo(MPI_COMM_WORLD);  
   
     if (mpi_info->rank == 0)  
     {  
         /* get file handle */  
         fileHandle_p = fopen(fname, "r");  
         if (fileHandle_p == NULL)  
         {  
             sprintf(error_msg, "Dudley_Mesh_read: Opening file %s for reading failed.", fname);  
             throw DudleyException(error_msg);  
         }  
24    
25          /* read header */  ElementFile* readElementFile(FILE* fileHandle, escript::JMPI mpiInfo)
26          sprintf(frm, "%%%d[^\n]", 1023);  {
27          scan_ret = fscanf(fileHandle_p, frm, name);      dim_t numEle;
28          if (scan_ret == EOF)      ElementTypeId typeID = Dudley_NoRef;
29              throw DudleyException("Mesh_read: Scan error while reading file");      char elementType[1024];
30          /* get the number of nodes */      int scan_ret;
         scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim, &numNodes);  
         if (scan_ret == EOF)  
             throw DudleyException("Mesh_read: Scan error while reading file");  
     }  
 #ifdef ESYS_MPI  
     /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs */  
     if (mpi_info->size > 1)  
     {  
         int temp1[3];  
         if (mpi_info->rank == 0)  
         {  
             temp1[0] = numDim;  
             temp1[1] = numNodes;  
             temp1[2] = strlen(name) + 1;  
         }  
         else  
         {  
             temp1[0] = 0;  
             temp1[1] = 0;  
             temp1[2] = 1;  
         }  
         MPI_Bcast(temp1, 3, MPI_INT, 0, mpi_info->comm);  
         numDim = temp1[0];  
         numNodes = temp1[1];  
         MPI_Bcast(name, temp1[2], MPI_CHAR, 0, mpi_info->comm);  
     }  
 #endif  
   
     /* allocate mesh */  
     mesh_p = Dudley_Mesh_alloc(name, numDim, mpi_info);  
   
     /* 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, nextCPU = 1;  
     int *tempInts = new  index_t[chunkSize * 3 + 1];        /* Stores the integer message data */  
     double *tempCoords = new  double[chunkSize * numDim];   /* 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]);  
                 } else 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]);  
                 } else 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]);  
                 }  
                 if (scan_ret == EOF)  
                     throw DudleyException("Mesh_read: Scan error while reading file");  
                 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)  
             {  
                 throw DudleyException("Dudley_Mesh_read: error reading chunks of mesh, data too large for message size");  
             }  
 #ifdef ESYS_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)  
             {  
                 tempInts[chunkSize * 3] = chunkNodes;       /* The message has one more int to send chunkNodes */  
                 MPI_Send(tempInts, chunkSize * 3 + 1, MPI_INT, nextCPU, 81720, mpi_info->comm);  
                 MPI_Send(tempCoords, chunkSize * numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);  
             }  
 #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 ESYS_MPI  
         /* Each worker receives two messages */  
         MPI_Status status;  
         MPI_Recv(tempInts, chunkSize * 3 + 1, MPI_INT, 0, 81720, mpi_info->comm, &status);  
         MPI_Recv(tempCoords, chunkSize * numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);  
         chunkNodes = tempInts[chunkSize * 3];       /* How many nodes are in this workers chunk? */  
 #endif  
     }                       /* Worker */  
   
     /* Copy node data from tempMem to mesh_p */  
     Dudley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);  
 #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];  
         }  
     }  
     delete[] tempInts;  
     delete[] tempCoords;  
31    
32      /****************************  read elements ***************************/      // Read the element typeID and number of elements
33      /* Read the element typeID */      if (mpiInfo->rank == 0) {
34      if (mpi_info->rank == 0)          scan_ret = fscanf(fileHandle, "%s %d\n", elementType, &numEle);
     {  
         scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
35          if (scan_ret == EOF)          if (scan_ret == EOF)
36              throw DudleyException("Mesh_read: Scan error while reading file");              throw IOError("Mesh::read: Scan error while reading file");
37          typeID = eltTypeFromString(element_type);          typeID = eltTypeFromString(elementType);
38      }      }
39  #ifdef ESYS_MPI  #ifdef ESYS_MPI
40      if (mpi_info->size > 1)      if (mpiInfo->size > 1) {
41      {          dim_t temp1[2];
42          int temp1[2], mpi_error;          temp1[0] = (dim_t)typeID;
         temp1[0] = (int)typeID;  
43          temp1[1] = numEle;          temp1[1] = numEle;
44          mpi_error = MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);          int mpiError = MPI_Bcast(temp1, 2, MPI_DIM_T, 0, mpiInfo->comm);
45          if (mpi_error != MPI_SUCCESS)          if (mpiError != MPI_SUCCESS) {
46          {              throw DudleyException("Mesh::read: broadcast of element typeID failed");
             throw DudleyException("Dudley_Mesh_read: broadcast of Element typeID failed");  
47          }          }
48          typeID = (Dudley_ElementTypeId) temp1[0];          typeID = static_cast<ElementTypeId>(temp1[0]);
49          numEle = temp1[1];          numEle = temp1[1];
50      }      }
51  #endif  #endif
52      if (typeID == Dudley_NoRef)      if (typeID == Dudley_NoRef) {
53      {          std::stringstream ss;
54          sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);          ss << "Mesh::read: Unidentified element type " << elementType;
55          throw DudleyException(error_msg);          throw IOError(ss.str());
56      }      }
57    
58      /* Allocate the ElementFile */      // Allocate the ElementFile
59      mesh_p->Elements = Dudley_ElementFile_alloc(typeID, mpi_info);      ElementFile* out = new ElementFile(typeID, mpiInfo);
60      numNodes = mesh_p->Elements->numNodes;  /* New meaning for numNodes: num nodes per element */      const int numNodes = out->numNodes;
61    
62      /********************** Read the element data **************************/      /********************** Read the element data **************************/
63      chunkSize = numEle / mpi_info->size + 1;      dim_t chunkSize = numEle / mpiInfo->size + 1;
64      int totalEle = 0;      dim_t totalEle = 0;
65      nextCPU = 1;      dim_t chunkEle = 0;
66      int chunkEle = 0;      int nextCPU = 1;
67      tempInts = new index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */      /// Store Id + Tag + node list (+ one int at end for chunkEle)
68      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */      index_t* tempInts = new index_t[chunkSize * (2 + numNodes) + 1];
69      if (mpi_info->rank == 0)        /* Master */      // Elements are specified as a list of integers...only need one message
70      {      // instead of two as with the nodes
71          for (;;)            /* Infinite loop */      if (mpiInfo->rank == 0) { // Master
72          {          for (;;) {            // Infinite loop
73  #pragma omp parallel for private (i0) schedule(static)  #pragma omp parallel for
74              for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)              for (index_t i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)
75                  tempInts[i0] = -1;                  tempInts[i0] = -1;
76              chunkEle = 0;              chunkEle = 0;
77              for (i0 = 0; i0 < chunkSize; i0++)              for (index_t i0 = 0; i0 < chunkSize; i0++) {
             {  
78                  if (totalEle >= numEle)                  if (totalEle >= numEle)
79                      break;  /* End inner loop */                      break; // End inner loop
80                  scan_ret =                  scan_ret = fscanf(fileHandle, "%d %d",
81                      fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],                                    &tempInts[i0 * (2 + numNodes) + 0],
82                             &tempInts[i0 * (2 + numNodes) + 1]);                                    &tempInts[i0 * (2 + numNodes) + 1]);
83                  if (scan_ret == EOF)                  if (scan_ret == EOF)
84                      throw DudleyException("Mesh_read: Scan error while reading file");                      throw IOError("Mesh::read: Scan error while reading file");
85                  for (i1 = 0; i1 < numNodes; i1++)                  for (int i1 = 0; i1 < numNodes; i1++) {
86                  {                      scan_ret = fscanf(fileHandle, " %d",
87                      scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);                                        &tempInts[i0 * (2 + numNodes) + 2 + i1]);
88                      if (scan_ret == EOF)                      if (scan_ret == EOF)
89                          throw DudleyException("Mesh_read: Scan error while reading file");                          throw IOError("Mesh::read: Scan error while reading file");
90                  }                  }
91                  scan_ret = fscanf(fileHandle_p, "\n");                  scan_ret = fscanf(fileHandle, "\n");
92                  if (scan_ret == EOF)                  if (scan_ret == EOF)
93                      throw DudleyException("Mesh_read: Scan error while reading file");                      throw IOError("Mesh::read: Scan error while reading file");
94                  totalEle++;                  totalEle++;
95                  chunkEle++;                  chunkEle++;
96              }              }
97  #ifdef ESYS_MPI  #ifdef ESYS_MPI
98              /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */              // Eventually we'll send chunk of elements to each CPU except 0
99              if (nextCPU < mpi_info->size)              // itself, here goes one of them
100              {              if (nextCPU < mpiInfo->size) {
101                  tempInts[chunkSize * (2 + numNodes)] = chunkEle;                  tempInts[chunkSize * (2 + numNodes)] = chunkEle;
102                  MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81722, mpi_info->comm);                  MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_DIM_T,
103                             nextCPU, 81722, mpiInfo->comm);
104              }              }
105  #endif  #endif
106              nextCPU++;              nextCPU++;
107              /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */              // Infinite loop ends when I've read a chunk for each of the worker
108              if (nextCPU > mpi_info->size)              // nodes plus one more chunk for the master
109                  break;      /* End infinite loop */              if (nextCPU > mpiInfo->size)
110          }                   /* Infinite loop */                  break; // End infinite loop
111      }                       /* End master */          } // Infinite loop
112      else      } // end master
113      {                       /* Worker */      else { // Worker
114  #ifdef ESYS_MPI  #ifdef ESYS_MPI
115          /* Each worker receives one message */          // Each worker receives one message
116          MPI_Status status;          MPI_Status status;
117          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81722, mpi_info->comm, &status);          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_DIM_T, 0,
118                     81722, mpiInfo->comm, &status);
119          chunkEle = tempInts[chunkSize * (2 + numNodes)];          chunkEle = tempInts[chunkSize * (2 + numNodes)];
120  #endif  #endif
121      }                       /* Worker */      } // Worker
122      Dudley_ElementFile_allocTable(mesh_p->Elements, chunkEle);  
123        out->allocTable(chunkEle);
124    
125      /* Copy Element data from tempInts to mesh_p */      // Copy Element data from tempInts to mesh
126      mesh_p->Elements->minColor = 0;      out->minColor = 0;
127      mesh_p->Elements->maxColor = chunkEle - 1;      out->maxColor = chunkEle - 1;
128  #pragma omp parallel for private (i0, i1) schedule(static)  #pragma omp parallel for
129      for (i0 = 0; i0 < chunkEle; i0++)      for (index_t i0 = 0; i0 < chunkEle; i0++) {
130      {          out->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];
131          mesh_p->Elements->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];          out->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];
132          mesh_p->Elements->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];          out->Owner[i0] = mpiInfo->rank;
133          mesh_p->Elements->Owner[i0] = mpi_info->rank;          out->Color[i0] = i0;
134          mesh_p->Elements->Color[i0] = i0;          for (int i1 = 0; i1 < numNodes; i1++) {
135          for (i1 = 0; i1 < numNodes; i1++)              out->Nodes[INDEX2(i1, i0, numNodes)] =
136          {                  tempInts[i0 * (2 + numNodes) + 2 + i1];
             mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];  
137          }          }
138      }      }
139      delete[] tempInts;      delete[] tempInts;
140        return out;
141    }
142    
143    } // anonymous
144    
145      /********************** read face elements ******************************/  namespace dudley {
146      /* Read the element typeID */  
147      if (mpi_info->rank == 0)  Mesh* Mesh::read(escript::JMPI mpiInfo, const std::string& filename,
148      {                   bool optimize)
149          scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  {
150          if (scan_ret == EOF)      dim_t numNodes;
151              throw DudleyException("Mesh_read: Scan error while reading file");      int numDim = 0;
152          typeID = eltTypeFromString(element_type);      char name[1024], frm[20];
153      }      FILE *fileHandle = NULL;
154  #ifdef ESYS_MPI      int scan_ret;
155      if (mpi_info->size > 1)  
156      {      if (mpiInfo->rank == 0) {
157          int temp1[2];          // open file
158          temp1[0] = (int)typeID;          fileHandle = fopen(filename.c_str(), "r");
159          temp1[1] = numEle;          if (!fileHandle) {
160          MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);              std::stringstream ss;
161          typeID = (Dudley_ElementTypeId) temp1[0];              ss << "Mesh::read: Opening file " << filename
162          numEle = temp1[1];                 << " for reading failed.";
163      }              throw DudleyException(ss.str());
 #endif  
     if (typeID == Dudley_NoRef)  
     {  
         sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);  
         throw DudleyException(error_msg);  
     }  
     /* Allocate the ElementFile */  
     mesh_p->FaceElements = Dudley_ElementFile_alloc(typeID, mpi_info);  
     numNodes = mesh_p->FaceElements->numNodes;  /* New meaning for numNodes: num nodes per element */  
   
     /********************** Read the face element data **********************/  
   
     chunkSize = numEle / mpi_info->size + 1;  
     totalEle = 0;  
     nextCPU = 1;  
     chunkEle = 0;  
     tempInts = new  index_t[chunkSize * (2 + numNodes) + 1];   /* 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]);  
                 if (scan_ret == EOF)  
                     throw DudleyException("Mesh_read: Scan error while reading file");  
                 for (i1 = 0; i1 < numNodes; i1++)  
                 {  
                     scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);  
                     if (scan_ret == EOF)  
                         throw DudleyException("Mesh_read: Scan error while reading file");  
                 }  
                 scan_ret = fscanf(fileHandle_p, "\n");  
                 if (scan_ret == EOF)  
                     throw DudleyException("Mesh_read: Scan error while reading file");  
                 totalEle++;  
                 chunkEle++;  
             }  
 #ifdef ESYS_MPI  
             /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */  
             if (nextCPU < mpi_info->size)  
             {  
                 tempInts[chunkSize * (2 + numNodes)] = chunkEle;  
                 MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81723, mpi_info->comm);  
             }  
 #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 ESYS_MPI  
         /* Each worker receives one message */  
         MPI_Status status;  
         MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81723, mpi_info->comm, &status);  
         chunkEle = tempInts[chunkSize * (2 + numNodes)];  
 #endif  
     }                       /* Worker */  
     Dudley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);  
     /* Copy Element data from tempInts to mesh_p */  
     mesh_p->FaceElements->minColor = 0;  
     mesh_p->FaceElements->maxColor = chunkEle - 1;  
 #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];  
164          }          }
     }  
     delete[] tempInts;  
165    
166      /************************* read nodal elements **************************/          // read header
167      // Read the element typeID          sprintf(frm, "%%%d[^\n]", 1023);
168      if (mpi_info->rank == 0)          scan_ret = fscanf(fileHandle, frm, name);
     {  
         scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
169          if (scan_ret == EOF)          if (scan_ret == EOF)
170              throw DudleyException("Mesh_read: Scan error while reading file");              throw IOError("Mesh::read: Scan error while reading file");
171          typeID = eltTypeFromString(element_type);  
172            // get the number of nodes
173            scan_ret = fscanf(fileHandle, "%1d%*s %d\n", &numDim, &numNodes);
174            if (scan_ret == EOF)
175                throw IOError("Mesh::read: Scan error while reading file");
176      }      }
177    
178  #ifdef ESYS_MPI  #ifdef ESYS_MPI
179      if (mpi_info->size > 1)      // MPI Broadcast numDim, numNodes, name if there are multiple MPI procs
180      {      if (mpiInfo->size > 1) {
181          int temp1[2];          dim_t temp1[3];
182          temp1[0] = (int)typeID;          if (mpiInfo->rank == 0) {
183          temp1[1] = numEle;              temp1[0] = numDim;
184          MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);              temp1[1] = numNodes;
185          typeID = (Dudley_ElementTypeId) temp1[0];              temp1[2] = strlen(name) + 1;
186          numEle = temp1[1];          } else {
187                temp1[0] = 0;
188                temp1[1] = 0;
189                temp1[2] = 1;
190            }
191            MPI_Bcast(temp1, 3, MPI_DIM_T, 0, mpiInfo->comm);
192            numDim = temp1[0];
193            numNodes = temp1[1];
194            MPI_Bcast(name, temp1[2], MPI_CHAR, 0, mpiInfo->comm);
195      }      }
196  #endif  #endif
197      if (typeID == Dudley_NoRef)  
198      {      // allocate mesh
199          sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);      Mesh* mesh = new Mesh(name, numDim, mpiInfo);
200          throw DudleyException(error_msg);  
201      }      // Each CPU will get at most chunkSize nodes so the message has to be
202      /* Allocate the ElementFile */      // sufficiently large
203      mesh_p->Points = Dudley_ElementFile_alloc(typeID, mpi_info);      dim_t chunkSize = numNodes / mpiInfo->size + 1;
204      numNodes = mesh_p->Points->numNodes;    /* New meaning for numNodes: num nodes per element */      dim_t totalNodes = 0;
205        dim_t chunkNodes = 0;
206      /******************  Read the nodal element data ************************/      int nextCPU = 1;
207      chunkSize = numEle / mpi_info->size + 1;      // Stores the integer message data
208      totalEle = 0;      index_t* tempInts = new index_t[chunkSize * 3 + 1];
209      nextCPU = 1;      // Stores the double message data
210      chunkEle = 0;      double* tempCoords = new double[chunkSize * numDim];
211      tempInts = new  index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */  
212      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */      // Read chunkSize nodes, send it in a chunk to worker CPU which copies
213      if (mpi_info->rank == 0)        /* Master */      // chunk into its local mesh.  It doesn't matter that a CPU has the wrong
214      {      // nodes for its elements, this is sorted out later. First chunk sent to
215          for (;;)            /* Infinite loop */      // CPU 1, second to CPU 2, ..., last chunk stays on CPU 0 (the master).
216          {      // The three columns of integers (Id, gDOF, Tag) are gathered into a single
217  #pragma omp parallel for private (i0) schedule(static)      // array tempInts and sent together in a single MPI message.
218              for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)      if (mpiInfo->rank == 0) { // Master
219            for (;;) {            // Infinite loop
220    #pragma omp parallel for
221                for (index_t i0 = 0; i0 < chunkSize * 3 + 1; i0++)
222                  tempInts[i0] = -1;                  tempInts[i0] = -1;
223              chunkEle = 0;  
224              for (i0 = 0; i0 < chunkSize; i0++)  #pragma omp parallel for
225              {              for (index_t i0 = 0; i0 < chunkSize * numDim; i0++)
226                  if (totalEle >= numEle)                  tempCoords[i0] = -1.0;
227                      break;  /* End inner loop */  
228                  scan_ret =              chunkNodes = 0;
229                      fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],              for (index_t i1 = 0; i1 < chunkSize; i1++) {
230                             &tempInts[i0 * (2 + numNodes) + 1]);                  if (totalNodes >= numNodes)
231                  if (scan_ret == EOF)                      break;  // End of inner loop
232                      throw DudleyException("Mesh_read: Scan error while reading file");                  if (1 == numDim) {
233                  for (i1 = 0; i1 < numNodes; i1++)                      scan_ret = fscanf(fileHandle, "%d %d %d %le\n",
234                  {                                        &tempInts[0 + i1],
235                      scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);                                        &tempInts[chunkSize + i1],
236                      if (scan_ret == EOF)                                        &tempInts[chunkSize * 2 + i1],
237                          throw DudleyException("Mesh_read: Scan error while reading file");                                        &tempCoords[i1 * numDim + 0]);
238                    } else if (2 == numDim) {
239                        scan_ret = fscanf(fileHandle, "%d %d %d %le %le\n",
240                                          &tempInts[0 + i1],
241                                          &tempInts[chunkSize + i1],
242                                          &tempInts[chunkSize * 2 + i1],
243                                          &tempCoords[i1 * numDim + 0],
244                                          &tempCoords[i1 * numDim + 1]);
245                    } else if (3 == numDim) {
246                        scan_ret = fscanf(fileHandle, "%d %d %d %le %le %le\n",
247                                          &tempInts[0 + i1],
248                                          &tempInts[chunkSize + i1],
249                                          &tempInts[chunkSize * 2 + i1],
250                                          &tempCoords[i1 * numDim + 0],
251                                          &tempCoords[i1 * numDim + 1],
252                                          &tempCoords[i1 * numDim + 2]);
253                  }                  }
                 scan_ret = fscanf(fileHandle_p, "\n");  
254                  if (scan_ret == EOF)                  if (scan_ret == EOF)
255                      throw DudleyException("Mesh_read: Scan error while reading file");                      throw IOError("Mesh::read: Scan error while reading file");
256                  totalEle++;                  totalNodes++; // When do we quit the infinite loop?
257                  chunkEle++;                  chunkNodes++; // How many nodes do we actually have in this chunk? It may be smaller than chunkSize.
258                }
259                if (chunkNodes > chunkSize) {
260                    throw DudleyException("Mesh::read: error reading chunks of mesh, data too large for message size");
261              }              }
262  #ifdef ESYS_MPI  #ifdef ESYS_MPI
263              /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */              // Eventually we'll send chunkSize nodes to each CPU numbered
264              if (nextCPU < mpi_info->size)              // 1 ... mpiInfo->size-1, here goes one of them
265              {              if (nextCPU < mpiInfo->size) {
266                  tempInts[chunkSize * (2 + numNodes)] = chunkEle;                  // The message has one more int to send chunkNodes
267                  MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81725, mpi_info->comm);                  tempInts[chunkSize * 3] = chunkNodes;
268                    MPI_Send(tempInts, chunkSize * 3 + 1, MPI_DIM_T, nextCPU, 81720, mpiInfo->comm);
269                    MPI_Send(tempCoords, chunkSize * numDim, MPI_DOUBLE, nextCPU, 81721, mpiInfo->comm);
270              }              }
271  #endif  #endif
272              nextCPU++;              nextCPU++;
273              /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */              // Infinite loop ends when I've read a chunk for each of the worker
274              if (nextCPU > mpi_info->size)              // nodes plus one more chunk for the master
275                  break;      /* End infinite loop */              if (nextCPU > mpiInfo->size)
276          }                   /* Infinite loop */                  break; // End infinite loop
277      }                       /* End master */          } // Infinite loop
278      else                    /* Worker */      } // End master
279      {      else { // Worker
280  #ifdef ESYS_MPI  #ifdef ESYS_MPI
281          /* Each worker receives one message */          // Each worker receives two messages
282          MPI_Status status;          MPI_Status status;
283          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81725, mpi_info->comm, &status);          MPI_Recv(tempInts, chunkSize * 3 + 1, MPI_DIM_T, 0, 81720, mpiInfo->comm, &status);
284          chunkEle = tempInts[chunkSize * (2 + numNodes)];          MPI_Recv(tempCoords, chunkSize * numDim, MPI_DOUBLE, 0, 81721, mpiInfo->comm, &status);
285  #endif          // How many nodes are in this worker's chunk?
286      }                       /* Worker */          chunkNodes = tempInts[chunkSize * 3];
287    #endif
288      /* Copy Element data from tempInts to mesh_p */      } // Worker
289      Dudley_ElementFile_allocTable(mesh_p->Points, chunkEle);  
290      mesh_p->Points->minColor = 0;      // Copy node data from tempMem to mesh
291      mesh_p->Points->maxColor = chunkEle - 1;      mesh->Nodes->allocTable(chunkNodes);
292  #pragma omp parallel for private (i0, i1) schedule(static)  
293      for (i0 = 0; i0 < chunkEle; i0++)  #pragma omp parallel for
294      {      for (index_t i0 = 0; i0 < chunkNodes; i0++) {
295          mesh_p->Points->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];          mesh->Nodes->Id[i0] = tempInts[0 + i0];
296          mesh_p->Points->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];          mesh->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize + i0];
297          mesh_p->Points->Owner[i0] = mpi_info->rank;          mesh->Nodes->Tag[i0] = tempInts[chunkSize * 2 + i0];
298          mesh_p->Points->Color[i0] = i0;          for (int i1 = 0; i1 < numDim; i1++) {
299          for (i1 = 0; i1 < numNodes; i1++)              mesh->Nodes->Coordinates[INDEX2(i1, i0, numDim)] = tempCoords[i0 * numDim + i1];
         {  
             mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];  
300          }          }
301      }      }
302      delete[] tempInts;      delete[] tempInts;
303        delete[] tempCoords;
304    
305      /******************  get the name tags **********************/      /*************************** read elements ******************************/
306      char *remainder = 0, *ptr;      mesh->Elements = readElementFile(fileHandle, mpiInfo);
307    
308        /************************ read face elements ****************************/
309        mesh->FaceElements = readElementFile(fileHandle, mpiInfo);
310    
311        /************************ read nodal elements ***************************/
312        mesh->Points = readElementFile(fileHandle, mpiInfo);
313    
314        /************************  get the name tags ****************************/
315        char *remainder = NULL, *ptr;
316      size_t len = 0;      size_t len = 0;
 #ifdef ESYS_MPI  
     int len_i;  
 #endif  
317      int tag_key;      int tag_key;
318      if (mpi_info->rank == 0)        /* Master */      if (mpiInfo->rank == 0) { // Master
319      {          // Read the word 'Tag'
320          /* Read the word 'Tag' */          if (!feof(fileHandle)) {
321          if (!feof(fileHandle_p))              scan_ret = fscanf(fileHandle, "%s\n", name);
         {  
             scan_ret = fscanf(fileHandle_p, "%s\n", name);  
322              if (scan_ret == EOF)              if (scan_ret == EOF)
323                  throw DudleyException("Mesh_read: Scan error while reading file");                  throw IOError("Mesh::read: Scan error while reading file");
324          }          }
325          // Read rest of file in one chunk, after using seek to find length          // Read rest of file in one chunk, after using seek to find length
326          long cur_pos, end_pos;          long cur_pos = ftell(fileHandle);
327          cur_pos = ftell(fileHandle_p);          fseek(fileHandle, 0L, SEEK_END);
328          fseek(fileHandle_p, 0L, SEEK_END);          long end_pos = ftell(fileHandle);
329          end_pos = ftell(fileHandle_p);          fseek(fileHandle, (long)cur_pos, SEEK_SET);
330          fseek(fileHandle_p, (long)cur_pos, SEEK_SET);          remainder = new char[end_pos - cur_pos + 1];
331          remainder = new  char[end_pos - cur_pos + 1];          if (!feof(fileHandle)) {
332          if (!feof(fileHandle_p))              scan_ret = fread(remainder, (size_t) end_pos - cur_pos,
333          {                               sizeof(char), fileHandle);
             scan_ret = fread(remainder, (size_t) end_pos - cur_pos, sizeof(char), fileHandle_p);  
334              if (scan_ret == EOF)              if (scan_ret == EOF)
335                  throw DudleyException("Mesh_read: Scan error while reading file");                  throw IOError("Mesh::read: Scan error while reading file");
336              remainder[end_pos - cur_pos] = 0;              remainder[end_pos - cur_pos] = 0;
337          }          }
338          len = strlen(remainder);          len = strlen(remainder);
339          while ((len > 1) && isspace(remainder[--len]))          while (len > 1 && isspace(remainder[--len])) {
         {  
340              remainder[len] = 0;              remainder[len] = 0;
341          }          }
342          len = strlen(remainder);          len = strlen(remainder);
343          // shrink the allocation unit      } // Master
 //          TMPMEMREALLOC(remainder, remainder, len + 1, char);  
     }                       /* Master */  
 #ifdef ESYS_MPI  
344    
345      len_i = (int)len;  #ifdef ESYS_MPI
346      MPI_Bcast(&len_i, 1, MPI_INT, 0, mpi_info->comm);      int len_i = static_cast<int>(len);
347      len = (size_t) len_i;      MPI_Bcast(&len_i, 1, MPI_INT, 0, mpiInfo->comm);
348      if (mpi_info->rank != 0)      len = static_cast<size_t>(len_i);
349      {      if (mpiInfo->rank != 0) {
350          remainder = new  char[len + 1];          remainder = new char[len + 1];
351          remainder[0] = 0;          remainder[0] = 0;
352      }      }
353      if (MPI_Bcast(remainder, len + 1, MPI_CHAR, 0, mpi_info->comm) != MPI_SUCCESS)      if (MPI_Bcast(remainder, len+1, MPI_CHAR, 0, mpiInfo->comm) != MPI_SUCCESS)
354          throw DudleyException("Dudley_Mesh_read: broadcast of remainder failed");          throw DudleyException("Mesh::read: broadcast of remainder failed");
355  #endif  #endif
356    
357      if (remainder[0]) {      if (remainder[0]) {
# Line 559  Dudley_Mesh *Dudley_Mesh_read(char *fnam Line 359  Dudley_Mesh *Dudley_Mesh_read(char *fnam
359          do {          do {
360              sscanf(ptr, "%s %d\n", name, &tag_key);              sscanf(ptr, "%s %d\n", name, &tag_key);
361              if (*name)              if (*name)
362                  Dudley_Mesh_addTagMap(mesh_p, name, tag_key);                  mesh->addTagMap(name, tag_key);
363              ptr++;              ptr++;
364          }          } while (NULL != (ptr = strchr(ptr, '\n')) && *ptr);
         while (NULL != (ptr = strchr(ptr, '\n')) && *ptr);  
365      }      }
366      if (remainder)      delete[] remainder;
         delete[] remainder;  
367    
368      if (mpi_info->rank == 0)      // close file
369          fclose(fileHandle_p);      if (mpiInfo->rank == 0)
370            fclose(fileHandle);
371    
372      Dudley_Mesh_resolveNodeIds(mesh_p);      mesh->resolveNodeIds();
373      Dudley_Mesh_prepare(mesh_p, optimize);      mesh->prepare(optimize);
374      return mesh_p;      return mesh;
375  }  }
376    
377  } // namespace dudley  } // namespace dudley

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26