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

  ViewVC Help
Powered by ViewVC 1.1.26