/[escript]/trunk/finley/src/Mesh_read.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_read.c

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26