/[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 1919 by phornby, Thu Oct 23 10:03:10 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;    FILE *fileHandle_p = NULL;
42    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    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    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,"Finley_Mesh_read: Opening file %s for reading failed.",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;
279             &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);    int scan_ret;
280    } else if (3 == numDim) {  
281        for (i0 = 0; i0 < numNodes; i0++)    Finley_resetError();
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],    if (mpi_info->rank == 0) {
284             &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],       /* get file handle */
285             &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],       fileHandle_p = fopen(fname, "r");
286             &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);       if (fileHandle_p==NULL) {
287    } /* if else else */         sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
288           Finley_setError(IO_ERROR,error_msg);
289    /* get the element type */         Paso_MPIInfo_free( mpi_info );
290           return NULL;
291    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);       }
292    typeID=Finley_RefElement_getTypeId(element_type);  
293    if (typeID==NoType) {       /* read header */
294      sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
295      Finley_setError(VALUE_ERROR,error_msg);       scan_ret = fscanf(fileHandle_p, frm, name);
296      return NULL;       FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
297    
298         /* get the number of nodes */
299         scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
300         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
301    }    }
302    /* read the elements */  
303  #ifndef PASO_MPI  #ifdef PASO_MPI
304    mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);    /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
305  #else    if (mpi_info->size > 1) {
306    /* TODO */      int temp1[3], error_code;
307  #endif      temp1[0] = numDim;
308    Finley_ElementFile_allocTable(mesh_p->Elements, numEle);      temp1[1] = numNodes;
309    mesh_p->Elements->minColor=0;      temp1[2] = strlen(name) + 1;
310    mesh_p->Elements->maxColor=numEle-1;      error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
311    for (i0 = 0; i0 < numEle; i0++) {      if (error_code != MPI_SUCCESS) {
312      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");
313      mesh_p->Elements->Color[i0]=i0;        return NULL;
314      for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {      }
315           fscanf(fileHandle_p, " %d",      numDim = temp1[0];
316              &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);      numNodes = temp1[1];
317      }   /* for i1 */      error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
318      fscanf(fileHandle_p, "\n");      if (error_code != MPI_SUCCESS) {
319    } /* for i0 */        Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
320          return NULL;
321    /* 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;  
322    }    }
 #ifndef PASO_MPI  
   mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);  
 #else  
   /* TODO */  
323  #endif  #endif
324    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
325    mesh_p->FaceElements->minColor=0;       /* allocate mesh */
326    mesh_p->FaceElements->maxColor=numEle-1;       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
327    for (i0 = 0; i0 < numEle; i0++) {       if (Finley_noError()) {
328      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 */
329      mesh_p->FaceElements->Color[i0]=i0;      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
330      for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
331           fscanf(fileHandle_p, " %d",      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
332              &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);  
333      }   /* for i1 */      /*
334      fscanf(fileHandle_p, "\n");        Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
335    } /* for i0 */        It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
336          First chunk sent to CPU 1, second to CPU 2, ...
337    /* get the Contact face element */        Last chunk stays on CPU 0 (the master)
338    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
339    contactTypeID=Finley_RefElement_getTypeId(element_type);      */
340    if (contactTypeID==NoType) {  
341      sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);      if (mpi_info->rank == 0) {  /* Master */
342      Finley_setError(VALUE_ERROR,error_msg);        for (;;) {            /* Infinite loop */
343      return NULL;          for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
344    }          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
345  #ifndef PASO_MPI          chunkNodes = 0;
346    mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);          for (i1=0; i1<chunkSize; i1++) {
347  #else            if (totalNodes >= numNodes) break;    /* End of inner loop */
348    /* TODO */                if (1 == numDim) {
349            scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
350              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
351              &tempCoords[i1*numDim+0]);
352            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
353              }
354                  if (2 == numDim) {
355            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
356              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
357              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
358            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
359              }
360                  if (3 == numDim) {
361            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
362              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
363              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
364            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
365              }
366              totalNodes++; /* When do we quit the infinite loop? */
367              chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
368            }
369            if (chunkNodes > chunkSize) {
370                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
371                  return NULL;
372            }
373    #ifdef PASO_MPI
374            /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
375            if (nextCPU < mpi_info->size) {
376                  int mpi_error;
377              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
378              mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
379              if ( mpi_error != MPI_SUCCESS ) {
380                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
381                    return NULL;
382              }
383              mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
384              if ( mpi_error != MPI_SUCCESS ) {
385                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
386                    return NULL;
387              }
388            }
389  #endif  #endif
390    Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);          nextCPU++;
391    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 */
392    mesh_p->ContactElements->maxColor=numEle-1;          if (nextCPU > mpi_info->size) break; /* End infinite loop */
393    for (i0 = 0; i0 < numEle; i0++) {        } /* Infinite loop */
394      fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);      }   /* End master */
395      mesh_p->ContactElements->Color[i0]=i0;      else {  /* Worker */
396      for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {  #ifdef PASO_MPI
397          fscanf(fileHandle_p, " %d",        /* Each worker receives two messages */
398             &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);        MPI_Status status;
399      }   /* for i1 */            int mpi_error;
400      fscanf(fileHandle_p, "\n");        mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
401    } /* for i0 */        if ( mpi_error != MPI_SUCCESS ) {
402                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
403    /* get the nodal element */              return NULL;
404    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        }
405    pointTypeID=Finley_RefElement_getTypeId(element_type);        mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
406    if (pointTypeID==NoType) {        if ( mpi_error != MPI_SUCCESS ) {
407      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");
408      Finley_setError(VALUE_ERROR,error_msg);              return NULL;
409      return NULL;        }
410    }        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 */  
411  #endif  #endif
412    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 */  
   
   
   /* close file */  
413    
414    fclose(fileHandle_p);      /* Copy node data from tempMem to mesh_p */
415            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
416            if (Finley_noError()) {
417          for (i0=0; i0<chunkNodes; i0++) {
418            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
419            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
420            mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
421            for (i1=0; i1<numDim; i1++) {
422              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
423            }
424              }
425            }
426    
427        TMPMEMFREE(tempInts);
428        TMPMEMFREE(tempCoords);
429    
430            /* read elements */
431    
432        /* Read the element typeID */
433            if (Finley_noError()) {
434          if (mpi_info->rank == 0) {
435                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
436            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
437                typeID=Finley_RefElement_getTypeId(element_type);
438          }
439    #ifdef PASO_MPI
440          if (mpi_info->size > 1) {
441            int temp1[2], mpi_error;
442            temp1[0] = (int) typeID;
443            temp1[1] = numEle;
444            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
445            if (mpi_error != MPI_SUCCESS) {
446              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
447              return NULL;
448            }
449            typeID = (ElementTypeId) temp1[0];
450            numEle = temp1[1];
451          }
452    #endif
453              if (typeID==NoType) {
454                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
455                Finley_setError(VALUE_ERROR, error_msg);
456              }
457        }
458    
459          /* Allocate the ElementFile */
460          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
461          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
462    
463          /* Read the element data */
464          if (Finley_noError()) {
465        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
466        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
467        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
468        if (mpi_info->rank == 0) {  /* Master */
469          for (;;) {            /* Infinite loop */
470            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
471            chunkEle = 0;
472            for (i0=0; i0<chunkSize; i0++) {
473              if (totalEle >= numEle) break; /* End inner loop */
474              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
475              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
476              for (i1 = 0; i1 < numNodes; i1++) {
477            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
478                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
479              }
480              scan_ret = fscanf(fileHandle_p, "\n");
481              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
482              totalEle++;
483              chunkEle++;
484            }
485    #ifdef PASO_MPI
486            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
487            if (nextCPU < mpi_info->size) {
488                  int mpi_error;
489              tempInts[chunkSize*(2+numNodes)] = chunkEle;
490              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
491              if ( mpi_error != MPI_SUCCESS ) {
492                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
493                    return NULL;
494              }
495            }
496    #endif
497            nextCPU++;
498            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
499            if (nextCPU > mpi_info->size) break; /* End infinite loop */
500          } /* Infinite loop */
501        }   /* End master */
502        else {  /* Worker */
503    #ifdef PASO_MPI
504          /* Each worker receives one message */
505          MPI_Status status;
506          int mpi_error;
507          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
508          if ( mpi_error != MPI_SUCCESS ) {
509                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
510                return NULL;
511          }
512          chunkEle = tempInts[chunkSize*(2+numNodes)];
513    #endif
514        }   /* Worker */
515    
516    /*   resolve id's : */      /* Copy Element data from tempInts to mesh_p */
517        Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
518            mesh_p->Elements->minColor=0;
519            mesh_p->Elements->maxColor=chunkEle-1;
520            if (Finley_noError()) {
521              #pragma omp parallel for private (i0, i1)
522          for (i0=0; i0<chunkEle; i0++) {
523            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
524            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
525                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
526                mesh_p->Elements->Color[i0] = i0;
527            for (i1 = 0; i1 < numNodes; i1++) {
528              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
529            }
530              }
531            }
532    
533        TMPMEMFREE(tempInts);
534          } /* end of Read the element data */
535    
536            /* read face elements */
537    
538        /* Read the element typeID */
539            if (Finley_noError()) {
540          if (mpi_info->rank == 0) {
541                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
542            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
543                typeID=Finley_RefElement_getTypeId(element_type);
544          }
545    #ifdef PASO_MPI
546          if (mpi_info->size > 1) {
547            int temp1[2], mpi_error;
548            temp1[0] = (int) typeID;
549            temp1[1] = numEle;
550            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
551            if (mpi_error != MPI_SUCCESS) {
552              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
553              return NULL;
554            }
555            typeID = (ElementTypeId) temp1[0];
556            numEle = temp1[1];
557          }
558    #endif
559              if (typeID==NoType) {
560                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
561                Finley_setError(VALUE_ERROR, error_msg);
562              }
563        }
564    
565          /* Allocate the ElementFile */
566          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
567          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
568    
569          /* Read the face element data */
570          if (Finley_noError()) {
571        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
572        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
573        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
574        if (mpi_info->rank == 0) {  /* Master */
575          for (;;) {            /* Infinite loop */
576            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
577            chunkEle = 0;
578            for (i0=0; i0<chunkSize; i0++) {
579              if (totalEle >= numEle) break; /* End inner loop */
580              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
581              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
582              for (i1 = 0; i1 < numNodes; i1++) {
583            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
584                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
585              }
586              scan_ret = fscanf(fileHandle_p, "\n");
587              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
588              totalEle++;
589              chunkEle++;
590            }
591    #ifdef PASO_MPI
592            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
593            if (nextCPU < mpi_info->size) {
594                  int mpi_error;
595              tempInts[chunkSize*(2+numNodes)] = chunkEle;
596              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
597              if ( mpi_error != MPI_SUCCESS ) {
598                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
599                    return NULL;
600              }
601            }
602    #endif
603            nextCPU++;
604            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
605            if (nextCPU > mpi_info->size) break; /* End infinite loop */
606          } /* Infinite loop */
607        }   /* End master */
608        else {  /* Worker */
609    #ifdef PASO_MPI
610          /* Each worker receives one message */
611          MPI_Status status;
612          int mpi_error;
613          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
614          if ( mpi_error != MPI_SUCCESS ) {
615                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
616                return NULL;
617          }
618          chunkEle = tempInts[chunkSize*(2+numNodes)];
619    #endif
620        }   /* Worker */
621    
622    Finley_Mesh_resolveNodeIds(mesh_p);      /* Copy Element data from tempInts to mesh_p */
623        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
624            mesh_p->FaceElements->minColor=0;
625            mesh_p->FaceElements->maxColor=chunkEle-1;
626            if (Finley_noError()) {
627              #pragma omp parallel for private (i0, i1)
628          for (i0=0; i0<chunkEle; i0++) {
629            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
630            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
631                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
632                mesh_p->FaceElements->Color[i0] = i0;
633            for (i1 = 0; i1 < numNodes; i1++) {
634              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
635            }
636              }
637            }
638    
639        TMPMEMFREE(tempInts);
640          } /* end of Read the face element data */
641    
642            /* read contact elements */
643    
644        /* Read the element typeID */
645            if (Finley_noError()) {
646          if (mpi_info->rank == 0) {
647                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
648            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
649                typeID=Finley_RefElement_getTypeId(element_type);
650          }
651    #ifdef PASO_MPI
652          if (mpi_info->size > 1) {
653            int temp1[2], mpi_error;
654            temp1[0] = (int) typeID;
655            temp1[1] = numEle;
656            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
657            if (mpi_error != MPI_SUCCESS) {
658              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
659              return NULL;
660            }
661            typeID = (ElementTypeId) temp1[0];
662            numEle = temp1[1];
663          }
664    #endif
665              if (typeID==NoType) {
666                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
667                Finley_setError(VALUE_ERROR, error_msg);
668              }
669        }
670    
671          /* Allocate the ElementFile */
672          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
673          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
674    
675          /* Read the contact element data */
676          if (Finley_noError()) {
677        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
678        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
679        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
680        if (mpi_info->rank == 0) {  /* Master */
681          for (;;) {            /* Infinite loop */
682            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
683            chunkEle = 0;
684            for (i0=0; i0<chunkSize; i0++) {
685              if (totalEle >= numEle) break; /* End inner loop */
686              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
687              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
688              for (i1 = 0; i1 < numNodes; i1++) {
689            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
690                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
691              }
692              scan_ret = fscanf(fileHandle_p, "\n");
693              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
694              totalEle++;
695              chunkEle++;
696            }
697    #ifdef PASO_MPI
698            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
699            if (nextCPU < mpi_info->size) {
700                  int mpi_error;
701              tempInts[chunkSize*(2+numNodes)] = chunkEle;
702              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
703              if ( mpi_error != MPI_SUCCESS ) {
704                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
705                    return NULL;
706              }
707            }
708    #endif
709            nextCPU++;
710            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
711            if (nextCPU > mpi_info->size) break; /* End infinite loop */
712          } /* Infinite loop */
713        }   /* End master */
714        else {  /* Worker */
715    #ifdef PASO_MPI
716          /* Each worker receives one message */
717          MPI_Status status;
718          int mpi_error;
719          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
720          if ( mpi_error != MPI_SUCCESS ) {
721                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
722                return NULL;
723          }
724          chunkEle = tempInts[chunkSize*(2+numNodes)];
725    #endif
726        }   /* Worker */
727    
728    /* rearrange elements: */      /* Copy Element data from tempInts to mesh_p */
729        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
730            mesh_p->ContactElements->minColor=0;
731            mesh_p->ContactElements->maxColor=chunkEle-1;
732            if (Finley_noError()) {
733              #pragma omp parallel for private (i0, i1)
734          for (i0=0; i0<chunkEle; i0++) {
735            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
736            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
737                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
738                mesh_p->ContactElements->Color[i0] = i0;
739            for (i1 = 0; i1 < numNodes; i1++) {
740              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
741            }
742              }
743            }
744    
745        TMPMEMFREE(tempInts);
746          } /* end of Read the contact element data */
747    
748            /* read nodal elements */
749    
750        /* Read the element typeID */
751            if (Finley_noError()) {
752          if (mpi_info->rank == 0) {
753                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
754            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
755                typeID=Finley_RefElement_getTypeId(element_type);
756          }
757    #ifdef PASO_MPI
758          if (mpi_info->size > 1) {
759            int temp1[2], mpi_error;
760            temp1[0] = (int) typeID;
761            temp1[1] = numEle;
762            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
763            if (mpi_error != MPI_SUCCESS) {
764              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
765              return NULL;
766            }
767            typeID = (ElementTypeId) temp1[0];
768            numEle = temp1[1];
769          }
770    #endif
771              if (typeID==NoType) {
772                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
773                Finley_setError(VALUE_ERROR, error_msg);
774              }
775        }
776    
777          /* Allocate the ElementFile */
778          mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
779          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
780    
781          /* Read the nodal element data */
782          if (Finley_noError()) {
783        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
784        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
785        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
786        if (mpi_info->rank == 0) {  /* Master */
787          for (;;) {            /* Infinite loop */
788            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
789            chunkEle = 0;
790            for (i0=0; i0<chunkSize; i0++) {
791              if (totalEle >= numEle) break; /* End inner loop */
792              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
793              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
794              for (i1 = 0; i1 < numNodes; i1++) {
795            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
796                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
797              }
798              scan_ret = fscanf(fileHandle_p, "\n");
799              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
800              totalEle++;
801              chunkEle++;
802            }
803    #ifdef PASO_MPI
804            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
805            if (nextCPU < mpi_info->size) {
806                  int mpi_error;
807              tempInts[chunkSize*(2+numNodes)] = chunkEle;
808              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
809              if ( mpi_error != MPI_SUCCESS ) {
810                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
811                    return NULL;
812              }
813            }
814    #endif
815            nextCPU++;
816            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
817            if (nextCPU > mpi_info->size) break; /* End infinite loop */
818          } /* Infinite loop */
819        }   /* End master */
820        else {  /* Worker */
821    #ifdef PASO_MPI
822          /* Each worker receives one message */
823          MPI_Status status;
824          int mpi_error;
825          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
826          if ( mpi_error != MPI_SUCCESS ) {
827                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
828                return NULL;
829          }
830          chunkEle = tempInts[chunkSize*(2+numNodes)];
831    #endif
832        }   /* Worker */
833    
834    Finley_Mesh_prepare(mesh_p);      /* Copy Element data from tempInts to mesh_p */
835        Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
836            mesh_p->Points->minColor=0;
837            mesh_p->Points->maxColor=chunkEle-1;
838            if (Finley_noError()) {
839              #pragma omp parallel for private (i0, i1)
840          for (i0=0; i0<chunkEle; i0++) {
841            mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
842            mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
843                mesh_p->Points->Owner[i0]  =mpi_info->rank;
844                mesh_p->Points->Color[i0] = i0;
845            for (i1 = 0; i1 < numNodes; i1++) {
846              mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
847            }
848              }
849            }
850    
851        TMPMEMFREE(tempInts);
852          } /* end of Read the nodal element data */
853    
854          /* get the name tags */
855          if (Finley_noError()) {
856            char *remainder, *ptr;
857            size_t len;
858            int tag_key;
859    
860    
861    
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    
869    #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
870    
871          remainder = NULL;
872              len=0;
873          while (1)
874              {
875                 size_t MALLOC_CHUNK = 1024;
876                 size_t buff_size = 0;
877                 int ch;
878    
879                 ch = fgetc(fileHandle_p);
880                 if( ch == '\r' )
881                 {
882                    continue;
883                 }
884                 if( len+1 > buff_size )
885                 {
886                    TMPMEMREALLOC(remainder,remainder,buff_size+MALLOC_CHUNK,char);
887                 }
888                 if( ch == EOF )
889                 {
890                    /* hit EOF */
891                    remainder[len] = (char)0;
892                    break;
893                 }
894                 remainder[len] = (char)ch;
895                 len++;
896              }
897    #else
898          /* Read rest of file in one chunk, after using seek to find length */
899              {
900                 long cur_pos, end_pos;
901    
902                 cur_pos = ftell(fileHandle_p);
903                 fseek(fileHandle_p, 0L, SEEK_END);
904                 end_pos = ftell(fileHandle_p);
905                 fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
906                 remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
907                 if (! feof(fileHandle_p))
908                 {
909                    scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
910                                     sizeof(char), fileHandle_p);
911    
912                    FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
913                    remainder[end_pos-cur_pos] = 0;
914                }
915              }
916    #endif
917          len = strlen(remainder);
918          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
919          len = strlen(remainder);
920              TMPMEMREALLOC(remainder,remainder,len+1,char);
921            }
922    #ifdef PASO_MPI
923    
924            if (MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm) != MPI_SUCCESS)
925            {
926               Finley_setError(PASO_MPI_ERROR,
927                               "Finley_Mesh_read: broadcast of tag len failed");
928               return NULL;
929            }
930        if (mpi_info->rank != 0) {
931          remainder = TMPMEMALLOC(len+1, char);
932          remainder[0] = 0;
933        }
934    
935            if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
936                MPI_SUCCESS)
937            {
938               Finley_setError(PASO_MPI_ERROR,
939                               "Finley_Mesh_read: broadcast of tags failed");
940               return NULL;
941            }
942    #endif
943        if (remainder[0]) {
944              ptr = remainder;
945              do {
946                sscanf(ptr, "%s %d\n", name, &tag_key);
947                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
948                ptr++;
949              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
950              TMPMEMFREE(remainder);
951        }
952          }
953    
954         }
955    
956         /* close file */
957         if (mpi_info->rank == 0) fclose(fileHandle_p);
958    
959         /*   resolve id's : */
960         /* rearrange elements: */
961    
962         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
963         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
964    
965         /* that's it */
966         #ifdef Finley_TRACE
967         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
968         #endif
969    
970    /* that's it */    /* that's it */
971    #ifdef Finley_TRACE    if (! Finley_noError()) {
972    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);  
973    }    }
974      Paso_MPIInfo_free( mpi_info );
975    return mesh_p;    return mesh_p;
976  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26