/[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 1402 by ksteube, Thu Jan 31 00:17: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;    index_t tag_key;
240      for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {    Finley_Mesh *mesh_p=NULL;
241           fscanf(fileHandle_p, " %d",    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
242              &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);    char error_msg[LenErrorMsg_MAX];
243      }   /* for i1 */    double time0=Finley_timer();
244      fscanf(fileHandle_p, "\n");    FILE *fileHandle_p = NULL;
245    } /* for i0 */    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
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 */
269    /* TODO */    if (mpi_info->size > 0) {
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 */  
288    
289         /* allocate mesh */
290         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
291         if (Finley_noError()) {
292        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1, mpi_error;
293        int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
294        double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
295    
296        /*
297          Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
298          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
299          First chunk sent to CPU 1, second to CPU 2, ...
300          Last chunk stays on CPU 0 (the master)
301          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
302        */
303    
304        if (mpi_info->rank == 0) {  /* Master */
305          for (;;) {            /* Infinite loop */
306            chunkNodes = 0;
307            for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
308            for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
309            for (i1=0; i1<chunkSize; i1++) {
310              if (totalNodes >= numNodes) break;
311                  if (1 == numDim)
312            fscanf(fileHandle_p, "%d %d %d %le\n",
313              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
314              &tempCoords[i1*numDim+0]);
315                  if (2 == numDim)
316            fscanf(fileHandle_p, "%d %d %d %le %le\n",
317              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
318              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
319                  if (3 == numDim)
320            fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
321              &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
322              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
323              totalNodes++;
324              chunkNodes++;
325            }
326            /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
327            if (nextCPU < mpi_info->size) {
328    #ifdef PASO_MPI
329              tempInts[numNodes*3] = chunkNodes;
330              /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
331              mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
332              if ( mpi_error != MPI_SUCCESS ) {
333                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
334                    return NULL;
335              }
336              mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
337              if ( mpi_error != MPI_SUCCESS ) {
338                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
339                    return NULL;
340              }
341    #endif
342              nextCPU++;
343            }
344            if (totalNodes >= numNodes) break;
345          } /* Infinite loop */
346        }   /* End master */
347        else {  /* Worker */
348    #ifdef PASO_MPI
349          /* Each worker receives two messages */
350          MPI_Status status;
351          mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
352          if ( mpi_error != MPI_SUCCESS ) {
353                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
354                return NULL;
355          }
356          mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
357          if ( mpi_error != MPI_SUCCESS ) {
358                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
359                return NULL;
360          }
361          chunkNodes = tempInts[numNodes*3];
362    #endif
363        }   /* Worker */
364    
365    #if 0
366            /* Display the temp mem for debugging */
367            printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
368            for (i0=0; i0<numNodes*3; i0++) {
369              printf(" %2d", tempInts[i0]);
370              if (i0%numNodes==numNodes-1) printf("\n");
371            }
372            printf("ksteube tempCoords:\n");
373            for (i0=0; i0<chunkNodes*numDim; i0++) {
374              printf(" %20.15e", tempCoords[i0]);
375              if (i0%numDim==numDim-1) printf("\n");
376            }
377    #endif
378    
379        printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
380        /* Copy node data from tempMem to mesh_p */
381            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
382            if (Finley_noError()) {
383          for (i0=0; i0<chunkNodes; i0++) {
384            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
385            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[numNodes+i0];
386            mesh_p->Nodes->Tag[i0]              = tempInts[numNodes*2+i0];
387            for (i1=0; i1<numDim; i1++) {
388              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
389            }
390              }
391            }
392    
393        TMPMEMFREE(tempInts);
394        TMPMEMFREE(tempCoords);
395    
396            /* read elements */
397    
398        /* Read the element typeID */
399            if (Finley_noError()) {
400          if (mpi_info->rank == 0) {
401                fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
402                typeID=Finley_RefElement_getTypeId(element_type);
403          }
404    #ifdef PASO_MPI
405          if (mpi_info->size > 0) {
406            int temp1[3];
407            temp1[0] = typeID;
408            temp1[1] = numEle;
409            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
410            if (mpi_error != MPI_SUCCESS) {
411              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
412              return NULL;
413            }
414            typeID = temp1[0];
415            numEle = temp1[1];
416          }
417    #endif
418              if (typeID==NoType) {
419                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
420                Finley_setError(VALUE_ERROR, error_msg);
421              }
422        }
423    
424          /* Read the element data */
425          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
426          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
427    
428          if (Finley_noError()) {
429        int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
430        int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1, mpi_error;
431        if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
432        if (mpi_info->rank == 0) {  /* Master */
433          for (;;) {            /* Infinite loop */
434            chunkEle = 0;
435            for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
436            for (i0=0; i0<chunkSize; i0++) {
437              if (totalEle >= numEle) break; /* End infinite loop */
438              fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
439              for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
440              fscanf(fileHandle_p, "\n");
441              totalEle++;
442              chunkEle++;
443            }
444            /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
445            if (nextCPU < mpi_info->size) {
446    #ifdef PASO_MPI
447              tempInts[numEle*(2+numNodes)] = chunkEle;
448    printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
449              mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
450              if ( mpi_error != MPI_SUCCESS ) {
451                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
452                    return NULL;
453              }
454    #endif
455              nextCPU++;
456            }
457            if (totalEle >= numEle) break; /* End infinite loop */
458          } /* Infinite loop */
459        }   /* End master */
460        else {  /* Worker */
461    #ifdef PASO_MPI
462          /* Each worker receives two messages */
463          MPI_Status status;
464    printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
465          mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
466          if ( mpi_error != MPI_SUCCESS ) {
467                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
468                return NULL;
469          }
470          chunkEle = tempInts[numEle*(2+numNodes)];
471    #endif
472        }   /* Worker */
473    #if 1
474        /* Display the temp mem for debugging */
475        printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
476        for (i0=0; i0<numEle*(numNodes+2); i0++) {
477          printf(" %2d", tempInts[i0]);
478          if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
479        }
480    #endif
481    
482        /* Copy Element data from tempInts to mesh_p */
483        Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
484            mesh_p->Elements->minColor=0;
485            mesh_p->Elements->maxColor=chunkEle-1;
486            if (Finley_noError()) {
487              #pragma omp parallel for private (i0, i1)
488          for (i0=0; i0<chunkEle; i0++) {
489            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
490            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
491            for (i1 = 0; i1 < numNodes; i1++) {
492              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
493            }
494              }
495            }
496    
497        TMPMEMFREE(tempInts);
498          }
499    
500    
501    #if 0 /* this is the original code for reading elements */
502            /* read elements */
503            if (Finley_noError()) {
504    
505               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
506               typeID=Finley_RefElement_getTypeId(element_type);
507               if (typeID==NoType) {
508                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
509                 Finley_setError(VALUE_ERROR,error_msg);
510               } else {
511                 /* read the elements */
512                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
513                 if (Finley_noError()) {
514                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
515                     mesh_p->Elements->minColor=0;
516                     mesh_p->Elements->maxColor=numEle-1;
517                     if (Finley_noError()) {
518                        for (i0 = 0; i0 < numEle; i0++) {
519                          fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
520                          mesh_p->Elements->Color[i0]=i0;
521                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
522                               fscanf(fileHandle_p, " %d",
523                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
524                          } /* for i1 */
525                          fscanf(fileHandle_p, "\n");
526                        } /* for i0 */
527                     }
528                 }
529              }
530            }
531    #endif
532    
533    printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
534    
535    /* close file */  #if 1
536    
537    fclose(fileHandle_p);    /* Define other structures to keep mesh_write from crashing */
538    
539    /*   resolve id's : */    mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
540      Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
541    
542    Finley_Mesh_resolveNodeIds(mesh_p);    mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
543      Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
544    
545    /* rearrange elements: */    mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
546      Finley_ElementFile_allocTable(mesh_p->Points, 0);
547    
548    Finley_Mesh_prepare(mesh_p);  #endif
549    
550    #if 0 /* comment out the rest of the un-implemented crap for now */
551            /* get the face elements */
552            if (Finley_noError()) {
553                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
554                 faceTypeID=Finley_RefElement_getTypeId(element_type);
555                 if (faceTypeID==NoType) {
556                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
557                   Finley_setError(VALUE_ERROR,error_msg);
558                 } else {
559                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
560                    if (Finley_noError()) {
561                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
562                       if (Finley_noError()) {
563                          mesh_p->FaceElements->minColor=0;
564                          mesh_p->FaceElements->maxColor=numEle-1;
565                          for (i0 = 0; i0 < numEle; i0++) {
566                            fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
567                            mesh_p->FaceElements->Color[i0]=i0;
568                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
569                                 fscanf(fileHandle_p, " %d",
570                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
571                            }   /* for i1 */
572                            fscanf(fileHandle_p, "\n");
573                          } /* for i0 */
574                       }
575                    }
576                 }
577            }
578            /* get the Contact face element */
579            if (Finley_noError()) {
580                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
581                 contactTypeID=Finley_RefElement_getTypeId(element_type);
582                 if (contactTypeID==NoType) {
583                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
584                   Finley_setError(VALUE_ERROR,error_msg);
585                 } else {
586                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
587                   if (Finley_noError()) {
588                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
589                       if (Finley_noError()) {
590                          mesh_p->ContactElements->minColor=0;
591                          mesh_p->ContactElements->maxColor=numEle-1;
592                          for (i0 = 0; i0 < numEle; i0++) {
593                            fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
594                            mesh_p->ContactElements->Color[i0]=i0;
595                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
596                                fscanf(fileHandle_p, " %d",
597                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
598                            }   /* for i1 */
599                            fscanf(fileHandle_p, "\n");
600                          } /* for i0 */
601                      }
602                   }
603                 }
604            }
605            /* get the nodal element */
606            if (Finley_noError()) {
607                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
608                 pointTypeID=Finley_RefElement_getTypeId(element_type);
609                 if (pointTypeID==NoType) {
610                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
611                   Finley_setError(VALUE_ERROR,error_msg);
612                 }
613                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
614                 if (Finley_noError()) {
615                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
616                    if (Finley_noError()) {
617                       mesh_p->Points->minColor=0;
618                       mesh_p->Points->maxColor=numEle-1;
619                       for (i0 = 0; i0 < numEle; i0++) {
620                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
621                         mesh_p->Points->Color[i0]=i0;
622                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
623                             fscanf(fileHandle_p, " %d",
624                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
625                         }  /* for i1 */
626                         fscanf(fileHandle_p, "\n");
627                       } /* for i0 */
628                    }
629                 }
630            }
631            /* get the name tags */
632            if (Finley_noError()) {
633               if (feof(fileHandle_p) == 0) {
634                  fscanf(fileHandle_p, "%s\n", name);
635                  while (feof(fileHandle_p) == 0) {
636                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
637                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
638                  }
639               }
640            }
641    #endif /* comment out the rest of the un-implemented crap for now */
642         }
643         /* close file */
644         if (mpi_info->rank == 0) fclose(fileHandle_p);
645    
646         /*   resolve id's : */
647         /* rearrange elements: */
648    
649         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
650         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
651         return mesh_p; /* ksteube temp return for debugging */
652    
653         /* that's it */
654         #ifdef Finley_TRACE
655         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
656         #endif
657    
658    /* that's it */    /* that's it */
659    #ifdef Finley_TRACE    if (! Finley_noError()) {
660    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);  
661    }    }
662      Paso_MPIInfo_free( mpi_info );
663    return mesh_p;    return mesh_p;
664  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26