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

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

  ViewVC Help
Powered by ViewVC 1.1.26