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

  ViewVC Help
Powered by ViewVC 1.1.26