/[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 1028 by gross, Wed Mar 14 00:15:24 2007 UTC revision 1725 by ksteube, Tue Aug 26 03:31:49 2008 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2  /**************************************************************/  /* $Id$ */
3    
4  /*   Finley: read mesh */  /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*   Author: gross@access.edu.au */  /*   Finley: read mesh */
 /*   Version: $Id$ */  
19    
20  /**************************************************************/  /**************************************************************/
21    
# Line 27  Line 25 
25    
26  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Finley file of name fname */
27    
28  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
29    
30    {
31    
32      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
33    dim_t numNodes, numDim, numEle, i0, i1;    dim_t numNodes, numDim, numEle, i0, i1;
34      index_t tag_key;
35    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
36    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
37    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
38    double time0=Finley_timer();    double time0=Finley_timer();
39    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
40    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
41      
42    Finley_resetError();    Finley_resetError();
43    
44    /* get file handle */    if (mpi_info->size > 1) {
45    fileHandle_p = fopen(fname, "r");       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");
46    if (fileHandle_p==NULL) {    } else {
47      sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);       /* get file handle */
48      Finley_setError(IO_ERROR,error_msg);       fileHandle_p = fopen(fname, "r");
49      return NULL;       if (fileHandle_p==NULL) {
50           sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
51           Finley_setError(IO_ERROR,error_msg);
52           Paso_MPIInfo_free( mpi_info );
53           return NULL;
54         }
55      
56         /* read header */
57         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58         fscanf(fileHandle_p, frm, name);
59      
60         /* get the nodes */
61      
62         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63         /* allocate mesh */
64         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
65         if (Finley_noError()) {
66      
67            /* read nodes */
68            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
69            if (Finley_noError()) {
70               if (1 == numDim) {
71                   for (i0 = 0; i0 < numNodes; i0++)
72                    fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
73                           &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
74                           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
75               } else if (2 == numDim) {
76                        for (i0 = 0; i0 < numNodes; i0++)
77                              fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
78                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
79                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
80                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
81               } else if (3 == numDim) {
82                        for (i0 = 0; i0 < numNodes; i0++)
83                              fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
84                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
85                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
86                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
87                                     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
88               } /* if else else */
89            }
90            /* read elements */
91            if (Finley_noError()) {
92      
93               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
94               typeID=Finley_RefElement_getTypeId(element_type);
95               if (typeID==NoType) {
96                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
97                 Finley_setError(VALUE_ERROR,error_msg);
98               } else {
99                 /* read the elements */
100                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
101                 if (Finley_noError()) {
102                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
103                     mesh_p->Elements->minColor=0;
104                     mesh_p->Elements->maxColor=numEle-1;
105                     if (Finley_noError()) {
106                        for (i0 = 0; i0 < numEle; i0++) {
107                          fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
108                          mesh_p->Elements->Color[i0]=i0;
109                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
110                               fscanf(fileHandle_p, " %d",
111                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
112                          } /* for i1 */
113                          fscanf(fileHandle_p, "\n");
114                        } /* for i0 */
115                     }
116                 }
117              }
118            }
119            /* get the face elements */
120            if (Finley_noError()) {
121                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
122                 faceTypeID=Finley_RefElement_getTypeId(element_type);
123                 if (faceTypeID==NoType) {
124                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
125                   Finley_setError(VALUE_ERROR,error_msg);
126                 } else {
127                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
128                    if (Finley_noError()) {
129                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
130                       if (Finley_noError()) {
131                          mesh_p->FaceElements->minColor=0;
132                          mesh_p->FaceElements->maxColor=numEle-1;
133                          for (i0 = 0; i0 < numEle; i0++) {
134                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
135                            mesh_p->FaceElements->Color[i0]=i0;
136                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
137                                 fscanf(fileHandle_p, " %d",
138                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
139                            }   /* for i1 */
140                            fscanf(fileHandle_p, "\n");
141                          } /* for i0 */
142                       }
143                    }
144                 }
145            }
146            /* get the Contact face element */
147            if (Finley_noError()) {
148                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
149                 contactTypeID=Finley_RefElement_getTypeId(element_type);
150                 if (contactTypeID==NoType) {
151                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
152                   Finley_setError(VALUE_ERROR,error_msg);
153                 } else {
154                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
155                   if (Finley_noError()) {
156                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
157                       if (Finley_noError()) {
158                          mesh_p->ContactElements->minColor=0;
159                          mesh_p->ContactElements->maxColor=numEle-1;
160                          for (i0 = 0; i0 < numEle; i0++) {
161                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
162                            mesh_p->ContactElements->Color[i0]=i0;
163                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
164                                fscanf(fileHandle_p, " %d",
165                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
166                            }   /* for i1 */
167                            fscanf(fileHandle_p, "\n");
168                          } /* for i0 */
169                      }
170                   }
171                 }
172            }  
173            /* get the nodal element */
174            if (Finley_noError()) {
175                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
176                 pointTypeID=Finley_RefElement_getTypeId(element_type);
177                 if (pointTypeID==NoType) {
178                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
179                   Finley_setError(VALUE_ERROR,error_msg);
180                 }
181                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
182                 if (Finley_noError()) {
183                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
184                    if (Finley_noError()) {
185                       mesh_p->Points->minColor=0;
186                       mesh_p->Points->maxColor=numEle-1;
187                       for (i0 = 0; i0 < numEle; i0++) {
188                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
189                         mesh_p->Points->Color[i0]=i0;
190                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
191                             fscanf(fileHandle_p, " %d",
192                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
193                         }  /* for i1 */
194                         fscanf(fileHandle_p, "\n");
195                       } /* for i0 */
196                    }
197                 }
198            }
199            /* get the name tags */
200            if (Finley_noError()) {
201               if (feof(fileHandle_p) == 0) {
202                  fscanf(fileHandle_p, "%s\n", name);
203                  while (feof(fileHandle_p) == 0) {
204                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
205                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
206                  }
207               }
208            }
209         }
210         /* close file */
211         fclose(fileHandle_p);
212      
213         /*   resolve id's : */
214         /* rearrange elements: */
215      
216         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
217         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
218      
219         /* that's it */
220         #ifdef Finley_TRACE
221         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
222         #endif
223    }    }
224    
225    /* read header */    /* that's it */
226    sprintf(frm,"%%%d[^\n]",LenString_MAX-1);    if (! Finley_noError()) {
227    fscanf(fileHandle_p, frm, name);         Finley_Mesh_free(mesh_p);
   
   /* get the nodes */  
   
   fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
   /* allocate mesh */  
 #ifndef PASO_MPI  
   mesh_p = Finley_Mesh_alloc(name,numDim,order);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
 #ifndef PASO_MPI  
   Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
   if (! Finley_noError()) return NULL;  
 #else  
   /* TODO */  
 #endif  
   
   if (1 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);  
   } else if (2 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);  
   } else if (3 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);  
   } /* if else else */  
   
   /* get the element type */  
   
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   typeID=Finley_RefElement_getTypeId(element_type);  
   if (typeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   /* read the elements */  
 #ifndef PASO_MPI  
   mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);  
   mesh_p->Elements->minColor=0;  
   mesh_p->Elements->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);  
     mesh_p->Elements->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {  
          fscanf(fileHandle_p, " %d",  
             &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   /* get the face elements */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   faceTypeID=Finley_RefElement_getTypeId(element_type);  
   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);  
     return NULL;  
228    }    }
229  #ifndef PASO_MPI    Paso_MPIInfo_free( mpi_info );
230    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);    return mesh_p;
231  #else  }
232    /* TODO */  
233  #endif  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
234    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
235    mesh_p->FaceElements->minColor=0;  {
236    mesh_p->FaceElements->maxColor=numEle-1;  
237    for (i0 = 0; i0 < numEle; i0++) {    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
238      fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);    dim_t numNodes, numDim, numEle, i0, i1;
239      mesh_p->FaceElements->Color[i0]=i0;    Finley_Mesh *mesh_p=NULL;
240      for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
241           fscanf(fileHandle_p, " %d",    char error_msg[LenErrorMsg_MAX];
242              &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);    double time0=Finley_timer();
243      }   /* for i1 */    FILE *fileHandle_p = NULL;
244      fscanf(fileHandle_p, "\n");    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
245    } /* for i0 */    index_t tag_key;
246    
247    /* get the Contact face element */    Finley_resetError();
248    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
249    contactTypeID=Finley_RefElement_getTypeId(element_type);    if (mpi_info->rank == 0) {
250    if (contactTypeID==NoType) {       /* get file handle */
251      sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);       fileHandle_p = fopen(fname, "r");
252      Finley_setError(VALUE_ERROR,error_msg);       if (fileHandle_p==NULL) {
253      return NULL;         sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
254           Finley_setError(IO_ERROR,error_msg);
255           Paso_MPIInfo_free( mpi_info );
256           return NULL;
257         }
258    
259         /* read header */
260         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
261         fscanf(fileHandle_p, frm, name);
262    
263         /* get the number of nodes */
264         fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
265    }    }
266  #ifndef PASO_MPI  
267    mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);  #ifdef PASO_MPI
268  #else    /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
269    /* TODO */    if (mpi_info->size > 1) {
270  #endif      int temp1[3], error_code;
271    Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);      temp1[0] = numDim;
272    mesh_p->ContactElements->minColor=0;      temp1[1] = numNodes;
273    mesh_p->ContactElements->maxColor=numEle-1;      temp1[2] = strlen(name) + 1;
274    for (i0 = 0; i0 < numEle; i0++) {      error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
275      fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);      if (error_code != MPI_SUCCESS) {
276      mesh_p->ContactElements->Color[i0]=i0;        Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
277      for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {        return NULL;
278          fscanf(fileHandle_p, " %d",      }
279             &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);      numDim = temp1[0];
280      }   /* for i1 */      numNodes = temp1[1];
281      fscanf(fileHandle_p, "\n");      error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
282    } /* for i0 */      if (error_code != MPI_SUCCESS) {
283          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
284    /* get the nodal element */        return NULL;
285    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);      }
   pointTypeID=Finley_RefElement_getTypeId(element_type);  
   if (pointTypeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
286    }    }
 #ifndef PASO_MPI  
   mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);  
 #else  
   /* TODO */  
287  #endif  #endif
   Finley_ElementFile_allocTable(mesh_p->Points, numEle);  
   mesh_p->Points->minColor=0;  
   mesh_p->Points->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
     mesh_p->Points->Color[i0]=i0;  
     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 */  
   
   
   /* close file */  
288    
289    fclose(fileHandle_p);       /* allocate mesh */
290         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
291         if (Finley_noError()) {
292        /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
293        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
294        int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
295        double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
296    
297        /*
298          Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
299          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
300          First chunk sent to CPU 1, second to CPU 2, ...
301          Last chunk stays on CPU 0 (the master)
302          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
303        */
304    
305        if (mpi_info->rank == 0) {  /* Master */
306          for (;;) {            /* Infinite loop */
307            for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
308            for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
309            chunkNodes = 0;
310            for (i1=0; i1<chunkSize; i1++) {
311              if (totalNodes >= numNodes) break;    /* Maybe end the infinite loop */
312                  if (1 == numDim)
313            fscanf(fileHandle_p, "%d %d %d %le\n",
314              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
315              &tempCoords[i1*numDim+0]);
316                  if (2 == numDim)
317            fscanf(fileHandle_p, "%d %d %d %le %le\n",
318              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
319              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
320                  if (3 == numDim)
321            fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
322              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
323              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
324              totalNodes++; /* When do we quit the infinite loop? */
325              chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
326            }
327            if (chunkNodes > chunkSize) {
328                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
329                  return NULL;
330            }
331    #ifdef PASO_MPI
332            /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
333            if (nextCPU < mpi_info->size) {
334                  int mpi_error;
335              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
336              mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
337              if ( mpi_error != MPI_SUCCESS ) {
338                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
339                    return NULL;
340              }
341              mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
342              if ( mpi_error != MPI_SUCCESS ) {
343                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
344                    return NULL;
345              }
346              nextCPU++;
347            }
348    #endif
349            if (totalNodes >= numNodes) break;  /* Maybe end the infinite loop */
350          } /* Infinite loop */
351        }   /* End master */
352        else {  /* Worker */
353    #ifdef PASO_MPI
354          /* Each worker receives two messages */
355          MPI_Status status;
356              int mpi_error;
357          mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
358          if ( mpi_error != MPI_SUCCESS ) {
359                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
360                return NULL;
361          }
362          mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
363          if ( mpi_error != MPI_SUCCESS ) {
364                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
365                return NULL;
366          }
367          chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
368    #endif
369        }   /* Worker */
370    
371    /*   resolve id's : */  #if 0
372            /* Display the temp mem for debugging */
373            printf("ksteube Nodes tempInts\n");
374            for (i0=0; i0<chunkSize*3; i0++) {
375              printf(" %2d", tempInts[i0]);
376              if (i0%chunkSize==chunkSize-1) printf("\n");
377            }
378            printf("ksteube tempCoords:\n");
379            for (i0=0; i0<chunkSize*numDim; i0++) {
380              printf(" %20.15e", tempCoords[i0]);
381              if (i0%numDim==numDim-1) printf("\n");
382            }
383                printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);
384    #endif
385    
386    Finley_Mesh_resolveNodeIds(mesh_p);      /* Copy node data from tempMem to mesh_p */
387            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
388            if (Finley_noError()) {
389          for (i0=0; i0<chunkNodes; i0++) {
390            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
391            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
392            mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
393            for (i1=0; i1<numDim; i1++) {
394              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
395            }
396              }
397            }
398    
399        TMPMEMFREE(tempInts);
400        TMPMEMFREE(tempCoords);
401    
402            /* read elements */
403    
404        /* Read the element typeID */
405            if (Finley_noError()) {
406          if (mpi_info->rank == 0) {
407                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
408                typeID=Finley_RefElement_getTypeId(element_type);
409          }
410    #ifdef PASO_MPI
411          if (mpi_info->size > 1) {
412            int temp1[2], mpi_error;
413            temp1[0] = (int) typeID;
414            temp1[1] = numEle;
415            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
416            if (mpi_error != MPI_SUCCESS) {
417              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
418              return NULL;
419            }
420            typeID = (ElementTypeId) temp1[0];
421            numEle = temp1[1];
422          }
423    #endif
424              if (typeID==NoType) {
425                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
426                Finley_setError(VALUE_ERROR, error_msg);
427              }
428        }
429    
430          /* Read the element data */
431          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
432          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
433    
434          if (Finley_noError()) {
435        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
436        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
437        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
438        if (mpi_info->rank == 0) {  /* Master */
439          for (;;) {            /* Infinite loop */
440            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
441            chunkEle = 0;
442            for (i0=0; i0<chunkSize; i0++) {
443              if (totalEle >= numEle) break; /* End infinite loop */
444              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
445              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
446              fscanf(fileHandle_p, "\n");
447              totalEle++;
448              chunkEle++;
449            }
450    #ifdef PASO_MPI
451            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
452            if (nextCPU < mpi_info->size) {
453                  int mpi_error;
454              tempInts[chunkSize*(2+numNodes)] = chunkEle;
455              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
456              if ( mpi_error != MPI_SUCCESS ) {
457                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
458                    return NULL;
459              }
460              nextCPU++;
461            }
462    #endif
463            if (totalEle >= numEle) break; /* End infinite loop */
464          } /* Infinite loop */
465        }   /* End master */
466        else {  /* Worker */
467    #ifdef PASO_MPI
468          /* Each worker receives one message */
469          MPI_Status status;
470          int mpi_error;
471          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
472          if ( mpi_error != MPI_SUCCESS ) {
473                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
474                return NULL;
475          }
476          chunkEle = tempInts[chunkSize*(2+numNodes)];
477    #endif
478        }   /* Worker */
479    
480    /* rearrange elements: */  #if 0
481        /* Display the temp mem for debugging */
482        for (i0=0; i0<chunkSize*(numNodes+2); i0++) {
483          printf(" %2d", tempInts[i0]);
484          if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
485        }
486    #endif
487    
488    Finley_Mesh_prepare(mesh_p);      /* Copy Element data from tempInts to mesh_p */
489        Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
490            mesh_p->Elements->minColor=0;
491            mesh_p->Elements->maxColor=chunkEle-1;
492            if (Finley_noError()) {
493              #pragma omp parallel for private (i0, i1)
494          for (i0=0; i0<chunkEle; i0++) {
495            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
496            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
497                mesh_p->Elements->Color[i0] = i0;
498            for (i1 = 0; i1 < numNodes; i1++) {
499              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
500            }
501              }
502            }
503    
504        TMPMEMFREE(tempInts);
505          }
506         }
507    
508         /* close file */
509         if (mpi_info->rank == 0) fclose(fileHandle_p);
510    
511         /*   resolve id's : */
512         /* rearrange elements: */
513    
514         return mesh_p; /* ksteube temp return for debugging */
515    
516         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
517         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
518    
519         /* that's it */
520         #ifdef Finley_TRACE
521         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
522         #endif
523    
524    /* that's it */    /* that's it */
525    #ifdef Finley_TRACE    if (! Finley_noError()) {
526    printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);         Finley_Mesh_free(mesh_p);
   #endif  
   if (Finley_noError()) {  
        if ( ! Finley_Mesh_isPrepared(mesh_p)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
        Finley_Mesh_dealloc(mesh_p);  
527    }    }
528      Paso_MPIInfo_free( mpi_info );
529    return mesh_p;    return mesh_p;
530  }  }

Legend:
Removed from v.1028  
changed lines
  Added in v.1725

  ViewVC Help
Powered by ViewVC 1.1.26