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

  ViewVC Help
Powered by ViewVC 1.1.26