/[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

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

Legend:
Removed from v.1387  
changed lines
  Added in v.2052

  ViewVC Help
Powered by ViewVC 1.1.26