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

revision 1713 by ksteube, Wed Aug 20 07:03:45 2008 UTC revision 2059 by gross, Mon Nov 17 22:57:50 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 19  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21    #include <ctype.h>
22  #include "Mesh.h"  #include "Mesh.h"
23    
24    #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Finley_setError(IO_ERROR,"scan error while reading finley file"); return NULL;} }
25    
26  /**************************************************************/  /**************************************************************/
27    
28  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Finley file of name fname */
# Line 35  Finley_Mesh* Finley_Mesh_read(char* fnam Line 37  Finley_Mesh* Finley_Mesh_read(char* fnam
37    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
38    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
39    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
40      #ifdef Finley_TRACE
41    double time0=Finley_timer();    double time0=Finley_timer();
42      #endif
43    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
44    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID=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 239  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 275  Finley_Mesh* Finley_Mesh_read_MPI(char*
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;    ElementTypeId typeID=NoType;
283    ElementTypeId faceTypeID, contactTypeID, pointTypeID;    int scan_ret;
   index_t tag_key;  
284    
285    Finley_resetError();    Finley_resetError();
286    
# Line 259  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 296  Finley_Mesh* Finley_Mesh_read_MPI(char*
296    
297       /* read header */       /* read header */
298       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
299       fscanf(fileHandle_p, frm, name);       scan_ret = fscanf(fileHandle_p, frm, name);
300         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
301    
302       /* get the number of nodes */       /* get the number of nodes */
303       fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);       scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
304         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
305    }    }
306    
307  #ifdef PASO_MPI  #ifdef PASO_MPI
# Line 290  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 329  Finley_Mesh* Finley_Mesh_read_MPI(char*
329       /* allocate mesh */       /* allocate mesh */
330       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
331       if (Finley_noError()) {       if (Finley_noError()) {
332      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large for that */      /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
333      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
334      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */      int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
335      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */      double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
# Line 309  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 348  Finley_Mesh* Finley_Mesh_read_MPI(char*
348          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;          for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
349          chunkNodes = 0;          chunkNodes = 0;
350          for (i1=0; i1<chunkSize; i1++) {          for (i1=0; i1<chunkSize; i1++) {
351            if (totalNodes >= numNodes) break;    /* Maybe end the infinite loop */            if (totalNodes >= numNodes) break;    /* End of inner loop */
352                if (1 == numDim)                if (1 == numDim) {
353          fscanf(fileHandle_p, "%d %d %d %le\n",          scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
354            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
355            &tempCoords[i1*numDim+0]);            &tempCoords[i1*numDim+0]);
356                if (2 == numDim)          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
357          fscanf(fileHandle_p, "%d %d %d %le %le\n",            }
358                  if (2 == numDim) {
359            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
360            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
361            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
362                if (3 == numDim)          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
363          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",            }
364                  if (3 == numDim) {
365            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
366            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],            &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
367            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
368            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
369              }
370            totalNodes++; /* When do we quit the infinite loop? */            totalNodes++; /* When do we quit the infinite loop? */
371            chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */            chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
372          }          }
# Line 344  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 389  Finley_Mesh* Finley_Mesh_read_MPI(char*
389                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
390                  return NULL;                  return NULL;
391            }            }
           nextCPU++;  
392          }          }
393  #endif  #endif
394          if (totalNodes >= numNodes) break;  /* Maybe end the infinite loop */          nextCPU++;
395            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
396            if (nextCPU > mpi_info->size) break; /* End infinite loop */
397        } /* Infinite loop */        } /* Infinite loop */
398      }   /* End master */      }   /* End master */
399      else {  /* Worker */      else {  /* Worker */
# Line 369  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 415  Finley_Mesh* Finley_Mesh_read_MPI(char*
415  #endif  #endif
416      }   /* Worker */      }   /* Worker */
417    
 #if 1  
         /* Display the temp mem for debugging */  
         printf("ksteube Nodes tempInts\n");  
         for (i0=0; i0<chunkSize*3; i0++) {  
           printf(" %2d", tempInts[i0]);  
           if (i0%chunkSize==chunkSize-1) printf("\n");  
         }  
         printf("ksteube tempCoords:\n");  
         for (i0=0; i0<chunkSize*numDim; i0++) {  
           printf(" %20.15e", tempCoords[i0]);  
           if (i0%numDim==numDim-1) printf("\n");  
         }  
 #endif  
   
         printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);  
   
418      /* Copy node data from tempMem to mesh_p */      /* Copy node data from tempMem to mesh_p */
419          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
420          if (Finley_noError()) {          if (Finley_noError()) {
# Line 406  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 436  Finley_Mesh* Finley_Mesh_read_MPI(char*
436      /* Read the element typeID */      /* Read the element typeID */
437          if (Finley_noError()) {          if (Finley_noError()) {
438        if (mpi_info->rank == 0) {        if (mpi_info->rank == 0) {
439              fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
440            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
441              typeID=Finley_RefElement_getTypeId(element_type);              typeID=Finley_RefElement_getTypeId(element_type);
442        }        }
443  #ifdef PASO_MPI  #ifdef PASO_MPI
444        if (mpi_info->size > 1) {        if (mpi_info->size > 1) {
445          int temp1[3], mpi_error;          int temp1[2], mpi_error;
446          temp1[0] = (int) typeID;          temp1[0] = (int) typeID;
447          temp1[1] = numEle;          temp1[1] = numEle;
448          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);          mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
# Line 429  Finley_Mesh* Finley_Mesh_read_MPI(char* Line 460  Finley_Mesh* Finley_Mesh_read_MPI(char*
460            }            }
461      }      }
462    
463        /* Read the element data */        /* Allocate the ElementFile */
464        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);        mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
465        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */        numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
466    
467          /* Read the element data */
468        if (Finley_noError()) {        if (Finley_noError()) {
469      int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
470      int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
471      if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
472      if (mpi_info->rank == 0) {  /* Master */      if (mpi_info->rank == 0) {  /* Master */
473        for (;;) {            /* Infinite loop */        for (;;) {            /* Infinite loop */
474            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
475          chunkEle = 0;          chunkEle = 0;
         for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;  
476          for (i0=0; i0<chunkSize; i0++) {          for (i0=0; i0<chunkSize; i0++) {
477            if (totalEle >= numEle) break; /* End infinite loop */            if (totalEle >= numEle) break; /* End inner loop */
478            fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);            scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
479            for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
480            fscanf(fileHandle_p, "\n");            for (i1 = 0; i1 < numNodes; i1++) {
481            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
482                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
483              }
484              scan_ret = fscanf(fileHandle_p, "\n");
485              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
486            totalEle++;            totalEle++;
487            chunkEle++;            chunkEle++;
488          }          }
         /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */  
         if (nextCPU < mpi_info->size) {  
489  #ifdef PASO_MPI  #ifdef PASO_MPI
490            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
491            if (nextCPU < mpi_info->size) {
492                int mpi_error;                int mpi_error;
493              tempInts[chunkSize*(2+numNodes)] = chunkEle;
494            tempInts[numEle*(2+numNodes)] = chunkEle;            mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);  
           mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);  
495            if ( mpi_error != MPI_SUCCESS ) {            if ( mpi_error != MPI_SUCCESS ) {
496                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
497                  return NULL;                  return NULL;
498            }            }
 #endif  
           nextCPU++;  
499          }          }
500          if (totalEle >= numEle) break; /* End infinite loop */  #endif
501            nextCPU++;
502            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
503            if (nextCPU > mpi_info->size) break; /* End infinite loop */
504        } /* Infinite loop */        } /* Infinite loop */
505      }   /* End master */      }   /* End master */
506      else {  /* Worker */      else {  /* Worker */
507  #ifdef PASO_MPI  #ifdef PASO_MPI
508        /* Each worker receives two messages */        /* Each worker receives one message */
509        MPI_Status status;        MPI_Status status;
510        int mpi_error;        int mpi_error;
511  printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);        mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
       mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);  
512        if ( mpi_error != MPI_SUCCESS ) {        if ( mpi_error != MPI_SUCCESS ) {
513              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
514              return NULL;              return NULL;
515        }        }
516        chunkEle = tempInts[numEle*(2+numNodes)];        chunkEle = tempInts[chunkSize*(2+numNodes)];
517  #endif  #endif
518      }   /* Worker */      }   /* Worker */
 #if 0  
     /* Display the temp mem for debugging */  
     printf("ksteube Elements tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);  
     for (i0=0; i0<numEle*(numNodes+2); i0++) {  
       printf(" %2d", tempInts[i0]);  
       if (i0%(numNodes+2)==numNodes+2-1) printf("\n");  
     }  
 #endif  
519    
520      /* Copy Element data from tempInts to mesh_p */      /* Copy Element data from tempInts to mesh_p */
521      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
# Line 499  printf("ksteube CPU=%d/%d recv on %d\n", Line 526  printf("ksteube CPU=%d/%d recv on %d\n",
526        for (i0=0; i0<chunkEle; i0++) {        for (i0=0; i0<chunkEle; i0++) {
527          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];          mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
528          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];          mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
529                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
530                mesh_p->Elements->Color[i0] = i0;
531          for (i1 = 0; i1 < numNodes; i1++) {          for (i1 = 0; i1 < numNodes; i1++) {
532            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];            mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
533          }          }
# Line 506  printf("ksteube CPU=%d/%d recv on %d\n", Line 535  printf("ksteube CPU=%d/%d recv on %d\n",
535          }          }
536    
537      TMPMEMFREE(tempInts);      TMPMEMFREE(tempInts);
538        }        } /* end of Read the element data */
539    
540            /* read face elements */
541    
542  #if 0 /* this is the original code for reading elements */      /* Read the element typeID */
         /* read elements */  
543          if (Finley_noError()) {          if (Finley_noError()) {
544          if (mpi_info->rank == 0) {
545                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
546            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
547                typeID=Finley_RefElement_getTypeId(element_type);
548          }
549    #ifdef PASO_MPI
550          if (mpi_info->size > 1) {
551            int temp1[2], mpi_error;
552            temp1[0] = (int) typeID;
553            temp1[1] = numEle;
554            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
555            if (mpi_error != MPI_SUCCESS) {
556              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
557              return NULL;
558            }
559            typeID = (ElementTypeId) temp1[0];
560            numEle = temp1[1];
561          }
562    #endif
563              if (typeID==NoType) {
564                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
565                Finley_setError(VALUE_ERROR, error_msg);
566              }
567        }
568    
569             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        /* Allocate the ElementFile */
570             typeID=Finley_RefElement_getTypeId(element_type);        mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
571             if (typeID==NoType) {        numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
572               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
573               Finley_setError(VALUE_ERROR,error_msg);        /* Read the face element data */
574             } else {        if (Finley_noError()) {
575               /* read the elements */      int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
576               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);      int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
577               if (Finley_noError()) {      /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
578                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);      if (mpi_info->rank == 0) {  /* Master */
579                   mesh_p->Elements->minColor=0;        for (;;) {            /* Infinite loop */
580                   mesh_p->Elements->maxColor=numEle-1;          for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
581                   if (Finley_noError()) {          chunkEle = 0;
582                      for (i0 = 0; i0 < numEle; i0++) {          for (i0=0; i0<chunkSize; i0++) {
583                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);            if (totalEle >= numEle) break; /* End inner loop */
584                        mesh_p->Elements->Color[i0]=i0;            scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
585                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
586                             fscanf(fileHandle_p, " %d",            for (i1 = 0; i1 < numNodes; i1++) {
587                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
588                        } /* for i1 */              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
589                        fscanf(fileHandle_p, "\n");            }
590                      } /* for i0 */            scan_ret = fscanf(fileHandle_p, "\n");
591                   }            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
592               }            totalEle++;
593              chunkEle++;
594            }
595    #ifdef PASO_MPI
596            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
597            if (nextCPU < mpi_info->size) {
598                  int mpi_error;
599              tempInts[chunkSize*(2+numNodes)] = chunkEle;
600              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
601              if ( mpi_error != MPI_SUCCESS ) {
602                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
603                    return NULL;
604              }
605            }
606    #endif
607            nextCPU++;
608            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
609            if (nextCPU > mpi_info->size) break; /* End infinite loop */
610          } /* Infinite loop */
611        }   /* End master */
612        else {  /* Worker */
613    #ifdef PASO_MPI
614          /* Each worker receives one message */
615          MPI_Status status;
616          int mpi_error;
617          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
618          if ( mpi_error != MPI_SUCCESS ) {
619                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
620                return NULL;
621          }
622          chunkEle = tempInts[chunkSize*(2+numNodes)];
623    #endif
624        }   /* Worker */
625    
626        /* Copy Element data from tempInts to mesh_p */
627        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
628            mesh_p->FaceElements->minColor=0;
629            mesh_p->FaceElements->maxColor=chunkEle-1;
630            if (Finley_noError()) {
631              #pragma omp parallel for private (i0, i1)
632          for (i0=0; i0<chunkEle; i0++) {
633            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
634            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
635                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
636                mesh_p->FaceElements->Color[i0] = i0;
637            for (i1 = 0; i1 < numNodes; i1++) {
638              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
639            }
640            }            }
641          }          }
642    
643        TMPMEMFREE(tempInts);
644          } /* end of Read the face element data */
645    
646            /* read contact elements */
647    
648        /* Read the element typeID */
649            if (Finley_noError()) {
650          if (mpi_info->rank == 0) {
651                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
652            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
653                typeID=Finley_RefElement_getTypeId(element_type);
654          }
655    #ifdef PASO_MPI
656          if (mpi_info->size > 1) {
657            int temp1[2], mpi_error;
658            temp1[0] = (int) typeID;
659            temp1[1] = numEle;
660            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
661            if (mpi_error != MPI_SUCCESS) {
662              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
663              return NULL;
664            }
665            typeID = (ElementTypeId) temp1[0];
666            numEle = temp1[1];
667          }
668  #endif  #endif
669              if (typeID==NoType) {
670                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
671                Finley_setError(VALUE_ERROR, error_msg);
672              }
673        }
674    
675          /* Allocate the ElementFile */
676          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
677          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
678    
679  printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);        /* Read the contact element data */
680          if (Finley_noError()) {
681        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
682        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
683        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
684        if (mpi_info->rank == 0) {  /* Master */
685          for (;;) {            /* Infinite loop */
686            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
687            chunkEle = 0;
688            for (i0=0; i0<chunkSize; i0++) {
689              if (totalEle >= numEle) break; /* End inner loop */
690              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
691              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
692              for (i1 = 0; i1 < numNodes; i1++) {
693            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
694                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
695              }
696              scan_ret = fscanf(fileHandle_p, "\n");
697              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
698              totalEle++;
699              chunkEle++;
700            }
701    #ifdef PASO_MPI
702            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
703            if (nextCPU < mpi_info->size) {
704                  int mpi_error;
705              tempInts[chunkSize*(2+numNodes)] = chunkEle;
706              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
707              if ( mpi_error != MPI_SUCCESS ) {
708                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
709                    return NULL;
710              }
711            }
712    #endif
713            nextCPU++;
714            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
715            if (nextCPU > mpi_info->size) break; /* End infinite loop */
716          } /* Infinite loop */
717        }   /* End master */
718        else {  /* Worker */
719    #ifdef PASO_MPI
720          /* Each worker receives one message */
721          MPI_Status status;
722          int mpi_error;
723          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
724          if ( mpi_error != MPI_SUCCESS ) {
725                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
726                return NULL;
727          }
728          chunkEle = tempInts[chunkSize*(2+numNodes)];
729    #endif
730        }   /* Worker */
731    
732  #if 1      /* Copy Element data from tempInts to mesh_p */
733        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
734            mesh_p->ContactElements->minColor=0;
735            mesh_p->ContactElements->maxColor=chunkEle-1;
736            if (Finley_noError()) {
737              #pragma omp parallel for private (i0, i1)
738          for (i0=0; i0<chunkEle; i0++) {
739            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
740            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
741                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
742                mesh_p->ContactElements->Color[i0] = i0;
743            for (i1 = 0; i1 < numNodes; i1++) {
744              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
745            }
746              }
747            }
748    
749    /* Define other structures to keep mesh_write from crashing */      TMPMEMFREE(tempInts);
750    /* Change the typeid from NoType later */        } /* end of Read the contact element data */
751    
752    mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);          /* read nodal elements */
   Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);  
753    
754    mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);      /* Read the element typeID */
755    Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);          if (Finley_noError()) {
756          if (mpi_info->rank == 0) {
757                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
758            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
759                typeID=Finley_RefElement_getTypeId(element_type);
760          }
761    #ifdef PASO_MPI
762          if (mpi_info->size > 1) {
763            int temp1[2], mpi_error;
764            temp1[0] = (int) typeID;
765            temp1[1] = numEle;
766            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
767            if (mpi_error != MPI_SUCCESS) {
768              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
769              return NULL;
770            }
771            typeID = (ElementTypeId) temp1[0];
772            numEle = temp1[1];
773          }
774    #endif
775              if (typeID==NoType) {
776                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
777                Finley_setError(VALUE_ERROR, error_msg);
778              }
779        }
780    
781    mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);        /* Allocate the ElementFile */
782    Finley_ElementFile_allocTable(mesh_p->Points, 0);        mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
783          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
784    
785          /* Read the nodal element data */
786          if (Finley_noError()) {
787        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
788        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
789        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
790        if (mpi_info->rank == 0) {  /* Master */
791          for (;;) {            /* Infinite loop */
792            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
793            chunkEle = 0;
794            for (i0=0; i0<chunkSize; i0++) {
795              if (totalEle >= numEle) break; /* End inner loop */
796              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
797              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
798              for (i1 = 0; i1 < numNodes; i1++) {
799            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
800                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
801              }
802              scan_ret = fscanf(fileHandle_p, "\n");
803              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
804              totalEle++;
805              chunkEle++;
806            }
807    #ifdef PASO_MPI
808            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
809            if (nextCPU < mpi_info->size) {
810                  int mpi_error;
811              tempInts[chunkSize*(2+numNodes)] = chunkEle;
812              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
813              if ( mpi_error != MPI_SUCCESS ) {
814                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
815                    return NULL;
816              }
817            }
818    #endif
819            nextCPU++;
820            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
821            if (nextCPU > mpi_info->size) break; /* End infinite loop */
822          } /* Infinite loop */
823        }   /* End master */
824        else {  /* Worker */
825    #ifdef PASO_MPI
826          /* Each worker receives one message */
827          MPI_Status status;
828          int mpi_error;
829          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
830          if ( mpi_error != MPI_SUCCESS ) {
831                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
832                return NULL;
833          }
834          chunkEle = tempInts[chunkSize*(2+numNodes)];
835  #endif  #endif
836        }   /* Worker */
837    
838  #if 0 /* comment out the rest of the un-implemented crap for now */      /* Copy Element data from tempInts to mesh_p */
839          /* get the face elements */      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
840          if (Finley_noError()) {          mesh_p->Points->minColor=0;
841               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          mesh_p->Points->maxColor=chunkEle-1;
              faceTypeID=Finley_RefElement_getTypeId(element_type);  
              if (faceTypeID==NoType) {  
                sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);  
                Finley_setError(VALUE_ERROR,error_msg);  
              } else {  
                 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
                 if (Finley_noError()) {  
                    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
                    if (Finley_noError()) {  
                       mesh_p->FaceElements->minColor=0;  
                       mesh_p->FaceElements->maxColor=numEle-1;  
                       for (i0 = 0; i0 < numEle; i0++) {  
                         fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);  
                         mesh_p->FaceElements->Color[i0]=i0;  
                         for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {  
                              fscanf(fileHandle_p, " %d",  
                                 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);  
                         }   /* for i1 */  
                         fscanf(fileHandle_p, "\n");  
                       } /* for i0 */  
                    }  
                 }  
              }  
         }  
         /* get the Contact face element */  
842          if (Finley_noError()) {          if (Finley_noError()) {
843               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);            #pragma omp parallel for private (i0, i1)
844               contactTypeID=Finley_RefElement_getTypeId(element_type);        for (i0=0; i0<chunkEle; i0++) {
845               if (contactTypeID==NoType) {          mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
846                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);          mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
847                 Finley_setError(VALUE_ERROR,error_msg);              mesh_p->Points->Owner[i0]  =mpi_info->rank;
848               } else {              mesh_p->Points->Color[i0] = i0;
849                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          for (i1 = 0; i1 < numNodes; i1++) {
850                 if (Finley_noError()) {            mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
851                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);          }
852                     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 */  
                   }  
                }  
              }  
853          }          }
854          /* get the nodal element */  
855          if (Finley_noError()) {      TMPMEMFREE(tempInts);
856               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);        } /* end of Read the nodal element data */
857               pointTypeID=Finley_RefElement_getTypeId(element_type);  
858               if (pointTypeID==NoType) {        /* get the name tags */
859                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);        if (Finley_noError()) {
860                 Finley_setError(VALUE_ERROR,error_msg);          char *remainder=0, *ptr;
861            size_t len;
862            int tag_key;
863    
864    
865    
866            if (mpi_info->rank == 0) {  /* Master */
867          /* Read the word 'Tag' */
868          if (! feof(fileHandle_p)) {
869            scan_ret = fscanf(fileHandle_p, "%s\n", name);
870            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
871          }
872    
873    #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
874    
875          remainder = NULL;
876              len=0;
877          while (1)
878              {
879                 size_t malloc_chunk = 1024;
880                 size_t buff_size = 0;
881                 int ch;
882    
883                 ch = fgetc(fileHandle_p);
884                 if( ch == '\r' )
885                 {
886                    continue;
887               }               }
888               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);               if( len+1 > buff_size )
889               if (Finley_noError()) {               {
890                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);                  TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
891                  if (Finley_noError()) {               }
892                     mesh_p->Points->minColor=0;               if( ch == EOF )
893                     mesh_p->Points->maxColor=numEle-1;               {
894                     for (i0 = 0; i0 < numEle; i0++) {                  /* hit EOF */
895                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);                  remainder[len] = (char)0;
896                       mesh_p->Points->Color[i0]=i0;                  break;
                      for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
                          fscanf(fileHandle_p, " %d",  
                             &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
                      }  /* for i1 */  
                      fscanf(fileHandle_p, "\n");  
                    } /* for i0 */  
                 }  
897               }               }
898                 remainder[len] = (char)ch;
899                 len++;
900              }
901    #else
902          /* Read rest of file in one chunk, after using seek to find length */
903              {
904                 long cur_pos, end_pos;
905    
906                 cur_pos = ftell(fileHandle_p);
907                 fseek(fileHandle_p, 0L, SEEK_END);
908                 end_pos = ftell(fileHandle_p);
909                 fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
910                 remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
911                 if (! feof(fileHandle_p))
912                 {
913                    scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
914                                     sizeof(char), fileHandle_p);
915    
916                    FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
917                    remainder[end_pos-cur_pos] = 0;
918                }
919              }
920    #endif
921          len = strlen(remainder);
922          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
923          len = strlen(remainder);
924              TMPMEMREALLOC(remainder,remainder,len+1,char);
925          }          }
926          /* get the name tags */  #ifdef PASO_MPI
927          if (Finley_noError()) {  
928             if (feof(fileHandle_p) == 0) {          if (MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm) != MPI_SUCCESS)
929                fscanf(fileHandle_p, "%s\n", name);          {
930                while (feof(fileHandle_p) == 0) {             Finley_setError(PASO_MPI_ERROR,
931                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);                             "Finley_Mesh_read: broadcast of tag len failed");
932                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);             return NULL;
933                }          }
934             }      if (mpi_info->rank != 0) {
935          remainder = TMPMEMALLOC(len+1, char);
936          remainder[0] = 0;
937        }
938    
939            if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
940                MPI_SUCCESS)
941            {
942               Finley_setError(PASO_MPI_ERROR,
943                               "Finley_Mesh_read: broadcast of tags failed");
944               return NULL;
945          }          }
946  #endif /* comment out the rest of the un-implemented crap for now */  #endif
947        if (remainder[0]) {
948              ptr = remainder;
949              do {
950                sscanf(ptr, "%s %d\n", name, &tag_key);
951                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
952                ptr++;
953              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
954              TMPMEMFREE(remainder);
955        }
956          }
957    
958       }       }
959    
960       /* close file */       /* close file */
961       if (mpi_info->rank == 0) fclose(fileHandle_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
962    
# Line 660  printf("ksteube CPU=%d/%d Element typeID Line 965  printf("ksteube CPU=%d/%d Element typeID
965    
966       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
967       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
      return mesh_p; /* ksteube temp return for debugging */  
968    
969       /* that's it */       /* that's it */
970       #ifdef Finley_TRACE       #ifdef Finley_TRACE

Legend:
Removed from v.1713  
changed lines
  Added in v.2059

  ViewVC Help
Powered by ViewVC 1.1.26