/[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 1028 by gross, Wed Mar 14 00:15:24 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) {  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;
36    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
37    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
38    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
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();
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 */  
 #ifndef PASO_MPI  
   mesh_p = Finley_Mesh_alloc(name,numDim,order);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
 #ifndef PASO_MPI  
   Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
   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 */  
 #ifndef PASO_MPI  
   mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   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;  
233    }    }
234  #ifndef PASO_MPI    Paso_MPIInfo_free( mpi_info );
235    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);    return mesh_p;
236  #else  }
237    /* TODO */  
238  #endif  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
239    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
240    mesh_p->FaceElements->minColor=0;  {
241    mesh_p->FaceElements->maxColor=numEle-1;  
242    for (i0 = 0; i0 < numEle; i0++) {    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
243      fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);    dim_t numNodes, numDim, numEle, i0, i1;
244      mesh_p->FaceElements->Color[i0]=i0;    Finley_Mesh *mesh_p=NULL;
245      for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
246           fscanf(fileHandle_p, " %d",    char error_msg[LenErrorMsg_MAX];
247              &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);    double time0=Finley_timer();
248      }   /* for i1 */    FILE *fileHandle_p = NULL;
249      fscanf(fileHandle_p, "\n");    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
250    } /* for i0 */    Finley_TagMap* tag_map;
251      index_t tag_key;
252    /* get the Contact face element */  
253    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);    Finley_resetError();
254    contactTypeID=Finley_RefElement_getTypeId(element_type);  
255    if (contactTypeID==NoType) {    if (mpi_info->rank == 0) {
256      sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);       /* get file handle */
257      Finley_setError(VALUE_ERROR,error_msg);       fileHandle_p = fopen(fname, "r");
258      return NULL;       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         }
264    
265         /* read header */
266         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
267         fscanf(fileHandle_p, frm, name);
268    
269         /* get the number of nodes */
270         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
271    }    }
272  #ifndef PASO_MPI  
273    mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);  #ifdef PASO_MPI
274  #else    /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
275    /* TODO */    if (mpi_info->size > 1) {
276  #endif      int temp1[3], error_code;
277    Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);      temp1[0] = numDim;
278    mesh_p->ContactElements->minColor=0;      temp1[1] = numNodes;
279    mesh_p->ContactElements->maxColor=numEle-1;      temp1[2] = strlen(name) + 1;
280    for (i0 = 0; i0 < numEle; i0++) {      error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
281      fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);      if (error_code != MPI_SUCCESS) {
282      mesh_p->ContactElements->Color[i0]=i0;        Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
283      for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {        return NULL;
284          fscanf(fileHandle_p, " %d",      }
285             &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);      numDim = temp1[0];
286      }   /* for i1 */      numNodes = temp1[1];
287      fscanf(fileHandle_p, "\n");      error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
288    } /* for i0 */      if (error_code != MPI_SUCCESS) {
289          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
290    /* get the nodal element */        return NULL;
291    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;  
292    }    }
 #ifndef PASO_MPI  
   mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);  
 #else  
   /* TODO */  
293  #endif  #endif
   Finley_ElementFile_allocTable(mesh_p->Points, numEle);  
   mesh_p->Points->minColor=0;  
   mesh_p->Points->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
     mesh_p->Points->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
         fscanf(fileHandle_p, " %d",  
            &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
294    
295    /* close file */       /* 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    fclose(fileHandle_p);      /* Copy node data from tempMem to mesh_p */
379            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    /*   resolve id's : */      /* 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    Finley_Mesh_resolveNodeIds(mesh_p);      /* Copy Element data from tempInts to mesh_p */
575        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    /* rearrange elements: */      /* 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    Finley_Mesh_prepare(mesh_p);      /* 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.1028  
changed lines
  Added in v.1771

  ViewVC Help
Powered by ViewVC 1.1.26