/[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 6008 by caltinay, Mon Feb 22 06:59:27 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 /************************************************************************************/  
   
 /*   Dudley: read mesh */  
   
 /************************************************************************************/  
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
   
 #include <ctype.h>  
17  #include "Mesh.h"  #include "Mesh.h"
18    
19  #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Dudley_setError(IO_ERROR,"scan error while reading dudley file"); return NULL;} }  namespace dudley {
20    
21  Dudley_Mesh *Dudley_Mesh_read(char *fname, index_t order, index_t reduced_order, bool optimize)  Dudley_Mesh *Dudley_Mesh_read(char *fname, index_t order, index_t reduced_order, bool optimize)
22  {  {
# Line 38  Dudley_Mesh *Dudley_Mesh_read(char *fnam Line 28  Dudley_Mesh *Dudley_Mesh_read(char *fnam
28      Dudley_ElementTypeId typeID = Dudley_NoRef;      Dudley_ElementTypeId typeID = Dudley_NoRef;
29      int scan_ret;      int scan_ret;
30    
     Dudley_resetError();  
31      /* No! Bad! take a parameter for this */      /* No! Bad! take a parameter for this */
32      esysUtils::JMPI mpi_info = esysUtils::makeInfo(MPI_COMM_WORLD);      escript::JMPI mpi_info = escript::makeInfo(MPI_COMM_WORLD);
33    
34      if (mpi_info->rank == 0)      if (mpi_info->rank == 0)
35      {      {
36      /* get file handle */          /* get file handle */
37      fileHandle_p = fopen(fname, "r");          fileHandle_p = fopen(fname, "r");
38      if (fileHandle_p == NULL)          if (fileHandle_p == NULL)
39      {          {
40          sprintf(error_msg, "Dudley_Mesh_read: Opening file %s for reading failed.", fname);              sprintf(error_msg, "Dudley_Mesh_read: Opening file %s for reading failed.", fname);
41          Dudley_setError(IO_ERROR, error_msg);              throw DudleyException(error_msg);
42          return NULL;          }
43      }  
44            /* read header */
45      /* read header */          sprintf(frm, "%%%d[^\n]", 1023);
46      sprintf(frm, "%%%d[^\n]", 1023);          scan_ret = fscanf(fileHandle_p, frm, name);
47      scan_ret = fscanf(fileHandle_p, frm, name);          if (scan_ret == EOF)
48      FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")              throw DudleyException("Mesh_read: Scan error while reading file");
49          /* get the number of nodes */          /* get the number of nodes */
50          scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim, &numNodes);          scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim, &numNodes);
51      FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")}          if (scan_ret == EOF)
52                throw DudleyException("Mesh_read: Scan error while reading file");
53        }
54  #ifdef ESYS_MPI  #ifdef ESYS_MPI
55      /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs */      /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs */
56      if (mpi_info->size > 1)      if (mpi_info->size > 1)
57      {      {
58      int temp1[3];          int temp1[3];
59      if (mpi_info->rank == 0)          if (mpi_info->rank == 0)
60      {          {
61          temp1[0] = numDim;              temp1[0] = numDim;
62          temp1[1] = numNodes;              temp1[1] = numNodes;
63          temp1[2] = strlen(name) + 1;              temp1[2] = strlen(name) + 1;
64      }          }
65      else          else
66      {          {
67          temp1[0] = 0;              temp1[0] = 0;
68          temp1[1] = 0;              temp1[1] = 0;
69          temp1[2] = 1;              temp1[2] = 1;
70      }          }
71      MPI_Bcast(temp1, 3, MPI_INT, 0, mpi_info->comm);          MPI_Bcast(temp1, 3, MPI_INT, 0, mpi_info->comm);
72      numDim = temp1[0];          numDim = temp1[0];
73      numNodes = temp1[1];          numNodes = temp1[1];
74      MPI_Bcast(name, temp1[2], MPI_CHAR, 0, mpi_info->comm);          MPI_Bcast(name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
75      }      }
76  #endif  #endif
77    
78      /* allocate mesh */      /* allocate mesh */
79      mesh_p = Dudley_Mesh_alloc(name, numDim, mpi_info);      mesh_p = Dudley_Mesh_alloc(name, numDim, mpi_info);
80    
81      if (Dudley_noError())      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
82        int chunkSize = numNodes / mpi_info->size + 1, totalNodes = 0, chunkNodes = 0, nextCPU = 1;
83        int *tempInts = new  index_t[chunkSize * 3 + 1];        /* Stores the integer message data */
84        double *tempCoords = new  double[chunkSize * numDim];   /* Stores the double message data */
85    
86        /*
87           Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
88           It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
89           First chunk sent to CPU 1, second to CPU 2, ...
90           Last chunk stays on CPU 0 (the master)
91           The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent
92           together in a single MPI message
93         */
94    
95        if (mpi_info->rank == 0)        /* Master */
96      {      {
97      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */          for (;;)            /* Infinite loop */
98      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 */  
         {  
99  #pragma omp parallel for private (i0) schedule(static)  #pragma omp parallel for private (i0) schedule(static)
100          for (i0 = 0; i0 < chunkSize * 3 + 1; i0++)              for (i0 = 0; i0 < chunkSize * 3 + 1; i0++)
101              tempInts[i0] = -1;                  tempInts[i0] = -1;
102    
103  #pragma omp parallel for private (i0) schedule(static)  #pragma omp parallel for private (i0) schedule(static)
104          for (i0 = 0; i0 < chunkSize * numDim; i0++)              for (i0 = 0; i0 < chunkSize * numDim; i0++)
105              tempCoords[i0] = -1.0;                  tempCoords[i0] = -1.0;
106    
107                chunkNodes = 0;
108                for (i1 = 0; i1 < chunkSize; i1++)
109                {
110                    if (totalNodes >= numNodes)
111                        break;  /* End of inner loop */
112                    if (1 == numDim) {
113                        scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
114                                          &tempInts[0 + i1], &tempInts[chunkSize + i1], &tempInts[chunkSize * 2 + i1],
115                                          &tempCoords[i1 * numDim + 0]);
116                    } else if (2 == numDim) {
117                        scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
118                                          &tempInts[0 + i1], &tempInts[chunkSize + i1], &tempInts[chunkSize * 2 + i1],
119                                          &tempCoords[i1 * numDim + 0], &tempCoords[i1 * numDim + 1]);
120                    } else if (3 == numDim) {
121                        scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
122                                          &tempInts[0 + i1], &tempInts[chunkSize + i1], &tempInts[chunkSize * 2 + i1],
123                                          &tempCoords[i1 * numDim + 0], &tempCoords[i1 * numDim + 1],
124                                          &tempCoords[i1 * numDim + 2]);
125                    }
126                    if (scan_ret == EOF)
127                        throw DudleyException("Mesh_read: Scan error while reading file");
128                    totalNodes++;       /* When do we quit the infinite loop? */
129                    chunkNodes++;       /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
130                }
131                if (chunkNodes > chunkSize)
132                {
133                    throw DudleyException("Dudley_Mesh_read: error reading chunks of mesh, data too large for message size");
134                }
135    #ifdef ESYS_MPI
136                /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
137                if (nextCPU < mpi_info->size)
138                {
139                    tempInts[chunkSize * 3] = chunkNodes;       /* The message has one more int to send chunkNodes */
140                    MPI_Send(tempInts, chunkSize * 3 + 1, MPI_INT, nextCPU, 81720, mpi_info->comm);
141                    MPI_Send(tempCoords, chunkSize * numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
142                }
143    #endif
144                nextCPU++;
145                /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
146                if (nextCPU > mpi_info->size)
147                    break;      /* End infinite loop */
148            }                   /* Infinite loop */
149        }                       /* End master */
150        else                    /* Worker */
151        {
152    #ifdef ESYS_MPI
153            /* Each worker receives two messages */
154            MPI_Status status;
155            MPI_Recv(tempInts, chunkSize * 3 + 1, MPI_INT, 0, 81720, mpi_info->comm, &status);
156            MPI_Recv(tempCoords, chunkSize * numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
157            chunkNodes = tempInts[chunkSize * 3];       /* How many nodes are in this workers chunk? */
158    #endif
159        }                       /* Worker */
160    
161          chunkNodes = 0;      /* Copy node data from tempMem to mesh_p */
162          for (i1 = 0; i1 < chunkSize; i1++)      Dudley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
         {  
             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, "Dudley_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, "Dudley_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, "Dudley_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)  
         {  
             Dudley_setError(ESYS_MPI_ERROR,  
                     "Dudley_Mesh_read: error reading chunks of mesh, data too large for message size");  
             return NULL;  
         }  
 #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);  
     if (Dudley_noError())  
     {  
163  #pragma omp parallel for private (i0, i1) schedule(static)  #pragma omp parallel for private (i0, i1) schedule(static)
164          for (i0 = 0; i0 < chunkNodes; i0++)      for (i0 = 0; i0 < chunkNodes; i0++)
165          {      {
166          mesh_p->Nodes->Id[i0] = tempInts[0 + i0];          mesh_p->Nodes->Id[i0] = tempInts[0 + i0];
167          mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize + i0];          mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize + i0];
168          mesh_p->Nodes->Tag[i0] = tempInts[chunkSize * 2 + i0];          mesh_p->Nodes->Tag[i0] = tempInts[chunkSize * 2 + i0];
169          for (i1 = 0; i1 < numDim; i1++)          for (i1 = 0; i1 < numDim; i1++)
170          {          {
171              mesh_p->Nodes->Coordinates[INDEX2(i1, i0, numDim)] = tempCoords[i0 * numDim + i1];              mesh_p->Nodes->Coordinates[INDEX2(i1, i0, numDim)] = tempCoords[i0 * numDim + i1];
172          }          }
         }  
     }  
     delete[] tempInts;  
     delete[] tempCoords;  
     }  
   
     /* ***********************************  read elements *************************************************************************************** */  
     if (Dudley_noError())  
     {  
     /* Read the element typeID */  
     if (mpi_info->rank == 0)  
     {  
         scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
         FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") typeID = eltTypeFromString(element_type);  
     }  
 #ifdef ESYS_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)  
         {  
         Dudley_setError(ESYS_MPI_ERROR, "Dudley_Mesh_read: broadcast of Element typeID failed");  
         return NULL;  
         }  
         typeID = (Dudley_ElementTypeId) temp1[0];  
         numEle = temp1[1];  
     }  
 #endif  
     if (typeID == Dudley_NoRef)  
     {  
         sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);  
         Dudley_setError(VALUE_ERROR, error_msg);  
     }  
173      }      }
174        delete[] tempInts;
175        delete[] tempCoords;
176    
177      /* Allocate the ElementFile */      /****************************  read elements ***************************/
178      if (Dudley_noError())      /* Read the element typeID */
179        if (mpi_info->rank == 0)
180        {
181            scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
182            if (scan_ret == EOF)
183                throw DudleyException("Mesh_read: Scan error while reading file");
184            typeID = eltTypeFromString(element_type);
185        }
186    #ifdef ESYS_MPI
187        if (mpi_info->size > 1)
188      {      {
189      mesh_p->Elements = Dudley_ElementFile_alloc(typeID, mpi_info);          int temp1[2], mpi_error;
190      numNodes = mesh_p->Elements->numNodes;  /* New meaning for numNodes: num nodes per element */          temp1[0] = (int)typeID;
191            temp1[1] = numEle;
192            mpi_error = MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);
193            if (mpi_error != MPI_SUCCESS)
194            {
195                throw DudleyException("Dudley_Mesh_read: broadcast of Element typeID failed");
196            }
197            typeID = (Dudley_ElementTypeId) temp1[0];
198            numEle = temp1[1];
199      }      }
200    #endif
201        if (typeID == Dudley_NoRef)
202        {
203            sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
204            throw DudleyException(error_msg);
205        }
206    
207        /* Allocate the ElementFile */
208        mesh_p->Elements = Dudley_ElementFile_alloc(typeID, mpi_info);
209        numNodes = mesh_p->Elements->numNodes;  /* New meaning for numNodes: num nodes per element */
210    
211      /* *************************** Read the element data **************************************************************************************** */      /********************** Read the element data **************************/
212      if (Dudley_noError())      chunkSize = numEle / mpi_info->size + 1;
213        int totalEle = 0;
214        nextCPU = 1;
215        int chunkEle = 0;
216        tempInts = new index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */
217        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
218        if (mpi_info->rank == 0)        /* Master */
219      {      {
220      int chunkSize = numEle / mpi_info->size + 1, totalEle = 0, nextCPU = 1, chunkEle = 0;          for (;;)            /* Infinite loop */
221      int *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 */  
         {  
222  #pragma omp parallel for private (i0) schedule(static)  #pragma omp parallel for private (i0) schedule(static)
223          for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)              for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)
224              tempInts[i0] = -1;                  tempInts[i0] = -1;
225          chunkEle = 0;              chunkEle = 0;
226          for (i0 = 0; i0 < chunkSize; i0++)              for (i0 = 0; i0 < chunkSize; i0++)
227          {              {
228              if (totalEle >= numEle)                  if (totalEle >= numEle)
229              break;  /* End inner loop */                      break;  /* End inner loop */
230              scan_ret =                  scan_ret =
231              fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],                      fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],
232                     &tempInts[i0 * (2 + numNodes) + 1]);                             &tempInts[i0 * (2 + numNodes) + 1]);
233              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") for (i1 = 0; i1 < numNodes; i1++)                  if (scan_ret == EOF)
234              {                      throw DudleyException("Mesh_read: Scan error while reading file");
235              scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);                  for (i1 = 0; i1 < numNodes; i1++)
236              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")}                  {
237              scan_ret = fscanf(fileHandle_p, "\n");                      scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);
238              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") totalEle++;                      if (scan_ret == EOF)
239              chunkEle++;                          throw DudleyException("Mesh_read: Scan error while reading file");
240          }                  }
241  #ifdef ESYS_MPI                  scan_ret = fscanf(fileHandle_p, "\n");
242          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */                  if (scan_ret == EOF)
243          if (nextCPU < mpi_info->size)                      throw DudleyException("Mesh_read: Scan error while reading file");
244          {                  totalEle++;
245              tempInts[chunkSize * (2 + numNodes)] = chunkEle;                  chunkEle++;
246              MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81722, mpi_info->comm);              }
247          }  #ifdef ESYS_MPI
248  #endif              /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
249          nextCPU++;              if (nextCPU < mpi_info->size)
250          /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */              {
251          if (nextCPU > mpi_info->size)                  tempInts[chunkSize * (2 + numNodes)] = chunkEle;
252              break;  /* End infinite loop */                  MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81722, mpi_info->comm);
253          }           /* Infinite loop */              }
254      }           /* End master */  #endif
255      else              nextCPU++;
256      {           /* Worker */              /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
257  #ifdef ESYS_MPI              if (nextCPU > mpi_info->size)
258          /* Each worker receives one message */                  break;      /* End infinite loop */
259          MPI_Status status;          }                   /* Infinite loop */
260          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81722, mpi_info->comm, &status);      }                       /* End master */
261          chunkEle = tempInts[chunkSize * (2 + numNodes)];      else
262  #endif      {                       /* Worker */
263      }           /* Worker */  #ifdef ESYS_MPI
264      Dudley_ElementFile_allocTable(mesh_p->Elements, chunkEle);          /* Each worker receives one message */
265            MPI_Status status;
266      /* Copy Element data from tempInts to mesh_p */          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81722, mpi_info->comm, &status);
267      if (Dudley_noError())          chunkEle = tempInts[chunkSize * (2 + numNodes)];
268      {  #endif
269          mesh_p->Elements->minColor = 0;      }                       /* Worker */
270          mesh_p->Elements->maxColor = chunkEle - 1;      Dudley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
271    
272        /* Copy Element data from tempInts to mesh_p */
273        mesh_p->Elements->minColor = 0;
274        mesh_p->Elements->maxColor = chunkEle - 1;
275  #pragma omp parallel for private (i0, i1) schedule(static)  #pragma omp parallel for private (i0, i1) schedule(static)
276          for (i0 = 0; i0 < chunkEle; i0++)      for (i0 = 0; i0 < chunkEle; i0++)
277          {      {
278          mesh_p->Elements->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];          mesh_p->Elements->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];
279          mesh_p->Elements->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];          mesh_p->Elements->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];
280          mesh_p->Elements->Owner[i0] = mpi_info->rank;          mesh_p->Elements->Owner[i0] = mpi_info->rank;
281          mesh_p->Elements->Color[i0] = i0;          mesh_p->Elements->Color[i0] = i0;
282          for (i1 = 0; i1 < numNodes; i1++)          for (i1 = 0; i1 < numNodes; i1++)
283          {          {
284              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];
285          }          }
286          }      }
287      }      delete[] tempInts;
288      delete[] tempInts;  
289      }  
290      /* ******************** end of Read the element data ***************************************************** */      /********************** read face elements ******************************/
291        /* Read the element typeID */
292      /* ********************* read face elements ************************************************************************************ */      if (mpi_info->rank == 0)
293      if (Dudley_noError())      {
294      {          scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
295      /* Read the element typeID */          if (scan_ret == EOF)
296      if (mpi_info->rank == 0)              throw DudleyException("Mesh_read: Scan error while reading file");
297      {          typeID = eltTypeFromString(element_type);
298          scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);      }
299          FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") typeID = eltTypeFromString(element_type);  #ifdef ESYS_MPI
300      }      if (mpi_info->size > 1)
301  #ifdef ESYS_MPI      {
302      if (mpi_info->size > 1)          int temp1[2];
303      {          temp1[0] = (int)typeID;
304          int temp1[2];          temp1[1] = numEle;
305          temp1[0] = (int)typeID;          MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);
306          temp1[1] = numEle;          typeID = (Dudley_ElementTypeId) temp1[0];
307          MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);          numEle = temp1[1];
308          typeID = (Dudley_ElementTypeId) temp1[0];      }
309          numEle = temp1[1];  #endif
310      }      if (typeID == Dudley_NoRef)
311  #endif      {
312      if (typeID == Dudley_NoRef)          sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
313      {          throw DudleyException(error_msg);
314          sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);      }
315          Dudley_setError(VALUE_ERROR, error_msg);      /* Allocate the ElementFile */
316      }      mesh_p->FaceElements = Dudley_ElementFile_alloc(typeID, mpi_info);
317      if (Dudley_noError())      numNodes = mesh_p->FaceElements->numNodes;  /* New meaning for numNodes: num nodes per element */
318      {  
319          /* Allocate the ElementFile */      /********************** Read the face element data **********************/
320          mesh_p->FaceElements = Dudley_ElementFile_alloc(typeID, mpi_info);  
321          numNodes = mesh_p->FaceElements->numNodes;  /* New meaning for numNodes: num nodes per element */      chunkSize = numEle / mpi_info->size + 1;
322      }      totalEle = 0;
323      }      nextCPU = 1;
324      /* ********************** Read the face element data ******************************************************************************* */      chunkEle = 0;
325        tempInts = new  index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */
326      if (Dudley_noError())      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
327      {      if (mpi_info->rank == 0)        /* Master */
328      int chunkSize = numEle / mpi_info->size + 1, totalEle = 0, nextCPU = 1, chunkEle = 0;      {
329      int *tempInts = new  index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */          for (;;)            /* Infinite loop */
330      /* 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 */  
         {  
331  #pragma omp parallel for private (i0) schedule(static)  #pragma omp parallel for private (i0) schedule(static)
332          for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)              for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)
333              tempInts[i0] = -1;                  tempInts[i0] = -1;
334          chunkEle = 0;              chunkEle = 0;
335          for (i0 = 0; i0 < chunkSize; i0++)              for (i0 = 0; i0 < chunkSize; i0++)
336          {              {
337              if (totalEle >= numEle)                  if (totalEle >= numEle)
338              break;  /* End inner loop */                      break;  /* End inner loop */
339              scan_ret =                  scan_ret =
340              fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],                      fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],
341                     &tempInts[i0 * (2 + numNodes) + 1]);                             &tempInts[i0 * (2 + numNodes) + 1]);
342              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") for (i1 = 0; i1 < numNodes; i1++)                  if (scan_ret == EOF)
343              {                      throw DudleyException("Mesh_read: Scan error while reading file");
344              scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);                  for (i1 = 0; i1 < numNodes; i1++)
345              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")}                  {
346              scan_ret = fscanf(fileHandle_p, "\n");                      scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);
347              FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") totalEle++;                      if (scan_ret == EOF)
348              chunkEle++;                          throw DudleyException("Mesh_read: Scan error while reading file");
349          }                  }
350  #ifdef ESYS_MPI                  scan_ret = fscanf(fileHandle_p, "\n");
351          /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */                  if (scan_ret == EOF)
352          if (nextCPU < mpi_info->size)                      throw DudleyException("Mesh_read: Scan error while reading file");
353          {                  totalEle++;
354              tempInts[chunkSize * (2 + numNodes)] = chunkEle;                  chunkEle++;
355              MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81723, mpi_info->comm);              }
356          }  #ifdef ESYS_MPI
357  #endif              /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
358          nextCPU++;              if (nextCPU < mpi_info->size)
359          /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */              {
360          if (nextCPU > mpi_info->size)                  tempInts[chunkSize * (2 + numNodes)] = chunkEle;
361              break;  /* End infinite loop */                  MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81723, mpi_info->comm);
362          }           /* Infinite loop */              }
363      }           /* End master */  #endif
364      else            /* Worker */              nextCPU++;
365      {              /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
366  #ifdef ESYS_MPI              if (nextCPU > mpi_info->size)
367          /* Each worker receives one message */                  break;      /* End infinite loop */
368          MPI_Status status;          }                   /* Infinite loop */
369          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81723, mpi_info->comm, &status);      }                       /* End master */
370          chunkEle = tempInts[chunkSize * (2 + numNodes)];      else                    /* Worker */
371  #endif      {
372      }           /* Worker */  #ifdef ESYS_MPI
373      Dudley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);          /* Each worker receives one message */
374      if (Dudley_noError())          MPI_Status status;
375      {          MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81723, mpi_info->comm, &status);
376          /* Copy Element data from tempInts to mesh_p */          chunkEle = tempInts[chunkSize * (2 + numNodes)];
377          mesh_p->FaceElements->minColor = 0;  #endif
378          mesh_p->FaceElements->maxColor = chunkEle - 1;      }                       /* Worker */
379        Dudley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
380        /* Copy Element data from tempInts to mesh_p */
381        mesh_p->FaceElements->minColor = 0;
382        mesh_p->FaceElements->maxColor = chunkEle - 1;
383  #pragma omp parallel for private (i0, i1)  #pragma omp parallel for private (i0, i1)
384          for (i0 = 0; i0 < chunkEle; i0++)      for (i0 = 0; i0 < chunkEle; i0++)
385          {      {
386          mesh_p->FaceElements->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];          mesh_p->FaceElements->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];
387          mesh_p->FaceElements->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];          mesh_p->FaceElements->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];
388          mesh_p->FaceElements->Owner[i0] = mpi_info->rank;          mesh_p->FaceElements->Owner[i0] = mpi_info->rank;
389          mesh_p->FaceElements->Color[i0] = i0;          mesh_p->FaceElements->Color[i0] = i0;
390          for (i1 = 0; i1 < numNodes; i1++)          for (i1 = 0; i1 < numNodes; i1++)
391          {          {
392              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];
393          }          }
         }  
     }  
   
     delete[] tempInts;  
     }  
     /* ************************************* end of Read the face element data *************************************** */  
   
     /* ********************************* read nodal elements ****************************************************** */  
     /* *******************************  Read the element typeID */  
     if (Dudley_noError())  
     {  
     if (mpi_info->rank == 0)  
     {  
         scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
         FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") typeID = eltTypeFromString(element_type);  
     }  
 #ifdef ESYS_MPI  
     if (mpi_info->size > 1)  
     {  
         int temp1[2];  
         temp1[0] = (int)typeID;  
         temp1[1] = numEle;  
         MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);  
         typeID = (Dudley_ElementTypeId) temp1[0];  
         numEle = temp1[1];  
     }  
 #endif  
     if (typeID == Dudley_NoRef)  
     {  
         sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);  
         Dudley_setError(VALUE_ERROR, error_msg);  
     }  
     }  
     if (Dudley_noError())  
     {  
     /* Allocate the ElementFile */  
     mesh_p->Points = Dudley_ElementFile_alloc(typeID, mpi_info);  
     numNodes = mesh_p->Points->numNodes;    /* New meaning for numNodes: num nodes per element */  
     }  
     /**********************************  Read the nodal element data **************************************************/  
     if (Dudley_noError())  
     {  
     int chunkSize = numEle / mpi_info->size + 1, totalEle = 0, nextCPU = 1, chunkEle = 0;  
     int *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]);  
             FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") for (i1 = 0; i1 < numNodes; i1++)  
             {  
             scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);  
             FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")}  
             scan_ret = fscanf(fileHandle_p, "\n");  
             FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") 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, 81725, 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, 81725, mpi_info->comm, &status);  
         chunkEle = tempInts[chunkSize * (2 + numNodes)];  
 #endif  
     }           /* Worker */  
   
     /* Copy Element data from tempInts to mesh_p */  
     Dudley_ElementFile_allocTable(mesh_p->Points, chunkEle);  
     if (Dudley_noError())  
     {  
         mesh_p->Points->minColor = 0;  
         mesh_p->Points->maxColor = chunkEle - 1;  
 #pragma omp parallel for private (i0, i1) schedule(static)  
         for (i0 = 0; i0 < chunkEle; i0++)  
         {  
         mesh_p->Points->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];  
         mesh_p->Points->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];  
         mesh_p->Points->Owner[i0] = mpi_info->rank;  
         mesh_p->Points->Color[i0] = i0;  
         for (i1 = 0; i1 < numNodes; i1++)  
         {  
             mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];  
         }  
         }  
     }  
   
     delete[] tempInts;  
     }  
     /* ******************************** end of Read the nodal element data ************************************ */  
  /******************  get the name tags *****************************************/  
     if (Dudley_noError())  
     {  
     char *remainder = 0, *ptr;  
     size_t len = 0;  
 #ifdef ESYS_MPI  
     int len_i;  
 #endif  
     int tag_key;  
     if (mpi_info->rank == 0)    /* Master */  
     {  
         /* Read the word 'Tag' */  
         if (!feof(fileHandle_p))  
         {  
         scan_ret = fscanf(fileHandle_p, "%s\n", name);  
         FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")}  
 #if defined(_WIN32)     /* windows ftell lies on unix formatted text files */  
         remainder = NULL;  
         len = 0;  
         while (1)  
         {  
         size_t malloc_chunk = 1024;  
         size_t buff_size = 0;  
         int ch;  
         ch = fgetc(fileHandle_p);  
         if (ch == '\r')  
         {  
             continue;  
         }  
         if (len + 1 > buff_size)  
         {  
             TMPMEMREALLOC(remainder, remainder, buff_size + malloc_chunk, char);  
         }  
         if (ch == EOF)  
         {  
             /* hit EOF */  
             remainder[len] = (char)0;  
             break;  
         }  
         remainder[len] = (char)ch;  
         len++;  
         }  
 #else  
         /* Read rest of file in one chunk, after using seek to find length */  
         {  
         long cur_pos, end_pos;  
         cur_pos = ftell(fileHandle_p);  
         fseek(fileHandle_p, 0L, SEEK_END);  
         end_pos = ftell(fileHandle_p);  
         fseek(fileHandle_p, (long)cur_pos, SEEK_SET);  
         remainder = new  char[end_pos - cur_pos + 1];  
         if (!feof(fileHandle_p))  
         {  
             scan_ret = fread(remainder, (size_t) end_pos - cur_pos, sizeof(char), fileHandle_p);  
             FSCANF_CHECK(scan_ret, "Dudley_Mesh_read") remainder[end_pos - cur_pos] = 0;  
         }  
         }  
 #endif  
         len = strlen(remainder);  
         while ((len > 1) && isspace(remainder[--len]))  
         {  
         remainder[len] = 0;  
         }  
         len = strlen(remainder);  
         // shrink the allocation unit  
 //      TMPMEMREALLOC(remainder, remainder, len + 1, char);  
     }           /* Master */  
 #ifdef ESYS_MPI  
   
     len_i = (int)len;  
     MPI_Bcast(&len_i, 1, MPI_INT, 0, mpi_info->comm);  
     len = (size_t) len_i;  
     if (mpi_info->rank != 0)  
     {  
         remainder = new  char[len + 1];  
         remainder[0] = 0;  
     }  
     if (MPI_Bcast(remainder, len + 1, MPI_CHAR, 0, mpi_info->comm) != MPI_SUCCESS)  
         Dudley_setError(ESYS_MPI_ERROR, "Dudley_Mesh_read: broadcast of remainder failed");  
 #endif  
   
     if (remainder[0])  
     {  
         ptr = remainder;  
         do  
         {  
         sscanf(ptr, "%s %d\n", name, &tag_key);  
         if (*name)  
             Dudley_Mesh_addTagMap(mesh_p, name, tag_key);  
         ptr++;  
         }  
         while (NULL != (ptr = strchr(ptr, '\n')) && *ptr);  
     }  
     if (remainder)  
         delete[] remainder;  
394      }      }
395        delete[] tempInts;
396    
397      /* close file */      /************************* read nodal elements **************************/
398        // Read the element typeID
399      if (mpi_info->rank == 0)      if (mpi_info->rank == 0)
400      fclose(fileHandle_p);      {
401            scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
402      /*   resolve id's : */          if (scan_ret == EOF)
403      /* rearrange elements: */              throw DudleyException("Mesh_read: Scan error while reading file");
404      if (Dudley_noError())          typeID = eltTypeFromString(element_type);
405      Dudley_Mesh_resolveNodeIds(mesh_p);      }
406      if (Dudley_noError())  #ifdef ESYS_MPI
407      Dudley_Mesh_prepare(mesh_p, optimize);      if (mpi_info->size > 1)
408        {
409            int temp1[2];
410            temp1[0] = (int)typeID;
411            temp1[1] = numEle;
412            MPI_Bcast(temp1, 2, MPI_INT, 0, mpi_info->comm);
413            typeID = (Dudley_ElementTypeId) temp1[0];
414            numEle = temp1[1];
415        }
416    #endif
417        if (typeID == Dudley_NoRef)
418        {
419            sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
420            throw DudleyException(error_msg);
421        }
422        /* Allocate the ElementFile */
423        mesh_p->Points = Dudley_ElementFile_alloc(typeID, mpi_info);
424        numNodes = mesh_p->Points->numNodes;    /* New meaning for numNodes: num nodes per element */
425    
426      /* that's it */      /******************  Read the nodal element data ************************/
427      if (!Dudley_noError())      chunkSize = numEle / mpi_info->size + 1;
428        totalEle = 0;
429        nextCPU = 1;
430        chunkEle = 0;
431        tempInts = new  index_t[chunkSize * (2 + numNodes) + 1];   /* Store Id + Tag + node list (+ one int at end for chunkEle) */
432        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
433        if (mpi_info->rank == 0)        /* Master */
434      {      {
435      Dudley_Mesh_free(mesh_p);          for (;;)            /* Infinite loop */
436            {
437    #pragma omp parallel for private (i0) schedule(static)
438                for (i0 = 0; i0 < chunkSize * (2 + numNodes) + 1; i0++)
439                    tempInts[i0] = -1;
440                chunkEle = 0;
441                for (i0 = 0; i0 < chunkSize; i0++)
442                {
443                    if (totalEle >= numEle)
444                        break;  /* End inner loop */
445                    scan_ret =
446                        fscanf(fileHandle_p, "%d %d", &tempInts[i0 * (2 + numNodes) + 0],
447                               &tempInts[i0 * (2 + numNodes) + 1]);
448                    if (scan_ret == EOF)
449                        throw DudleyException("Mesh_read: Scan error while reading file");
450                    for (i1 = 0; i1 < numNodes; i1++)
451                    {
452                        scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0 * (2 + numNodes) + 2 + i1]);
453                        if (scan_ret == EOF)
454                            throw DudleyException("Mesh_read: Scan error while reading file");
455                    }
456                    scan_ret = fscanf(fileHandle_p, "\n");
457                    if (scan_ret == EOF)
458                        throw DudleyException("Mesh_read: Scan error while reading file");
459                    totalEle++;
460                    chunkEle++;
461                }
462    #ifdef ESYS_MPI
463                /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
464                if (nextCPU < mpi_info->size)
465                {
466                    tempInts[chunkSize * (2 + numNodes)] = chunkEle;
467                    MPI_Send(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, nextCPU, 81725, mpi_info->comm);
468                }
469    #endif
470                nextCPU++;
471                /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
472                if (nextCPU > mpi_info->size)
473                    break;      /* End infinite loop */
474            }                   /* Infinite loop */
475        }                       /* End master */
476        else                    /* Worker */
477        {
478    #ifdef ESYS_MPI
479            /* Each worker receives one message */
480            MPI_Status status;
481            MPI_Recv(tempInts, chunkSize * (2 + numNodes) + 1, MPI_INT, 0, 81725, mpi_info->comm, &status);
482            chunkEle = tempInts[chunkSize * (2 + numNodes)];
483    #endif
484        }                       /* Worker */
485    
486        /* Copy Element data from tempInts to mesh_p */
487        Dudley_ElementFile_allocTable(mesh_p->Points, chunkEle);
488        mesh_p->Points->minColor = 0;
489        mesh_p->Points->maxColor = chunkEle - 1;
490    #pragma omp parallel for private (i0, i1) schedule(static)
491        for (i0 = 0; i0 < chunkEle; i0++)
492        {
493            mesh_p->Points->Id[i0] = tempInts[i0 * (2 + numNodes) + 0];
494            mesh_p->Points->Tag[i0] = tempInts[i0 * (2 + numNodes) + 1];
495            mesh_p->Points->Owner[i0] = mpi_info->rank;
496            mesh_p->Points->Color[i0] = i0;
497            for (i1 = 0; i1 < numNodes; i1++)
498            {
499                mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0 * (2 + numNodes) + 2 + i1];
500            }
501        }
502        delete[] tempInts;
503    
504        /******************  get the name tags **********************/
505        char *remainder = 0, *ptr;
506        size_t len = 0;
507    #ifdef ESYS_MPI
508        int len_i;
509    #endif
510        int tag_key;
511        if (mpi_info->rank == 0)        /* Master */
512        {
513            /* Read the word 'Tag' */
514            if (!feof(fileHandle_p))
515            {
516                scan_ret = fscanf(fileHandle_p, "%s\n", name);
517                if (scan_ret == EOF)
518                    throw DudleyException("Mesh_read: Scan error while reading file");
519            }
520            // Read rest of file in one chunk, after using seek to find length
521            long cur_pos, end_pos;
522            cur_pos = ftell(fileHandle_p);
523            fseek(fileHandle_p, 0L, SEEK_END);
524            end_pos = ftell(fileHandle_p);
525            fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
526            remainder = new  char[end_pos - cur_pos + 1];
527            if (!feof(fileHandle_p))
528            {
529                scan_ret = fread(remainder, (size_t) end_pos - cur_pos, sizeof(char), fileHandle_p);
530                if (scan_ret == EOF)
531                    throw DudleyException("Mesh_read: Scan error while reading file");
532                remainder[end_pos - cur_pos] = 0;
533            }
534            len = strlen(remainder);
535            while ((len > 1) && isspace(remainder[--len]))
536            {
537                remainder[len] = 0;
538            }
539            len = strlen(remainder);
540            // shrink the allocation unit
541    //          TMPMEMREALLOC(remainder, remainder, len + 1, char);
542        }                       /* Master */
543    #ifdef ESYS_MPI
544    
545        len_i = (int)len;
546        MPI_Bcast(&len_i, 1, MPI_INT, 0, mpi_info->comm);
547        len = (size_t) len_i;
548        if (mpi_info->rank != 0)
549        {
550            remainder = new  char[len + 1];
551            remainder[0] = 0;
552        }
553        if (MPI_Bcast(remainder, len + 1, MPI_CHAR, 0, mpi_info->comm) != MPI_SUCCESS)
554            throw DudleyException("Dudley_Mesh_read: broadcast of remainder failed");
555    #endif
556    
557        if (remainder[0]) {
558            ptr = remainder;
559            do {
560                sscanf(ptr, "%s %d\n", name, &tag_key);
561                if (*name)
562                    Dudley_Mesh_addTagMap(mesh_p, name, tag_key);
563                ptr++;
564            }
565            while (NULL != (ptr = strchr(ptr, '\n')) && *ptr);
566      }      }
567      /* free up memory */      if (remainder)
568            delete[] remainder;
569    
570        if (mpi_info->rank == 0)
571            fclose(fileHandle_p);
572    
573        Dudley_Mesh_resolveNodeIds(mesh_p);
574        Dudley_Mesh_prepare(mesh_p, optimize);
575      return mesh_p;      return mesh_p;
576  }  }
577    
578    } // namespace dudley
579    

Legend:
Removed from v.6008  
changed lines
  Added in v.6009

  ViewVC Help
Powered by ViewVC 1.1.26