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

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

  ViewVC Help
Powered by ViewVC 1.1.26