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

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

  ViewVC Help
Powered by ViewVC 1.1.26