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

Legend:
Removed from v.1628  
changed lines
  Added in v.2052

  ViewVC Help
Powered by ViewVC 1.1.26