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

Legend:
Removed from v.1384  
changed lines
  Added in v.1887

  ViewVC Help
Powered by ViewVC 1.1.26