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

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

  ViewVC Help
Powered by ViewVC 1.1.26