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

Legend:
Removed from v.964  
changed lines
  Added in v.1907

  ViewVC Help
Powered by ViewVC 1.1.26