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

Legend:
Removed from v.1028  
changed lines
  Added in v.2311

  ViewVC Help
Powered by ViewVC 1.1.26