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

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

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

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

Legend:
Removed from v.1156  
changed lines
  Added in v.1771

  ViewVC Help
Powered by ViewVC 1.1.26